Updated BootstrapParallel.R demo to use mxData()
[openmx:openmx.git] / demo / BootstrapParallel.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 require(OpenMx)
17
18 # parameters for the simulation: lambda = factor loadings,
19 # specifics = specific variances
20 lambda <- matrix(c(.8, .5, .7, 0), 4, 1)
21 nObs <- 500
22 nReps <- 10
23 nVar <- nrow(lambda)
24 specifics <- diag(nVar)
25 chl <- chol(lambda %*% t(lambda) + specifics)
26
27 # indices for parameters and hessian estimate in results
28 pStrt <- 3
29 pEnd <- pStrt + 2*nVar - 1
30 hStrt <- pEnd + 1
31 hEnd <- hStrt + 2*nVar - 1
32
33 # dimension names for OpenMx
34 dn <- list()
35 dn[[1]] <- paste("Var", 1:4, sep="")
36 dn[[2]] <- dn[[1]]
37
38 # function to get a covariance matrix
39 randomCov <- function(nObs, nVar, chl, dn) {
40   x <- matrix(rnorm(nObs*nVar), nObs, nVar)
41   x <- x %*% chl
42   thisCov <- cov(x)
43   dimnames(thisCov) <- dn
44   return(thisCov)  
45 }
46
47 createNewModel <- function(index, prefix, model) {
48         modelname <- paste(prefix, index, sep='')
49         data <- mxData(randomCov(nObs, nVar, chl, dn), type="cov", numObs=nObs)
50         model <- mxModel(model, data)
51         model <- mxModel(model, name = modelname)
52         return(model)
53 }
54
55 getStats <- function(model) {
56         retval <- c(model@output$status[[1]],
57                 max(abs(model@output$gradient)),
58                 model@output$estimate,
59                 sqrt(diag(solve(model@output$hessian))))
60         return(retval)
61 }
62
63
64 # initialize obsCov for MxModel
65 obsCov <- randomCov(nObs, nVar, chl, dn)
66
67 # results matrix: get results for each simulation
68 results <- matrix(0, nReps, hEnd)
69 dnr <- c("inform", "maxAbsG", paste("lambda", 1:nVar, sep=""),
70          paste("specifics", 1:nVar, sep=""),
71          paste("hessLambda", 1:nVar, sep=""),
72          paste("hessSpecifics", 1:nVar, sep=""))
73 dimnames(results)[[2]] <- dnr
74
75 # instantiate MxModel
76 template <- mxModel(name="stErrSim",
77                        mxMatrix(name="lambda", type="Full", nrow=4, ncol=1,
78                                 free=TRUE, values=c(.8, .5, .7, 0)),
79                        mxMatrix(name="specifics", type="Diag", nrow=4,
80                                 free=TRUE, values=rep(1, 4)),
81                        mxAlgebra(lambda %*% t(lambda) + specifics,
82                                  name="preCov", dimnames=dn),
83                        mxData(observed=obsCov, type="cov", numObs=nObs),
84                        mxMLObjective(covariance='preCov'),
85                        independent = TRUE)
86
87 topModel <- mxModel(name = 'container')
88
89 submodels <- lapply(1:nReps, createNewModel, 'stErrSim', template)
90
91 names(submodels) <- imxExtractNames(submodels)
92 topModel@submodels <- submodels
93
94 modelResults <- mxRun(topModel, silent=TRUE, suppressWarnings=TRUE)
95
96 results <- t(omxSapply(modelResults@submodels, getStats))
97
98 # get rid of bad covergence results
99 results2 <- data.frame(results[which(results[,1] <= 1),])
100
101 # summarize the results
102 means <- mean(results2)
103 stdevs <- sd(results2)
104 sumResults <- data.frame(matrix(dnr[pStrt:pEnd], 2*nVar, 1,
105                                 dimnames=list(NULL, "Parameter")))
106 sumResults$mean <- means[pStrt:pEnd]
107 sumResults$obsStDev <- stdevs[pStrt:pEnd]
108 sumResults$meanHessEst <- means[hStrt:hEnd]
109 sumResults$sqrt2meanHessEst <- sqrt(2) * sumResults$meanHessEst
110
111 # print results
112 print(sumResults)
113
114