Multigroup two-tier
[openmx:openmx.git] / models / nightly / ifa-cai2009.R
1 # This data is from an email:
2 #
3 # Date: Wed, 06 Feb 2013 19:49:24 -0800
4 # From: Li Cai <lcai@ucla.edu>
5 # To: Joshua N Pritikin <jpritikin@pobox.com>
6 # Subject: Re: how did you control item bias in Cai (2010, p. 592) ?
7
8 #options(error = utils::recover)
9 library(OpenMx)
10 library(rpf)
11 library(ggplot2)
12 library(stringr)
13
14 correct.LL <- 29995.30
15 data.raw <- suppressWarnings(try(read.csv("models/nightly/data/cai2009.csv"), silent=TRUE))
16 if (is(data.raw, "try-error")) data.raw <- read.csv("data/cai2009.csv")
17 data.g1 <- as.data.frame(data.raw[data.raw$G==1, 2:13])
18 data.g2 <- as.data.frame(data.raw[data.raw$G==2, 2:17])
19 if (0) {
20   # for flexmirt
21   write.table(data.g1, "cai2009-g1.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)
22   write.table(data.g2, "cai2009-g2.csv", quote=FALSE, row.names=FALSE, col.names=FALSE)  
23   fm <- read.flexmirt("~/irt/cai2009/cai2009-prm.txt")
24   fm$G1$spec <- NULL
25   fm$G2$spec <- NULL
26 }
27
28 fm <- structure(list(G1 = structure(list(param = structure(c(0.992675,  0.646717, 0, 0, 0.876469, 1.41764, 1.25402, 0, 0, 0.0826927,  1.76547, 1.20309, 0, 0, -0.346706, 2.1951, 0.844399, 0, 0, -0.978301,  1.37774, 0, 1.06694, 0, 0.992373, 1.80365, 0, 0.814109, 0, 0.213559,  2.15718, 0, 1.58086, 0, -0.418129, 1.18201, 0, 1.56533, 0, -1.24173,  1.80474, 0, 0, 1.0774, 0.810718, 2.60754, 0, 0, 1.23507, 0.0598008,  1.01874, 0, 0, 0.724402, -0.294029, 1.68916, 0, 0, 1.37546, -1.13333 ), .Dim = c(5L, 12L), .Dimnames = list(NULL, c("i1", "i2", "i3",  "i4", "i5", "i6", "i7", "i8", "i9", "i10", "i11", "i12"))), mean = structure(c(0.822622,  -0.290462, 0.19672, 0.733993), .Names = c("X6", "X7", "X8", "X9" )), cov = structure(c(0.826046, 0, 0, 0, 0, 1.656, 0, 0, 0, 0,  1.11263, 0, 0, 0, 0, 1.07878), .Dim = c(4L, 4L))), .Names = c("param",  "mean", "cov")),
29                      G2 = structure(list(param = structure(c(0.992675,  0.646717, 0, 0, 0, 0.876469, 1.41764, 1.25402, 0, 0, 0, 0.0826927,  1.76547, 1.20309, 0, 0, 0, -0.346706, 2.1951, 0.844399, 0, 0,  0, -0.978301, 1.37774, 0, 1.06694, 0, 0, 0.992373, 1.80365, 0,  0.814109, 0, 0, 0.213559, 2.15718, 0, 1.58086, 0, 0, -0.418129,  1.18201, 0, 1.56533, 0, 0, -1.24173, 1.80474, 0, 0, 1.0774, 0,  0.810718, 2.60754, 0, 0, 1.23507, 0, 0.0598008, 1.01874, 0, 0,  0.724402, 0, -0.294029, 1.68916, 0, 0, 1.37546, 0, -1.13333,  1.75531, 0, 0, 0, 1.20652, 0.875564, 1.26308, 0, 0, 0, 1.25013,  0.196607, 1.44526, 0, 0, 0, 0.990354, -0.351181, 1.89461, 0,  0, 0, 0.85611, -1.09382), .Dim = c(6L, 16L), .Dimnames = list(     NULL, c("i1", "i2", "i3", "i4", "i5", "i6", "i7", "i8", "i9",      "i10", "i11", "i12", "i13", "i14", "i15", "i16"))), mean = structure(c(0,  0, 0, 0, 0), .Names = c("X6", "X7", "X8", "X9", "X10")), cov = structure(c(1,  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,  0, 0, 1), .Dim = c(5L, 5L))), .Names = c("param", "mean", "cov" ))), .Names = c("G1", "G2"))
30
31 for (col in colnames(data.g1)) data.g1[[col]] <- ordered(data.g1[[col]])
32 for (col in colnames(data.g2)) data.g2[[col]] <- ordered(data.g2[[col]])
33
34 mk.model <- function(model.name, data, latent.free) {
35 #   name <- "g1"
36 #   data <- data.g1
37 #   latent.free <- TRUE
38 #   
39   numItems <- dim(data)[2]
40   numPersons <- dim(data)[1]
41   spec <- list()
42   spec[1:numItems] <- rpf.grm(factors = 2)
43   
44   dims <- (1 + numItems/4)
45   design <- matrix(c(rep(1,numItems),
46                      kronecker(2:dims,rep(1,4))), byrow=TRUE, ncol=numItems)
47   
48   ispec <- mxMatrix(name="ItemSpec", nrow=3, ncol=numItems,
49            values=sapply(spec, function(m) slot(m,'spec')),
50            free=FALSE)
51
52   ip.mat <- mxMatrix(name="ItemParam", nrow=3, ncol=numItems,
53                      values=c(1.4,1,0),
54                      free=c(TRUE,TRUE,TRUE),
55                      lbound=c(1e-5, 1e-5, NA))
56   ip.mat@free.group <- 'param'
57   
58   for (ix in 1:numItems) {
59     for (px in 1:3) {
60       name <- paste(c('p',ix,',',px), collapse='')
61       ip.mat@labels[px,ix] <- name
62     }
63   }
64   eip.mat <- mxAlgebra(ItemParam, name="EItemParam", fixed=TRUE)
65
66   m.mat <- mxMatrix(name="mean", nrow=1, ncol=dims, values=0, free=latent.free)
67   cov.mat <- mxMatrix(name="cov", nrow=dims, ncol=dims, values=diag(dims),
68                       free=latent.free)
69   
70   m1 <- mxModel(model=model.name, ip.mat, eip.mat, ispec, m.mat, cov.mat,
71                 mxMatrix(name="Design", nrow=dim(design)[1], ncol=numItems, values=design),
72                 mxData(observed=data, type="raw"),
73                 mxExpectationBA81(
74                   ItemSpec="ItemSpec",
75                   Design="Design",
76                   EItemParam="EItemParam",
77                   mean="mean", cov="cov",
78                   qpoints=21, qwidth=5),
79                 mxFitFunctionBA81(ItemParam="ItemParam", rescale=FALSE))
80   m1
81 }
82
83 groups <- paste("g", 1:2, sep="")
84
85 if (1) {
86   g1 <- mk.model("g1", data.g1, TRUE)
87   g2 <- mk.model("g2", data.g2, FALSE)
88   g1@matrices$ItemParam@values <-
89     rbind(fm$G1$param[1,], apply(fm$G1$param[2:4,], 2, sum), fm$G1$param[5,])
90   g1@matrices$mean@values <- t(fm$G1$mean)
91   g1@matrices$cov@values <- fm$G1$cov
92   g2@matrices$ItemParam@values <-
93     rbind(fm$G2$param[1,], apply(fm$G2$param[2:5,], 2, sum), fm$G2$param[6,])
94   
95   cModel <- mxModel(model="cModel", g1,g2,
96                     mxComputeOnce(paste(groups, 'expectation', sep='.'), context='E'))
97 #  cModel <- mxOption(cModel, "Number of Threads", 1)
98   for (grp in groups) cModel@submodels[[grp]]@expectation@scores <- 'full'
99   cModel.eap <- mxRun(cModel)
100
101   fm.sco.g1 <- suppressWarnings(try(read.table("models/nightly/data/cai2009-sco-g1.txt"), silent=TRUE))
102   if (is(fm.sco.g1, "try-error")) fm.sco.g1 <- read.table("data/cai2009-sco-g1.txt")
103   fm.sco.g2 <- suppressWarnings(try(read.table("models/nightly/data/cai2009-sco-g2.txt"), silent=TRUE))
104   if (is(fm.sco.g2, "try-error")) fm.sco.g2 <- read.table("data/cai2009-sco-g2.txt")
105   colnames(fm.sco.g1) <- c("grp","id",colnames(cModel.eap@submodels$g1@expectation@scores.out))
106   colnames(fm.sco.g2) <- c("grp","id",colnames(cModel.eap@submodels$g2@expectation@scores.out))
107   
108   scores.g1 <- cModel.eap@submodels$g1@expectation@scores.out
109   omxCheckCloseEnough(as.matrix(fm.sco.g1[,-1:-2]),
110                       scores.g1, 1e-3)
111   omxCheckCloseEnough(as.matrix(fm.sco.g2[,-1:-2]),
112                       cModel.eap@submodels$g2@expectation@scores.out, 1e-3)
113   
114     cModel <- mxModel(cModel,
115                       mxFitFunctionMultigroup(paste(groups, "fitfunction", sep=".")),
116                       mxComputeSequence(steps=list(
117                         mxComputeOnce(paste(groups, 'expectation', sep="."), context='M'),
118                         mxComputeOnce('fitfunction'))))
119     cModel <- mxRun(cModel)
120     omxCheckCloseEnough(cModel@output$minimum, correct.LL, .01)
121 }
122
123 if(1) {
124   g1 <- mk.model("g1", data.g1, TRUE)
125   g2 <- mk.model("g2", data.g2, FALSE)
126   grpModel <- mxModel(model="groupModel", g1, g2,
127                       mxFitFunctionMultigroup(paste(groups, "fitfunction", sep=".")),
128                       mxComputeIterate(steps=list(
129                         mxComputeOnce(paste(groups, "EItemParam", sep=".")),
130                         mxComputeOnce(paste(groups, 'expectation', sep='.'), context='E'),
131                         #                      mxComputeGradientDescent(free.group='param'),
132                         mxComputeNewtonRaphson(free.group='param'),
133                         mxComputeOnce(paste(groups, 'expectation', sep="."), context='M'),
134                         mxComputeOnce('fitfunction')
135                       ), verbose=TRUE))
136   
137   #grpModel <- mxOption(grpModel, "Number of Threads", 1)
138   
139   # NPSOL options:
140   grpModel <- mxOption(grpModel, "Analytic Gradients", 'Yes')
141   grpModel <- mxOption(grpModel, "Verify level", '-1')
142   grpModel <- mxOption(grpModel, "Function precision", '1.0E-7')
143   
144   grpModel <- mxRun(grpModel)
145     
146   omxCheckCloseEnough(grpModel@fitfunction@result, correct.LL, .01)
147   omxCheckCloseEnough(grpModel@submodels$g1@matrices$mean@values, t(fm$G1$mean), .01)
148   omxCheckCloseEnough(grpModel@submodels$g1@matrices$cov@values, fm$G1$cov, .01)
149   
150   #qplot(c(0, 3), stat="function", fun=function (x) dlnorm(x, sdlog=.25), geom="line") + ylim(0,1.5)
151   if (0) {
152     var.type <- apply(ip.mat@free,1,sum)
153     bias <- calc.bias(bank[sapply(bank, function (b) b$ghp == 40)])
154     df <- data.frame(x=t(correct)[t(ip.mat@free)],
155                      var=c(rep('a', var.type[1]),
156                            rep('b', var.type[2])),
157                      bias=t(bias)[t(ip.mat@free)])
158     df$var <- factor(df$var)
159     ggplot(df, aes(x, bias, color=var)) + geom_point() + xlab("true parameter value")
160   }
161 }