Updated copyright to 2013 for R/ demo/ models/passing and src/ folders, and also...
[openmx:openmx.git] / models / passing / UnivariateRandomInterceptWide.R
1 #
2 #   Copyright 2007-2013 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16 # ---------------------------------------------------------------------
17 # Program: UniRandomIntTest-120815.R
18 #  Author: Steve Boker
19 #    Date: Wed Aug 15 10:50:12 CEST 2012
20 #
21 # This program simulates some univariate multilevel data with random intercepts only,
22 #   fits it with lme(), fits a naive wide format multilevel OpenMx model and
23 #   checks the results
24 #
25 # ---------------------------------------------------------------------
26 # Revision History
27 #   Steve Boker -- Wed Aug 15 10:50:14 CEST 2012
28 #      Created UniRandomIntTest-120815.R
29 #
30 # ---------------------------------------------------------------------
31
32 # ----------------------------------
33 # Read libraries and set options.
34
35 options(width=110)
36 library(nlme)
37 library(OpenMx)
38
39 # ----------------------------------
40 # Set constants.
41
42 sdLevelOneE <- sqrt(.2)
43 sdIntercepts <- sqrt(.5)
44 sdX <- sqrt(1)
45
46 N <- 400    # number of participants
47 P <- 100    # number of observations per participant
48 b0 <- .5    # Fixed effect intercept
49 b1 <- .8    # Fixed effect slope
50
51 set.seed(1)
52
53 # ----------------------------------
54 # Simulate the data.
55
56 X <- rnorm(N*P, 0, sd=sdX)
57 ID <- rep(1:N, each=P)
58 b0i <- b0 + rnorm(N, 0, sd=sdIntercepts)
59 Y <- rep(b0i, each=P) + b1*X + rnorm(N*P, 0, sd=sdLevelOneE)
60
61 SimUniRandomIntFrame <- data.frame(ID, X, Y)
62
63 # ----------------------------------
64 # Test with lme().
65
66 lmeOut <- summary(lme(Y ~ X, random= list(~ 1 | ID), data=SimUniRandomIntFrame))
67
68 # ----------------------------------
69 # Set constants.
70
71 theIDs <- unique(SimUniRandomIntFrame$ID)
72 totalN <- length(theIDs)
73 totalVars <- 2
74
75 maxP <- 0
76 for (tID in theIDs) {
77     tLen <- length(SimUniRandomIntFrame$ID[SimUniRandomIntFrame$ID==tID])
78     if (tLen > maxP) 
79         maxP <- tLen
80 }
81
82 # ----------------------------------
83 # Wide-format the data frame from tall format.
84
85 wideMatrix <- matrix(NA, nrow=totalN, ncol=1 + (maxP*totalVars))
86 colnames(wideMatrix) <- c("ID", paste("Y",1:maxP, sep=""), paste("X",1:maxP, sep=""))
87 i <- 1
88 for (tID in theIDs) {
89     wideMatrix[i, 1] <- tID
90     tY <- SimUniRandomIntFrame$Y[SimUniRandomIntFrame$ID==tID]
91     wideMatrix[i, 2:(length(tY)+1)] <- tY
92     tX <- SimUniRandomIntFrame$X[SimUniRandomIntFrame$ID==tID]
93     wideMatrix[i, (2+maxP):(length(tY)+1+maxP)] <- tX
94     i <- i + 1
95 }
96 wideFrame <- data.frame(wideMatrix)
97
98 manifestNames <- colnames(wideFrame)[2:dim(wideFrame)[2]]
99 xNames <- paste("X",1:maxP, sep="")
100 yNames <- paste("Y",1:maxP, sep="")
101 latentNames <- c("b0i")
102
103 # ----------------------------------
104 # Build the OpenMx wide model.
105
106 OpenMxModelUniRandomIntModel1 <- mxModel("OpenMxModelUniRandomIntModel1",
107         type="RAM", 
108         manifestVars=manifestNames,
109     latentVars=latentNames,
110     mxPath(from=xNames, to=yNames, connect="single", arrows=1, free=TRUE, values=.2, labels="b1"),
111     mxPath(from=xNames, to=xNames, connect="single", arrows=2, free=TRUE, values=.8, labels="vX"),
112     mxPath(from=yNames, to=yNames, connect="single", arrows=2, free=TRUE, values=.8, labels="eY"),
113     mxPath(from=latentNames, to=yNames, arrows=1, free=FALSE, values=1),
114     mxPath(from=latentNames, to=latentNames, connect="single", arrows=2, free=TRUE, values=.8, labels="vb0i"),
115     mxPath(from="one", to=c(xNames), arrows=1, free=TRUE, values=1, labels="mX"),
116     mxPath(from="one", to=c(latentNames), arrows=1, free=TRUE, values=1, labels="mb0i"),
117     mxData(observed=wideFrame, type="raw")
118 )
119
120 # ----------------------------------
121 # Fit the model and examine the summary results.
122
123 OpenMxModelUniRandomIntModel1Fit <- mxRun(OpenMxModelUniRandomIntModel1)
124
125 summary(OpenMxModelUniRandomIntModel1Fit)
126
127 omxCheckCloseEnough(lmeOut$coefficients$fixed[1], mxEval(mb0i, model=OpenMxModelUniRandomIntModel1Fit), 0.001)
128
129 omxCheckCloseEnough(lmeOut$coefficients$fixed[2], mxEval(b1, model=OpenMxModelUniRandomIntModel1Fit), 0.001)
130
131 omxCheckCloseEnough(lmeOut$sigma, mxEval(sqrt(eY), model=OpenMxModelUniRandomIntModel1Fit), 0.001)
132
133 omxCheckCloseEnough(sd(c(lmeOut$coefficients$random$ID)), mxEval(sqrt(vb0i), model=OpenMxModelUniRandomIntModel1Fit), 0.001)
134