ifa: Gradients
[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 setClassUnion("MxOptionalFunction", c("NULL", "function"))
18
19 setClass(Class = "MxExpectationBA81",
20          representation = representation(
21            ItemSpec = "MxCharOrNumber",
22            Design = "MxOptionalCharOrNumber",
23            ItemParam = "MxCharOrNumber",
24            RPF = "MxOptionalFunction",
25            ItemPrior = "MxCharOrNumber",
26            GHpoints = "numeric",
27            GHarea = "numeric",
28            cache = "logical"),
29          contains = "MxBaseExpectation")
30
31 setMethod("initialize", "MxExpectationBA81",
32           function(.Object, ItemSpec, ItemParam, ItemPrior, Design, RPF, gh, cache, name = 'expectation') {
33             .Object@name <- name
34             .Object@ItemSpec <- ItemSpec
35             .Object@Design <- Design
36             .Object@ItemParam <- ItemParam
37             .Object@RPF <- RPF
38             .Object@GHpoints <- gh[[1]]
39             .Object@GHarea <- gh[[2]]
40             .Object@cache <- cache
41             .Object@ItemPrior <- ItemPrior
42             .Object@data <- as.integer(NA)
43             return(.Object)
44           }
45 )
46
47 setMethod("genericExpDependencies", signature("MxExpectationBA81"),
48           function(.Object, dependencies) {
49                   sources <- c(.Object@ItemPrior, .Object@ItemParam,
50                                .Object@ItemSpec, .Object@Design)
51                   dependencies <- imxAddDependency(sources, .Object@name, dependencies)
52                   return(dependencies)
53           })
54
55 setMethod("genericExpFunConvert", signature("MxExpectationBA81"), 
56           function(.Object, flatModel, model, labelsData, defVars, dependencies) {
57                   modelname <- imxReverseIdentifier(model, .Object@name)[[1]]
58                   if(is.na(.Object@data)) {
59                           msg <- paste(typeof(.Object),
60                                        "does not have a dataset associated with it in model",
61                                        omxQuotes(modelname))
62                           stop(msg, call.=FALSE)
63                   }
64                   name <- .Object@name
65                   for (s in c("data", "ItemParam", "ItemPrior", "ItemSpec", "Design")) {
66                           if (is.null(slot(.Object, s))) next;
67                           slot(.Object, s) <-
68                             imxLocateIndex(flatModel, slot(.Object, s), name)
69                   }
70
71                                         # How to get the data object?
72                   ## if (.Object@data@type != 'raw') {
73                   ##   error(paste(typeof(.Object), "only supports raw data"));
74                   ## }
75                   return(.Object)
76           })
77
78 setMethod("genericExpFunNamespace", signature("MxExpectationBA81"), 
79         function(.Object, modelname, namespace) {
80                 .Object@name <- imxIdentifier(modelname, .Object@name)
81                 for (s in c("ItemParam", "ItemPrior", "ItemSpec", "Design")) {
82                         if (is.null(slot(.Object, s))) next;
83                         slot(.Object, s) <-
84                           imxConvertIdentifier(slot(.Object, s), modelname, namespace)
85                 }
86                 return(.Object)
87 })
88
89 setMethod("genericExpRename", signature("MxExpectationBA81"),
90         function(.Object, oldname, newname) {
91           # not sure what goes here yet
92                 return(.Object)
93 })
94
95 ##' Create a Bock & Aitkin (1981) expectation
96 ##'
97 ##' The RPF function is an optional argument for experimenting with
98 ##' novel response probability functions. The model ID in the ItemSpec
99 ##' is ignored if you provide your own function. There is no provision
100 ##' to provide analytic gradients from R. If you want to provide
101 ##' analytic gradients then you need to implement them in C.
102 ##'
103 ##' The standard Normal distribution of the quadrature acts like a
104 ##' prior distribution for difficulty. It is not necessary to impose
105 ##' any additional Bayesian prior on difficulty estimates (Baker &
106 ##' Kim, 2004, p. 196).
107 ##' 
108 ##' The cache more than halves CPU time, especially if you have many
109 ##' specific dimensions. With the cache enabled, memory usage is about
110 ##' 8 * (numItems + numUnique * (2 + numThreads * (numSpecific+1)) +
111 ##' numUnique * product(quadrature grids) * numSpecific). With the
112 ##' cache disabled, memory use is limited to about 8 * (numItems + 2 *
113 ##' numUnique + (numUnique+1) * (numThreads * (numSpecific+1)) bytes.
114 ##' For example, suppose your model has 20 items, 500 unique response
115 ##' patterns, 2 primary dimensions, 3 specific dimensions, a 40 point
116 ##' quadrature and your computer offers 8 way multiprocessing. With
117 ##' the cache enabled, 8 * (20 + 500 * (2+8 * (3+1)) + 500 * 40^3 * 3)
118 ##' = 733 MiB. With the cache disabled, 8 * (20 + 2 * 500 + (500 + 1) *
119 ##' (8 * (3 + 1))) = 134 KiB.
120 ##'
121 ##' @param ItemSpec one column for each item containing the model ID
122 ##' (see \code{\link{mxLookupIRTItemModelID}}), numDimensions, and
123 ##' numOutcomes (3 rows)
124 ##' @param ItemParam one column for each item with parameters starting
125 ##' at row 1 and extra rows filled with NA
126 ##' @param Design one column per item assignment of person abilities
127 ##' to item dimensions (optional)
128 ##' @param RPF a function to turn item parameters into a response
129 ##' probability matrix (optional)
130 ##' @param GHpoints number of points to use for Gauss-Hermite quadrature integrations (default 10)
131 ##' See Seong (1990) for some considerations on specifying this parameter.
132 ##' @param cache whether to cache part of the expectation calculation
133 ##' (enabled by default).
134 ##' @references
135 ##' Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item
136 ##' parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.
137 ##'
138 ##' Cai, L. (2010). A two-tier full-information item factor analysis
139 ##' model with applications. Psychometrika, 75, 581-612.
140 ##'
141 ##' Seong, T. J. (1990). Sensitivity of marginal maximum likelihood
142 ##' estimation of item and ability parameters to the characteristics
143 ##' of the prior ability distributions. Applied Psychological
144 ##' Measurement, 14(3), 299-311.
145
146 mxExpectationBA81 <- function(ItemSpec, ItemParam, ItemPrior, Design=NULL, RPF=NULL, GHpoints=10, cache=TRUE) {
147   if (GHpoints < 3) {
148     error("GHpoints should be 3 or greater")
149   }
150   gh <- imxGaussHermiteData(GHpoints)
151   
152   return(new("MxExpectationBA81", ItemSpec, ItemParam, ItemPrior, Design, RPF, gh, cache))
153 }
154
155 # Not used, just informational
156 imxAdhocGaussHermiteData <- function(n, width) {
157   x <- seq(-width, width, length.out=n)
158   list(x=x, w=dnorm(x)/sum(dnorm(x)))
159 }
160
161 imxGaussHermiteData <- function(n) .Call(omxGaussHermiteData, n)
162
163 ##' Convert an IRT item model name to an ID
164 ##'
165 ##' drm1 is the standard 3PL. drm is the multidimensional version of
166 ##' the 3PL. gpcm1 is Generalized Partial Credit Model.
167 ##'
168 ##' @param name name of the item model (string)
169 ##' @return the integer ID assigned to the given model
170 mxLookupIRTItemModelID <- function(name) {
171         models <- .Call(omx_get_rpf_names)
172         match(name, models) - 1;
173 }