Move all fitfunction args to expectation
[openmx:openmx.git] / models / nightly / RegimeSwitching_MatrixRawNoCholDifferentTransitions.R
1 #\r
2 #   Copyright 2007-2010 The OpenMx Project\r
3 #\r
4 #   Licensed under the Apache License, Version 2.0 (the "License");\r
5 #   you may not use this file except in compliance with the License.\r
6 #   You may obtain a copy of the License at\r
7\r
8 #        http://www.apache.org/licenses/LICENSE-2.0\r
9\r
10 #   Unless required by applicable law or agreed to in writing, software\r
11 #   distributed under the License is distributed on an "AS IS" BASIS,\r
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\r
13 #   See the License for the specific language governing permissions and\r
14 #   limitations under the License.\r
15 \r
16 # -----------------------------------------------------------------------------\r
17 # Program: RegimeSwitching_MatrixRaw.R  \r
18 # Author: Michael Neale\r
19 # Date: 2010.09.17 \r
20 #\r
21 # ModelType: Growth Mixture with Switching\r
22 # DataType: Continuous\r
23 # Field: None\r
24 #\r
25 # Purpose: \r
26 #       Regime Switching Growth Mixture Model\r
27 #       Matrix style model input - Raw data input\r
28 #       Borrows from seminal work by Conor Dolan et al 2005 in SEM Journal\r
29 #\r
30 # RevisionHistory:\r
31 # -----------------------------------------------------------------------------\r
32 \r
33 # Load Libraries\r
34 require(OpenMx)\r
35 # -----------------------------------------------------------------------------\r
36 \r
37 # Prepare Data\r
38 dolan2005Data <- read.table(file='data/sel.txt',header=FALSE,na.strings='-9.00',col.names=c(paste("time",1:4,sep="")))\r
39 # -----------------------------------------------------------------------------\r
40 \r
41 # Number of occasions, factors & regimes\r
42 nocc <- 4\r
43 nfac <- 2\r
44 nregime <- 3\r
45 nfacxnregime <- nfac * nregime\r
46 # -----------------------------------------------------------------------------\r
47 \r
48 # Set up labels for the models; these populate larger vectors\r
49 baseMeanLabels<-c("meanIntercept","meanSlope")\r
50 baseCovarianceLabels<-c("varIntercept","covInterceptSlope","covInterceptSlope","varSlope")\r
51 # -----------------------------------------------------------------------------\r
52 \r
53 # Set up symmetric factor covariance matrix - same thing for all nregime^nocc models\r
54 factorCovariances <- mxMatrix(type="Symm",nfacxnregime,nfacxnregime,name="factorCovariances")\r
55 \r
56 j <- 1\r
57 for (i in 1:nregime)\r
58\r
59     factorCovariances@labels[j:(j+1),j:(j+1)] <- paste(baseCovarianceLabels,i,sep="")\r
60     j <- j + nfac\r
61 }\r
62 factorCovariances@free[2,2] <- TRUE\r
63 factorCovariances@free[6,6] <- TRUE\r
64 factorCovariances@lbound[2,2] <- .0001\r
65 factorCovariances@lbound[6,6] <- .0001\r
66 factorCovariances@values[2,2] <- .1\r
67 factorCovariances@values[6,6] <- .05\r
68 #\r
69 # -----------------------------------------------------------------------------\r
70 \r
71 # Set up base factor loading matrix; its values change with each regime-switched models\r
72 factorLoadings <- mxMatrix(type="Full", nrow=nocc, ncol=nfacxnregime, values=c(rep(1,nocc),0:(nocc-1),rep(0,(nocc*(nfacxnregime-2)))), free=F, name="factorLoadings")\r
73 # -----------------------------------------------------------------------------\r
74 \r
75 # Factor Mean matrix; same for all models.  Funky method of getting cartesian product of baseLabels and 1:nregime\r
76 factorMeans <- mxMatrix(type="Full", nrow=nfacxnregime, ncol=1, free=c(T,T,T,F,T,T), values=c(3.35,0.45,10.1,0,1.17,0.22), name="factorMeans")\r
77 factorMeans@labels[] <- apply(expand.grid(baseMeanLabels,1:nregime),1,function(x) paste(x,collapse=""))\r
78 # -----------------------------------------------------------------------------\r
79 \r
80 # Set up base residual variance matrix; its parameter labels change with each regime-switched model\r
81 residuals <- mxMatrix(type="Diag", nrow=nocc, ncol=nocc, free=T, values=1, lbound=.1, labels="residual1", name="residuals")\r
82 # -----------------------------------------------------------------------------\r
83 \r
84 \r
85 # Create a generic MxModel object for Regime 1, allowing for nregime * nfac latent variables which influence nocc observed variables\r
86 # It is done this way to allow for possible covariances between Regime latent growth curve factors\r
87 regime1 <- mxModel("regime1", factorCovariances, residuals, factorMeans, factorLoadings,\r
88     mxAlgebra(expression=factorLoadings %&% (factorCovariances) + residuals, name="expectedCovariances"),\r
89     mxAlgebra(expression=t(factorLoadings %*% factorMeans), name="expectedMeans"),\r
90     mxExpectationNormal(covariance="expectedCovariances", means="expectedMeans", dimnames = names(dolan2005Data)),\r
91                    mxFitFunctionML(vector=TRUE)\r
92                  )\r
93 # -----------------------------------------------------------------------------\r
94 \r
95 # Create a list of MxModel objects and their names for regimeswitches 1...nregime^nocc by overwriting the parameter labels \r
96 #   of the factor means and covariances of class 1's model\r
97 modelList <- vector("list",(nregime^nocc))\r
98 modelNames <- vector("list",(nregime^nocc))\r
99 # -----------------------------------------------------------------------------\r
100 \r
101 # Construct a nregime^nocc by nocc matrix with the relevant switching patterns in it\r
102 switchPatterns<-matrix(0,nregime^nocc,0)\r
103 for (j in 1:nocc) {\r
104     conor <- matrix(1,nregime,nocc)\r
105     conor[,j]<-1:nregime\r
106     tmp1 <- conor[,1]\r
107     for (k in 2:nocc) {\r
108         tmp1 <- tmp1 %x% conor[,k]\r
109     }\r
110     switchPatterns<-cbind(switchPatterns,tmp1)\r
111     }\r
112 # -----------------------------------------------------------------------------\r
113 \r
114 # Use the switching patterns to modify the models in the list\r
115 # First set up a couple matrices to hold parameter label numbers and starting values\r
116 resNumber <- matrix(c(rep((1:nregime),nocc)),nrow=nregime,ncol=nocc)\r
117 resValues <- matrix(rep(c(2,5,3)),nrow=nregime,ncol=nocc)\r
118 # Then another one for the basis functions (factor loadings)\r
119 ncomponents <- 2\r
120 # NB, lots of things need to change, including line below, if ncomponents >2, particularly factorCovariances and factorLoadings\r
121 basisFunctions <- matrix( c(rep(1,nocc),0:(nocc-1)),nocc,ncomponents)\r
122 # This loop makes use of a patternVector z to describe the particular combination of regimes that a person may be in\r
123 ii <- 0\r
124 for (i in 1:(nregime^nocc))\r
125 {\r
126     temp <- regime1\r
127     patternVector <- t(switchPatterns[i,])\r
128     for (j in 1:nocc) \r
129         {\r
130         ii <- ii+1\r
131         z <- matrix(0,1,nregime)\r
132         z[patternVector[j]] <- 1\r
133         temp$residuals@labels[j,j] <- paste("residual",resNumber[patternVector[j],j],sep="")\r
134         temp$residuals@values[j,j] <- resValues[patternVector[j],j]\r
135         temp$factorLoadings@values[j,] <- z %x% basisFunctions[j,]\r
136         }\r
137     temp@name <- paste("Regime",i,sep="")\r
138     modelNames[i] <- paste("Regime",i,sep="")\r
139     modelList[i] <- temp\r
140 }\r
141 # -----------------------------------------------------------------------------\r
142 # Now we construct the weight matrix algebra\r
143 # in Classic Mx it was:\r
144 # g unit nr 1                   \r
145 # p full nr 1 fr                                ! probs initial\r
146 # q full nr nr fr                       ! transition probs \r
147 # end matrices;\r
148 # !\r
149 # begin algebra;                \r
150 # d=\m2v(q');                           !  transition probs into vector\r
151 # ! mixing proportions \r
152 # ! a=(p@g@g).(d@g).(g@d)                                       ! nt=3\r
153 # a=(p@g@g@g).(d@g@g).(g@d@g).(g@g@d);          ! nt=4\r
154 # and I am not going to f with making it general just now...\r
155 \r
156 # To avoid constraints, we express initial and transition matrices as proportions by bounding them to be positive and expressing as proportions (TBD)\r
157 rawInitialProbs <- mxMatrix(type="Full", nrow=nregime, ncol=1, free=c(F,rep(T,(nregime-1))), values=c(.35,.15,.50), lbound=.00001, name="rawInitialProbs")\r
158 rawTransitionProbs1 <- mxMatrix(type="Full", nrow=nregime, ncol=nregime, free=c(F,rep(T,(nregime-1))), values=c(1,.2,.1,1,.7,.05,1,.05,.65), lbound=.00001, name="rawTransitionProbs1")\r
159 initialProbs <- mxAlgebra(rawInitialProbs %x% (1/sum(rawInitialProbs)), name="initialProbs")\r
160 # was\r
161 #transitionProbs1 <- mxAlgebra(rawTransitionProbs1 %x% (1/sum(rawTransitionProbs1)), name="transitionProbs1")\r
162 transitionProbs1 <- mxAlgebra(rawTransitionProbs1 / (unitNregime %x% (t(unitNregime) %*% rawTransitionProbs1)), name="transitionProbs1")\r
163 transitionVector1 <- mxAlgebra(cvectorize(transitionProbs1), name="transitionVector1")\r
164 \r
165 rawTransitionProbs2 <- mxMatrix(type="Full", nrow=nregime, ncol=nregime, free=c(F,rep(T,(nregime-1))), values=c(1,.2,.1,1,.7,.05,1,.05,.65), lbound=.00001, name="rawTransitionProbs2")\r
166 # was\r
167 #transitionProbs2 <- mxAlgebra(rawTransitionProbs2 %x% (1/sum(rawTransitionProbs2)), name="transitionProbs2")\r
168 transitionProbs2 <- mxAlgebra(rawTransitionProbs2 / (unitNregime %x% (t(unitNregime) %*% rawTransitionProbs2)), name="transitionProbs2")\r
169 transitionVector2 <- mxAlgebra(cvectorize(transitionProbs2), name="transitionVector2")\r
170 \r
171 rawTransitionProbs3 <- mxMatrix(type="Full", nrow=nregime, ncol=nregime, free=c(F,rep(T,(nregime-1))), values=c(1,.2,.1,1,.7,.05,1,.05,.65), lbound=.00001, name="rawTransitionProbs3")\r
172 # was\r
173 #transitionProbs3 <- mxAlgebra(rawTransitionProbs3 %x% (1/sum(rawTransitionProbs3)), name="transitionProbs3")\r
174 transitionProbs3 <- mxAlgebra(rawTransitionProbs3 / (unitNregime %x% (t(unitNregime) %*% rawTransitionProbs3)), name="transitionProbs3")\r
175 transitionVector3 <- mxAlgebra(cvectorize(transitionProbs3), name="transitionVector3")\r
176 \r
177 \r
178 unitNregime <- mxMatrix(type="Unit",nrow=nregime, ncol=1, name="unitNregime")\r
179 weights <- mxAlgebra((initialProbs %x% unitNregime %x% unitNregime %x% unitNregime) * (transitionVector1 %x% unitNregime %x% unitNregime) * (unitNregime %x% transitionVector2 %x% unitNregime) * (unitNregime %x% unitNregime %x% transitionVector3), name="weights")\r
180 # -----------------------------------------------------------------------------\r
181 \r
182 # Build up objective function using generous doses of paste\r
183 objectives <- paste(modelNames, 'objective', sep=".")\r
184 modelnumbers <- 1:(nregime^nocc)\r
185 \r
186 components <- paste("weights[", modelnumbers, ",1]", " %x% ", objectives, sep = '')\r
187 componentsSum <- paste(components, collapse = " + ")\r
188 algebraString <- paste("mxAlgebra(-2*sum(log(", componentsSum, ")), name='mixtureObj')", sep = "")\r
189 algebraObjective <- eval(parse(text=algebraString)[[1]])\r
190 obj <- mxFitFunctionAlgebra("mixtureObj")\r
191 # -----------------------------------------------------------------------------\r
192       \r
193 # Construct MxModel from all the various pieces; add options and then run it and summarize results\r
194 rsgcmmDiffT <- mxModel("Regime Switching Growth Curve Model", mxData(observed=dolan2005Data, type="raw"), modelList, rawInitialProbs, initialProbs, rawTransitionProbs1, transitionProbs1, transitionVector1, rawTransitionProbs2, transitionProbs2, transitionVector2, rawTransitionProbs3, transitionProbs3, transitionVector3, unitNregime, weights, algebraObjective, obj)      \r
195 rsgcmmDiffT <- mxOption(rsgcmmDiffT, 'Checkpoint Count', 1)\r
196 #rsgcmm <- mxOption(rsgcmm, "Standard Errors", "No")\r
197 #rsgcmm <- mxOption(rsgcmm, "Calculate Hessian", "No")\r
198 \r
199 rsgcmmDiffTFit <- mxRun(rsgcmmDiffT, unsafe=TRUE)\r
200 \r
201 summary(rsgcmmDiffTFit)\r
202 \r
203 # -----------------------------------------------------------------------------\r
204 \r
205 # Check results to see if they are within specified bounds\r
206 #omxCheckCloseEnough(rsgcmmDiffTFit@output$Minus2LogLikelihood, 9918.364, 0.1)\r
207 #omxCheckCloseEnough(max(mxEval(rsgcmmFit@output$transitionProbs, rsgcmmFit)), 0.4790, 0.01)\r
208 #omxCheckCloseEnough(min(mxEval(rsgcmmFit@output$transitionProbs, rsgcmmFit)), 0.1608, 0.01)\r
209 ## -----------------------------------------------------------------------------\r