Zhiyong Zhang's Psychometric Website
 
Home Control Panel Login Register Search Submit Logout About me Contact me
My Research
Google
Web
This Site
Home > My Research > R > Using R and Mplus for growth curve power analysis based on likelihood ratio test
Using R and Mplus for growth curve power analysis based on likelihood ratio test
2010-02-28    Zhang, Z., von Oertzen, T., & Volkle, V.       Read: 11369 times
Cite this page: Zhang, Z., von Oertzen, T., & Volkle, V. (2010). Using R and Mplus for growth curve power analysis based on likelihood ratio test. Retrieved December 16, 2018, from http://www.psychstat.org/us/article.php/87.htm.
Using R and Mplus for growth curve power analysis based on likelihood ratio test

A briefly annotated version can be downloaded and viewed at http://www.psychstat.org/us/showimg.php?iid=66.

## Test if the chisq difference test applys when missing data are available
## Generate data in R and analyze data in Mplus
## The chisq value is read in for density plot.

library(mvtnorm)

gen.model<-function(T){
## generate mplus scripts for data analysis
## Mplus input scripts
mplus.null<-'m1null.inp'
mplus.alt<-'m2alt.inp'

## for null model
cat('TITLE: Null model\n',file=mplus.null);
cat('DATA:\n',file=mplus.null, append=T);
cat('  FILE=data.txt;\n',file=mplus.null, append=T);
cat('VARIABLE: \n',file=mplus.null, append=T);
cat('    NAMES ARE y1-y',T, ' x;\n', sep='', file=mplus.null, append=T);
cat('    USEVARIABLES ARE y1-y',T, ' x;\n', sep='', file=mplus.null, append=T);
cat('    MISSING = ALL(999); \n',file=mplus.null, append=T);

cat('ANALYSIS: \n',file=mplus.null, append=T);
cat('    COVERAGE=.01; \n',file=mplus.null, append=T);

cat('MODEL: \n',file=mplus.null, append=T);
cat('    i s q| ',file=mplus.null, append=T);
for (i in 1:T){
  cat('y',i,'@',i-1,' ',sep='',file=mplus.null,append=T)
}
cat(';\n',file=mplus.null,append=T)

cat('    [i * 50, s *, q * ]; \n',file=mplus.null, append=T);
cat('    i * 100, s *, q * ; \n',file=mplus.null, append=T);
cat('    q on x; \n',file=mplus.null, append=T);
cat('    i with s * 0; \n',file=mplus.null, append=T);
cat('    i with q * 0; \n',file=mplus.null, append=T);
cat('    s with q * 0; \n',file=mplus.null, append=T);
cat('    y1-y',T, ' (10);\n', sep='', file=mplus.null, append=T);
cat('    [y1-y',T, '@0];\n', sep='', file=mplus.null, append=T);

cat('SAVEDATA:\n', sep='', file=mplus.null, append=T);
cat('    results=res.txt;\n', sep='', file=mplus.null, append=T);

## for alternative model
cat('TITLE: Null model\n',file=mplus.alt);
cat('DATA:\n',file=mplus.alt, append=T);
cat('  FILE=data.txt;\n',file=mplus.alt, append=T);
cat('VARIABLE: \n',file=mplus.alt, append=T);
cat('    NAMES ARE y1-y',T, ' x;\n', sep='', file=mplus.alt, append=T);
cat('    USEVARIABLES ARE y1-y',T, ' x;\n', sep='', file=mplus.alt, append=T);
cat('    MISSING = ALL(999); \n',file=mplus.alt, append=T);

cat('ANALYSIS: \n',file=mplus.alt, append=T);
cat('    COVERAGE=.01; \n',file=mplus.alt, append=T);

cat('MODEL: \n',file=mplus.alt, append=T);
cat('    i s q| ',file=mplus.alt, append=T);
for (i in 1:T){
  cat('y',i,'@',i-1,' ',sep='',file=mplus.alt,append=T)
}
cat(';\n',file=mplus.alt,append=T)

cat('    [i * 50, s *, q *]; \n',file=mplus.alt, append=T);
cat('    i * 100, s *, q *; \n',file=mplus.alt, append=T);
cat('    q on x @ 0; \n',file=mplus.alt, append=T);
cat('    i with s * 0; \n',file=mplus.alt, append=T);
cat('    i with q * 0; \n',file=mplus.alt, append=T);
cat('    q with s * 0; \n',file=mplus.alt, append=T);
cat('    y1-y',T, ' (10);\n', sep='', file=mplus.alt, append=T);
cat('    [y1-y',T, '@0];\n', sep='', file=mplus.alt, append=T);

cat('SAVEDATA:\n', sep='', file=mplus.alt, append=T);
cat('    results=res.txt;\n', sep='', file=mplus.alt, append=T);
}

## Data generation

## quadratic model with covariate
grm.quadratic<-function(N, T, R){
## Constants and parameters
mL<-50
vL<-100
mS<-5
vS<-25
vLS<-0
mQ<--1
vQ<-9
vLQ<-0
vSQ<-0

bQ<-0

vE<-25

sigma<-array(c(vL,vLS,vLQ, vLS,vS,vSQ, vLQ,vSQ,vQ), dim=c(3,3))


y<-array(NA, dim=c(N, T))
x<-rep(c(0,1),N/2)

for (i in 1:N){ 
  mu<-c(mL,mS,mQ+x[i]*bQ)
  LSQ<-rmvnorm(1, mu, sigma)
  for (j in 1:T){
    y[i, j] <- LSQ[1] + LSQ[2]*(j-1)/(T-1) + LSQ[3]*((j-1)/(T-1))^2 + rnorm(1, 0, sqrt(vE))
  }
}

data.file<-paste('Q-N-',N,'-T-',T,'-',R,'.txt',sep='')
write.table(cbind(y,x), data.file, row.names=F, col.names=F)

}


for (i in 1:10000){
  grm.quadratic(500,5,i)
}

## process data
data.gen<-function(i,N,T, r=0.5, alpha=0){
##read in data 
data.file<-paste('Q-N-500-T-5-',i,'.txt',sep='')
  y<-read.table(data.file)
##missing data or complete data
  if (r>0){
     for (i in 2:T){
       ## find which element is mssing for y[,t-1]
       ind.c<-NULL
       ind.m<-NULL
       for (j in 1:N){
           if (!is.na(y[j,i-1])){
              ind.c<-c(ind.c,j)
           }else{
              ind.m<-c(ind.m,j)}
       }

       if (length(ind.m)!=0){
          y[ind.m,i]<-NA }

       cut.r<-quantile(y[1:N,i-1],r,na.rm=T)

       for (k in 1:N){
          if ( !is.na(y[k,i])){
             if (y[k,i-1]>cut.r){
                if (runif(1)<r*alpha){y[k,i]<-NA}
             }else{
                if (runif(1)<1-alpha+r*alpha){y[k,i]<-NA}
             }
          }
       }
     }
  }
  write.table(cbind(y[1:N,1:T],y[1:N,6]),'data.txt',row.names=F,col.names=F,na='999') 
}

## Run Mplus and get the chi-square statistics
ptm <- proc.time()
R<-10000

N.v<-c(500,400, 300, 200, 100)
T.v<-4:5
alpha.v<-c(0, 0.5, 1)
r.v<-c(.3,.1, 0)

for (T in T.v){
  gen.model(T)
  for (N in N.v){
    for (r in r.v){

       for (alpha in alpha.v){
            res.file<-paste('res-N',N,'-T',T,'-alpha',alpha,'-R',r,'-a.txt',sep='')
            for (i in 1:R){
                data.gen(i,N,T,r,alpha)
                ## Model 1
                system("c:\\progra~1\\mplus\\mplus.exe m1null.inp",show.output.on.console=F)
                if (file.exists('res.txt')){
                   m1.out<-scan('res.txt')
                   unlink('res.txt')
                }else{ m1.out<-NA}

                system("c:\\progra~1\\mplus\\mplus.exe m2alt.inp",show.output.on.console=F)
                if (file.exists('res.txt')){
                   m2.out<-scan('res.txt')
                   unlink('res.txt')
                   unlink('data.txt')
                }else{ m2.out<-NA}

                cat(c(i,m1.out,m2.out), file=res.file,append=T)
                cat("\n", file=res.file,append=T)
             }

          }
       }
    }
}
proc.time()-ptm

Submitted by: johnny
Add a comment View comment Add to my favorite Email to a friend Print
If you want to rate, please login first, or click here to register. Or you can use USERNAME: guest and PASSWORD: guest to log in.
Average score 0, based on 0 comments
1 2 3 4 5 6 7 8 9 10
Copyright © 2003-13 Zhiyong Zhang \'s Psychometric Website
Maintained by Zhiyong Zhang