Pass item models in a list instead of an MxMatrix
[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 = "list",
20            Design = "MxOptionalCharOrNumber",
21            EItemParam = "MxOptionalCharOrNumber",
22            qpoints = "numeric",
23            qwidth = "numeric",
24            cache = "logical",
25            scores = "character",
26            scores.out = "matrix",
27            mean = "MxCharOrNumber",
28            cov = "MxCharOrNumber"),
29          contains = "MxBaseExpectation")
30
31 setMethod("initialize", "MxExpectationBA81",
32           function(.Object, ItemSpec, EItemParam, Design,
33                    qpoints, qwidth, cache, mean, cov, scores,
34                    name = 'expectation') {
35             .Object@name <- name
36             .Object@ItemSpec <- ItemSpec
37             .Object@Design <- Design
38             .Object@EItemParam <- EItemParam
39             .Object@qpoints <- qpoints
40             .Object@qwidth <- qwidth
41             .Object@cache <- cache
42             .Object@scores <- scores
43             .Object@data <- as.integer(NA)
44             .Object@mean <- mean
45             .Object@cov <- cov
46             .Object@scores.out <- matrix()
47             return(.Object)
48           }
49 )
50
51 setMethod("genericExpDependencies", signature("MxExpectationBA81"),
52           function(.Object, dependencies) {
53                   sources <- c(.Object@EItemParam,
54                                .Object@mean, .Object@cov,
55                                .Object@Design)
56                   dependencies <- imxAddDependency(sources, .Object@name, dependencies)
57                   return(dependencies)
58           })
59
60 setMethod("genericExpFunConvert", signature("MxExpectationBA81"), 
61           function(.Object, flatModel, model, labelsData, defVars, dependencies) {
62                   modelname <- imxReverseIdentifier(model, .Object@name)[[1]]
63                   if(is.na(.Object@data)) {
64                           msg <- paste(typeof(.Object),
65                                        "does not have a dataset associated with it in model",
66                                        omxQuotes(modelname))
67                           stop(msg, call.=FALSE)
68                   }
69                   name <- .Object@name
70                   for (s in c("data", "EItemParam",
71                               "Design", "mean", "cov")) {
72                           if (is.null(slot(.Object, s))) next;
73                           slot(.Object, s) <-
74                             imxLocateIndex(flatModel, slot(.Object, s), name)
75                   }
76
77                                         # How to get the data object?
78                   ## if (.Object@data@type != 'raw') {
79                   ##   stop(paste(typeof(.Object), "only supports raw data"));
80                   ## }
81                   return(.Object)
82           })
83
84 setMethod("qualifyNames", signature("MxExpectationBA81"), 
85         function(.Object, modelname, namespace) {
86                 .Object@name <- imxIdentifier(modelname, .Object@name)
87                 for (s in c("EItemParam",
88                             "Design", "mean", "cov")) {
89                         if (is.null(slot(.Object, s))) next;
90                         slot(.Object, s) <-
91                           imxConvertIdentifier(slot(.Object, s), modelname, namespace)
92                 }
93                 return(.Object)
94 })
95
96 setMethod("genericExpRename", signature("MxExpectationBA81"),
97         function(.Object, oldname, newname) {
98           # not sure what goes here yet
99                 return(.Object)
100 })
101
102 ##' Create a Bock & Aitkin (1981) expectation
103 ##'
104 ##' The standard Normal distribution of the quadrature acts like a
105 ##' prior distribution for difficulty. It is not necessary to impose
106 ##' any additional Bayesian prior on difficulty estimates (Baker &
107 ##' Kim, 2004, p. 196).
108 ##' 
109 ##' @param ItemParam one column for each item with parameters starting
110 ##' at row 1 and extra rows filled with NA
111 ##' @param Design one column per item assignment of person abilities
112 ##' to item dimensions (optional)
113 ##' @param qpoints number of points to use for rectangular quadrature integrations (default 49)
114 ##' See Seong (1990) for some considerations on specifying this parameter.
115 ##' @param cache whether to cache part of the expectation calculation
116 ##' (enabled by default).
117 ##' @references
118 ##' Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item
119 ##' parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.
120 ##'
121 ##' Cai, L. (2010). A two-tier full-information item factor analysis
122 ##' model with applications. Psychometrika, 75, 581-612.
123 ##'
124 ##' Seong, T. J. (1990). Sensitivity of marginal maximum likelihood
125 ##' estimation of item and ability parameters to the characteristics
126 ##' of the prior ability distributions. Applied Psychological
127 ##' Measurement, 14(3), 299-311.
128
129 mxExpectationBA81 <- function(ItemSpec, EItemParam, Design=NULL,
130                               qpoints=NULL, qwidth=6.0, cache=TRUE, mean=NULL, cov=NULL,
131                               scores="omit") {
132
133         if (missing(qpoints)) qpoints <- 49
134         if (qpoints < 3) {
135                 stop("qpoints should be 3 or greater")
136         }
137         if (qpoints %% 2 == 0) {
138                 warning("An even number of qpoints can obtain a better than true fit; Use an odd number of qpoints")
139         }
140         if (missing(qwidth)) qwidth <- 6
141         if (qwidth <= 0) {
142                 stop("qwidth must be positive")
143         }
144   
145         score.options <- c("omit", "unique", "full")
146         if (!match(scores, score.options)) {
147                 stop(paste("Valid score options are", deparse(score.options)))
148         }
149
150         return(new("MxExpectationBA81", ItemSpec, EItemParam, Design,
151                    qpoints, qwidth, cache, mean, cov, scores))
152 }