mxSimplify2Array doc
[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            item = "MxCharOrNumber",
21            EstepItem = "MxOptionalMatrix",
22            qpoints = "numeric",
23            qwidth = "numeric",
24            mean = "MxOptionalCharOrNumber",
25            cov = "MxOptionalCharOrNumber",
26              debugInternal="logical",
27              dataColumns="integer",
28            dims = "character",
29            numStats = "numeric",
30            verbose = "integer",
31              output = "list",
32              debug = "list",
33              naAction = "character",
34              minItemsPerScore = "integer",
35              weightColumn = "MxCharOrNumber"),
36          contains = "MxBaseExpectation")
37
38 setMethod("initialize", "MxExpectationBA81",
39           function(.Object, ItemSpec, item, EstepItem,
40                    qpoints, qwidth, mean, cov, verbose, debugInternal,
41                    naAction, minItemsPerScore, weightColumn, name = 'expectation') {
42             .Object@name <- name
43             .Object@ItemSpec <- ItemSpec
44             .Object@item <- item
45             .Object@EstepItem <- EstepItem
46             .Object@qpoints <- qpoints
47             .Object@qwidth <- qwidth
48             .Object@data <- as.integer(NA)
49             .Object@mean <- mean
50             .Object@cov <- cov
51             .Object@verbose <- verbose
52             .Object@debugInternal <- debugInternal
53             .Object@naAction <- naAction
54             .Object@minItemsPerScore <- minItemsPerScore
55             .Object@weightColumn <- weightColumn
56             return(.Object)
57           }
58 )
59
60 setMethod("genericExpDependencies", signature("MxExpectationBA81"),
61           function(.Object, dependencies) {
62                   sources <- c(.Object@mean, .Object@cov, .Object@item)
63                   dependencies <- imxAddDependency(sources, .Object@name, dependencies)
64                   return(dependencies)
65           })
66
67 setMethod("genericExpFunConvert", signature("MxExpectationBA81"), 
68           function(.Object, flatModel, model, labelsData, defVars, dependencies) {
69                   modelname <- imxReverseIdentifier(model, .Object@name)[[1]]
70                   if(is.na(.Object@data)) {
71                           msg <- paste(typeof(.Object),
72                                        "does not have a dataset associated with it in model",
73                                        omxQuotes(modelname))
74                           stop(msg, call.=FALSE)
75                   }
76
77                   verifyMvnNames(.Object@cov, .Object@mean, "prior", flatModel, model@name, class(.Object))
78
79                   fnames <- dimnames(flatModel[[.Object@cov]])[[1]]
80
81                   item <- flatModel@matrices[[.Object@item]]
82                   if (is.null(rownames(item)) || any(rownames(item)[1:length(fnames)] != fnames)) {
83                           msg <- paste("The first", length(fnames), "rownames of item",
84                                        "must be", omxQuotes(fnames))
85                           stop(msg)
86                   }
87
88                   name <- .Object@name
89                   for (s in c("data", "item")) {
90                           if (is.null(slot(.Object, s))) next
91                           slot(.Object, s) <-
92                             imxLocateIndex(flatModel, slot(.Object, s), name)
93                   }
94                   for (s in c("mean", "cov")) {
95                           name <- slot(.Object, s)
96                           matrixNumber <- match(name, names(flatModel@matrices))
97                           if (is.na(matrixNumber)) {
98                                   slot(.Object, s) <- NULL
99                           } else {
100                                   slot(.Object, s) <- -matrixNumber
101                           }
102                   }
103
104                   .Object@dims <- colnames(flatModel@datasets[[.Object@data + 1]]@observed) # for summary
105
106                   if (is.null(colnames(item))) {
107                           stop(paste(class(.Object),
108                                      ": the item matrix must have column names",
109                                      "to match against the observed data column names."))
110                   } else {
111                           dc <- match(colnames(item), .Object@dims) - 1L
112                           if (any(is.na(dc))) {
113                                   msg <- paste("Some items are not found among the observed data columns:",
114                                                omxQuotes(colnames(item)[is.na(dc)]))
115                                   stop(msg, call.=FALSE)
116                           }
117                           .Object@dataColumns <- dc
118                           if (!is.na(.Object@weightColumn)) {
119                                   wc <- match(.Object@weightColumn, .Object@dims) - 1L
120                                   if (is.na(wc)) {
121                                           msg <- paste("Weight column",.Object@weightColumn,
122                                                        "not found in observed data columns")
123                                           stop(msg, call.=FALSE)
124                                   }
125                                   .Object@weightColumn <- wc
126                           }
127                   }
128                   return(.Object)
129           })
130
131 setMethod("qualifyNames", signature("MxExpectationBA81"), 
132         function(.Object, modelname, namespace) {
133                 .Object@name <- imxIdentifier(modelname, .Object@name)
134                 for (s in c("item", "mean", "cov")) {
135                         if (is.null(slot(.Object, s))) next;
136                         slot(.Object, s) <-
137                           imxConvertIdentifier(slot(.Object, s), modelname, namespace)
138                 }
139                 return(.Object)
140 })
141
142 setMethod("genericExpRename", signature("MxExpectationBA81"),
143         function(.Object, oldname, newname) {
144           # not sure what goes here yet
145                 return(.Object)
146 })
147
148 ##' Create a Bock & Aitkin (1981) expectation
149 ##'
150 ##' When a two-tier covariance matrix is recognized, this expectation
151 ##' automatically enables analytic dimension reduction (Cai, 2010).
152 ##' 
153 ##' The standard Normal distribution of the quadrature acts like a
154 ##' prior distribution for difficulty. It is not necessary to impose
155 ##' any additional Bayesian prior on difficulty estimates (Baker &
156 ##' Kim, 2004, p. 196).
157 ##'
158 ##' @param ItemSpec a single item model (to replicate) or a list of
159 ##' item models in the same order as the column of \code{ItemParam}
160 ##' @param item the name of the mxMatrix holding item parameters
161 ##' with one column for each item model with parameters starting at
162 ##' row 1 and extra rows filled with NA
163 ##' @param ...  Not used.  Forces remaining arguments to be specified by name.
164 ##' @param qpoints number of points to use for equal interval quadrature integration (default 49L)
165 ##' @param qwidth the width of the quadrature as a positive Z score (default 6.0)
166 ##' @param mean the name of the mxMatrix holding the mean vector
167 ##' @param cov the name of the mxMatrix holding the covariance matrix
168 ##' @param verbose the level of runtime diagnostics (default 0L)
169 ##' @param minItemsPerScore scores with fewer than this many items are considered NA (default 1L)
170 ##' @param naAction how to deal with NA scores. If "fail" then an
171 ##' error will be thrown. If "pass" then NAs will be excluded from the
172 ##' estimation but faithfully filled into the EAP scores. (default "fail")
173 ##' @param weightColumn the name of the column in the data containing the row weights (default NA)
174 ##' @param EstepItem a simple matrix of item parameters for the
175 ##' E-step. This option is mainly of use for debugging derivatives.
176 ##' @param debugInternal when enabled, some of the internal tables are
177 ##' returned in $debug. This is mainly of use to developers.
178 ##' @seealso \href{http://cran.r-project.org/web/packages/rpf/index.html}{RPF}
179 ##' @references
180 ##' Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item
181 ##' parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.
182 ##'
183 ##' Cai, L. (2010). A two-tier full-information item factor analysis
184 ##' model with applications. Psychometrika, 75, 581-612.
185 ##'
186 ##' Seong, T. J. (1990). Sensitivity of marginal maximum likelihood
187 ##' estimation of item and ability parameters to the characteristics
188 ##' of the prior ability distributions. Applied Psychological
189 ##' Measurement, 14(3), 299-311.
190
191 mxExpectationBA81 <- function(ItemSpec, item="item", ...,
192                               qpoints=49L, qwidth=6.0, mean="mean", cov="cov",
193                               verbose=0L, naAction="fail", minItemsPerScore=1L, weightColumn=NA_integer_,
194                               EstepItem=NULL, debugInternal=FALSE) {
195
196         if (length(list(...)) > 0) {
197                 stop(paste("Remaining parameters must be passed by name", deparse(list(...))))
198         }
199
200         if (packageVersion("rpf") < "0.28") stop("Please install 'rpf' version 0.28 or newer")
201         if (qpoints < 3) {
202                 stop("qpoints must be 3 or greater")
203         }
204         if (missing(qwidth)) qwidth <- 6
205         if (qwidth <= 0) {
206                 stop("qwidth must be positive")
207         }
208   
209         if (!is.list(ItemSpec)) ItemSpec <- list(ItemSpec)
210
211         actions <- c("pass", "fail")
212         if (!match(naAction, actions)) stop(paste("Valid choices for naAction are", omxQuotes(actions)))
213
214         if (minItemsPerScore < 0L) stop("minItemsPerScore must be non-negative")
215
216         if (is.na(weightColumn)) weightColumn <- as.integer(weightColumn)
217
218         return(new("MxExpectationBA81", ItemSpec, item, EstepItem,
219                    qpoints, qwidth, mean, cov, verbose, debugInternal,
220                    naAction, minItemsPerScore, weightColumn))
221 }
222
223 ##' Like simplify2array but works with vectors of different lengths
224 ##'
225 ##' Vectors are filled column-by-column into a matrix. Shorter vectors
226 ##' are padded with NAs to fill whole columns.
227 ##' @param x a list of vectors
228 ##' @param higher whether to produce a higher rank array (defaults to FALSE)
229 ##' @examples
230 ##' v1 <- 1:3
231 ##' v2 <- 4:5
232 ##' v3 <- 6:10
233 ##' mxSimplify2Array(list(v1,v2,v3))
234 ##'
235 ##' #     [,1] [,2] [,3]
236 ##' # [1,]    1    4    6
237 ##' # [2,]    2    5    7
238 ##' # [3,]    3   NA    8
239 ##' # [4,]   NA   NA    9
240 ##' # [5,]   NA   NA   10
241
242 mxSimplify2Array <- function(x, higher=FALSE) {
243         if (higher) {
244                 stop("higher=TRUE is not implemented. Consider using simplify2array")
245         }
246         par <- x
247   len <- sapply(par, length)
248   biggest <- which(len == max(len))[1]
249   out <- matrix(NA, nrow=max(len), ncol=length(par))
250   for (x in 1:length(par)) {
251     out[1:len[x],x] <- par[[x]]
252   }
253   rownames(out) <- names(par[[biggest]])
254   out
255 }