===== Power analysis =====
possible journals:

Methodology: http://www.psycontent.com/content/120235/

==== FIML ====
[[R-Mplus-MNAR-Simulation]]

[[FIML-MNAR-0.1]]

[[FIML-MNAR-0.2]]

[[FIML-MNAR-0.4]]
==== Scripts and codes ====
  * [[R codes for plot two-sided Q-Q plots]]
  * [[Generate complete data for linear model]]

==== Output parameter bias in a table ====
<code>

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')
  }
}

</code>

==== R codes for process the results for estimate bias ====
<code>## 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")
          }
       }
    }
 }
</code>
==== R codes to extract the chi-square differences ====
<code>
## 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)           }
       }
    }
 }


</code>
 ==== R codes for Q-Q plots ====
<code>
## 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)
       }
    }
 }
</code>

==== R codes to calculate Type I error ====
<code>
## 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')
  }
}
</code>
==== R codes to calculate power ====
<code>
## 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')
  }
}
</code>
==== Q-Q plots with double y-axis ====

<code>
## 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)
       }
    }
 }
</code>