====== Program scripts for growth rate analysis ======
===== Quadratic growth rate model =====
## 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
===== Exponential growth rate model =====
## 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