## R function for quadratic growth rate model
grm<-function(i,T){
    ## write the model in Mplus language
    filename<-'model.inp';
    cat('TITLE: Growth rate model;\n',file=filename)
    cat('DATA: FILE IS data.txt;\n',file=filename,append=T)
    cat('VARIABLE:\n',file=filename,append=T)
    cat('  NAMES ARE math1-math10 bpi gender;\n',file=filename,append=T)
    cat('  USEVARIABLES ARE\n',file=filename,append=T)
    cat('    math1-math10 bpi gender bgender;\n',file=filename,append=T)
    cat('  MISSING=.;\n',file=filename,append=T)
    cat('DEFINE: bgender=bpi*gender;\n',file=filename,append=T)
    cat('ANALYSIS: ITERATIONS = 2000;\n',file=filename,append=T)
    cat('  COVERAGE = .00;\n',file=filename,append=T)
    cat('MODEL:\n',file=filename,append=T)
    for (j in 1:T){
       cat(paste('  level BY math',j,'@1;\n',sep=''),file=filename,append=T)
       cat(paste('  slope BY math',j,'@',j-5,';\n',sep=''),file=filename,append=T)
       cat(paste('  qslope BY math',j,'@',(j-5)^2-2*(j-5)*(i-5),';\n',sep=''),
           file=filename,append=T)
     }
    cat('  [math1-math10@0];\n',file=filename,append=T)
    cat('  [level slope qslope];\n',file=filename,append=T)
    cat('  math1-math10;\n',file=filename,append=T)
    cat('  level qslope slope;\n',file=filename,append=T)
    cat('  level on bpi gender bgender;\n',file=filename,append=T)
    cat('  slope on bpi gender bgender;\n',file=filename,append=T)
    cat('  qslope on bpi gender bgender;\n',file=filename,append=T)
    cat('SAVEDATA:\n',file=filename,append=T)
    cat('  RESULTS=parest.txt;\n',file=filename,append=T)
    cat('  FORMAT is F10.5;\n',file=filename,append=T)
    ## run Mplus
    system("c:/programs/mplus/mplus model.inp",invisible = TRUE)
}
## obtain the rate of growth
res<-NULL
for (i in 1:10){
  grm(i,10)
  out<-scan('parest.txt')
  res<-rbind(res,out)
}
## print the regression coefficients
rate.reg<-res[,c(17,45,18,46,19,47)]
colnames(rate.reg)<-c("BPI","s.e.","Gender","s.e.","interaction","s.e.")
rownames(rate.reg)<-paste('Age',6:15)
rate.reg
## print the model fit statistics
fit.stat<-t(t(res[1,c(57:61,65:72)]))
colnames(fit.stat)<-"Fit statistics"
rownames(fit.stat)<-c("Chi-square", "d.f.", "p-value","CFI","TLI",
         "AIC","BIC","a.BIC","RMSEA","Low CI","High CI","probability","SRMR")
fit.stat

## Run mplus in R for exponential growth rate analysis

grm<-function(i,T){
    ## write the model in Mplus language
    filename<-'model.inp';
    cat('TITLE: Exponential growth rate model;\n',file=filename)
    cat('DATA: FILE IS wide2.txt;\n',file=filename,append=T)
    cat('VARIABLE:\n',file=filename,append=T)
    cat('  NAMES ARE math1-math10 bpi gender;\n',file=filename,append=T)
    cat('  USEVARIABLES ARE\n',file=filename,append=T)
    cat('    math1-math10 bpi gender bsex;\n',file=filename,append=T)
    cat('  MISSING=.;\n',file=filename,append=T)
    cat('DEFINE: bsex=bpi*gender;\n',file=filename,append=T)
    cat('ANALYSIS: ITERATIONS = 20000;\n',file=filename,append=T)
    cat('  COVERAGE = .00;\n',file=filename,append=T)

    cat('MODEL:\n',file=filename,append=T)
    for (j in 1:T){
       cat(paste('  level BY math',j,'@1;\n',sep=''),file=filename,append=T)
     }
    cat('  slope BY math1* (L1);\n',file=filename,append=T) 
    for (j in 2:T){
       cat(paste('  slope BY math',j,'(L',j,');\n',sep=''),file=filename,append=T)
     }   
    cat('  [math1-math10@0];\n',file=filename,append=T)
    cat('  [level*8 slope*-6];\n',file=filename,append=T)
    cat('  math1-math10;\n',file=filename,append=T)
    cat('  level on bpi gender bsex;\n',file=filename,append=T)
    cat('  slope on bpi gender bsex;\n',file=filename,append=T)
    cat('MODEL CONSTRAINT:\n',file=filename,append=T)
    cat('  new (p1*.156);\n',file=filename,append=T)
    #cat('  p1>0;\n',file=filename,append=T)
    for (j in 1:T){
       cat(paste('  L',j,'=-exp(',i-j,'*p1)/p1;\n',sep=''),file=filename,append=T)
     }
    cat('SAVEDATA:\n',file=filename,append=T)
    cat('  RESULTS=parest.txt;\n',file=filename,append=T)
    cat('  FORMAT is F10.5;\n',file=filename,append=T)

    ## run Mplus
    system("c:/programs/mplus/mplus model.inp",invisible = TRUE)
}

## obtain the rate of growth
res<-NULL
for (i in 1:10){
  grm(i,10)
  out<-scan('parest')
  res<-rbind(res,out)
}

## print the regression coefficients 
rate.reg<-res[,c(26,58,27,59,28,60)] 
colnames(rate.reg)<-c("BPI","s.e.","Gender","s.e.","interaction","s.e.") 
rownames(rate.reg)<-paste("age",6:15) 
rate.reg 
## print the model fit statistics 
fit.stat<-t(t(res[1,c(57:61,65:72)])) 
colnames(fit.stat)<-"Fit statistics" 
rownames(fit.stat)<-c("Chi-square", "d.f.", "p-value","CFI","TLI", 
"AIC","BIC","a.BIC","RMSEA","Low CI","High CI","probability","SRMR") 
fit.stat