lab:projects:02robust_growth_curve_modeling_using_student-t_distribution:r_codes_for_final_results_process
Table of Contents
## 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] & temp[,4]>par[i+1] )/N
cvg1.l<-sum(temp[,3]>par[i+1])/N
cvg1.u<-sum(temp[,4]<par[i+1])/N
cvg2<-sum( temp[,5]<par[i+1] & temp[,6]>par[i+1] )/N
cvg2.l<-sum(temp[,5]>par[i+1])/N
cvg2.u<-sum(temp[,6]<par[i+1])/N
power1<-sum( temp[,3]>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] & temp[,4]>par[i+1] )/N
cvg1.l<-sum(temp[,3]>par[i+1])/N
cvg1.u<-sum(temp[,4]<par[i+1])/N
cvg2<-sum( temp[,5]<par[i+1] & temp[,6]>par[i+1] )/N
cvg2.l<-sum(temp[,5]>par[i+1])/N
cvg2.u<-sum(temp[,6]<par[i+1])/N
power1<-sum( temp[,3]>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] & temp[,4]>par[i+1] )/N
cvg1.l<-sum(temp[,3]>par[i+1])/N
cvg1.u<-sum(temp[,4]<par[i+1])/N
cvg2<-sum( temp[,5]<par[i+1] & temp[,6]>par[i+1] )/N
cvg2.l<-sum(temp[,5]>par[i+1])/N
cvg2.u<-sum(temp[,6]<par[i+1])/N
power1<-sum( temp[,3]>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] & temp[,4]>par[i+1] )/N
cvg2<-sum( temp[,5]<par[i+1] & temp[,6]>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] & (temp[,1]+temp[,2]*1.96) >par[i+1] )/N
cvg1<-sum( temp[,3]<par[i+1] & temp[,4]>par[i+1] )/N
cvg2<-sum( temp[,5]<par[i+1] & temp[,6]>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] & (temp[,1]+temp[,2]*1.96) >par[i+1] )/N
cvg1<-sum( temp[,3]<par[i+1] & temp[,4]>par[i+1] )/N
cvg2<-sum( temp[,5]<par[i+1] & temp[,6]>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]<results[,i]+1.96*results[,i+9] )/n
power[i]<-sum( results[,i]-1.96*results[,i+9]>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)