Fix regression in GrowthMixtureModelRandomStarts
[openmx:openmx.git] / models / nightly / GrowthMixtureModelRandomStarts.R
1 #
2 #   Copyright 2007-2012 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: GrowthMixtureModelRandomStarts.R  
18 # Author: Unknown
19 # Date: 9999.99.99 
20 #
21 # ModelType: Growth Mixture
22 # DataType: Continuous
23 # Field: None
24 #
25 # Purpose: 
26 #      Growth mixture model with random starts
27
28 #
29 # RevisionHistory:
30 #      Ross Gore -- 2011.06.16 updated & reformatted
31 # -----------------------------------------------------------------------
32
33 library(OpenMx)
34 demo(GrowthMixtureModel_MatrixRaw)
35
36
37
38
39 trials <- 10
40 # how many trials?
41 # -------------------------------------
42
43
44 parNames <- names(omxGetParameters(gmm))
45 # place all of the parameter 
46 # names in a vector
47 # -------------------------------------
48
49
50 input <- matrix(NA, trials, length(parNames))
51 dimnames(input) <- list(c(1: trials), c(parNames))
52
53 output <- matrix(NA, trials, length(parNames))
54 dimnames(output) <- list(c(1: trials), c(parNames))
55
56 fit <- matrix(NA, trials, 5)
57 dimnames(fit) <- list(c(1: trials), c("Minus2LL", "Status", "Iterations", "pclass1", "time"))
58 # make matrices to hold everything
59 # -------------------------------------
60
61 input[,"p1"] <- runif(trials, 0.1, 0.9)
62 input[,"p1"] <- input[,"p1"]/(1-input[,"p1"])
63 # populate the class probabilities
64 # -------------------------------------
65
66
67 v <- c("vari1", "vars1", "vari2", "vars2", "residual")
68 input[,v] <- runif(trials*5, 0, 10)
69 # populate the variances
70 # -------------------------------------
71
72
73 m <- c("meani1", "means1", "meani2", "means2")
74 input[,m] <- runif(trials*4, -5, 5)
75 # populate the means
76 # -------------------------------------
77
78
79 r <- runif(trials*2, -0.9, 0.9)
80 scale <- c(
81     sqrt(input[,"vari1"]*input[,"vars1"]),
82     sqrt(input[,"vari2"]*input[,"vars2"]))
83 input[,c("cov1", "cov2")] <- r * scale
84 # populate the covariances
85 # -------------------------------------
86
87 for (i in 1: trials){
88         temp1 <- omxSetParameters(gmm,
89                 labels=parNames,
90                 values=input[i,]
91                 )
92                 
93         temp1@name <- paste("Starting Values Set", i)
94                 
95         temp2 <- mxRun(temp1, unsafe=TRUE, suppressWarnings=TRUE)
96         
97         output[i,] <- omxGetParameters(temp2)
98         fit[i,] <- c(
99                 temp2@output$Minus2LogLikelihood,
100                 temp2@output$status[[1]],
101                 temp2@output$iterations,
102                 round(temp2$classProbs@result[1,1], 4),
103                 temp2@output$wallTime
104                 )
105         }
106         
107 fit
108 table(round(fit[,1], 3), fit[,2])
109
110 # Serial Respecification
111 # -----------------------------------------------------------------------------
112
113
114
115
116 #require(snowfall)
117 #sfInit(parallel=TRUE, cpus=4)
118 #sfLibrary(OpenMx)
119
120 topModel <- mxModel("Top")      
121
122 makeModel <- function(modelNumber){
123         temp <- mxModel(gmm, 
124                 independent=TRUE,
125                 name=paste("Iteration", modelNumber, sep=""))
126         temp <- omxSetParameters(temp,
127                 labels=parNames,
128                 values=input[modelNumber,])
129         return(temp)
130 }
131         
132 mySubs <- lapply(1:trials, makeModel)
133         
134 topModel@submodels <- mySubs
135
136 results <- mxRun(topModel, suppressWarnings=TRUE)
137
138 fitStats <- function(model){
139         retval <- c(
140                 model@output$Minus2LogLikelihood,
141                 model@output$status[[1]],
142                 model@output$iterations,
143                 round(model$classProbs@result[1,1], 4)
144                 )       
145         return(retval)
146 }
147
148 resultsFit <- t(omxSapply(results@submodels, fitStats))
149 #sfStop()
150
151 # Parallel Respecification
152 # -----------------------------------------------------------------------------
153
154
155 sum(fit[,5])
156 # Serial Time, in seconds (ignoring the overhead caused by the for loop)
157 # -------------------------------------
158
159 results@output$wallTime
160 # Parallel Time, in seconds
161 # -------------------------------------
162
163 # Compare Model Estimation Times
164 # -----------------------------------------------------------------------------
165
166
167 library(OpenMx)
168
169 # Reload OpenMx to offload the snowfall OpenMx library
170 # -----------------------------------------------------------------------------