add MG support

philchalmers · Nov 13, 2024 · 9d7ec9a · 9d7ec9a
1 parent 5df35d5
commit 9d7ec9a
Show file tree

Hide file tree

Showing 7 changed files with 78 additions and 51 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,5 +1,5 @@
 Package: mirt
-Version: 1.42.6
+Version: 1.42.7
 Type: Package
 Title: Multidimensional Item Response Theory
 Authors@R: c( person("Phil", family="Chalmers", email =

diff --git a/NEWS.md b/NEWS.md
@@ -16,10 +16,10 @@
   slopes fixed to 1 and estimation of the latent variance term, mimicking the 
   `Rasch` itemtype family
 
-- `mirt()` gains `itemdesign` and `item.formula` arguments to fit fixed 
-  item design characteristics (e.g. LLTMs; Fischer, 1983) to all or a subset
-  of items. Arguments are similar to those in `mixedmirt()`, though generally 
-  not as flexible
+- `mirt()` and `multipleGroup()` gain `itemdesign` and `item.formula` arguments 
+  to fit fixed item design characteristics (e.g. LLTMs; Fischer, 1983) 
+  to all or a subset of items. Arguments are similar to those in `mixedmirt()`, 
+  though currently not as flexible
 
 - Partially-compensatory family of `itemtypes` now behave more consistently
   when loading structures specified where trace lines products are only 

diff --git a/R/mirt.R b/R/mirt.R
@@ -1438,50 +1438,8 @@ mirt <- function(data, model = 1, itemtype = NULL, guess = 0, upper = 1, SE = FA
                                               dentype=dentype, formula=formula, method=method)
     if(!is.null(latent.regression$data))
         data <- latent.regression$data
-    mixed.design <- NULL
-    if(!is.null(itemdesign)){
-        stopifnot('itemdesign only supported for dichotmous item tests' =
-                      all(apply(data, 2, \(x) length(na.omit(unique(x)))) == 2))
-        if(nrow(itemdesign) < ncol(data)){
-            has_idesign <- colnames(data) %in% rownames(itemdesign)
-            if(!any(has_idesign))
-                stop('No rownames in itemdesign match colnames(data)', call.=FALSE)
-            dummy <- as.data.frame(matrix(NA, sum(!has_idesign), ncol(itemdesign)))
-            colnames(dummy) <- colnames(itemdesign)
-            itemdesign <- rbind(dummy, itemdesign)
-        } else {
-            has_idesign <- rep(TRUE, nrow(itemdesign))
-            rownames(itemdesign) <- colnames(data)
-        }
-        itemdesignold <- itemdesign
-        if(is.list(item.formula)){
-            mf <- lapply(item.formula, \(x){
-                if(length(x) == 3){
-                    ghost <- x[[2]]
-                    itemdesign[[ghost]] <- 1
-                }
-                model.frame(x, itemdesign, na.action=NULL)
-            })
-            mm <- lapply(1:length(mf), \(i){
-                ret <- model.matrix(item.formula[[i]], mf[[i]])
-                ret[rowSums(is.na(ret)) > 0, ] <- NA
-                ret
-            })
-            names(mm) <- do.call(c, lapply(item.formula,
-                                           \(x) if(length(x) == 3) as.character(x[[2]]) else ""))
-            for(i in 1:length(mm))
-                if(names(mm)[i] != "")
-                    colnames(mm[[i]]) <- paste0(names(mm)[i], '.', colnames(mm[[i]]))
-            mm <- do.call(cbind, mm)
-        } else {
-            mf <- model.frame(item.formula, itemdesign, na.action = NULL)
-            mm <- model.matrix(item.formula, mf)
-        }
-
-        mixed.design <- list(random=NULL, fixed=mm, from='mirt',
-                             lr.random=NULL, lr.fixed=NULL, has_idesign=has_idesign)
-        attr(mixed.design, 'itemdesign') <- itemdesignold
-    }
+    mixed.design <- make.mixed.design(item.formula=item.formula,
+                                      item.design=item.design, data=data)
     mod <- ESTIMATION(data=data, model=model, group=rep('all', nrow(data)),
                       itemtype=itemtype, guess=guess, upper=upper, grsm.block=grsm.block,
                       pars=pars, method=method, constrain=constrain, SE=SE, TOL=TOL,

diff --git a/R/multipleGroup.R b/R/multipleGroup.R
@@ -52,6 +52,8 @@
 #'   differential item functioning (DIF) across groups
 #' @param method a character object that is either \code{'EM'}, \code{'QMCEM'}, or \code{'MHRM'}
 #'   (default is \code{'EM'}). See \code{\link{mirt}} for details
+#' @param itemdesign see \code{\link{mirt}} for details
+#' @param item.formula see \code{\link{mirt}} for details
 #' @param dentype type of density form to use for the latent trait parameters. Current options include
 #'   all of the methods described in \code{\link{mirt}}, as well as
 #'
@@ -534,10 +536,12 @@
 #' }
 multipleGroup <- function(data, model = 1, group, itemtype = NULL,
                           invariance = '', method = 'EM',
-                          dentype = 'Gaussian', ...)
+                          dentype = 'Gaussian', itemdesign=NULL, item.formula = NULL, ...)
 {
     Call <- match.call()
     dots <- list(...)
+    mixed.design <- make.mixed.design(item.formula=item.formula,
+                                      item.design=item.design, data=data)
     if(is.character(model)) model <- mirt.model(model)
     if(!is.null(dots$formula))
         stop('latent regression models not supported for multiple group yet', call.=FALSE) #TODO
@@ -565,7 +569,7 @@ multipleGroup <- function(data, model = 1, group, itemtype = NULL,
     }
     if(grepl('mixture', dentype)) group <- rep('full', nrow(data))
     mod <- ESTIMATION(data=data, model=model, group=group, invariance=invariance, method=method,
-                      itemtype=itemtype, dentype=dentype, ...)
+                      itemtype=itemtype, dentype=dentype, mixed.design=mixed.design, ...)
     if(is(mod, 'MultipleGroupClass') || is(mod, 'MixtureClass'))
         mod@Call <- Call
     return(mod)

diff --git a/R/utils.R b/R/utils.R
@@ -1903,6 +1903,54 @@ loadESTIMATEinfo <- function(info, ESTIMATE, constrain, warn){
     return(ESTIMATE)
 }
 
+make.mixed.design <- function(item.formula, item.design, data){
+    mixed.design <- NULL
+    if(!is.null(itemdesign)){
+        stopifnot('itemdesign only supported for dichotmous item tests' =
+                      all(apply(data, 2, \(x) length(na.omit(unique(x)))) == 2))
+        if(nrow(itemdesign) < ncol(data)){
+            has_idesign <- colnames(data) %in% rownames(itemdesign)
+            if(!any(has_idesign))
+                stop('No rownames in itemdesign match colnames(data)', call.=FALSE)
+            dummy <- as.data.frame(matrix(NA, sum(!has_idesign), ncol(itemdesign)))
+            colnames(dummy) <- colnames(itemdesign)
+            itemdesign <- rbind(dummy, itemdesign)
+        } else {
+            has_idesign <- rep(TRUE, nrow(itemdesign))
+            rownames(itemdesign) <- colnames(data)
+        }
+        itemdesignold <- itemdesign
+        if(is.list(item.formula)){
+            mf <- lapply(item.formula, \(x){
+                if(length(x) == 3){
+                    ghost <- x[[2]]
+                    itemdesign[[ghost]] <- 1
+                }
+                model.frame(x, itemdesign, na.action=NULL)
+            })
+            mm <- lapply(1:length(mf), \(i){
+                ret <- model.matrix(item.formula[[i]], mf[[i]])
+                ret[rowSums(is.na(ret)) > 0, ] <- NA
+                ret
+            })
+            names(mm) <- do.call(c, lapply(item.formula,
+                                           \(x) if(length(x) == 3) as.character(x[[2]]) else ""))
+            for(i in 1:length(mm))
+                if(names(mm)[i] != "")
+                    colnames(mm[[i]]) <- paste0(names(mm)[i], '.', colnames(mm[[i]]))
+            mm <- do.call(cbind, mm)
+        } else {
+            mf <- model.frame(item.formula, itemdesign, na.action = NULL)
+            mm <- model.matrix(item.formula, mf)
+        }
+
+        mixed.design <- list(random=NULL, fixed=mm, from='mirt',
+                             lr.random=NULL, lr.fixed=NULL, has_idesign=has_idesign)
+        attr(mixed.design, 'itemdesign') <- itemdesignold
+    }
+    mixed.design
+}
+
 make.randomdesign <- function(random, longdata, covnames, itemdesign, N, LR=FALSE){
     ret <- vector('list', length(random))
     for(i in seq_len(length(random))){

diff --git a/man/multipleGroup.Rd b/man/multipleGroup.Rd
diff --git a/tests/testthat/test-18-LLTM.R b/tests/testthat/test-18-LLTM.R
@@ -33,6 +33,17 @@ test_that('LLTM', {
     m2 <- M2(lltm.4)
     expect_equal(m2$TLI, 1.001862, tolerance=1e-2)
 
+    group <- factor(rep(c('G1', 'G2'), each=500))
+    lltm.G <- multipleGroup(dat, group=group, itemtype = 'Rasch', SE=TRUE, verbose=FALSE,
+                            item.formula = ~ 0 + difficulty, itemdesign=itemdesign)
+    cfsG1 <- coef(lltm.G, simplify=TRUE)$G1$items
+    cfsG2 <- coef(lltm.G, simplify=TRUE)$G2$items
+    expect_equal(as.vector(cfsG1[1, 1:3]), c(1.009664, 0.000000, 0.000000), tolerance=1e-2)
+    expect_equal(as.vector(cfsG2[1, 1:3]), c(0.9074375, 0.000000, 0.000000), tolerance=1e-2)
+    fs <- fscores(lltm.G)
+    expect_equal(fs[1:3], c(1.0170938, -0.2457873, -0.2457873), tolerance=1e-2)
+
+
 })
 
 test_that('MLTM', {