Move all fitfunction args to expectation
[openmx:openmx.git] / models / nightly / PPML_100days_cov.R
1 # One hundred days model
2 # Each manifest represents a single variable as measured on a different day, days 1-100
3 # Three predicting latents. Over the trial, one represents a constant part,
4 # one represents a linear part, and one an exponential part
5
6 require(OpenMx)
7 require(MASS)
8
9 # Random covariance matrix
10
11 manifests <- unlist(lapply(1:100, function(n) { paste("Day",n, sep="") } ))
12 latents <- c('Const', 'Lin', 'Exp')
13
14 amtT <- 100
15 amtStep <- amtT / 99
16 ConstLoadings <- rep(1,100)
17 LinLoadings <- unlist(lapply(1:100, function (n) { n*amtStep/(amtT-1) } ))
18 EXPFACTOR <- 0.22
19 ExpLoadings <- unlist(lapply(1:100, function (n) { -exp(-EXPFACTOR*n*amtStep) } ))
20
21 lambda <- cbind(ConstLoadings, LinLoadings, ExpLoadings)
22 dataLatents <- matrix(rnorm(3*3), 3, 3)
23 dataTest <- lambda %*% dataLatents %*% t(dataLatents) %*% t(lambda) + diag(rep(1,100))
24 dataTest <- (dataTest + t(dataTest)) / 2
25
26 dataRaw <- mvrnorm(n=5000, rep(0,100), dataTest)
27 colnames(dataRaw) <- manifests
28
29 colnames(dataTest) <- manifests
30 rownames(dataTest) <- manifests
31
32 factorModel <- mxModel("One Hundred Days Model",
33       type="RAM",
34           # Vars
35       manifestVars = manifests,
36       latentVars = latents,
37       # A
38           mxPath(from='Const',  to=manifests,value=ConstLoadings,       free=FALSE),
39           mxPath(from='Lin',    to=manifests,value=LinLoadings,         free=FALSE),
40           mxPath(from='Exp',    to=manifests,value=ExpLoadings,         free=FALSE),
41       # S
42           mxPath(from=manifests, arrows=2,value=rep(1.0, 100), labels=rep("Res", 100)),
43           #mxPath(from=latents, arrows=2,values=1.0),
44           mxPath(from="Const", to=latents, arrows=2,values=1.0),
45           mxPath(from="Lin", to=latents, arrows=2,values=1.0),
46           mxPath(from="Exp", to=latents, arrows=2,values=1.0),
47           
48           # Means
49 #         mxPath(from="one", to=latents, values=0, free=TRUE),
50           
51           # Data
52       mxData(dataTest, type="cov", numObs=100)
53         )
54
55 factorModelOut <- mxRun(imxPPML(factorModel, TRUE))