testcases/case1.R

In brandmaier/semper: Structural Equation Modeling, power, effect size, and reliability

#
# Five occasion LGCM, no latent correlation, no missingness
# Power curves for different mean slope sizes
#
# see also: Figure 3 in http://nd.psychstat.org/_media/research/zhangwang2009.pdf
require("semper")
#
# This model specification was automatically generated by Onyx
# 
require("OpenMx");

manifests<-c("x1","x2","x3","x4","x5")
latents<-c("icept","slope")
model <- mxModel("Unnamed_Model", 
                 type="RAM",
                 manifestVars = manifests,
                 latentVars = latents,
                 mxPath(from="icept",to=c("x1","x2","x3","x4","x5"), free=c(FALSE,FALSE,FALSE,FALSE,FALSE), value=c(1.0,1.0,1.0,1.0,1.0) , arrows=1, label=c("icept__x1","icept__x2","icept__x3","icept__x4","icept__x5") ),
                 mxPath(from="slope",to=c("x2","x3","x4","x5"), free=c(FALSE,FALSE,FALSE,FALSE), value=c(1.0,2.0,3.0,4.0) , arrows=1, label=c("slope__x2","slope__x3","slope__x4","slope__x5") ),
                 mxPath(from="one",to=c("icept","slope"), 
                        free=c(TRUE,TRUE), value=c(10.0,0.0) , arrows=1, label=c("__mu_L","__mu_S") ),
                 mxPath(from="icept",to=c("icept"), free=c(TRUE), value=c(4.0) , arrows=2, label=c("__sigma_2_L") ),
                 mxPath(from="slope",to=c("slope"), free=c(TRUE), 
                        value=c(1.0) , arrows=2, label=c("__sigma_2_S") ),
                 mxPath(from="x1",to=c("x1"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("e") ),
                 mxPath(from="x2",to=c("x2"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("e") ),
                 mxPath(from="x3",to=c("x3"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("e") ),
                 mxPath(from="x4",to=c("x4"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("e") ),
                 mxPath(from="x5",to=c("x5"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("e") ),
                 mxPath(from="one",to=c("x1","x2","x3","x4","x5"), free=F, value=0, arrows=1),
                 mxData(data.frame(x1=0,x2=0,x3=0,x4=0,x5=0), type = "raw")
);

model <- omxSetParameters(model, labels="__mu_S", values=.1)
simp1 <- simPowerZeroRestriction(trueModel=model, restrictions="__mu_S",N=c(10,10000),repetitions=1000)
model <- omxSetParameters(model, labels="__mu_S", values=.2)
simp2 <- simPowerZeroRestriction(trueModel=model, restrictions="__mu_S",N=c(10,10000),repetitions=1000)
model <- omxSetParameters(model, labels="__mu_S", values=.3)
simp3 <- simPowerZeroRestriction(trueModel=model, restrictions="__mu_S",N=c(10,10000),repetitions=1000)

prng <- c(seq(0.0001,0.99999,0.01),0.99999)
plot.simPower(simp3, xlim=c(0,1000), prng=prng)
plot(simp2,add=T,lty=2, prng=prng)
plot(simp1,add=T,lty=3, prng=prng)