ifa: Snapshot
[openmx:openmx.git] / models / failing / irt-cai2009.R
1 options(error = utils::recover)
2 library(OpenMx)
3 library(rpf)
4 library(ggplot2)
5 library(stringr)
6
7 qwidth <- 5
8 quad <- imxEqualIntervalQuadratureData(21, qwidth)
9
10 data.raw <- read.csv("/opt/OpenMx/models/failing/data/cai-2009.csv")
11 data.g1 <- data.raw[data.raw$G==1, 2:13]
12 data.g2 <- data.raw[data.raw$G==2, 2:17]
13 colnames(data.g2) <- str_replace(colnames(data.g2), 'I', 'J')
14 data <- as.data.frame(cbind(data.g1, data.g2))
15 for (col in colnames(data)) data[[col]] <- ordered(data[[col]])
16
17 correct.g1 <- matrix(c(0.99, 1.42, 1.77, 2.19,  1.38, 1.8, 2.16, 1.18, 1.8, 2.61, 1.02, 1.69,
18                        0.65,1.26, 1.21, 0.85, 1.06, 0.81, 1.58, 1.56, 1.08, 1.24, 0.73, 1.39,
19                        0.88, 0.08, -0.35, -0.98, 0.99, 0.21, -0.42, -1.24, 0.81, 0.06, -0.29, -1.14,
20                        rep(0,12)), byrow=TRUE, nrow=4)
21 correct.g2 <- matrix(c(0.99, 1.42, 1.77, 2.19, 1.38, 1.8,  2.16, 1.18, 1.8,  2.61, 1.02, 1.69, 1.76, 1.26, 1.45, 1.9,
22                        0.65, 1.26, 1.21, 0.85, 1.06, 0.81, 1.58, 1.56, 1.08, 1.24, 0.73, 1.39, 1.21, 1.25, 0.99, 0.85,
23                        0.88, 0.08,-0.35,-0.98, 0.99, 0.21,-0.42,-1.24, 0.81, 0.06,-0.29,-1.14, 0.88, 0.2, -0.35,-1.09,
24                        rep(0,16)), byrow=TRUE, nrow=4)
25
26
27 numItems <- dim(data)[2]
28 numPersons <- dim(data)[1]
29
30 spec <- list()
31 spec[1:numItems] <- rpf.drm(dimensions = 2)
32
33 fit.g1 <- rpf.ot2000.chisq(spec[1:12], correct.g1, correct.g1!=0, data[1:12], quad)
34 fit.g2 <- rpf.ot2000.chisq(spec[13:numItems], correct.g2, correct.g2!=0, data[13:numItems], quad)
35 print(sum(sapply(fit.g1, function (f) f$statistic), sapply(fit.g2, function (f) f$statistic)))
36
37 design <- matrix(c(rep(c(1,2),4),
38                    rep(c(1,3),4),
39                    rep(c(1,4),4),
40                    rep(c(1,5),4),
41                    rep(c(1,6),4),
42                    rep(c(1,7),4),
43                    rep(c(1,8),4)), ncol=12+16)
44
45 ip.mat <- mxMatrix(name="ItemParam", nrow=4, ncol=numItems,
46                    values=c(1.4,1,0,0),
47                    free=c(TRUE,TRUE,TRUE,FALSE),
48                    lbound=c(1e-5, 1e-5, -qwidth, 0),
49                    ubound=c(10, 10, qwidth, 0))
50
51 for (ix in 1:12) {
52   for (px in 1:3) {
53     name <- paste(c('p',px,',',ix), collapse='')
54     ip.mat@labels[px,ix] <- name
55     ip.mat@labels[px,12+ix] <- name
56   }
57 }
58
59 fit1 <- function() {  
60   m1 <- mxModel(model="drm1", ip.mat,
61               mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
62                        values=sapply(spec, function(m) slot(m,'spec')),
63                        free=FALSE, byrow=TRUE),
64               mxMatrix(name="Design", nrow=dim(design)[1], ncol=numItems, values=design),
65               mxData(observed=data, type="raw"),
66               mxExpectationBA81(
67                 ItemSpec="ItemSpec",
68                 ItemParam="ItemParam",
69                 Design="Design",
70                         quadrature=quad),
71               mxFitFunctionBA81())
72
73   if (1) {
74     m1 <- mxOption(m1, "Analytic Gradients", 'yes')
75     m1 <- mxOption(m1, "Verify level", '-1')
76   } else {
77     m1 <- mxOption(m1, "Analytic Gradients", 'no')
78   }
79   m1 <- mxOption(m1, "Function precision", '1.0E-5')
80   m1 <- mxOption(m1, "Calculate Hessian", "No")
81   m1 <- mxOption(m1, "Standard Errors", "No")
82   m1 <- mxRun(m1)
83   m1
84 }
85
86 calc.bias <- function (bank) {
87   bias <- matrix(0, nrow=dim(correct)[1], ncol=dim(correct)[2])
88   for (sx in 1:length(bank)) {
89     bias <- bias + bank[[sx]]$param
90   }
91   bias <- (bias / length(bank)) - correct
92   bias
93 }
94
95 rda <- "cai2009-fit.rda"
96 if(0) {
97   fit <- fit1()
98   save(fit, file=rda)
99 } else {
100   load(rda)
101 }
102
103
104 openmx.itemparam <- fit@matrices$ItemParam@values
105 openmx.g1 <- openmx.itemparam[,1:12]
106 openmx.g2 <- openmx.itemparam[,13:28]
107
108 fit.g1 <- rpf.ot2000.chisq(spec[1:12], openmx.g1, openmx.g1!=0, data[1:12], quad)
109 fit.g2 <- rpf.ot2000.chisq(spec[13:numItems], openmx.g2, openmx.g2!=0, data[13:numItems], quad)
110 print(sum(sapply(fit.g1, function (f) f$statistic), sapply(fit.g2, function (f) f$statistic)))
111 # 881
112
113 #abs(calc.bias(bank[filter])[ip.mat@free])
114
115 if(0) {
116   bygh <- bank[sapply(bank, function (b) b$seed == 1)]
117   df <- as.data.frame(t(sapply(bygh, function(b) c(points=b$ghp, LL=b$LL))))
118   ggplot(df, aes(x=points, y=LL)) + geom_line()
119 }
120
121 #qplot(c(0, 3), stat="function", fun=function (x) dlnorm(x, sdlog=.25), geom="line") + ylim(0,1.5)
122 if (0) {
123   var.type <- apply(ip.mat@free,1,sum)
124   bias <- calc.bias(bank[sapply(bank, function (b) b$ghp == 40)])
125   df <- data.frame(x=t(correct)[t(ip.mat@free)],
126                    var=c(rep('a', var.type[1]),
127                          rep('b', var.type[2])),
128                    bias=t(bias)[t(ip.mat@free)])
129   df$var <- factor(df$var)
130   ggplot(df, aes(x, bias, color=var)) + geom_point() + xlab("true parameter value")
131 }