ifa: Store latent mean & cov in regular MxMatrix
[openmx:openmx.git] / models / failing / ExtraOrdinalVars.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 #Ordinal Data test, based on poly3dz.mx
19
20 # Data
21 data <- read.table("../passing/data/mddndzf.dat", na.string=".", col.names=c("t1neur1", "t1mddd4l", "t2neur1", "t2mddd4l"))
22
23 nthresh1 <- 1
24 nthresh2 <- 12
25 diff <- nthresh2-nthresh1
26 nvar <- 4
27
28 Mx1Threshold <- rbind(
29 c(-1.9209, 0.3935, -1.9209, 0.3935),
30 c(-0.5880, NA    , -0.5880, NA    ),
31 c(-0.0612, NA    , -0.0612, NA    ),
32 c( 0.3239, NA    ,  0.3239, NA    ),
33 c( 0.6936, NA    ,  0.6936, NA    ),
34 c( 0.8856, NA    ,  0.8856, NA    ),
35 c( 1.0995, NA    ,  1.0995, NA    ),
36 c( 1.3637, NA    ,  1.3637, NA    ),
37 c( 1.5031, NA    ,  1.5031, NA    ),
38 c( 1.7498, NA    ,  1.7498, NA    ),
39 c( 2.0733, NA    ,  2.0733, NA    ),
40 c( 2.3768, NA    ,  2.3768, NA    ))
41
42 Mx1R <- rbind(
43     c(1.0000,  0.2955,  0.1268, 0.0760),
44     c(0.2955,  1.0000, -0.0011, 0.1869),
45     c(0.1268, -0.0011,  1.0000, 0.4377),
46     c(0.0760,  0.1869,  0.4377, 1.0000))
47
48 nameList <- names(data)
49 data$breaking <- rep(1, length(data$t1neur1))  ### <-- This line causes the model to break.
50
51 # Define the model
52 model <- mxModel('model')
53 model <- mxModel(model, mxMatrix("Stand", name = "R", nrow = nvar, ncol = nvar, free=TRUE))
54 model <- mxModel(model, mxMatrix("Zero", name = "M", nrow = 1, ncol = nvar, free=FALSE))
55 model <- mxModel(model, mxMatrix("Full", 
56             name="thresh", 
57             values=cbind(
58                     seq(-1.9, 1.9, by=(3.8/(nthresh2-1))),  # t1Neur1: 12 thresholds evenly spaced from -1.9 to 1.9
59                    c(rep(1, nthresh1),rep(NA, diff)),       # t1mddd4l: 1 threshold at 1
60                     seq(-1.9, 1.9, by=(3.8/(nthresh2-1))),  # t2Neur1: 12 thresholds same as t1Neur1
61                    c(rep(1, nthresh1),rep(NA, diff))        # t2mddd4l: 1 threshold same as t1mddd4l
62                     ),
63             free = c(rep(c( rep(TRUE, nthresh2), 
64                             rep(TRUE, nthresh1), rep(FALSE, diff)
65                             ), 2)), 
66             dimnames = list(c(), nameList), 
67             labels = rep(c(paste("neur", 1:nthresh2, sep=""),
68                         paste("mddd4l", 1:nthresh1, sep=""), rep(NA, diff))
69                         )))
70
71 # Define the objective function
72 objective <- mxFIMLObjective(covariance="R", means="M", dimnames=nameList, thresholds="thresh")
73
74 # Define the observed covariance matrix
75 dataMatrix <- mxData(data, type='raw')
76
77 # Add the objective function and the data to the model
78 model <- mxModel(model, objective, dataMatrix)
79
80 # Run the job
81 model <- mxRun(model)
82
83 estimates <- model@output$estimate
84
85 # Results from old Mx:
86 omxCheckCloseEnough(mxEval(thresh, model)[,1], Mx1Threshold[,1], 0.01)
87 omxCheckCloseEnough(mxEval(thresh, model)[1,2], Mx1Threshold[1,2], 0.01)
88 omxCheckCloseEnough(mxEval(R, model), Mx1R, 0.01)
89 omxCheckCloseEnough(model@output$Minus2LogLikelihood, 4081.48, 0.02)