models/passing/univSatM.R

   1 #
   2 #   Copyright 2007-2010 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 require(OpenMx)
  18
  19 #Simulate Data
  20 set.seed(100)
  21 x <- rnorm (1000, 0, 1)
  22 testData <- as.matrix(x)
  23 selVars <- c("X")
  24 dimnames(testData) <- list(NULL, selVars)
  25 summary(testData)
  26 mean(testData)
  27 var(testData)
  28
  29 #example 1: Saturated Model with Cov Matrices and Paths input
  30 univSatModel1 <- mxModel("univSat1",
  31     manifestVars= selVars,
  32     mxPath(
  33         from=c("X"),
  34         arrows=2,
  35         free=T,
  36         values=1,
  37         lbound=.01,
  38         labels="vX"
  39     ),
  40     mxData(
  41         observed=var(testData),
  42         type="cov",
  43         numObs=1000
  44     ),
  45     type="RAM"
  46     )
  47 univSatFit1 <- mxRun(univSatModel1)
  48 EC1 <- mxEval(S, univSatFit1)
  49 LL1 <- mxEval(objective, univSatFit1)
  50 SL1 <- univSatFit1@output$SaturatedLikelihood
  51 Chi1 <- LL1-SL1
  52
  53 #example 1m: Saturated Model with Cov Matrices & Means and Paths input
  54 univSatModel1m <- mxModel("univSat1m",
  55     manifestVars= selVars,
  56     mxPath(
  57         from=c("X"),
  58         arrows=2,
  59         free=T,
  60         values=1,
  61         lbound=.01,
  62         labels="vX"
  63     ),
  64     mxPath(
  65         from="one",
  66         to="X",
  67         arrows=1,
  68         free=T,
  69         values=0,
  70         labels="mX"
  71     ),
  72     mxData(
  73         observed=var(testData),
  74         type="cov",
  75         numObs=1000,
  76         means=mean(testData)
  77     ),
  78     type="RAM"
  79     )
  80 univSatFit1m <- mxRun(univSatModel1m)
  81 EM1m <- mxEval(M, univSatFit1m)
  82 EC1m <- mxEval(S, univSatFit1m)
  83 LL1m <- mxEval(objective,univSatFit1m);
  84 SL1m <- univSatFit1m@output$SaturatedLikelihood
  85 Chi1m <- LL1m-SL1m
  86
  87 #example 2: Saturated Model with Raw Data and Path input
  88 univSatModel2 <- mxModel("univSat2",
  89     manifestVars= selVars,
  90     mxPath(
  91         from=c("X"),
  92         arrows=2,
  93         free=T,
  94         values=1,
  95         lbound=.01,
  96         labels="vX"
  97     ),
  98     mxData(
  99         observed=testData,
 100         type="raw",
 101     ),
 102     type="RAM"
 103     )
 104 univSatFit2 <- mxRun(univSatModel2)
 105 EM2 <- mxEval(M, univSatFit2)
 106 EC2 <- mxEval(S, univSatFit2)
 107 LL2 <- mxEval(objective,univSatFit2);
 108
 109 #example 2s: Saturated Model with Raw Data and Path input built upon Cov/Means version
 110 univSatModel2s <- mxModel(univSatModel1,
 111     mxData(
 112         observed=testData,
 113         type="raw"
 114     ),
 115     type="RAM"
 116     )
 117 univSatFit2s <- mxRun(univSatModel2s)
 118 EM2s <- mxEval(M, univSatFit2s)
 119 EC2s <- mxEval(S, univSatFit2s)
 120 LL2s <- mxEval(objective,univSatFit2s);
 121
 122
 123 #example 3: Saturated Model with Cov Matrices and Matrices input
 124 univSatModel3 <- mxModel("univSat3",
 125     mxMatrix(
 126         type="Symm",
 127         nrow=1,
 128         ncol=1,
 129         free=T,
 130         values=1,
 131         dimnames=list(selVars,selVars),
 132         name="expCov"
 133     ),
 134     mxData(
 135         observed=var(testData),
 136         type="cov",
 137         numObs=1000
 138     ),
 139     mxMLObjective(
 140         covariance="expCov"
 141     )
 142     )
 143 univSatFit3 <- mxRun(univSatModel3)
 144 EC3 <- mxEval(expCov, univSatFit3)
 145 LL3 <- mxEval(objective, univSatFit3)
 146 SL3 <- univSatFit3@output$SaturatedLikelihood
 147 Chi3 <- LL3-SL3
 148
 149 #example 3m: Saturated Model with Cov Matrices & Means and Matrices input
 150 univSatModel3m <- mxModel("univSat3m",
 151     mxMatrix(
 152         type="Symm",
 153         nrow=1,
 154         ncol=1,
 155         free=T,
 156         values=1,
 157         dimnames=list(selVars,selVars),
 158         name="expCov"
 159     ),
 160     mxMatrix(
 161         type="Full",
 162         nrow=1,
 163         ncol=1,
 164         free=T,
 165         values=0,
 166         dimnames=list(NULL,selVars),
 167         name="expMean"
 168     ),
 169     mxData(
 170         observed=var(testData),
 171         type="cov",
 172         numObs=1000,
 173         means=mean(testData)
 174     ),
 175     mxMLObjective(
 176         covariance="expCov",
 177         means="expMean"
 178     )
 179     )
 180 univSatFit3m <- mxRun(univSatModel3m)
 181 EM3m <- mxEval(expMean, univSatFit3m)
 182 EC3m <- mxEval(expCov, univSatFit3m)
 183 LL3m <- mxEval(objective, univSatFit3m);
 184 SL3m <- univSatFit3m@output$SaturatedLikelihood
 185 Chi3m <- LL3m-SL3m
 186
 187 #examples 4: Saturated Model with Raw Data and Matrices input
 188 univSatModel4 <- mxModel("univSat4",
 189     mxMatrix(
 190         type="Symm",
 191         nrow=1,
 192         ncol=1,
 193         free=T,
 194         values=1,
 195         dimnames=list(selVars,selVars),
 196         name="expCov"
 197     ),
 198     mxMatrix(
 199         type="Full",
 200         nrow=1,
 201         ncol=1,
 202         free=T,
 203         values=0,
 204         dimnames=list(NULL,selVars),
 205         name="expMean"
 206     ),
 207     mxData(
 208         observed=testData,
 209         type="raw",
 210     ),
 211     mxFIMLObjective(
 212         covariance="expCov",
 213         means="expMean"
 214     )
 215     )
 216 univSatFit4 <- mxRun(univSatModel4)
 217 EM4 <- mxEval(expMean, univSatFit4)
 218 EC4 <- mxEval(expCov, univSatFit4)
 219 LL4 <- mxEval(objective, univSatFit4);
 220
 221
 222 #Mx answers hard-coded
 223 #example Mx..1: Saturated Model with Cov Matrices
 224 Mx.EC1 <-  1.062112
 225 Mx.LL1 <- -1.474434e-17
 226
 227 #example Mx..1m: Saturated Model with Cov Matrices & Means
 228 Mx.EM1m <- 0.01680509
 229 Mx.EC1m <- 1.062112
 230 Mx.LL1m <- -1.108815e-13
 231
 232 Mx.EM2 <- 0.01680516
 233 Mx.EC2 <- 1.061050
 234 Mx.LL2 <- 2897.135
 235
 236
 237 #OpenMx summary
 238 cov <- rbind(cbind(EC1,EC1m,EC2),cbind(EC3,EC3m,EC4))
 239 mean <- rbind(cbind(EM1m, EM2),cbind(EM3m,EM4))
 240 like <- rbind(cbind(LL1,LL1m,LL2),cbind(LL3,LL3m,LL4))
 241 cov; mean; like
 242
 243 #old Mx summary
 244 Mx.cov <- cbind(Mx.EC1,Mx.EC1m,Mx.EC2)
 245 Mx.mean <- cbind(Mx.EM1m,Mx.EM2)
 246 Mx.like <- cbind(Mx.LL1,Mx.LL1m,Mx.LL2)
 247 Mx.cov; Mx.mean; Mx.like
 248
 249
 250 #Compare OpenMx results to Mx results (LL: likelihood; EC: expected covariance, EM: expected means)
 251 #1:CovPat
 252 omxCheckCloseEnough(Chi1,Mx.LL1,.001)
 253 omxCheckCloseEnough(EC1,Mx.EC1,.001)
 254 #1m:CovMPat
 255 omxCheckCloseEnough(Chi1m,Mx.LL1m,.001)
 256 omxCheckCloseEnough(EC1m,Mx.EC1m,.001)
 257 omxCheckCloseEnough(EM1m,Mx.EM1m,.001)
 258 #2:RawPat
 259 omxCheckCloseEnough(LL2,Mx.LL2,.001)
 260 omxCheckCloseEnough(EC2,Mx.EC2,.001)
 261 omxCheckCloseEnough(EM2,Mx.EM2,.001)
 262 #2:RawSPat
 263 omxCheckCloseEnough(LL2s,Mx.LL2,.001)
 264 omxCheckCloseEnough(EC2s,Mx.EC2,.001)
 265 omxCheckCloseEnough(EM2s,Mx.EM2,.001)
 266 #3:CovMat
 267 omxCheckCloseEnough(Chi3,Mx.LL1,.001)
 268 omxCheckCloseEnough(EC3,Mx.EC1,.001)
 269 #3m:CovMPat
 270 omxCheckCloseEnough(Chi3m,Mx.LL1m,.001)
 271 omxCheckCloseEnough(EC3m,Mx.EC1m,.001)
 272 omxCheckCloseEnough(EM3m,Mx.EM1m,.001)
 273 #4:RawMat
 274 omxCheckCloseEnough(LL4,Mx.LL2,.001)
 275 omxCheckCloseEnough(EC4,Mx.EC2,.001)
 276 omxCheckCloseEnough(EM4,Mx.EM2,.001)
 277