ifa: Add options for EAP scores
[openmx:openmx.git] / R / MxExpectationBA81.R
1 #
2 #   Copyright 2012 The OpenMx Project
3 #
4 #   Licensed under the Apache License, Version 2.0 (the "License");
5 #   you may not use this file except in compliance with the License.
6 #   You may obtain a copy of the License at
7
8 #        http://www.apache.org/licenses/LICENSE-2.0
9
10 #   Unless required by applicable law or agreed to in writing, software
11 #   distributed under the License is distributed on an "AS IS" BASIS,
12 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 #   See the License for the specific language governing permissions and
14 #   limitations under the License.
15
16
17 setClass(Class = "MxExpectationBA81",
18          representation = representation(
19            ItemSpec = "MxCharOrNumber",
20            Design = "MxOptionalCharOrNumber",
21            ItemParam = "MxCharOrNumber",
22            EItemParam = "MxOptionalCharOrNumber",
23            CustomPrior = "MxOptionalCharOrNumber",
24            qpoints = "numeric",
25            qwidth = "numeric",
26            cache = "logical",
27            rescale = "logical",
28            scores = "character",
29            free.mean = "MxOptionalLogical",
30            free.cov = "MxOptionalMatrix"),
31          contains = "MxBaseExpectation")
32
33 setMethod("initialize", "MxExpectationBA81",
34           function(.Object, ItemSpec, ItemParam, EItemParam, CustomPrior, Design,
35                    qpoints, qwidth, cache, free.mean, free.cov, rescale, scores,
36                    name = 'expectation') {
37             .Object@name <- name
38             .Object@ItemSpec <- ItemSpec
39             .Object@Design <- Design
40             .Object@ItemParam <- ItemParam
41             .Object@EItemParam <- EItemParam
42             .Object@qpoints <- qpoints
43             .Object@qwidth <- qwidth
44             .Object@cache <- cache
45             .Object@rescale <- rescale
46             .Object@scores <- scores
47             .Object@CustomPrior <- CustomPrior
48             .Object@data <- as.integer(NA)
49             .Object@free.mean <- free.mean
50             .Object@free.cov <- free.cov
51             return(.Object)
52           }
53 )
54
55 setMethod("genericExpDependencies", signature("MxExpectationBA81"),
56           function(.Object, dependencies) {
57                   sources <- c(.Object@CustomPrior, .Object@ItemParam,
58                                .Object@EItemParam,
59                                .Object@ItemSpec, .Object@Design)
60                   dependencies <- imxAddDependency(sources, .Object@name, dependencies)
61                   return(dependencies)
62           })
63
64 setMethod("genericExpFunConvert", signature("MxExpectationBA81"), 
65           function(.Object, flatModel, model, labelsData, defVars, dependencies) {
66                   modelname <- imxReverseIdentifier(model, .Object@name)[[1]]
67                   if(is.na(.Object@data)) {
68                           msg <- paste(typeof(.Object),
69                                        "does not have a dataset associated with it in model",
70                                        omxQuotes(modelname))
71                           stop(msg, call.=FALSE)
72                   }
73                   name <- .Object@name
74                   for (s in c("data", "ItemParam", "EItemParam", "CustomPrior", "ItemSpec", "Design")) {
75                           if (is.null(slot(.Object, s))) next;
76                           slot(.Object, s) <-
77                             imxLocateIndex(flatModel, slot(.Object, s), name)
78                   }
79
80                                         # How to get the data object?
81                   ## if (.Object@data@type != 'raw') {
82                   ##   stop(paste(typeof(.Object), "only supports raw data"));
83                   ## }
84                   return(.Object)
85           })
86
87 setMethod("genericExpFunNamespace", signature("MxExpectationBA81"), 
88         function(.Object, modelname, namespace) {
89                 .Object@name <- imxIdentifier(modelname, .Object@name)
90                 for (s in c("ItemParam", "EItemParam", "CustomPrior", "ItemSpec", "Design")) {
91                         if (is.null(slot(.Object, s))) next;
92                         slot(.Object, s) <-
93                           imxConvertIdentifier(slot(.Object, s), modelname, namespace)
94                 }
95                 return(.Object)
96 })
97
98 setMethod("genericExpRename", signature("MxExpectationBA81"),
99         function(.Object, oldname, newname) {
100           # not sure what goes here yet
101                 return(.Object)
102 })
103
104 ##' Create a Bock & Aitkin (1981) expectation
105 ##'
106 ##' The standard Normal distribution of the quadrature acts like a
107 ##' prior distribution for difficulty. It is not necessary to impose
108 ##' any additional Bayesian prior on difficulty estimates (Baker &
109 ##' Kim, 2004, p. 196).
110 ##' 
111 ##' @param ItemSpec one column for each item containing the model ID
112 ##' (see \code{\link{mxLookupIRTItemModelID}}), numDimensions, and
113 ##' numOutcomes (3 rows)
114 ##' @param ItemParam one column for each item with parameters starting
115 ##' at row 1 and extra rows filled with NA
116 ##' @param Design one column per item assignment of person abilities
117 ##' to item dimensions (optional)
118 ##' @param qpoints number of points to use for rectangular quadrature integrations (default 49)
119 ##' See Seong (1990) for some considerations on specifying this parameter.
120 ##' @param cache whether to cache part of the expectation calculation
121 ##' (enabled by default).
122 ##' @references
123 ##' Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item
124 ##' parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.
125 ##'
126 ##' Cai, L. (2010). A two-tier full-information item factor analysis
127 ##' model with applications. Psychometrika, 75, 581-612.
128 ##'
129 ##' Seong, T. J. (1990). Sensitivity of marginal maximum likelihood
130 ##' estimation of item and ability parameters to the characteristics
131 ##' of the prior ability distributions. Applied Psychological
132 ##' Measurement, 14(3), 299-311.
133
134 mxExpectationBA81 <- function(ItemSpec, ItemParam, EItemParam=NULL, CustomPrior=NULL, Design=NULL,
135                               qpoints=NULL, qwidth=6.0, cache=TRUE, free.mean=NULL, free.cov=NULL,
136                               rescale=TRUE, scores="omit") {
137
138         if (missing(qpoints)) qpoints <- 49
139         if (qpoints < 3) {
140                 stop("qpoints should be 3 or greater")
141         }
142         if (missing(qwidth)) qwidth <- 6
143         if (qwidth <= 0) {
144                 stop("qwidth must be positive")
145         }
146   
147         if (!is.null(free.cov) && any(free.cov != t(free.cov))) {
148                 stop("free.cov must be symmetric")
149         }
150
151         score.options <- c("omit", "unique", "full")
152         if (!match(scores, score.options)) {
153                 stop(paste("Valid score options are", deparse(score.options)))
154         }
155
156         return(new("MxExpectationBA81", ItemSpec, ItemParam, EItemParam, CustomPrior, Design,
157                    qpoints, qwidth, cache, free.mean, free.cov, rescale, scores))
158 }
159
160 imxEqualIntervalQuadratureData <- function(n, width) {
161   x <- seq(-width, width, length.out=n)
162   list(x=x, w=dnorm(x)/sum(dnorm(x)))
163 }