## Calculate bias, cvg, type1 error, and power from simulation results ## For student t analysis proc.res<-function(resname){ res<-read.table(resname,header=T,row.names=NULL) N<-dim(res)[1] par.name<-c('mL','vL','mS','vS','vLS','vE','k') par<-c(6,3,2,1,0,1,3) for (i in 0:6){ ind<-c(i+2, i+9, i+16, i+23, 2*i+30, 2*i+31) temp<-res[,ind] avg<-mean(temp[,1]) if (par[i+1]==0){ bias<-(avg-par[i+1])*100 }else{ bias<-(avg-par[i+1])/par[i+1]*100 } esd<-sd(temp[,1]) msd<-mean(temp[,2]) cvg1<-sum( temp[,3]par[i+1] )/N cvg1.l<-sum(temp[,3]>par[i+1])/N cvg1.u<-sum(temp[,4]par[i+1] )/N cvg2.l<-sum(temp[,5]>par[i+1])/N cvg2.u<-sum(temp[,6]0 | temp[,4]<0 )/N power1.l<-sum(temp[,3]>0)/N power1.u<-sum(temp[,4]<0)/N power2<-sum( temp[,5]>0 | temp[,6]<0 )/N power2.l<-sum(temp[,5]>0)/N power2.u<-sum(temp[,6]<0)/N cat(par.name[i+1], c(avg, bias, msd, esd, cvg1, cvg1.l, cvg1.u, power1, power1.l, power1.u, cvg2, cvg2.l, cvg2.u, power2, power2.l, power2.u, N), '\n') } } proc.res('Ndata-Tdata-100-0.1.txt') ## for normal data analysis proc.res<-function(resname){ res<-read.table(resname,header=T,row.names=NULL) N<-dim(res)[1] par.name<-c('mL','vL','mS','vS','vLS','vE') par<-c(6,3, 2,1,0,1) for (i in 0:5){ ind<-c(i+1, i+7, i+13, i+19, 2*i+25, 2*i+26) temp<-res[,ind] avg<-mean(temp[,1]) if (par[i+1]==0){ bias<-(avg-par[i+1])*100 }else{ bias<-(avg-par[i+1])/par[i+1]*100 } esd<-sd(temp[,1]) msd<-mean(temp[,2]) cvg1<-sum( temp[,3]par[i+1] )/N cvg1.l<-sum(temp[,3]>par[i+1])/N cvg1.u<-sum(temp[,4]par[i+1] )/N cvg2.l<-sum(temp[,5]>par[i+1])/N cvg2.u<-sum(temp[,6]0 | temp[,4]<0 )/N power1.l<-sum(temp[,3]>0)/N power1.u<-sum(temp[,4]<0)/N power2<-sum( temp[,5]>0 | temp[,6]<0 )/N power2.l<-sum(temp[,5]>0)/N power2.u<-sum(temp[,6]<0)/N cat(par.name[i+1], c(avg, bias, msd, esd, cvg1, cvg1.l, cvg1.u, power1, power1.l, power1.u, cvg2, cvg2.l, cvg2.u, power2, power2.l, power2.u, N), '\n') } } ## for the latent basis model results process ## Calculate bias, cvg, type1 error, and power from simulation results ## For student t analysis proc.res<-function(resname){ res<-read.table(resname,header=T,row.names=NULL) res<-res[res[,62]<2,] N<-dim(res)[1] par.name<-c('mL','vL','mS','vS','vLS','vE','k','a[2]','a[3]','a[4]') par<-c(6,1,2,.5,0,1.5, 3, .25, .5, .75) for (i in 0:9){ ind<-c(i+2, i+12, i+22, i+32, 2*i+42, 2*i+43) temp<-res[,ind] avg<-mean(temp[,1]) if (par[i+1]==0){ bias<-(avg-par[i+1])*100 }else{ bias<-(avg-par[i+1])/par[i+1]*100 } esd<-sd(temp[,1]) msd<-mean(temp[,2]) cvg1<-sum( temp[,3]par[i+1] )/N cvg1.l<-sum(temp[,3]>par[i+1])/N cvg1.u<-sum(temp[,4]par[i+1] )/N cvg2.l<-sum(temp[,5]>par[i+1])/N cvg2.u<-sum(temp[,6]0 | temp[,4]<0 )/N power1.l<-sum(temp[,3]>0)/N power1.u<-sum(temp[,4]<0)/N power2<-sum( temp[,5]>0 | temp[,6]<0 )/N power2.l<-sum(temp[,5]>0)/N power2.u<-sum(temp[,6]<0)/N cat(par.name[i+1], c(avg, bias, msd, esd, cvg1, cvg1.l, cvg1.u, power1, power1.l, power1.u, cvg2, cvg2.l, cvg2.u, power2, power2.l, power2.u, N), '\n') } } ==== Normal model ==== res<-read.table('res-Tdata-N-N-500-Uniform.txt',row.names=NULL) head(res) N<-dim(res)[1] par<-c(6,1,2,.5,0,1.5) for (i in 0:5){ ind<-c(i+2, i+8, i+14, i+20, 2*i+26, 2*i+27) temp<-res[,ind] avg<-mean(temp[,1]) if (par[i+1]==0) { bias<-(avg-par[i+1])*100 }else{ bias<-(avg-par[i+1])/par[i+1]*100 } msd<-mean(temp[,2]) cvg1<-sum( temp[,3]par[i+1] )/N cvg2<-sum( temp[,5]par[i+1] )/N power1<-sum( temp[,3]>0 | temp[,4]<0 )/N power2<-sum( temp[,5]>0 | temp[,6]<0 )/N cat(c(i+1, avg, bias, msd, cvg1, power1, cvg2, power2), '\n') } ==== Latent basis models ==== proc.res<-function(resname){ res<-read.table(resname, row.names=NULL) res<-res[res[,62]<2,] N<-dim(res)[1] par.name<-c('a[2]','a[3]','a[4]','mL','mS','vL','vS','vLS','vE','k') par<-c(6,1,2,.5,0,1.5, 3, .25, .5, .75) allres<-array(NA, dim=c(10, 10)) for (i in 0:9){ ind<-c(i+2, i+12, i+22, i+32, 2*i+42, 2*i+43) temp<-res[,ind] avg<-mean(temp[,1]) if (par[i+1]==0) { bias<-(avg-par[i+1])*100 }else{ bias<-(avg-par[i+1])/par[i+1]*100 } esd<-sd(temp[,1]) msd<-mean(temp[,2]) cvg0<-sum( (temp[,1]-temp[,2]*1.96) par[i+1] )/N cvg1<-sum( temp[,3]par[i+1] )/N cvg2<-sum( temp[,5]par[i+1] )/N power0<-sum( (temp[,1]-temp[,2]*1.96) >0 | (temp[,1]+temp[,2]*1.96) <0 )/N power1<-sum( temp[,3]>0 | temp[,4]<0 )/N power2<-sum( temp[,5]>0 | temp[,6]<0 )/N allres[i+1, ]<-c(avg, bias, esd, msd, cvg0, power0, cvg1, power1, cvg2, power2) } allres<-allres[c(8:10,1:7),] rownames(allres)<-par.name round(allres,3) } proc.res('res-Tdata-T-N-100-Wishart-LBM.txt') proc.res('res-Tdata-T-N-200-Wishart-LBM.txt') proc.res('res-Tdata-T-N-300-Wishart-LBM.txt') proc.res('res-Tdata-T-N-400-Wishart-LBM.txt') proc.res('res-Tdata-T-N-500-Wishart-LBM.txt') proc.res('res-Ndata-T-N-100-Wishart-LBM.txt') proc.res('res-Ndata-T-N-200-Wishart-LBM.txt') proc.res('res-Ndata-T-N-300-Wishart-LBM.txt') proc.res('res-Ndata-T-N-400-Wishart-LBM.txt') proc.res('res-Ndata-T-N-500-Wishart-LBM.txt') ## for the normal latent basis model proc.res<-function(resname){ res<-read.table(resname, row.names=NULL) res<-res[res[,56]<2,] N<-dim(res)[1] par.name<-c('a[2]','a[3]','a[4]','mL','mS','vL','vS','vLS','vE') par<-c(6,1,2,.5,0,1.5,.25, .5, .75) allres<-array(NA, dim=c(9, 10)) for (i in 0:8){ ind<-c(i+2, i+11, i+20, i+29, 2*i+38, 2*i+39) temp<-res[,ind] avg<-mean(temp[,1]) if (par[i+1]==0) { bias<-(avg-par[i+1])*100 }else{ bias<-(avg-par[i+1])/par[i+1]*100 } esd<-sd(temp[,1]) msd<-mean(temp[,2]) cvg0<-sum( (temp[,1]-temp[,2]*1.96) par[i+1] )/N cvg1<-sum( temp[,3]par[i+1] )/N cvg2<-sum( temp[,5]par[i+1] )/N power0<-sum( (temp[,1]-temp[,2]*1.96) >0 | (temp[,1]+temp[,2]*1.96) <0 )/N power1<-sum( temp[,3]>0 | temp[,4]<0 )/N power2<-sum( temp[,5]>0 | temp[,6]<0 )/N allres[i+1, ]<-c(avg, bias, esd, msd, cvg0, power0, cvg1, power1, cvg2, power2) } allres<-allres[c(7:9,1:6),] rownames(allres)<-par.name round(allres,3) } ==== Mplus ==== ## Simulated t data y<-read.table('../data-generated-student-coda-LBM.txt') ymat<-array(y[, 2], dim=c(1000,5, 500)) # t(ymat[2,,])->temp plot(1:5, temp[1,1:5], ylim=c(30, 150), type='l') for (i in 2:200){ lines(1:5, temp[i,1:5]) } ## write data into txt files for (i in 1001:2000){ t(ymat[i-1000,,1:100])->temp filename<-paste('c:/zzy/student-t/temp/TdataRN500I',i,'.txt', sep='') write.table(temp, filename, row.names=F, col.names=F) } source("http://nd.psychstat.org/_media/classdata/rtobugs.r") for (i in 1001:2000){ t(ymat[i-1000,,])->temp filename<-paste('../Tdata/TdataN500I',i,'bugs.txt', sep='') rtobugs(list(y=temp), filename) } ## Simulated t data y<-read.table('N:\\Private\\zzy\\research\\Bayesian analysis of longitudinal data\\Student t distribution\\Simulation\\data-generated-student-coda-300.txt') ymat<-array(y[, 2], dim=c(1000,5, 300)) ## write data into txt files for (i in 1001:2000){ t(ymat[i-1000,,])->temp filename<-paste('TdataN500I',i,'.txt', sep='') write.table(temp, filename, row.names=F, col.names=F) } source("http://nd.psychstat.org/_media/classdata/rtobugs.r") for (i in 1001:2000){ t(ymat[i-1000,,])->temp filename<-paste('TdataN500I',i,'bugs.txt', sep='') rtobugs(list(y=temp), filename) } for (i in 1001:2000){ #filename<-paste('TdataN500I',i,'.txt', sep='') filename<-paste('N:\\Private\\zzy\\research\\Bayesian analysis of longitudinal data\\Student t distribution\\Simulation\\Tdata\\TdataN500I',i,'.txt', sep='') data<-read.table(filename) temp<-data[1:500,] outname<-paste('c:/zzy/student-t/temp/TdataN500I',i,'.txt', sep='') write.table(temp, outname, row.names=F, col.names=F) } ## process results from Mplus analysis res<-scan('c:\\mplus\\all_results.txt') res<-scan('c:\\mplus\\T-N-400-LBM.txt') results<-matrix(res, byrow=T, ncol=34) results<-results[results[,1]<10, ] results<-results[results[,1]>-2, ] results<-results[results[,2]<6, ] results<-results[results[,3]<6, ] results<-results[results[,4]<6, ] results<-results[results[,7]<15, ] results<-results[results[,8]>-6, ] results<-results[results[,9]<20, ] n<-dim(results)[1] t.est<-c(.25, .5, .75, 1.5,6,1,2,0,.5) avg<-apply(results[,1:9],2,mean) bias<-avg-t.est sd.est<-apply(results[,1:9],2,sd) msd<-apply(results[,10:18],2,mean) cvg<-rep(0,9) power<-rep(0,9) for (i in 1:9){ cvg[i]<-sum( t.est[i]>results[,i]-1.96*results[,i+9] & t.est[i]0 | results[,i]+1.96*results[,i+9]<0 )/n } cbind(t.est, avg, bias,sd.est, msd, cvg, power)[c(1,2,3,5,6,7,9,8,4),] ## y<-read.table('N:/Private/zzy/research/Bayesian analysis of longitudinal data/Student t distribution/real-data-analysis/scale.txt', na.string='.') N<-dim(y)[1] y1<-y[sample(1:N, 338), ] boxplot(round(y,1)) write.table(round(y1,1), 'N:/Private/zzy/research/Bayesian analysis of longitudinal data/Student t distribution/real-data-analysis/scale.txt', col.names=F, row.names=F, na='.') y2<-as.matrix(y) bugsdata<-list(N=dim(y2)[1], y=round(y2,1)) rtobugs(bugsdata, 'N:/Private/zzy/research/Bayesian analysis of longitudinal data/Student t distribution/real-data-analysis/tempbugsdata.txt') y<-y[y[,1]<10,] y<-y[y[,2]<10,] y<-y[y[,3]<10,] y<-y[y[,4]<10,] y<-y[y[,5]<10,] y<-y[y[,1]>0,] y<-y[y[,2]>0,] y<-y[y[,3]>0,] y<-y[y[,4]>0,] y<-y[y[,5]>0,] y<-y[y[,2]>y[,1],] y<-y[y[,3]>y[,2],] y<-y[y[,4]>y[,3],] y<-y[y[,5]>y[,4],] boxplot(y[,1],y[,2],y[,3],y[,4],y[,5]) anscombe.test(y[,1]) y<-read.table('N:/Private/zzy/research/Bayesian analysis of longitudinal data/Student t distribution/real-data-analysis/TdataN500I1705.txt') y<-round(y, 1)