fix bugs in meansTest() when random effects are specified

NAMESPACE

@@ -52,7 +52,8 @@ importFrom("irlba",
 	"irlba")
 
 importFrom("nlme",
-	"lme")
+	"lme",
+	"lme.formula")
 
 importFrom("parallel",
 	"nextRNGStream")

NEWS

@@ -1,4 +1,11 @@
 
+CHANGES IN VERSION 3.6.3 [2024-6-26]
+-----------------------------------
+
+BUG FIXES
+
+    o Bug fixes for 'meansTest()' when random effects are specified
+
 CHANGES IN VERSION 3.6.2 [2024-6-10]
 -----------------------------------
 

R/stats-meansTest.R

@@ -52,14 +52,19 @@ setMethod("meansTest", "ANY",
 	datalist <- apply(y, 1L, function(yi)
 		{
 			data[[response]] <- yi
+			ok <- vars %in% names(data)
+			if ( !all(ok) ) {
+				labs <- paste0(vars[!ok], collapse=", ")
+				stop("couldn't find variable: ", labs)
+			}
 			as.data.frame(data[vars])
 		})
 	# fit models
 	if ( verbose ) {
 		lab <- if (n != 1L) "models" else "model"
 		message("fitting ", n, " ", lab)
 	}
-	FIT <- .lm_fit_fun(fixed, random)
+	FIT <- .lmFit_fun(fixed, random)
 	models <- chunkLapply(datalist, FIT,
 		nchunks=nchunks, verbose=verbose,
 		BPPARAM=BPPARAM, ...)
@@ -69,12 +74,14 @@ setMethod("meansTest", "ANY",
 		lab <- if (n != 1L) "models" else "model"
 		message("testing ", n, " ", lab)
 	}
-	TEST <- .lm_test_fun(reduced)
+	TEST <- .lmTest_fun(reduced, random)
 	tests <- chunkMapply(TEST, models, datalist,
 		nchunks=nchunks, verbose=verbose,
 		BPPARAM=BPPARAM, ...)
 	tests <- DataFrame(do.call(rbind, tests))
 	# return results
+	if ( anyNA(tests$statistic) )
+		warning(sum(is.na(tests$statistic)), " tests could not be performed")
 	if ( is.null(random) ) {
 		mcols <- DataFrame(fixed=deparse1(fixed), tests)
 	} else {
@@ -84,46 +91,45 @@ setMethod("meansTest", "ANY",
 	as(ResultsList(models, mcols=mcols), "MeansTest")
 })
 
-.lm_fit_fun <- function(fixed, random)
+.lmFit_fun <- function(fixed, random)
 {
 	FIT <- function(data, ...)
 	{
 		model <- NULL
 		if ( is.null(random) ) {
-			try(fit <- lm(fixed, data=data, ...), silent=TRUE)
-			try(model <- update(fit, . ~ ., data=data), silent=TRUE)
+			model <- try(lm(fixed, data=data, ...), silent=TRUE)
 		} else {
-			try(fit <- lme(fixed, data=data, random=random, ...), silent=TRUE)
-			try(model <- update(fit, . ~ ., data=data, method="ML"), silent=TRUE)
+			model <- try(lme(fixed, data=data,
+				random=random, method="ML", ...), silent=TRUE)
 		}
-		if ( !is.null(model) )
+		if ( !inherits(model, "try-error") )
 		{
-			model$model <- NULL
+			model <- update(model, . ~ .)
 			model$data <- data
 		}
 		model
 	}
 	FIT
 }
 
-.lm_test_fun <- function(reduced)
+.lmTest_fun <- function(reduced, random)
 {
 	TEST <- function(model, data)
 	{
-		if ( is.null(model) ) {
-			return(c(LR=NA, PValue=NA))
+		if ( inherits(model, "try-error") ) {
+			return(c(statistic=NA, pvalue=NA))
 		} else {
 			full <- model
 		}
 		if ( inherits(model, "lm") ) {
-			null <- update(full, reduced, data=data)
+			null <- update(full, reduced)
 			num <- as.numeric(logLik(null))
 			den <- as.numeric(logLik(full))
 			df <- abs(null$df.residual - full$df.residual)
 			LR <- -2 * (num - den)
 			PValue <- pchisq(LR, df, lower.tail=FALSE)
 		} else if ( inherits(model, "lme") ) {
-			null <- update(full, reduced, data=data, method="ML")
+			null <- update(full, reduced)
 			aov <- anova(null, full)
 			df <- abs(diff(aov[,"df"]))
 			LR <- aov[2L,"L.Ratio"]
@@ -275,7 +281,7 @@ setMethod("meansTest", "SpatialDGMM",
 		lab <- if (n != 1L) "models" else "model"
 		message("fitting ", n, " ", lab)
 	}
-	FIT <- .lm_fit_fun(fixed, random)
+	FIT <- .lmFit_fun(fixed, random)
 	models <- chunkLapply(datalist, FIT,
 		nchunks=nchunks, verbose=verbose,
 		BPPARAM=BPPARAM, ...)
@@ -285,7 +291,7 @@ setMethod("meansTest", "SpatialDGMM",
 		lab <- if (n != 1L) "models" else "model"
 		message("testing ", n, " ", lab)
 	}
-	TEST <- .lm_test_fun(reduced)
+	TEST <- .lmTest_fun(reduced, random)
 	tests <- chunkMapply(TEST, models, datalist,
 		nchunks=nchunks, verbose=verbose,
 		BPPARAM=BPPARAM, ...)

man/meansTest.Rd

@@ -60,7 +60,7 @@ segmentationTest(x, fixed, random, samples = run(x),
 
     \item{response}{The name of to assign the response variable in the fitted models.}
 
-    \item{reduced}{A one-sided formula specifying the reduced model for the null hypothesis. The default is an intercept-only model.}
+    \item{reduced}{A one-sided formula specifying the fixed effects in the reduced model for the null hypothesis. The default is an intercept-only model. Random effects are retained.}
 
     \item{byrow}{For the default method, are the rows or the columns the \code{x} .}
 
@@ -87,6 +87,10 @@ segmentationTest(x, fixed, random, samples = run(x),
     \item{layout}{A vector of the form \code{c(nrow, ncol)} specifying the number of rows and columns in the facet grid.}
 }
 
+\details{
+    Likelihood ratio tests are used for the hypothesis testing for consistency between fixed-effects models and mixed-effects models.
+}
+
 \value{
     An object of class \code{MeansTest} derived from \code{ResultsList}, where each element contains a linear model.
 }

tests/testthat/test-stats.R

@@ -279,7 +279,7 @@ test_that("spatialDGMM", {
 test_that("meansTest", {
 
 	set.seed(1)
-	s <- simulateImage(preset=4, dim=c(10L, 10L), nrun=3,
+	s <- simulateImage(preset=4, dim=c(10L, 10L), nrun=4,
 		representation="centroid")
 	s$truecondition <- ifelse(s$circleA | s$circleB, s$condition, NA)
 	s$truecondition <- factor(s$truecondition)
@@ -304,6 +304,18 @@ test_that("meansTest", {
 	expect_is(topf, "DataFrame")
 	expect_false(is.unsorted(rev(topf$statistic)))
 
+	s$random <- rep.int(NA_character_, length(s))
+	s$random <- replace(s$random, run(s) %in% c("runA1", "runA2"), "X")
+	s$random <- replace(s$random, run(s) %in% c("runA3", "runA4"), "Y")
+	s$random <- replace(s$random, run(s) %in% c("runB1", "runB2"), "U")
+	s$random <- replace(s$random, run(s) %in% c("runB3", "runB4"), "V")
+
+	mtr <- meansTest(s, fixed=~condition, random=~1|random)
+
+	expect_true(validObject(mtr))
+	expect_true(all(mcols(mtr)$statistic > 0))
+	expect_true(all(mcols(mtr)$pvalue > 0))
+
 	set.seed(2)
 	s2 <- simulateImage(preset=4, dim=c(10L, 10L), nrun=1,
 		representation="centroid")