#multiple covariates measurement error model for NHIS/NHANES to predict BMI: #Structural ME approach: #Two ME (reported BMI and CHOLES from NHIS & one covariate W_ij(SMOKE) from NHANES ): rm(list=ls(all=TRUE)) library(sas7bdat) library(MASS) getwd() setwd() dataNHIS <- read.sas7bdat("samadult_nhis2014.sas7bdat") #n=35928 attach(dataNHIS) summary(dataNHIS) dim(dataNHIS) sumt=35928 m=50 BMI.x<- rep(0,sumt) CHOLES.x<-rep(0,sumt) smoke.tem<-rep(0,sumt) #----------------- a<-1 to<-rep(0,m) for(i in 1:2){ for(j in 1:5){ for(k in 1:5){ for(l in 1:35928){ if(SEX[l]==i & agecat[l]==j & RIDRETH1[l]==k){ to[a]<-to[a]+1 } } a<-a+1 } } } #----------------- cum<-rep(0,m) cum[1]<-to[1] for(i in 2:m){ cum[i]<-cum[i-1]+to[i] } a<-1 b<-to[1] county<-rep(0,sumt) for(i in 1:m){ county[a:b]<-rep(i,to[i]) a<-b+1 b<-a+to[i+1]-1 } #----------------- jj=1 for(i in 1:2){ for(j in 1:5){ for(k in 1:5){ for(l in 1:35928){ if(SEX[l]==i & agecat[l]==j & RIDRETH1[l]==k){ BMI.x[jj]<-BMI[l] CHOLES.x[jj]<-CHLEV_R[l] #131 NA smoke.tem[jj]<-currentsmoker[l] #199 NA jj<-jj+1 } } } } } #----------- CHOLES1.x<-CHOLES.x CHOLES.x.tem1 <- CHOLES1.x[!is.na(CHOLES1.x)] #to remove a vector with NA's summary(CHOLES.x.tem1) CHOLES1.x[is.na(CHOLES1.x)]<-0 #to replace NA's with 0 summary(CHOLES1.x) smoke1.tem<-smoke.tem smoke1.tem1 <- smoke1.tem[!is.na(smoke1.tem)] #to remove a vector with NA's summary(smoke1.tem1) smoke1.tem[is.na(smoke1.tem)]<-0 #to replace NA's with 0 summary(smoke1.tem) BMI.x.mean<-rep(0,m) CHOLES.x.mean<-rep(0,m) smoke.tem.mean<-rep(0,m) a<-1 b<-to[1] for(i in 1:m){ BMI.x.mean[i]<-mean(BMI.x[a:b]) CHOLES.x.mean[i]<-mean(CHOLES1.x[a:b]) smoke.tem.mean[i]<-mean(smoke1.tem[a:b]) a<-b+1 b<-a+to[i+1]-1 } ########################################## dataNHANES <- read.sas7bdat("nhanes_bmx201314.sas7bdat") #n=5588 attach(dataNHANES) summary(dataNHANES) dim(dataNHANES) no<-rep(0,m) n<-length(SEX) BMI.y<-rep(0,n) #waist.w<-rep(0,n) smoker.w<-rep(0,n) #BPsys.w<-rep(0,n) a<-1 b1<-1 for(i in 1:2){ for(j in 1:5){ for(k in 1:5){ for(l in 1:5588){ if(SEX[l]==i & agecat[l]==j & RIDRETH1[l]==k){ no[a]<-no[a]+1 BMI.y[b1]<-BMXBMI[l] #some NA #waist.w[b1]<-BMXWAIST[l] #some NA smoker.w[b1]<-currentsmoker[l] #2 NA's #BPsys.w[b1]<-BPXSY23[l] #200 NA's b1<-b1+1 } } a<-a+1 } } } #----------------- BMI1.y<-BMI.y BMI.y.tem1 <- BMI1.y[!is.na(BMI1.y)] #to remove a vector with NA's summary(BMI.y.tem1) BMI1.y[is.na(BMI1.y)]<-29.1 #to replace NA's with mean of other elements of the vector summary(BMI1.y) #waist1.w<-waist.w #waist.w.tem1 <- waist1.w[!is.na(waist1.w)] #to remove a vector with NA's #summary(waist.w.tem1) #waist1.w[is.na(waist1.w)]<-99.19 #to replace NA's with mean of other elements of the vector #summary(waist1.w) smoker1.w<-smoker.w smoker.w.tem1 <- smoker1.w[!is.na(smoker1.w)] #to remove a vector with NA's summary(smoker.w.tem1) smoker1.w[is.na(smoker1.w)]<-1 #to replace NA's with 1 summary(smoker1.w) #BPsys1.w<-BPsys.w #BPsys.w.tem1 <- BPsys1.w[!is.na(BPsys1.w)] #to remove a vector with NA's #summary(BPsys.w.tem1) #BPsys1.w[is.na(BPsys1.w)]<-123.2 #to replace NA's with mean of other elements of the vector #summary(BPsys1.w) #------------------ #waist.w.mean<-rep(0,m) #BPsys.w.mean<-rep(0,m) smoker.w.mean<-rep(0,m) a<-1 b<-no[1] alpha<-0 n.tem<-0 sum.xt<-0 for(i in 1:m){ # waist.w.mean[i]<-mean(waist1.w[a:b]) # BPsys.w.mean[i]<-mean(BPsys1.w[a:b]) smoker.w.mean[i]<-mean(smoker1.w[a:b]) a<-b+1 b<-a+no[i+1]-1 alpha<-alpha+(no[i]*smoker.w.mean[i])/n n.tem<-n.tem+no[i]*(no[i]-1) sum.xt<-sum.xt+(no[i]*smoker.w.mean[i])^2 } a.hat<-((n.tem*alpha*(1-alpha))/(sum.xt-n.tem*alpha*alpha-n*alpha))-1 #func.a<-function(a) #{ans<-(n.tem*alpha*(1-alpha))/(a+1)+n.tem*alpha*alpha+n*alpha # return(ans) #} #fit <- optim(a, func.a, method="L-BFGS-B", lower=c(0.01), # upper=c(20), hessian = TRUE) #, control=list(fnscale=1)) #, factr=0.1)) #control = list(maxit=1000),trace=1 #a.hat<-fit$par smok.pop.mean<-rep(0,m) for(i in 1:m){ smok.pop.mean[i]<-(no[i]*smoker.w.mean[i]+a.hat*alpha)/(no[i]+a.hat) } ############################################### #response var: BMI.y #covariates: smoker.w, BMI.x.mean, CHOLES.x.mean y.ran<-BMI1.y p1=1 w.ij.ran<-rbind(smoker1.w) #(BPsys1.w ,waist1.w, smoker1.w) w.i.bar.ran<-rbind(smoker.w.mean) w.bar.ran<-rep(0,p1) for(i in 1:m){ w.bar.ran<-w.bar.ran+(no[i]/n)*w.i.bar.ran[,i] } p2=2 X.ran<-matrix(c(BMI.x,CHOLES1.x),p2,sumt,byrow=T) X1.ran<-X.ran[1,] X2.ran<-X.ran[2,] ################################ I=m+1 s<-rep(0,I) st<-rep(0,I) s[1]<-1 st[1]<-1 for(i in 2:I){ s[i]<-s[i-1]+no[i-1] st[i]<-st[i-1]+to[i-1] } s st #---------------------------------------------- X.X<-matrix(0,p1+p2+1,p1+p2+1) weight<-rep(0,m) for(i in 1:m){ a<-st[i] b<-st[i+1]-1 a1<-s[i] b1<-s[i+1]-1 X.X<-X.X+c(1,w.i.bar.ran[1,i], mean(X1.ran[a:b]), mean(X2.ran[a:b]) )%*%t(c(1,w.i.bar.ran[1,i], mean(X1.ran[a:b]), mean(X2.ran[a:b]))) } #when p2=2 #print("E") for(i in 1:m){ a<-st[i] b<-st[i+1]-1 #a1<-s[i] #b1<-s[i+1]-1 temp<-c(1,w.i.bar.ran[1,i], mean(X1.ran[a:b]), mean(X2.ran[a:b])) #when p2=2 weight[i]<-1-t(temp)%*%solve(X.X)%*%temp} #print("F") ############################### ybar.ran<-rep(0,m) Xbar.ran<-matrix(0,p2,m) sum.y<-0 sum.X<-rep(0,p2) XBARR<-rep(0,p2) for(i in 1:m){ a<-s[i] b<-s[i+1]-1 at<-st[i] bt<-st[i+1]-1 ybar.ran[i]<-mean(y.ran[a:b]) Xbar.ran[1,i]<-mean(X1.ran[at:bt]) #when p2=2 Xbar.ran[2,i]<-mean(X2.ran[at:bt]) sum.y<-sum.y+no[i]*ybar.ran[i] sum.X[1]<-sum.X[1]+to[i]*Xbar.ran[1,i] #when p2=2 sum.X[2]<-sum.X[2]+to[i]*Xbar.ran[2,i] } ybarr<-sum.y/n XBARR[1]<-sum.X[1]/sumt #estimate of mu_x (FME) XBARR[2]<-sum.X[2]/sumt #MUx.p[r,]<-XBARR ##################################### sum.XX<-rep(0,p2) mXX<-rep(0,p2) for(i in 1:m){ a<-st[i] for(k in 1:to[i]){ sum.XX[1]<-sum.XX[1]+(X1.ran[a+k-1]-Xbar.ran[1,i])^2 sum.XX[2]<-sum.XX[2]+(X2.ran[a+k-1]-Xbar.ran[2,i])^2 } } mXX[1]<-sum.XX[1]/(sumt-m) mXX[2]<-sum.XX[2]/(sumt-m) sum.yy<-0 sum.ww<-matrix(0,p1,p1) msb.w<-matrix(0,p1,p1) for(i in 1:m){ a<-s[i] for(k in 1:no[i]){ sum.yy<-sum.yy+(y.ran[a+k-1]-ybar.ran[i])^2 sum.ww<-sum.ww+(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]) } msb.w<-msb.w+(no[i]/(m-1))*(w.i.bar.ran[,i]-w.bar.ran)%*%t(w.i.bar.ran[,i]-w.bar.ran) } myy<-sum.yy/(n-m) mWW<-sum.ww/(n-m) #################################### vareta.p<-mXX #estimate of var(eta) when p2=2 #################################### varx.p<-matrix(0,p2,p2) sum.X<-matrix(0,p2,p2) #sum.A<-matrix(0,p2,p2) tem.g1<-0 for(i in 1:m){ sum.X<-sum.X+to[i]*(Xbar.ran[,i]-XBARR)%*%t(Xbar.ran[,i]-XBARR) #sum.A<-sum.A+diag(mXX,p2)/to[i] tem.g1<-tem.g1+to[i]^2 } g1m<-sumt-tem.g1/sumt varx.p<-((1/(m-1))*sum.X-mXX)*(m-1)/g1m #(1/m)*sum.A if(varx.p[1,1]<=0 | varx.p[2,2]<=0){ varx.p.check<-varx.p varx.p.check[1,1]<-ifelse(varx.p[1,1]<=0,mXX[1]/m^{2/3},varx.p[1,1]) varx.p.check[2,2]<-ifelse(varx.p[2,2]<=0,mXX[2]/m^{2/3},varx.p[2,2]) sigma.check.av<-(1/2)*(varx.p.check[1,1]+varx.p.check[2,2]) if(eigen(varx.p.check)$values[1]<0){ a1<-(2*sigma.check.av)/((eigen(varx.p.check)$values[2])+sigma.check.av/(sqrt(m))) eig1<-a1*sigma.check.av/(sqrt(m)) eig2<-a1*(eigen(varx.p.check)$values[2]) } if(eigen(varx.p.check)$values[2]<0){ a2<-(2*sigma.check.av)/((eigen(varx.p.check)$values[1])+sigma.check.av/(sqrt(m))) eig2<-a2*sigma.check.av/(sqrt(m)) eig1<-a2*(eigen(varx.p.check)$values[1]) } varx.p<-eig1*(eigen(varx.p.check)$vectors[1,])%*%t(eigen(varx.p.check)$vectors[1,])+ eig2*(eigen(varx.p.check)$vectors[2,])%*%t(eigen(varx.p.check)$vectors[2,]) } ################################### tem.w<-matrix(0,p1,p1) tem.y<-rep(0,p1) for(i in 1:m){ a<-s[i] for(j in 1:no[i]){ tem.w<-tem.w+(w.ij.ran[,a+j-1]-w.bar.ran)%*%t(w.ij.ran[,a+j-1]) tem.y<-tem.y+(w.ij.ran[,a+j-1]-w.bar.ran)*y.ran[a+j-1] } } beta1.p<-solve(tem.w)%*%tem.y ################################# sum.b2<-rep(0,p2) tem.g1<-0 sum.b22<-0 tem.g2<-0 temm<-0 for(i in 1:m){ a<-s[i] for(j in 1:no[i]){ temm<-temm+(to[i]*(1-to[i]/sumt)*t(w.i.bar.ran[,i])%*%solve(tem.w)%*%(w.ij.ran[,a+j-1]-w.bar.ran)) } sum.b2[1]<-sum.b2[1]+to[i]*ybar.ran[i]*(Xbar.ran[1,i]-XBARR[1]) sum.b2[2]<-sum.b2[2]+to[i]*ybar.ran[i]*(Xbar.ran[2,i]-XBARR[2]) tem.g1<-tem.g1+to[i]^2 tem.g2<-tem.g2+no[i]^2 sum.b22<-sum.b22+no[i]*(ybar.ran[i]-ybarr)^2 } g1m<-sumt-tem.g1/sumt g2m<-n-tem.g2/n beta2.p<-solve(sum.X/(m-1)-mXX)%*%(sum.b2/(m-1)) #g1m-temm)) beta0.p<-ybarr-t(beta1.p)%*%w.bar.ran-t(beta2.p)%*%XBARR ################################# tem.w.bar<-0 tem.w1<-0 for(i in 1:m){ a<-s[i] for(j in 1:no[i]){ tem.w.bar<-tem.w.bar+(t(beta1.p)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]))^2 tem.w1<-tem.w1+no[i]*(t(beta1.p)%*%(w.i.bar.ran[,i]-w.bar.ran))^2 } } varv.p<-(((sum.b22/(m-1)) - myy- t(beta1.p)%*% (msb.w-mWW) %*%beta1.p )*(m-1)/g2m) - t(beta2.p)%*%(sum.b2/(g1m-0)) #-(1/g2m)*(tem.w1-((m-1)/(n-m))*tem.w.bar) #varv.p<-ifelse(varv.p<=0, t(beta2.p)%*%(sum.b2/(g1m-temm)), varv.p ) if(varv.p<=0)varv.p<-0.0001 #estimate of var(u) # varv.p<-1.49 # taken from the funcional ME #--------------------- tem.yw<-0 tem.Aw1<-matrix(0,p1,p1) tem.Bw1<-matrix(0,p1,p1) A.w<-0 B.w<-0 for(i in 1:m){ a<-s[i] for(j in 1:no[i]){ tem.yw<-tem.yw+((y.ran[a+j-1]-ybar.ran[i])-t(beta1.p)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]))^2 tem.Aw1<-tem.Aw1+(w.ij.ran[,a+j-1]-w.bar.ran)%*%t(w.ij.ran[,a+j-1]) } tem.Bw1<-tem.Bw1+(no[i]^2)*(w.i.bar.ran[,i]-w.bar.ran)%*%t(w.i.bar.ran[,i]-w.bar.ran) } tem.Aw2<-0 for(i in 1:m){ a<-s[i] for(j in 1:no[i]){ tem.Aw2<-tem.Aw2+t(w.ij.ran[,a+j-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]) B.w<-B.w+t(w.ij.ran[,a+j-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%tem.Bw1%*%solve(tem.Aw1)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]) } } A.w<-(n-m)-tem.Aw2 vare.p<-tem.yw/(n-m-p1) #(1/A.w)*(tem.yw-(varv.p+t(beta2.p)%*%(sum.b2/(g1m-temm)))*B.w) ############################################################################### muhat.p<-rep(0,m) # mu.b.p<-rep(0,m) W.bar<-matrix(c(smok.pop.mean),1,m) #for now, we take the population mean of smoke from the NHIS(bigger survey): for(i in 1:m){ fi<-1 #(1-(no[i]/to[i])) #Ai<-vareta0/to[i] #Hi<-diag(p2)+Ai%*%solve(varx0) #Di.b<-vare0/(vare0+no[i]*(varv0+beta2%*%solve(Hi)%*%Ai%*%beta2)) #mu.b.p[i]<-(1-fi*Di.b)*(ybar.ran[i]+t(beta1)%*%(W.bar[,i]-w.i.bar.ran[,i]))+fi*Di.b*(beta0+t(beta1)%*%W.bar[,i]+t(beta2)%*%mux0)+fi*Di.b*(t(beta2)%*%varx0%*%solve(varx0+Ai)%*%(Xbar.ran[,i]-mux0) ) Ai.p<-diag(c(mXX[1]/to[i],mXX[2]/to[i])) varx.p<-diag(c(varx.p[1,1],varx.p[2,2])) Hi.p<-diag(p2)+Ai.p%*%solve(varx.p) Di.p<-vare.p/(vare.p+no[i]*(varv.p+t(beta2.p)%*%solve(Hi.p)%*%Ai.p%*%beta2.p)) muhat.p[i]<-(1-fi*Di.p)*(ybar.ran[i]+t(beta1.p)%*%(W.bar[,i]-w.i.bar.ran[,i]))+fi*Di.p*(beta0.p+t(beta1.p)%*%W.bar[,i]+t(beta2.p)%*%XBARR)+fi*Di.p*(t(beta2.p)%*%varx.p%*%solve(varx.p+Ai.p)%*%(Xbar.ran[,i]-XBARR) ) } #mu.b.p muhat.p #----------------------------------------------------------------- g1.p<-rep(0,m) for(i in 1:m){ Ai.p<-diag(c(mXX[1]/to[i],mXX[2]/to[i])) Hi.p<-diag(p2)+Ai.p%*%solve(varx.p) Di.p<-vare.p/no[i] Mi.p<-solve(Ai.p)+solve(varx.p) g1.p[i]<-(Di.p*(varv.p+t(beta2.p)%*%solve(Mi.p)%*%beta2.p))/(Di.p+varv.p+t(beta2.p)%*%solve(Mi.p)%*%beta2.p) } g1.p ##########################Jackknife method:####################### beta0.p.j<-rep(0,m) beta1.p.j<-matrix(0,p1,m) beta2.p.j<-matrix(0,p2,m) vareta.p.j<-matrix(0,p2,m) vare.p.j<-rep(0,m) varv.p.j<-rep(0,m) mux.p.j<-matrix(0,p2,m) varx.p.j<-array(0,c(p2,p2,m)) XBARR.j<-matrix(0,p2,m) ybarr.j<-rep(0,m) w.bar.j<-matrix(0,p1,m) #print(r) for(j in 1:m){ sum.y.j<-0 sum.X.j<-rep(0,p2) sum.w.j<-rep(0,p1) for(i in 1:m){if (i !=j){ sum.y.j<-sum.y.j+no[i]*ybar.ran[i] sum.X.j[1]<-sum.X.j[1]+to[i]*Xbar.ran[1,i] sum.X.j[2]<-sum.X.j[2]+to[i]*Xbar.ran[2,i] sum.w.j<-sum.w.j+no[i]*w.i.bar.ran[,i] } } ybarr.j[j]<-sum.y.j/(n-no[j]) XBARR.j[1,j]<-sum.X.j[1]/(sumt-to[j]) XBARR.j[2,j]<-sum.X.j[2]/(sumt-to[j]) w.bar.j[,j]<-sum.w.j/(n-no[j]) sum.XX.j<-rep(0,p2) for(i in 1:m){if(i !=j){ a<-st[i] for(k in 1:to[i]){ sum.XX.j[1]<-sum.XX.j[1]+(X1.ran[a+k-1]-Xbar.ran[1,i])^2 sum.XX.j[2]<-sum.XX.j[2]+(X2.ran[a+k-1]-Xbar.ran[2,i])^2} } } mXX.j<-rep(0,p2) mXX.j[1]<-sum.XX.j[1]/((sumt-to[j])-(m-1)) mXX.j[2]<-sum.XX.j[2]/((sumt-to[j])-(m-1)) sum.yy.j<-0 sum.ww.j<-matrix(0,p1,p1) msb.w.j<-matrix(0,p1,p1) for(i in 1:m){if(i !=j){ a<-s[i] for(k in 1:no[i]){ sum.yy.j<-sum.yy.j+(y.ran[a+k-1]-ybar.ran[i])^2 sum.ww.j<-sum.ww.j+(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]) } msb.w.j<-msb.w.j+(no[i]/(m-1-1))*(w.i.bar.ran[,i]-w.bar.j[,j])%*%t(w.i.bar.ran[,i]-w.bar.j[,j]) } } myy.j<-sum.yy.j/((n-no[j])-(m-1)) mWW.j<-sum.ww.j/((n-no[j])-(m-1)) #print(r) mux.p.j[,j]<-XBARR.j[,j] vareta.p.j[,j]<-mXX.j ################################### sum.X<-matrix(0,p2,p2) sum.A<-matrix(0,p2,p2) tem.g1<-0 for(i in 1:m){if(i !=j){ sum.X<-sum.X+to[i]*(Xbar.ran[,i]-XBARR.j[,j])%*%t(Xbar.ran[,i]-XBARR.j[,j]) sum.A<-sum.A+diag(c(mXX.j[1]/to[i],mXX.j[2]/to[i])) tem.g1<-tem.g1+to[i]^2 }} g1m.j<-(sumt-to[j])-tem.g1/(sumt-to[j]) varx.p.j[,,j]<-((1/((m-1)-1))*sum.X-mXX.j)*(m-1-1)/g1m.j #(1/(m-1))*sum.A if(varx.p.j[1,1,j]<=0 | varx.p.j[2,2,j]<=0){ varx.p.check<-varx.p.j[,,j] varx.p.check[1,1]<-ifelse(varx.p.j[1,1,j]<=0,mXX.j[1]/(m-1)^{2/3},varx.p.j[1,1,j]) varx.p.check[2,2]<-ifelse(varx.p.j[2,2,j]<=0,mXX.j[2]/(m-1)^{2/3},varx.p.j[2,2,j]) sigma.check.av<-(1/2)*(varx.p.check[1,1]+varx.p.check[2,2]) if(eigen(varx.p.check)$values[1]<0){ aa<-(2*sigma.check.av)/((eigen(varx.p.check)$values[2])+sigma.check.av/(sqrt(m-1))) eig1<-aa*sigma.check.av/(sqrt(m-1)) eig2<-aa*(eigen(varx.p.check)$values[2]) } if(eigen(varx.p.check)$values[2]<0){ ab<-(2*sigma.check.av)/((eigen(varx.p.check)$values[1])+sigma.check.av/(sqrt(m-1))) eig2<-ab*sigma.check.av/(sqrt(m-1)) eig1<-ab*(eigen(varx.p.check)$values[1]) } varx.p.j[,,j]<-eig1*(eigen(varx.p.check)$vectors[1,])%*%t(eigen(varx.p.check)$vectors[1,])+ eig2*(eigen(varx.p.check)$vectors[2,])%*%t(eigen(varx.p.check)$vectors[2,]) } #################################### tem.w<-matrix(0,p1,p1) tem.y<-rep(0,p1) for(i in 1:m){if(i !=j){ a<-s[i] for(k in 1:no[i]){ tem.w<-tem.w+(w.ij.ran[,a+k-1]-w.bar.j[,j])%*%t(w.ij.ran[,a+k-1]) tem.y<-tem.y+(w.ij.ran[,a+k-1]-w.bar.j[,j])*y.ran[a+k-1] } }} beta1.p.j[,j]<-solve(tem.w)%*%tem.y ################################# sum.b2.j<-rep(0,p2) # tem.g1<-0 sum.b22.j<-0 tem.g2<-0 temm<-0 for(i in 1:m){if(i !=j){ a<-s[i] for(k in 1:no[i]){ temm<-temm+(to[i]*(1-to[i]/(sumt-to[j]))*t(w.i.bar.ran[,i])%*%solve(tem.w)%*%(w.ij.ran[,a+k-1]-w.bar.j[,j])) } sum.b2.j[1]<-sum.b2.j[1]+to[i]*ybar.ran[i]*(Xbar.ran[1,i]-XBARR.j[1,j]) sum.b2.j[2]<-sum.b2.j[2]+to[i]*ybar.ran[i]*(Xbar.ran[2,i]-XBARR.j[2,j]) # tem.g1<-tem.g1+to[i]^2 tem.g2<-tem.g2+no[i]^2 sum.b22.j<-sum.b22.j+no[i]*(ybar.ran[i]-ybarr.j[j])^2 }} # g1m.j<-(sumt-to[j])-tem.g1/(sumt-to[j]) g2m.j<-(n-no[j])-tem.g2/(n-no[j]) beta2.p.j[,j]<-solve(sum.X/(m-1-1)-mXX.j)%*%(sum.b2.j/(m-1-1)) #g1m.j-temm)) beta0.p.j[j]<-ybarr.j[j]-t(beta1.p.j[,j])%*%w.bar.j[,j]-t(beta2.p.j[,j])%*%XBARR.j[,j] ################################# tem.w.bar<-0 tem.w1<-0 for(i in 1:m){if(i !=j){ a<-s[i] for(k in 1:no[i]){ tem.w.bar<-tem.w.bar+(t(beta1.p.j[,j])%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]))^2 tem.w1<-tem.w1+no[i]*(t(beta1.p.j[,j])%*%(w.i.bar.ran[,i]-w.bar.j[,j]))^2 } }} varv.p.j[j]<-(((sum.b22.j/((m-1)-1)) - myy.j- t(beta1.p.j[,j])%*% (msb.w.j-mWW.j) %*%beta1.p.j[,j] )*((m-1)-1)/g2m.j) - t(beta2.p.j[,j])%*%(sum.b2.j/(g1m.j-0)) #- (1/g2m.j)*(tem.w1-(((m-1)-1)/((n-no[j])-(m-1)))*tem.w.bar) if(varv.p.j[j]<=0)varv.p.j[j]<-0.0001 #estimate of var(u) # varv.p.j[j]<-1.49 #taken from the functional ME ################################## tem.yw<-0 tem.Aw1<-matrix(0,p1,p1) tem.Bw1<-matrix(0,p1,p1) A.w<-0 B.w<-0 for(i in 1:m){if(i !=j){ a<-s[i] for(k in 1:no[i]){ tem.yw<-tem.yw+((y.ran[a+k-1]-ybar.ran[i])-t(beta1.p.j[,j])%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]))^2 tem.Aw1<-tem.Aw1+(w.ij.ran[,a+k-1]-w.bar.j[,j])%*%t(w.ij.ran[,a+k-1]) } tem.Bw1<-tem.Bw1+(no[i]^2)*(w.i.bar.ran[,i]-w.bar.j[,j])%*%t(w.i.bar.ran[,i]-w.bar.j[,j]) }} tem.Aw2<-0 for(i in 1:m){if(i !=j){ a<-s[i] for(k in 1:no[i]){ tem.Aw2<-tem.Aw2+t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]) B.w<-B.w+t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%tem.Bw1%*%solve(tem.Aw1)%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]) } }} A.w<-((n-no[j])-(m-1))-tem.Aw2 vare.p.j[j]<-tem.yw/((n-no[j])-(m-1)-p1) #(1/A.w)*(tem.yw-(varv.p.j[j]+t(beta2.p.j[,j]) %*%(sum.b2.j/(g1m.j-temm)) )*B.w) ###### }# loop of j ####################### g.hat.p.j<-rep(0,m) g1.hat.p.j<-rep(0,m) g.hat.wp.j<-rep(0,m) g1.hat.wp.j<-rep(0,m) c1<-matrix(0,m,m) mse.p.j<-rep(0,m) mse1.p.j<-rep(0,m) mse2.p.j<-rep(0,m) mse.wp.j<-rep(0,m) mse1.wp.j<-rep(0,m) mse2.wp.j<-rep(0,m) for(z in 1:m){ fi<-1 #((NO[z]-no[z])/NO[z]) for(i in 1:m){ #print("H1") Ai.p.j<-diag(c(mXX.j[1]/to[z],mXX.j[2]/to[z])) varx.p.j[,,i]<-diag(c(varx.p.j[1,1,i], varx.p.j[2,2,i])) Hi.p.j<-diag(p2)+Ai.p.j%*%solve(varx.p.j[,,i]) Di.p.j<-vare.p.j[i]/(vare.p.j[i]+no[z]*(varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Hi.p.j)%*%Ai.p.j%*%beta2.p.j[,i])) Di1.p.j<-vare.p.j[i]/no[z] Mi.p.j<-solve(Ai.p.j)+solve(varx.p.j[,,i]) g1.hat.p.j[z]<-g1.hat.p.j[z]+(Di1.p.j*(varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]))/(Di1.p.j+varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]) - g1.p[z] g1.hat.wp.j[z]<-g1.hat.wp.j[z]+ weight[i]*( (Di1.p.j*(varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]))/(Di1.p.j+varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]) - g1.p[z] ) muhat.p.j<-(1-fi*Di.p.j)*(ybar.ran[z]+t(beta1.p.j[,i])%*%(W.bar[,z]-w.i.bar.ran[,z]))+fi*Di.p.j*(beta0.p.j[i]+t(beta1.p.j[,i])%*%W.bar[,z]+t(beta2.p.j[,i])%*%XBARR.j[,i])+fi*Di.p.j*(t(beta2.p.j[,i])%*%varx.p.j[,,i]%*%solve(varx.p.j[,,i]+Ai.p.j)%*%(Xbar.ran[,z]-XBARR.j[,i]) ) g.hat.p.j[z]<-g.hat.p.j[z]+ (muhat.p.j -muhat.p[z])^2 g.hat.wp.j[z]<-g.hat.wp.j[z]+weight[i]*((muhat.p.j -muhat.p[z] )^2) c1[i,z]<- muhat.p.j } #print("H2") } for (z in 1:m){ #print("H3") mse.p.j[z]<-((m-1)/m)*g.hat.p.j[z]+g1.p[z]-((m-1)/m)*g1.hat.p.j[z] mse1.p.j[z]<-g1.p[z]-((m-1)/m)*g1.hat.p.j[z] mse2.p.j[z]<-((m-1)/m)*g.hat.p.j[z] mse.wp.j[z]<-g.hat.wp.j[z]+g1.p[z]-g1.hat.wp.j[z] mse1.wp.j[z]<-g1.p[z]-g1.hat.wp.j[z] mse2.wp.j[z]<-g.hat.wp.j[z] } #print("H4") ############################################################################## muhat.p mse.p.j mse1.p.j mse2.p.j mse.wp.j mse1.wp.j mse2.wp.j area<-c(1:m) result<-cbind(area,no,to,muhat.p, g1.p, mse.p.j,mse1.p.j,mse2.p.j, mse.wp.j, mse1.wp.j, mse2.wp.j) result #########################To produce Figure 1####################################### BMI<-muhat.p quartz() boxplot(BMI,cex.lab=1.3, cex.axis=1.5) box(lwd=3) ##########################Figure 2################################################# Weighted<-mse.wp.j Unweighted<-mse.p.j Res<-c(Unweighted, Weighted) region<-c(1:m) boxplot(Res,cex.lab=1.3, cex.axis=1.5, ylab="MSPE estimation") box(lwd=3) ####################################################################################