ifa: snapshot with irt-2d working
[openmx:openmx.git] / models / failing / irt-drm1-bias.R
1 #options(error = utils::recover)
2 library(OpenMx)
3 library(rpf)
4 library(ggplot2)
5
6 numItems <- 20
7 numPersons <- 500
8
9 i1 <- rpf.drm(numChoices=4, a.prior.sdlog=.5, multidimensional=TRUE)
10 #i1 <- rpf.drm(numChoices=4, a.prior.sdlog=.5)
11 items <- list()
12 items[1:numItems] <- i1
13
14 if (0) {
15   correct <- sapply(items, rpf.rparam)
16   correct[2,] <- scale(correct[2,])
17   correct[3,] <- 0
18   cat(deparse(correct))
19 } else {
20   correct <- structure(c(1.54250842425733, 1.29865606045677, 0, 0.582738899753006,  0.439074405155082, 0, 0.479997481961601, 1.6288972150727, 0,  1.69444870619279, 0.415803213429429, 0, 0.617838919096744, 0.0339637225500039,  0, 0.89104263140429, -0.310656987834084, 0, 0.856862341680302,  -1.39777061130952, 0, 0.399785394294146, 0.452539004701956, 0,  1.13750982818114, -1.12551515354692, 0, 0.551920087175635, 0.234962883698946,  0, 1.317228172309584, -0.767585189231164, 0, 1.06680521649664,  -0.17744848799824, 0, 0.820181809806619, 0.697978914209079, 0,  1.2446674744156, 0.579329622962913, 0, 3.17834154396065, 0.0422351712197786,  0, 1.29268733199003, -1.50421669038342, 0, 0.620813738032935,  -1.72835375813844, 0, 0.57191713624839, 0.748756461527565, 0,  1.44045131430693, 1.69557125401783, 0, 0.462792889649201, -0.956221050560265,  0), .Dim = c(3L, 20L), .Dimnames = list(c("a", "b", "c"), NULL))
21 }
22
23 ip.mat <- mxMatrix(name="ItemParam", nrow=3, ncol=numItems,
24                    values=c(1,0,0),
25                    free=c(TRUE,TRUE,FALSE),
26                    lbound=c(1e-5, -1e6, 0))
27
28 fit1 <- function(seed=5, ghp=11) {
29   result <- list(seed=seed, ghp=ghp, sdlog=i1@a.prior.sdlog, mdim=(class(i1) == "rpf.mdim.drm"))
30   
31   set.seed(seed)
32   ability <- rnorm(numPersons)
33   gen.param <- correct
34   if (class(i1) == "rpf.mdim.drm") {
35     gen.param[2,] <- gen.param[2,] * -gen.param[1,]
36   }
37   data <- rpf.sample(ability, items, gen.param)
38
39   m1 <- mxModel(model="drm1", ip.mat,
40               mxMatrix(name="ItemSpec", nrow=6, ncol=numItems,
41                        values=sapply(items, function(m) slot(m,'spec')),
42                        free=FALSE, byrow=TRUE),
43               mxData(observed=data, type="raw"),
44               mxExpectationBA81(
45                 ItemSpec="ItemSpec",
46                 ItemParam="ItemParam",
47                         GHpoints=ghp),
48               mxFitFunctionBA81())
49
50   if (1) {
51     result$gradient <- 1
52     m1 <- mxOption(m1, "Analytic Gradients", 'yes')
53     m1 <- mxOption(m1, "Verify level", '-1')
54   } else {
55     result$gradient <- 0
56     m1 <- mxOption(m1, "Analytic Gradients", 'no')
57   }
58   m1 <- mxOption(m1, "Function precision", '1.0E-5')
59   m1 <- mxOption(m1, "Calculate Hessian", "No")
60   m1 <- mxOption(m1, "Standard Errors", "No")
61   m1 <- mxRun(m1)
62
63   result$cpuTime <- m1@output$cpuTime
64   result$LL <- m1@output$Minus2LogLikelihood
65   result$param <- m1@matrices$ItemParam@values
66   if (class(i1) == "rpf.mdim.drm") {
67     result$param[2,] <- result$param[2,] / -result$param[1,]
68   }
69   result
70 }
71
72 calc.bias <- function (bank) {
73   bias <- matrix(0, nrow=dim(correct)[1], ncol=dim(correct)[2])
74   for (sx in 1:length(bank)) {
75     bias <- bias + bank[[sx]]$param
76   }
77   bias <- (bias / length(bank)) - correct
78   bias
79 }
80
81 if (0) {
82   bank <- list()
83 }
84 setwd("/opt/OpenMx")
85 rda <- "irt-drm1-bias.rda"
86 load(rda)
87 for (cnt in 1:500) {
88   if (length(bank)) {
89     filter <- sapply(bank, function (b) b$mdim==(class(i1) == "rpf.mdim.drm") & b$sdlog == .5 & b$ghp==17)
90     if (any(sapply(bank, function (b) b$seed==cnt) & filter)) next;
91   }
92   bi <- length(bank) + 1
93   bank[[bi]] <- fit1(ghp=17, seed=cnt)
94   save(bank, file=rda)
95   
96   print(cor(bank[[bi]]$param[ip.mat@free], correct[ip.mat@free]))
97 #  print(cor(c(calc.bias(bank[filter])[ip.mat@free]), c(correct[ip.mat@free])))
98 }
99
100 for (cnt in 7:50) {
101   filter <- sapply(bank, function (b) b$mdim==(class(i1) == "rpf.mdim.drm") & b$sdlog == .5)
102   if (any(sapply(bank, function (b) b$ghp==cnt) & filter)) next;
103   bi <- length(bank) + 1
104   bank[[bi]] <- fit1(ghp=cnt, seed=1)
105   save(bank, file=rda)
106   
107   print(cor(bank[[bi]]$param[ip.mat@free], correct[ip.mat@free]))
108   #  print(cor(c(calc.bias(bank[filter])[ip.mat@free]), c(correct[ip.mat@free])))
109 }
110
111 #abs(calc.bias(bank[filter])[ip.mat@free])
112
113 if(1) {
114   df <- list()
115   for (mdim in c(TRUE,FALSE)) {
116     x <- correct[ip.mat@free]
117     bias <- calc.bias(bank[sapply(bank, function (b) b$mdim == mdim & b$ghp==17)])[ip.mat@free]
118     df <- rbind(df, data.frame(mdim=mdim, bias=bias, x=x))
119   }
120   df$mdim <- factor(df$mdim)
121   ggplot(df, aes(x=x, y=bias, color=mdim)) + geom_point()
122 }
123
124 if(1) {
125   bygh <- bank[sapply(bank, function (b) b$seed == 1)]
126   df <- as.data.frame(t(sapply(bygh, function(b) c(points=b$ghp, LL=b$LL, mdim=b$mdim))))
127   df$mdim <- factor(df$mdim)
128   ggplot(df, aes(x=points, y=LL, color=mdim)) + geom_line()
129 }
130
131 if(1) {
132   bygh <- bank[sapply(bank, function (b) b$seed == 1)]
133   df <- as.data.frame(t(sapply(bygh, function(b) c(mdim=b$mdim, points=b$ghp,
134                                                    time=as.numeric(b$cpuTime)))))
135   df$mdim <- factor(df$mdim)
136   ggplot(df, aes(x=points, y=time, color=mdim)) + geom_line()
137 }
138
139 #qplot(c(0, 3), stat="function", fun=function (x) dlnorm(x, sdlog=.25), geom="line") + ylim(0,1.5)
140 if (0) {
141   plot(correct[ip.mat@free],
142        calc.bias(bank[sapply(bank, function (b) b$mdim==TRUE & b$sdlog==.5)])[ip.mat@free])
143 }