Fix for #131

R/LmWrapper.R

@@ -52,15 +52,17 @@ setMethod('initialize', 'LMlike', function(.Object, ..., design=data.frame()){
 
 setMethod('coef', signature=c(object='LMlike'), function(object, which, singular=TRUE, ...){
     stopifnot(which %in% c('C', 'D'))
-    co <- object@defaultCoef
-    if(which=='C' & object@fitted['C']){
-        co[names(coef(object@fitC))] <- coef(object@fitC)
+    co = object@defaultCoef
+    if(which == 'C' & object@fitted['C']){
+        some_co = coef(object@fitC)
+    } else if(which == 'D' & object@fitted['D']){
+        some_co = coef(object@fitD)
+    } else{
+        #!fitted
+        some_co = co
     }
-
-    if(which=='D' & object@fitted['D']){
-        co <- coef(object@fitD)
-    }
-    if(!singular) co <- co[!is.na(co)]
+    co[names(some_co)] = some_co
+    if(!singular) co = co[!is.na(co)]
     co
 })
 
@@ -79,16 +81,22 @@ setMethod('summary', signature=c(object='LMlike'), function(object){
 ##' @describeIn LMlike update the formula or design from which the \code{model.matrix} is constructed
 ##' @param formula. \code{formula}
 ##' @param design something coercible to a \code{data.frame}
+##' @param keepDefaultCoef \code{logical}. Should the coefficient names be preserved from \code{object} or updated if the model matrix has changed?
 ##' @param ... passed to \code{model.matrix}
 ##' @importMethodsFrom stats4 update coef vcov summary
-setMethod('update', signature=c(object='LMlike'), function(object, formula., design, ...){
+setMethod('update', signature=c(object='LMlike'), function(object, formula., design, keepDefaultCoef = FALSE, ...){
+    o_old = object
     if(!missing(formula.)){
         object@formula <- update.formula(object@formula, formula.)
     }
     if(!missing(design)){
         object@design <- as(design, 'data.frame')
     }
     model.matrix(object) <- model.matrix(object@formula, object@design, ...)
+    if(keepDefaultCoef){
+        object@defaultCoef = o_old@defaultCoef
+        object@defaultVcov = o_old@defaultVcov
+    }
     object@fitC <- object@fitD <- numeric(0)
     object@fitted <- c(C=FALSE, D=FALSE)
     object

R/ZlmFit-bootstrap.R

@@ -56,8 +56,8 @@ bootVcov1 <- function(zlmfit, R=99, boot_index = NULL){
 
     manyvc <- laply(boot_index, function(s){
         newsca <- sca[,s]
-        LMlike <- update(LMlike, design=colData(newsca))
-        zlm(sca=newsca, LMlike=LMlike, onlyCoef=TRUE)
+        LMlike <- update(LMlike, design=colData(newsca), keepDefaultCoef = TRUE)
+        zlm(sca=newsca, LMlike = LMlike, onlyCoef = TRUE)
     })
 
     manyvc

R/lmWrapper-glmer.R

@@ -45,7 +45,9 @@ toAdditiveFormula <- function(string){
 ##' @describeIn LMERlike update the formula or design matrix
 ##' @param formula. \code{formula}
 ##' @param design  something coercible to a \code{data.frame}
-setMethod('update', signature=c(object='LMERlike'), function(object, formula., design, ...){
+##' @param  keepDefaultCoef \code{logical}. Should the coefficient names be preserved from \code{object} or updated if the model matrix has changed?
+setMethod('update', signature=c(object='LMERlike'), function(object, formula., design, keepDefaultCoef = FALSE, ...){
+    o_old = object
     if(!missing(formula.)){
         object@formula <- update.formula(object@formula, formula.)
     }
@@ -54,6 +56,10 @@ setMethod('update', signature=c(object='LMERlike'), function(object, formula., d
         object@design <- as(design, 'data.frame')
     }
     model.matrix(object) <- model.matrix(as.formula(paste0('~', reComponents$FEform)), object@design, ...)
+    if(keepDefaultCoef){
+      object@defaultCoef = o_old@defaultCoef
+      object@defaultVcov = o_old@defaultVcov
+    }
     object@fitC <- object@fitD <- numeric(0)
     object@fitted <- c(C=FALSE, D=FALSE)
     object

inst/NEWS

@@ -37,3 +37,7 @@ Transition to use `SingleCellExperiment` and manage log counts slot.
 
 - deviance_residuals_hook had dependencies on lm.influence internals that abruptly changed in R 4.0 with no obvious fix.  This function is now defunct.
 - Defunct functions `cData`, `fData`, `zlm.SingleCellAssay`, `exprs`, `combine` have been removed.
+
+## 1.17.2
+
+- Fix issue in bootstrapping in which some factor levels aren't sampled in a the bag

tests/testthat/test-bootstrap.R

@@ -103,24 +103,25 @@ test_that('Bootstrap recovers covariance', {
 
 parallel::stopCluster(cl)
 
-# context("Nearly singular designs")
-# 
-# test_that('Bootstrap results are padded appropriately', {
-#     N <- 12
-#     m <- 20
-#     p <- 3
-#     X <- getX(p, N/p, N)
-#     beta <- rbind(2, matrix(0, nrow = p-1, ncol = m))
-#     Y <- simYs(m, X, beta, rho=0, sigma=1, p=.7)
-#     cData <- data.frame(group = attr(X, 'group'))
-#     cData$group = factor(cData$group)
-#     sca <- suppressMessages(suppressWarnings(FromMatrix(t(Y$Y), cData=cData)))
-#     zfit <- suppressWarnings(zlm(~group, sca=sca))
-#     # Only fit on groupA/groupB samples
-#     boot <- bootVcov1(zfit,R = NULL,boot_index  = list(c(rep(1, 4), 1:8)))
-#     expect_equal(colnames(boot), colnames(coef(zfit, 'D')))
-# })
-
-
+context("Nearly singular designs")
+
+test_that('Bootstrap results are padded appropriately', {
+    N <- 12
+    m <- 20
+    p <- 3
+    X <- getX(p, N/p, N)
+    beta <- rbind(2, matrix(0, nrow = p-1, ncol = m))
+    Y <- simYs(m, X, beta, rho=0, sigma=1, p=.7)
+    cData <- data.frame(group = attr(X, 'group'))
+    cData$group = factor(cData$group)
+    sca <- suppressMessages(suppressWarnings(FromMatrix(t(Y$Y), cData=cData)))
+    zfit <- suppressWarnings(zlm(~group, sca=sca))
+    # Only fit on groupA/groupB samples
+    boot <- bootVcov1(zfit,R = NULL,boot_index  = list(c(rep(1, 4), 1:8)))
+    expect_equal(colnames(boot), colnames(coef(zfit, 'D')))
+    expect_true(all(is.na(boot[,'groupC',])))
+})
+
+