Generic ComputeEM with Ramsay (1975) acceleration
[openmx:openmx.git] / models / nightly / ifa-2d-mg.R
1 # echo 0 > /proc/self/coredump_filter  # normally 023
2 # R --vanilla --no-save -f models/failing/bock-aitkin-1981.R
3 # R -d gdb --vanilla --no-save -f models/failing/bock-aitkin-1981.R
4
5 #options(error = browser)
6 require(OpenMx)
7 require(rpf)
8 #library(mvtnorm)
9
10 set.seed(7)
11 correct.LL <- 48939.35
12
13 numItems <- 30
14 numPeople <- 500
15 maxDim <- 2
16
17 items <- vector("list", numItems)
18 correct <- vector("list", numItems)
19 for (ix in 1:numItems) {
20         items[[ix]] <- rpf.grm(factors=2)
21         correct[[ix]] <- rpf.rparam(items[[ix]])
22 }
23 correct.mat <- simplify2array(correct)
24 correct.mat[1,1:10] <- 0
25 correct.mat[2,20:30] <- 0
26
27 maxParam <- max(vapply(items, function(i) rpf.numParam(i), 0))
28 maxOutcomes <- max(vapply(items, function(i) i@outcomes, 0))
29
30 data.g1 <- rpf.sample(numPeople, items, correct.mat)
31 data.g2 <- rpf.sample(numPeople, items, correct.mat, mean=c(-.5,.8), cov=matrix(c(2,.5,.5,2),nrow=2))
32 data.g3 <- rpf.sample(numPeople, items, correct.mat, mean=c(-.1,-.8), cov=matrix(c(.9,-.5,-.5,.9),nrow=2))
33
34 if (0) {
35   # for flexMIRT, write CSV
36   write.table(sapply(data.g1, unclass)-1, file="2d-mg-g1.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
37   write.table(sapply(data.g2, unclass)-1, file="2d-mg-g2.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
38   write.table(sapply(data.g3, unclass)-1, file="2d-mg-g3.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
39   fm <- read.flexmirt("/home/joshua/irt/ifa-2d-mg/2d-mg-prm.txt")
40 }
41
42 mkgroup <- function(model.name, data, latent.free) {  
43   ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam, ncol=numItems,
44                      values=c(1, 1.4, 0), free=TRUE)
45   ip.mat@values[correct.mat==0] <- 0
46   ip.mat@free[correct.mat==0] <- FALSE
47   
48   for (ix in 1:numItems) {
49     for (px in 1:3) {
50       name <- paste(c('p',ix,',',px), collapse='')
51       ip.mat@labels[px,ix] <- name
52     }
53   }
54   
55   m.mat <- mxMatrix(name="mean", nrow=1, ncol=2, values=0, free=latent.free)
56   cov.mat <- mxMatrix(name="cov", nrow=2, ncol=2, values=diag(2), free=latent.free)
57   if (latent.free) {
58     for (c1 in 1:2) {
59       for (c2 in 1:c1) {
60         cov.mat@labels[c1,c2] <- paste(model.name, paste(c1, c2, sep=","),sep="-")
61         cov.mat@labels[c2,c1] <- paste(model.name, paste(c1, c2, sep=","),sep="-")
62       }
63     }  
64   }
65   
66   m1 <- mxModel(model=model.name,
67                 ip.mat, m.mat, cov.mat,
68                 mxData(observed=data, type="raw"),
69                 mxExpectationBA81(mean="mean", cov="cov",
70                                   ItemSpec=items,
71                                   ItemParam="ItemParam",
72                                   qpoints=21, qwidth=5),
73                 mxFitFunctionML())
74   m1
75 }
76
77 groups <- paste("g", 1:3, sep="")
78
79 if (1) {
80   fm <- structure(list(G1 = structure(list(param = structure(c(0, 0.547418,  -0.746423, 0, 0.675093, -1.01778, 0, 1.00217, 0.066893, 0, 1.21546,  2.68894, 0, 1.08965, 1.93943, 0, 0.73088, -0.235478, 0, 1.69934,  0.816274, 0, 1.99294, -1.34435, 0, 1.13576, 0.805195, 0, 0.500186,  -0.281224, 0.717699, 1.43936, 0.172262, 0.978424, 0.756688, -0.643423,  1.52223, 0.709748, -0.239081, 1.14345, 2.11652, -0.719694, 0.82372,  0.57575, -0.485341, 0.820995, 2.08929, 0.575516, 0.954974, 1.58002,  -0.28948, 0.837872, 1.13202, 1.51745, 1.27315, 1.32316, -1.56062,  1.20298, 0, -1.05121, 0.971849, 0, 0.49741, 1.31535, 0, 1.01287,  1.63251, 0, 1.24117, 1.48944, 0, -0.943851, 0.664333, 0, -1.71981,  0.46719, 0, -0.518989, 0.848514, 0, 0.371468, 0.924849, 0, -0.365851,  0.641203, 0, -0.303372, 0.884883, 0, 0.600367), .Dim = c(3L,  30L), .Dimnames = list(NULL, c("i1", "i2", "i3", "i4", "i5",  "i6", "i7", "i8", "i9", "i10", "i11", "i12", "i13", "i14", "i15",  "i16", "i17", "i18", "i19", "i20", "i21", "i22", "i23", "i24",  "i25", "i26", "i27", "i28", "i29", "i30"))), mean = structure(c(0,  0), .Names = c("X6", "X7")), cov = structure(c(1, 0, 0, 1), .Dim = c(2L,  2L))), .Names = c("param", "mean", "cov")),
81                        G2 = structure(list(param = structure(c(0, 0.547418, -0.746423, 0, 0.675093,      -1.01778, 0, 1.00217, 0.066893, 0, 1.21546, 2.68894, 0, 1.08965,      1.93943, 0, 0.73088, -0.235478, 0, 1.69934, 0.816274, 0,      1.99294, -1.34435, 0, 1.13576, 0.805195, 0, 0.500186, -0.281224,      0.717699, 1.43936, 0.172262, 0.978424, 0.756688, -0.643423,      1.52223, 0.709748, -0.239081, 1.14345, 2.11652, -0.719694,      0.82372, 0.57575, -0.485341, 0.820995, 2.08929, 0.575516,      0.954974, 1.58002, -0.28948, 0.837872, 1.13202, 1.51745,      1.27315, 1.32316, -1.56062, 1.20298, 0, -1.05121, 0.971849,      0, 0.49741, 1.31535, 0, 1.01287, 1.63251, 0, 1.24117, 1.48944,      0, -0.943851, 0.664333, 0, -1.71981, 0.46719, 0, -0.518989,      0.848514, 0, 0.371468, 0.924849, 0, -0.365851, 0.641203,      0, -0.303372, 0.884883, 0, 0.600367), .Dim = c(3L, 30L), .Dimnames = list(         NULL, c("i1", "i2", "i3", "i4", "i5", "i6", "i7", "i8",          "i9", "i10", "i11", "i12", "i13", "i14", "i15", "i16",          "i17", "i18", "i19", "i20", "i21", "i22", "i23", "i24",          "i25", "i26", "i27", "i28", "i29", "i30"))), mean = structure(c(-0.507038,      0.703947), .Names = c("X6", "X7")), cov = structure(c(1.87718,      0.491825, 0.491825, 2.05385), .Dim = c(2L, 2L))), .Names = c("param",  "mean", "cov")),
82                        G3 = structure(list(param = structure(c(0, 0.547418,  -0.746423, 0, 0.675093, -1.01778, 0, 1.00217, 0.066893, 0, 1.21546,  2.68894, 0, 1.08965, 1.93943, 0, 0.73088, -0.235478, 0, 1.69934,  0.816274, 0, 1.99294, -1.34435, 0, 1.13576, 0.805195, 0, 0.500186,  -0.281224, 0.717699, 1.43936, 0.172262, 0.978424, 0.756688, -0.643423,  1.52223, 0.709748, -0.239081, 1.14345, 2.11652, -0.719694, 0.82372,  0.57575, -0.485341, 0.820995, 2.08929, 0.575516, 0.954974, 1.58002,  -0.28948, 0.837872, 1.13202, 1.51745, 1.27315, 1.32316, -1.56062,  1.20298, 0, -1.05121, 0.971849, 0, 0.49741, 1.31535, 0, 1.01287,  1.63251, 0, 1.24117, 1.48944, 0, -0.943851, 0.664333, 0, -1.71981,  0.46719, 0, -0.518989, 0.848514, 0, 0.371468, 0.924849, 0, -0.365851,  0.641203, 0, -0.303372, 0.884883, 0, 0.600367), .Dim = c(3L,  30L), .Dimnames = list(NULL, c("i1", "i2", "i3", "i4", "i5",  "i6", "i7", "i8", "i9", "i10", "i11", "i12", "i13", "i14", "i15",  "i16", "i17", "i18", "i19", "i20", "i21", "i22", "i23", "i24",  "i25", "i26", "i27", "i28", "i29", "i30"))), mean = structure(c(-0.0282586,  -0.827363), .Names = c("X6", "X7")), cov = structure(c(0.891858,  -0.422099, -0.422099, 0.835799), .Dim = c(2L, 2L))), .Names = c("param",  "mean", "cov"))), .Names = c("G1", "G2", "G3"))
83   
84   g1 <- mkgroup("g1", data.g1, FALSE)
85   g2 <- mkgroup("g2", data.g2, TRUE)
86   g3 <- mkgroup("g3", data.g3, TRUE)
87   
88   g1@matrices$ItemParam@values <- fm$G1$param
89   g2@matrices$ItemParam@values <- fm$G2$param
90   g3@matrices$ItemParam@values <- fm$G3$param
91   g2@matrices$mean@values <- t(fm$G2$mean)
92   g3@matrices$mean@values <- t(fm$G3$mean)
93   g2@matrices$cov@values <- fm$G2$cov
94   g3@matrices$cov@values <- fm$G3$cov
95   
96   cModel <- mxModel(model="cModel", g1,g2,g3,
97                     mxComputeOnce(paste(groups, 'expectation', sep='.'), context='EM'))
98   for (grp in groups) cModel@submodels[[grp]]@expectation@scores <- 'full'
99   cModel.eap <- mxRun(cModel)
100   
101   fm.sco <- suppressWarnings(try(read.table("models/nightly/data/2d-mg-sco.txt"), silent=TRUE))
102   if (is(fm.sco, "try-error")) fm.sco <- read.table("data/2d-mg-sco.txt")
103   colnames(fm.sco) <- c("group", "row", "s1", "s2", "se1", "se2", paste("cov",1:3,sep=""))
104   
105   omxCheckCloseEnough(as.matrix(fm.sco[fm.sco$group==1,-1:-2]),
106                       cModel.eap@submodels$g1@expectation@scores.out, 1e-3)
107   omxCheckCloseEnough(as.matrix(fm.sco[fm.sco$group==2,-1:-2]),
108                       cModel.eap@submodels$g2@expectation@scores.out, 1e-3)
109   omxCheckCloseEnough(as.matrix(fm.sco[fm.sco$group==3,-1:-2]),
110                       cModel.eap@submodels$g3@expectation@scores.out, 1e-3)
111
112   cModel <- mxModel(cModel,
113                     mxFitFunctionMultigroup(paste(groups, "fitfunction", sep=".")),
114                     mxComputeSequence(steps=list(
115                       mxComputeOnce(paste(groups, 'expectation', sep='.')),
116                       mxComputeOnce('fitfunction', fit=TRUE,
117                                     free.set=apply(expand.grid(groups, c('mean','cov')), 1, paste, collapse='.')))))
118   for (grp in groups) cModel@submodels[[grp]]@expectation@scores <- 'omit'
119   cModel.fit <- mxRun(cModel)
120   omxCheckCloseEnough(cModel.fit@output$minimum, correct.LL, .01)
121 }
122
123 if (1) {
124   g1 <- mkgroup("g1", data.g1, FALSE)
125   g2 <- mkgroup("g2", data.g2, TRUE)
126   g3 <- mkgroup("g3", data.g3, TRUE)
127   
128   grpModel <- mxModel(model="groupModel", g1, g2, g3,
129                       mxFitFunctionMultigroup(paste(groups, "fitfunction", sep=".")),
130                       mxComputeEM(paste(groups, 'expectation', sep='.'),
131                                   mxComputeNewtonRaphson(free.set=paste(groups, 'ItemParam', sep=".")),
132                                   mxComputeOnce('fitfunction', fit=TRUE,
133                                                 free.set=apply(expand.grid(groups, c('mean','cov')), 1, paste, collapse='.'))))
134
135   grpModel <- mxRun(grpModel, silent=TRUE)
136   omxCheckCloseEnough(grpModel@output$minimum, correct.LL, .01)
137   omxCheckCloseEnough(c(grpModel@submodels$g2@matrices$mean@values), c(-.507, .703), .01)
138   omxCheckCloseEnough(c(grpModel@submodels$g2@matrices$cov@values), c(1.877, .491, .491, 2.05), .01)
139   omxCheckCloseEnough(c(grpModel@submodels$g3@matrices$mean@values), c(-.028, -.827), .01)
140   omxCheckCloseEnough(c(grpModel@submodels$g3@matrices$cov@values), c(.892, -.422, -.422, .836), .01)
141 }