lab:zhang:power_analysis_with_timo_manuel
Table of Contents
Power analysis
possible journals:
Methodology: http://www.psycontent.com/content/120235/
FIML
Scripts and codes
Output parameter bias in a table
linear.res<-function(filename)
{
##process the results
# res.file<-paste('res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
## filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\linear\\power\\res-N100-T3-alpha0-R0-a.txt'
res.comp<-read.table(filename,col.names=c(1:55),fill=T)
## data in the file
## Model with freely estimated change score
## 1 id 2 var_E 3 mu_I 4 mu_S 5 psi_I 6 psi_IS 7 psi_S
## s.e. 8 9 10 11 12 13
## 14 chi-value 15 chi-df 16 chi-p 17 CFI 18 TLI 19 H0ll 20 H1ll
## 21 para 22 AIC 23 BIC 24 A-BIC 25 RMSEA 26 RMSEA.l 27 RMSEA.h 28 RMSEA.p 29 SRMR
## model with zero change score
## 30 var_E 31 mu_I 32 psi_I 33 psi_IS 34 psi_S
## 35 36 37 38 39
## 40 chi-value 41 chi-df 42 chi-p 43 CFI 44 TLI 45 H0ll 46 H1ll
## 47 para 48 AIC 49 BIC 50 A-BIC 51 RMSEA 52 RMSEA.l 53 RMSEA.h 54 RMSEA.p 55 SRMR
res.comp<-res.comp[!is.na(res.comp[,55]),]
chi.diff<-(res.comp[,40]-res.comp[,14])
res.comp<-res.comp[chi.diff>0,]
mu_S <- mean(res.comp[,4],na.rm=T)
mu_S
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
r.v<-c(.1,.3)
M.M<-c('MCAR', 'MAR', 'MNAR')
alpha.v<-c(1,0,0)
a<-c('-a','-a', '')
for (T in T.v){
for (N in N.v){
cat(T, ' ', N, ' ')
filename<-paste('C:\\zzy\\LongPower\\Results\\MCAR\\linear\\power\\res-N',N,'-T',T,'-alpha',alpha,'-R',0,'-a.txt',sep='')
cat(linear.res(filename))
cat(' ')
for (r in r.v){
for (i in 1:3){
filename<-paste('C:\\zzy\\LongPower\\Results\\',M.M[i],'\\linear\\power\\res-N',N,'-T',T,'-alpha',alpha.v[i],'-R',r, a[i] ,'.txt',sep='')
cat(linear.res(filename))
cat(' ')
}
}
cat('\n')
}
}
quad.res<-function(filename)
{
##process the results
#filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\quadratic\\power\\res-N100-T4-alpha0-R0-a.txt'
res.comp<-read.table(filename,col.names=c(1:83),fill=T)
## data in the file
## Model with freely estimated change score
## 1 id 2 var_E 3 mu_I 4 mu_S 5 mu_Q 6 Beta 7 psi_I 8 psi_IS 9 psi_S 10 PSI_IQ 11 psi_QS 12 psi_Q 13 psi_XI 14 psi_XQ
## s.e. 15 16 17 18 19 20 21 22 23 24 25 26 27
## 28 chi-value 29 chi-df 30 chi-p 31 CFI 32 TLI 33 H0ll 34 H1ll
## 35 para 36 AIC 37 BIC 38 A-BIC 39 RMSEA 40 RMSEA.l 41 RMSEA.h 42 RMSEA.p 43 SRMR
## model with zero change score
## 44 var_E 45 mu_I 46 mu_S 47 mu_Q 48 Beta 49 psi_I 50 psi_IS 51 psi_S 52 PSI_IQ 53 psi_QS 54 psi_Q 55 psi_XI 56 psi_XQ
## s.e. 57 58 59 60 61 62 63 64 65 66 67
## 68 chi-value 69 chi-df 70 chi-p 71 CFI 72 TLI 73 H0ll 74 H1ll
## 75 para 76 AIC 77 BIC 78 A-BIC 79 RMSEA 80 RMSEA.l 81 RMSEA.h 82 RMSEA.p 83 SRMR
res.comp<-res.comp[!is.na(res.comp[,83]),]
chi.diff1<-(res.comp[,68]-res.comp[,28])
res.comp<-res.comp[chi.diff1>0,]
Beta <- mean(res.comp[,6])
Beta
}
N.v<-c(100,200,300,400,500)
T.v<-4:5
r.v<-c(.1,.3)
M.M<-c('MCAR', 'MAR', 'MNAR')
alpha.v<-c(1,0,0)
a<-c('-a','-a', '-a')
for (T in T.v){
for (N in N.v){
cat(T, ' ', N, ' ')
filename<-paste('C:\\zzy\\LongPower\\Results\\MCAR\\quadratic\\power\\res-N',N,'-T',T,'-alpha',alpha,'-R',0,'-a.txt',sep='')
cat(quad.res(filename))
cat(' ')
for (r in r.v){
for (i in 1:3){
filename<-paste('C:\\zzy\\LongPower\\Results\\',M.M[i],'\\quadratic\\power\\res-N',N,'-T',T,'-alpha',alpha.v[i],'-R',r, a[i] ,'.txt',sep='')
cat(quad.res(filename))
cat(' ')
}
}
cat('\n')
}
}
R codes for process the results for estimate bias
## Processing results for the power analysis
## Function for the linear model
linear.res<-function(filename)
{
##process the results
# res.file<-paste('res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
## filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\linear\\power\\res-N100-T3-alpha0-R0-a.txt'
res.comp<-read.table(filename,col.names=c(1:55),fill=T)
## data in the file
## Model with freely estimated change score
## 1 id 2 var_E 3 mu_I 4 mu_S 5 psi_I 6 psi_IS 7 psi_S
## s.e. 8 9 10 11 12 13
## 14 chi-value 15 chi-df 16 chi-p 17 CFI 18 TLI 19 H0ll 20 H1ll
## 21 para 22 AIC 23 BIC 24 A-BIC 25 RMSEA 26 RMSEA.l 27 RMSEA.h 28 RMSEA.p 29 SRMR
## model with zero change score
## 30 var_E 31 mu_I 32 psi_I 33 psi_IS 34 psi_S
## 35 36 37 38 39
## 40 chi-value 41 chi-df 42 chi-p 43 CFI 44 TLI 45 H0ll 46 H1ll
## 47 para 48 AIC 49 BIC 50 A-BIC 51 RMSEA 52 RMSEA.l 53 RMSEA.h 54 RMSEA.p 55 SRMR
res.comp<-res.comp[!is.na(res.comp[,55]),]
chi.diff<-(res.comp[,40]-res.comp[,14])
res.comp<-res.comp[chi.diff>0,]
chi.diff1<-chi.diff[chi.diff>0]
nn<-length(chi.diff1)
pow1<-sum(chi.diff1>3.841459)/nn
var_E <- mean(res.comp[,2],na.rm=T)
mu_I <- mean(res.comp[,3],na.rm=T)
mu_S <- mean(res.comp[,4],na.rm=T)
psi_I <- mean(res.comp[,5],na.rm=T)
psi_IS <- mean(res.comp[,6],na.rm=T)
psi_S <- mean(res.comp[,7],na.rm=T)
CI<-cbind(res.comp[,4] -1.96*res.comp[,10],res.comp[,4]+1.96*res.comp[,10])
cvg<-sum(CI[,1]<1 & CI[,2]>1)/nn
cng<-nn/10000
c(var_E, mu_I, mu_S, psi_I, psi_IS, psi_S, pow1, cvg, cng)
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
alpha.v<-c(0)
r.v<-c(0, .1,.3)
for (T in T.v){
for (N in N.v){
for (r in r.v){
for (alpha in alpha.v){
## Power
filename<-paste('C:\\zzy\\research\\LongPower\\Results\\MNAR\\linear\\power\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'.txt',sep='')
## Type1
#filename<-paste('C:\\zzy\\research\\LongPower\\Results\\MNAR\\linear\\type1\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'.txt',sep='')
res<-linear.res(filename)
cat(c(T, N, r, alpha, res))
cat("\n")
}
}
}
}
##################################################
## 2. parameter bias quadratic model
quad.res<-function(filename)
{
##process the results
#filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\quadratic\\power\\res-N100-T4-alpha0-R0-a.txt'
res.comp<-read.table(filename,col.names=c(1:83),fill=T)
## data in the file
## Model with freely estimated change score
## 1 id 2 var_E 3 mu_I 4 mu_S 5 mu_Q 6 Beta 7 psi_I 8 psi_IS 9 psi_S 10 PSI_IQ 11 psi_QS 12 psi_Q 13 psi_XI 14 psi_XQ
## s.e. 15 16 17 18 19 20 21 22 23 24 25 26 27
## 28 chi-value 29 chi-df 30 chi-p 31 CFI 32 TLI 33 H0ll 34 H1ll
## 35 para 36 AIC 37 BIC 38 A-BIC 39 RMSEA 40 RMSEA.l 41 RMSEA.h 42 RMSEA.p 43 SRMR
## model with zero change score
## 44 var_E 45 mu_I 46 mu_S 47 mu_Q 48 Beta 49 psi_I 50 psi_IS 51 psi_S 52 PSI_IQ 53 psi_QS 54 psi_Q 55 psi_XI 56 psi_XQ
## s.e. 57 58 59 60 61 62 63 64 65 66 67
## 68 chi-value 69 chi-df 70 chi-p 71 CFI 72 TLI 73 H0ll 74 H1ll
## 75 para 76 AIC 77 BIC 78 A-BIC 79 RMSEA 80 RMSEA.l 81 RMSEA.h 82 RMSEA.p 83 SRMR
res.comp<-res.comp[!is.na(res.comp[,83]),]
chi.diff1<-(res.comp[,68]-res.comp[,28])
res.comp<-res.comp[chi.diff1>0,]
chi.diff1<-chi.diff1[chi.diff1>0]
nn<-length(chi.diff1)
pow1<-sum(chi.diff1>3.84)/nn
var_E <- mean(res.comp[,2])
mu_I <- mean(res.comp[,3])
mu_S <- mean(res.comp[,4])
mu_Q <- mean(res.comp[,5])
Beta <- mean(res.comp[,6])
psi_I <- mean(res.comp[,7])
psi_IS <- mean(res.comp[,8])
psi_S <- mean(res.comp[,9])
psi_IQ <- mean(res.comp[,10])
psi_QS <- mean(res.comp[,11])
psi_Q <- mean(res.comp[,12])
#psi_XI <- mean(res.comp[,13])
#psi_XS <- mean(res.comp[,14])
CI<-cbind(res.comp[,6] -1.96*res.comp[,19],res.comp[,6]+1.96*res.comp[,19])
cvg<-sum(CI[,1]< -1 & CI[,2]> -1)/nn
cng<-nn/10000
c(var_E, mu_I, mu_S, mu_Q, Beta, psi_I, psi_IS, psi_S, psi_IQ, psi_QS, psi_Q, pow1, cvg, cng)
#c(mu_I, mu_S, mu_Q, Beta, nn)
}
M<-c('power', 'type1')
N.v<-c(100,200,300,400,500)
T.v<-4:5
alpha.v<-c(0)
r.v<-c(0, .1,.3)
for (T in T.v){
for (N in N.v){
for (r in r.v){
for (alpha in alpha.v){
filename<-paste('C:\\zzy\\research\\LongPower\\Results\\MNAR\\quadratic\\type1\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
# filename<-paste('C:\\zzy\\research\\LongPower\\Results\\MCAR\\quadratic\\power\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
res<-quad.res(filename)
cat(c(T, N, r, alpha, res))
cat("\n")
}
}
}
}
R codes to extract the chi-square differences
## Processing results for the power analysis
## Function for the linear model
linear.res<-function(filename)
{
##process the results
# res.file<-paste('res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
## filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\linear\\power\\res-N100-T3-alpha0-R0-a.txt'
res.comp<-read.table(filename,col.names=c(1:55),fill=T)
## data in the file
## Model with freely estimated change score
## 1 id 2 var_E 3 mu_I 4 mu_S 5 psi_I 6 psi_IS 7 psi_S
## s.e. 8 9 10 11 12 13
## 14 chi-value 15 chi-df 16 chi-p 17 CFI 18 TLI 19 H0ll 20 H1ll
## 21 para 22 AIC 23 BIC 24 A-BIC 25 RMSEA 26 RMSEA.l 27 RMSEA.h 28 RMSEA.p 29 SRMR
## model with zero change score
## 30 var_E 31 mu_I 32 psi_I 33 psi_IS 34 psi_S
## 35 36 37 38 39
## 40 chi-value 41 chi-df 42 chi-p 43 CFI 44 TLI 45 H0ll 46 H1ll
## 47 para 48 AIC 49 BIC 50 A-BIC 51 RMSEA 52 RMSEA.l 53 RMSEA.h 54 RMSEA.p 55 SRMR
res.comp<-res.comp[!is.na(res.comp[,55]),]
chi.diff<-(res.comp[,40]-res.comp[,14])
res.comp<-res.comp[chi.diff>0,]
chi.diff1<-chi.diff[chi.diff>0]
chi.diff1
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
alpha.v<-c(1)
r.v<-c(0, .1,.3)
MM<-'MCAR'
for (T in T.v){
for (N in N.v){
for (r in r.v){
for (alpha in alpha.v){
## Power
#filename<-paste('C:\\zzy\\research\\LongPower\\Results\\MNAR\\linear\\power\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'.txt',sep='')
## Type1
filename<-paste('C:\\zzy\\research\\LongPower\\Results\\',MM,'\\linear\\type1\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'.txt',sep='')
res<-linear.res(filename)
fileout<-paste('C:\\zzy\\research\\LongPower\\Results\\qqplot\\linear-',MM,'T-',T,'-N',N,'-R',r,'.txt',sep='')
write.table(res,fileout, row.names=F, col.names=F)
}
}
}
}
##################################################
## 2. parameter bias quadratic model
quad.res<-function(filename)
{
##process the results
#filename<-'C:\\zzy\\research\\LongPower\\Results\\MAR\\quadratic\\power\\res-N100-T4-alpha0-R0-a.txt'
res.comp<-read.table(filename,col.names=c(1:83),fill=T)
## data in the file
## Model with freely estimated change score
## 1 id 2 var_E 3 mu_I 4 mu_S 5 mu_Q 6 Beta 7 psi_I 8 psi_IS 9 psi_S 10 PSI_IQ 11 psi_QS 12 psi_Q 13 psi_XI 14 psi_XQ
## s.e. 15 16 17 18 19 20 21 22 23 24 25 26 27
## 28 chi-value 29 chi-df 30 chi-p 31 CFI 32 TLI 33 H0ll 34 H1ll
## 35 para 36 AIC 37 BIC 38 A-BIC 39 RMSEA 40 RMSEA.l 41 RMSEA.h 42 RMSEA.p 43 SRMR
## model with zero change score
## 44 var_E 45 mu_I 46 mu_S 47 mu_Q 48 Beta 49 psi_I 50 psi_IS 51 psi_S 52 PSI_IQ 53 psi_QS 54 psi_Q 55 psi_XI 56 psi_XQ
## s.e. 57 58 59 60 61 62 63 64 65 66 67
## 68 chi-value 69 chi-df 70 chi-p 71 CFI 72 TLI 73 H0ll 74 H1ll
## 75 para 76 AIC 77 BIC 78 A-BIC 79 RMSEA 80 RMSEA.l 81 RMSEA.h 82 RMSEA.p 83 SRMR
res.comp<-res.comp[!is.na(res.comp[,83]),]
chi.diff1<-(res.comp[,68]-res.comp[,28])
res.comp<-res.comp[chi.diff1>0,]
chi.diff1<-chi.diff1[chi.diff1>0]
chi.diff1
}
N.v<-c(100,200,300,400,500)
T.v<-4:5
alpha.v<-c(0)
r.v<-c(0, .1,.3)
MM<-'MAR'
for (T in T.v){
for (N in N.v){
for (r in r.v){
for (alpha in alpha.v){
filename<-paste('C:\\zzy\\research\\LongPower\\Results\\',MM, '\\quadratic\\type1\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
# filename<-paste('C:\\zzy\\research\\LongPower\\Results\\MCAR\\quadratic\\power\\res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
res<-quad.res(filename)
fileout<-paste('C:\\zzy\\research\\LongPower\\Results\\qqplot\\quad-',MM,'T-',T,'-N',N,'-R',r,'.txt',sep='')
write.table(res,fileout, row.names=F, col.names=F) }
}
}
}
R codes for Q-Q plots
## qqplot for testing the distributions
## for linear model
qq.linear<-function(T, N, r){
prob<-seq(.01,.99,.01)
plotname<-paste('QQ-Linear-', T, '-N-', N, '-R-', r, '.pdf', sep='')
pdf(file=plotname, width=5, height=5)
## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
qqplot(qchisq(prob, 1), quantile(y.comp, prob=prob), type='l', lwd=3, ylab='Sample quantiles', xlab='Theoretical quantiles')
abline(0,1)
## MCAR data
file.mcar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mcar<-scan(file.mcar)
lines(qchisq(prob, 1), quantile(y.mcar, prob=prob), lty=2, lwd=3)
## MAR data
file.mar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mar<-scan(file.mar)
lines(qchisq(prob, 1), quantile(y.mar, prob=prob), lty=3, lwd=3)
## MNAR data
file.mnar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MNAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mnar<-scan(file.mnar)
lines(qchisq(prob, 1), quantile(y.mnar, prob=prob), lty=4, lwd=3)
legend('bottomright', c('Complete','MCAR', 'MAR', 'MNAR'), lty=c(1,2,3,4), lwd=c(3,3,3,3))
dev.off()
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
r.v<-c(.1,.3)
for (T in T.v){
for (N in N.v){
for (r in r.v){
qq.linear(T, N, r)
}
}
}
## Q-Q plot for quadratic model
qq.quad<-function(T, N, r){
prob<-seq(.01,.99,.01)
plotname<-paste('QQ-Quad-', T, '-N-', N, '-R-', r, '.pdf', sep='')
pdf(file=plotname, width=5, height=5)
## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
qqplot(qchisq(prob, 1), quantile(y.comp, prob=prob), type='l', lwd=3, ylab='Sample quantiles', xlab='Theoretical quantiles')
abline(0,1)
## MCAR data
file.mcar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MCAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mcar<-scan(file.mcar)
lines(qchisq(prob, 1), quantile(y.mcar, prob=prob), lty=2, lwd=3)
## MAR data
file.mar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mar<-scan(file.mar)
lines(qchisq(prob, 1), quantile(y.mar, prob=prob), lty=3, lwd=3)
## MNAR data
file.mnar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MNAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mnar<-scan(file.mnar)
lines(qchisq(prob, 1), quantile(y.mnar, prob=prob), lty=4, lwd=3)
legend('bottomright', c('Complete','MCAR', 'MAR', 'MNAR'), lty=c(1,2,3,4), lwd=c(3,3,3,3))
dev.off()
}
N.v<-c(100,200,300,400,500)
T.v<-4:5
r.v<-c(.1,.3)
for (T in T.v){
for (N in N.v){
for (r in r.v){
qq.quad(T, N, r)
}
}
}
R codes to calculate Type I error
## qqplot for testing the distributions
## for linear model
qq.linear<-function(MM, T, N, r){
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-', MM,'T-',T,'-N',N,'-R',r,'.txt',sep='')
y.comp<-scan(file.comp, quiet=T)
100*sum(y.comp>3.84)/length(y.comp)
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
r.v<-c(.1,.3)
M.M<-c('MCAR', 'MAR', 'MNAR')
for (T in T.v){
for (N in N.v){
cat(T, ' ', N, ' ')
cat(qq.linear('MCAR',T, N, 0))
cat(' ')
for (r in r.v){
for (i in 1:3){
cat(qq.linear(M.M[i], T, N, r))
cat(' ')
}
}
cat('\n')
}
}
## for linear model
qq.quad<-function(MM, T, N, r){
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-', MM,'T-',T,'-N',N,'-R',r,'.txt',sep='')
y.comp<-scan(file.comp, quiet=T)
100*sum(y.comp>3.84)/length(y.comp)
}
N.v<-c(100,200,300,400,500)
T.v<-4:5
r.v<-c(.1,.3)
M.M<-c('MCAR', 'MAR', 'MNAR')
for (T in T.v){
for (N in N.v){
cat(T, ' ', N, ' ')
cat(qq.quad('MCAR',T, N, 0))
cat(' ')
for (r in r.v){
for (i in 1:3){
cat(qq.quad(M.M[i], T, N, r))
cat(' ')
}
}
cat('\n')
}
}
R codes to calculate power
## qqplot for testing the distributions
## for linear model
qq.linear<-function(MM, T, N, r){
file.comp<-paste('C:\\zzy\\LongPower\\Results\\power\\linear-', MM,'T-',T,'-N',N,'-R',r,'.txt',sep='')
y.comp<-scan(file.comp, quiet=T)
100*sum(y.comp>3.84)/length(y.comp)
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
r.v<-c(.1,.3)
M.M<-c('MCAR', 'MAR', 'MNAR')
for (T in T.v){
for (N in N.v){
cat(T, ' ', N, ' ')
cat(qq.linear('MCAR',T, N, 0))
cat(' ')
for (r in r.v){
for (i in 1:3){
cat(qq.linear(M.M[i], T, N, r))
cat(' ')
}
}
cat('\n')
}
}
## for linear model
qq.quad<-function(MM, T, N, r){
file.comp<-paste('C:\\zzy\\LongPower\\Results\\power\\quad-', MM,'T-',T,'-N',N,'-R',r,'.txt',sep='')
y.comp<-scan(file.comp, quiet=T)
100*sum(y.comp>3.84)/length(y.comp)
}
N.v<-c(100,200,300,400,500)
T.v<-4:5
r.v<-c(.1,.3)
M.M<-c('MCAR', 'MAR', 'MNAR')
for (T in T.v){
for (N in N.v){
cat(T, ' ', N, ' ')
cat(qq.quad('MCAR',T, N, 0))
cat(' ')
for (r in r.v){
for (i in 1:3){
cat(qq.quad(M.M[i], T, N, r))
cat(' ')
}
}
cat('\n')
}
}
Q-Q plots with double y-axis
## qqplot for testing the distributions
## for linear model
qq.linear<-function(T, N, r){
T=4
N=300
r=0.1
prob<-seq(.01,.99,.01)
plotname<-paste('QQ-Linear-', T, '-N-', N, '-R-', r, '.pdf', sep='')
pdf(file=plotname, width=5, height=5)
par(mar=c(5, 4, 4, 5) + 0.1)
## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
plot(qchisq(prob, 1), quantile(y.comp, prob=prob), ylim=c(0,7), xlim=c(0,7),type='l', lwd=2, ylab='Sample quantiles', xlab='Theoretical quantiles')
#abline(0,1)
## MCAR data
file.mcar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mcar<-scan(file.mcar)
lines(qchisq(prob, 1), quantile(y.mcar, prob=prob), lty=2, lwd=2)
## MAR data
file.mar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mar<-scan(file.mar)
lines(qchisq(prob, 1), quantile(y.mar, prob=prob), lty=3, lwd=2)
par(new=T, mar=c(5, 4, 4, 5) + 0.1)
plot(1,1,type='n',ylim=c(0,55), xlim=c(0,7), axes=F, ylab='', xlab='')
## MNAR data
file.mnar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MNAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mnar<-scan(file.mnar)
quantile(y.mnar, .99)
lines(qchisq(prob, 1), quantile(y.mnar, prob=prob), lty=4, lwd=2, ylim=c(0,160), xlim=c(0,7))
axis(4)
mtext('MNAR sample quantitles', 4, line=3)
legend('bottomright', c('Complete','MCAR', 'MAR', 'MNAR'), lty=c(1,2,3,4), lwd=c(3,3,3,3))
dev.off()
}
N.v<-c(100,200,300,400,500)
T.v<-3:5
r.v<-c(.1,.3)
for (T in T.v){
for (N in N.v){
for (r in r.v){
qq.linear(T, N, r)
}
}
}
## Q-Q plot for quadratic model
qq.quad<-function(T, N, r){
T=4
N=300
r=0.3
prob<-seq(.01,.99,.01)
plotname<-paste('QQ-Quad-', T, '-N-', N, '-R-', r, '.pdf', sep='')
pdf(file=plotname, width=5, height=5)
par(mar=c(5, 4, 4, 5) + 0.1)
## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
qqplot(qchisq(prob, 1), quantile(y.comp, prob=prob), ylim=c(0,7), xlim=c(0,7),type='l', lwd=2, ylab='Sample quantiles', xlab='Theoretical quantiles')
#abline(0,1)
## MCAR data
file.mcar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MCAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mcar<-scan(file.mcar)
lines(qchisq(prob, 1), quantile(y.mcar, prob=prob), lty=2, lwd=2)
## MAR data
file.mar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mar<-scan(file.mar)
lines(qchisq(prob, 1), quantile(y.mar, prob=prob), lty=3, lwd=2)
## MNAR data
file.mnar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MNAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mnar<-scan(file.mnar)
lines(qchisq(prob, 1), quantile(y.mnar, prob=prob), lty=4, lwd=2)
legend('bottomright', c('Complete','MCAR', 'MAR', 'MNAR'), lty=c(1,2,3,4), lwd=c(3,3,3,3))
dev.off()
}
N.v<-c(100,200,300,400,500)
T.v<-4:5
r.v<-c(.1,.3)
for (T in T.v){
for (N in N.v){
for (r in r.v){
qq.quad(T, N, r)
}
}
}