ifa: Snapshot
[openmx:openmx.git] / models / failing / irt-cai2010-1.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 = utils::recover)
6 library(ggplot2)
7 library(OpenMx)
8 library(rpf)
9 library(mvtnorm)
10
11 numPersons <- 500
12
13 items.3d <- list()
14 items.3d[1:20] <- rpf.drm(dimensions=3)
15 #lapply(items, function (i) i@spec[4] <- 100)  # flatten prior on slope
16 numItems <- length(items.3d)
17 maxParam.3d <- max(vapply(items.3d, function(i) i@numParam, 0))
18
19 if(0) {
20   correct <- sapply(items.3d, rpf.rparam)
21   correct[5,] <- 0
22   correct[1,11:20] <- 0
23   correct[2,1:10] <- 0
24   correct[3,c(3:5,8:10,13:15,18:20)] <- 0
25   for (d in c(1,6,11,16)) correct[2,d+1] <- correct[2,d]
26 } else {
27   # Cai (2010, p. 592)
28   correct <- matrix(c(1.7,1.5,2.06,2.56,1.96,2.82,1,2.6,1.15,2.76,rep(0,10),
29                       rep(0,10),2.62,2.3,1.94,1.43,1.95,1.67,2.67,1.98,2.22,2.19,
30                       1.63,1.64,rep(0,3),1.56,1.56,rep(0,3),1.89,1.89,rep(0,3),2,2,rep(0,3),
31                       .91,.61,.02,.04,.02,.21,-.4,-.24,-.74,.3,1.85,
32                       1.12,.58,-1.26,.07,1.05,2.05,.16,-.44,.3,
33                       rep(0,20)), byrow=TRUE, ncol=numItems)
34 }
35
36 design.3d <- matrix(c(rep(1,10),rep(NA,10),
37                       rep(NA,10),rep(2,10),
38                       3,3,rep(NA,3),4,4,rep(NA,3),5,5,rep(NA,3),6,6,rep(NA,3)),
39                     byrow=TRUE, nrow=3)
40
41 cov <- matrix(c(1, .68, .68, 1), nrow=2)
42
43 ip.labels <- matrix(NA, nrow=maxParam.3d, ncol=numItems)
44 ip.labels[3,1:2] <- "s1"
45 ip.labels[3,6:7] <- "s2"
46 ip.labels[3,11:12] <- "s3"
47 ip.labels[3,16:17] <- "s4"
48
49 ip.free <- correct != 0
50 ip.mat <- mxMatrix(name="ItemParam", nrow=maxParam.3d, ncol=numItems,
51                    values=c(1.414, 1.414, 1, 0, 0),
52                    free=ip.free, labels=ip.labels,
53                    lbound=c(1*10^-9, 1*10^-9, 1*10^-9, -10^9,0))
54 ip.mat@values[!ip.free] <- 0
55
56 fit1 <- function(seed=5, ghp=11) {
57   result <- list(seed=seed, ghp=ghp)
58   
59   set.seed(seed)
60   ability <- cbind(rmvnorm(numPersons, sigma=cov),
61                    matrix(rnorm(numPersons * 4), ncol=4))
62   design.full <- design.3d
63   design.full[is.na(design.3d)] <- 1  # ignored because of item parameters
64   data.3d <- rpf.sample(ability, items.3d, correct, design.full)
65   
66   m1 <- mxModel(model="cai1",
67                 mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
68                          values=sapply(items.3d, function(m) slot(m,'spec')),
69                          free=FALSE, byrow=TRUE),
70                 mxMatrix(name="Design", nrow=dim(design.3d)[1], ncol=numItems, values=design.3d),
71                 ip.mat,
72                 mxData(observed=data.3d, type="raw"),
73                 mxExpectationBA81(ItemSpec="ItemSpec", Design="Design", ItemParam="ItemParam",
74                                   GHpoints=ghp),
75                 mxFitFunctionBA81()
76   )
77   
78   m1 <- mxOption(m1, "Analytic Gradients", 'no')
79   if (1) {
80     m1 <- mxOption(m1, "Analytic Gradients", 'yes')
81     m1 <- mxOption(m1, "Verify level", '-1')
82   }
83   m1 <- mxOption(m1, "Function precision", '1.0E-5')
84   m1 <- mxOption(m1, "Calculate Hessian", "No")
85   m1 <- mxOption(m1, "Standard Errors", "No")
86   
87   m1 <- mxRun(m1, silent=TRUE)
88   
89   result$cpuTime <- m1@output$cpuTime
90   result$LL <- m1@output$Minus2LogLikelihood
91   result$param <- m1@matrices$ItemParam@values
92   result
93 }
94
95 calc.bias <- function (bank) {
96   bias <- matrix(0, nrow=dim(correct)[1], ncol=dim(correct)[2])
97   for (sx in 1:length(bank)) {
98     bias <- bias + bank[[sx]]$param
99   }
100   bias <- (bias / length(bank)) - correct
101   bias
102 }
103
104 if (0) {
105   bank <- list()
106 }
107 setwd("/opt/OpenMx")
108 rda <- "irt-cai2010.rda"
109 load(rda)
110 for (seed in 1:500) {
111   if (length(bank)) {
112     if (any(seed == sapply(bank, function (b) b$seed))) next;
113   }
114   bi <- length(bank) + 1
115   bank[[bi]] <- fit1(ghp=13, seed=seed)
116   save(bank, file=rda)
117
118   print(cor(bank[[bi]]$param[ip.mat@free], correct[ip.mat@free]))
119 }
120
121 calc.bias(bank[sapply(bank, function (b) b$ghp == 13)])
122
123 if (1) {
124   var.type <- apply(ip.mat@free,1,sum)
125   bias <- calc.bias(bank[sapply(bank, function (b) b$ghp == 13)])
126   df <- data.frame(x=t(correct)[t(ip.mat@free)],
127                    var=c(rep('a1', var.type[1]),
128                          rep('a2', var.type[2]),
129                          rep('a3', var.type[3]),
130                          rep('b', var.type[4])),
131                    bias=t(bias)[t(ip.mat@free)])
132   df$var <- factor(df$var)
133   ggplot(df, aes(x, bias, color=var)) + geom_point() + xlab("true parameter value")
134 }
135
136 if(0) {
137   bygh <- bank #bank[sapply(bank, function (b) b$seed == 1)]
138   df <- data.frame(time=sapply(bygh, function(b) c(as.numeric(b$cpuTime))))
139   ggplot(df, aes(time)) + geom_histogram(binwidth=5)
140 }