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