From 0fefc9111160cac31747d881821dea7b7bb9bc2c Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Fri, 25 Oct 2024 13:02:49 -0400 Subject: [PATCH] Fix methods and add tests --- DESCRIPTION | 4 +- R/mixedModelTools.R | 30 ++++--- tests/testthat/test_random_effects.R | 120 +++++++++++++++++++++++++++ 3 files changed, 140 insertions(+), 14 deletions(-) create mode 100644 tests/testthat/test_random_effects.R diff --git a/DESCRIPTION b/DESCRIPTION index b0b3814..c604cd4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.4.3.9004 -Date: 2024-09-19 +Version: 1.4.3.9005 +Date: 2024-10-25 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( diff --git a/R/mixedModelTools.R b/R/mixedModelTools.R index abecc96..b249bac 100644 --- a/R/mixedModelTools.R +++ b/R/mixedModelTools.R @@ -28,7 +28,7 @@ get_Z <- function(formula, data, newdata=NULL){ # so d d d e e e --> d e d e d e # this should only change things for factors out <- lapply(Zt_list, function(x){ - id <- ave(rownames(x) == rownames(x), rownames(x), FUN=cumsum) + id <- stats::ave(rownames(x) == rownames(x), rownames(x), FUN=cumsum) x[order(id),,drop=FALSE] }) Zt <- do.call(rbind, out) @@ -89,8 +89,11 @@ check_formula <- function(formula){ if(form_has_correlated_effects(formula)){ stop("Correlated random slopes and intercepts are not supported. Replace | with || in your formula(s).", call.=FALSE) } + + rhs <- lapply(rand, function(x) x[[3]]) + rhs_all <- lapply(rhs, function(x) paste(deparse(x), collapse="")) - char <- paste(formula, collapse=" ") + char <- paste(rhs_all, collapse=" ") if(grepl(":|/", char)){ stop("Nested random effects (using / and :) are not supported", call.=FALSE) @@ -156,11 +159,12 @@ get_randvar_info <- function(tmb_report, type, formula, data){ sigma_est <- tmb_report$par.fixed[sigma_ind] sigma_cov <- as.matrix(tmb_report$cov.fixed[sigma_ind,sigma_ind]) re <- get_reTrms(formula, data) + Z <- get_Z(formula, data) # inefficient!! list(names=sigma_names(formula, data), estimates=sigma_est, covMat=sigma_cov, invlink="exp", invlinkGrad="exp", n_obs=nrow(data), n_levels=lapply(re$flist, function(x) length(levels(x))), cnms=re$cnms, - levels=rownames(re$Zt)) + levels=colnames(Z)) } get_fixed_names <- function(tmb_report){ @@ -175,7 +179,8 @@ print_randvar_info <- function(object){ val <- do.call(object$invlink, list(object$estimates)) - disp <- data.frame(Groups=names(object$cnms), Name=unlist(object$cnms), + disp <- data.frame(Groups=rep(names(object$cnms), sapply(object$cnms, length)), + Name=unlist(object$cnms), Variance=round(val^2,3), Std.Dev.=round(val,3)) cat("Random effects:\n") print(disp, row.names=FALSE) @@ -246,7 +251,7 @@ setMethod("sigma", "unmarkedEstimate", function(object, level=0.95, ...){ ses <- sqrt(diag(rinf$covMat)) lower <- vals - z*ses upper <- vals + z*ses - Groups <- names(rinf$cnms) + Groups <- rep(names(rinf$cnms), sapply(rinf$cnms, length)) Name <- unlist(rinf$cnms) data.frame(Model=object@short.name, Groups=Groups, Name=Name, sigma=exp(vals), lower=exp(lower), upper=exp(upper)) @@ -280,16 +285,17 @@ setMethod("randomTerms", "unmarkedEstimate", function(object, level=0.95, ...){ stop("No random effects in this submodel", call.=FALSE) } - Groups <- lapply(1:length(rv$cnms), function(x){ - gn <- names(rv$cnms)[x] - rep(gn, rv$n_levels[[gn]]) + gname <- rep(names(rv$cnms), sapply(rv$cnms, length)) + + Groups <- lapply(1:length(gname), function(x){ + rep(gname[x], rv$n_levels[[gname[x]]]) }) Groups <- do.call(c, Groups) - Name <- lapply(1:length(rv$cnms), function(x){ - gn <- names(rv$cnms)[x] - var <- rv$cnms[[x]] - rep(var, rv$n_levels[[gn]]) + Name <- unlist(rv$cnms) + + Name <- lapply(1:length(Name), function(x){ + rep(Name[x], rv$n_levels[[gname[x]]]) }) Name <- do.call(c, Name) diff --git a/tests/testthat/test_random_effects.R b/tests/testthat/test_random_effects.R new file mode 100644 index 0000000..3cacfc9 --- /dev/null +++ b/tests/testthat/test_random_effects.R @@ -0,0 +1,120 @@ +context("random effects tools") +skip_on_cran() + +set.seed(123) +M <- 10 +J <- 3 + +sc <- data.frame(x1 = rnorm(M), + x2 = factor(sample(letters[1:2], M, replace=TRUE)), + x3 = factor(sample(letters[3:4], M, replace=TRUE)), + g = factor(sample(letters[5:7], M, replace=TRUE))) + +test_that("random factor slopes work", { + + test1 <- get_Z(~x1 + x2 + (x1 + x2 || g), sc) + expect_is(test1, "dgCMatrix") + + test1 <- as.matrix(test1) + expect_equal(colnames(test1), rep(letters[5:7], 4)) + + expect_equal(test1, structure(c(0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, +0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1.55870831414912, +0, 0.129287735160946, 1.71506498688328, 0, 0, 0, -0.445661970099958, +0, -0.23017748948328, 0, 0, 0, 0, 0.460916205989202, 0, 0, 0, +-0.560475646552213, 0, 0, 0.070508391424576, 0, 0, 0, -1.26506123460653, +-0.686852851893526, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 1, 0, 0), dim = c(10L, 12L), dimnames = list(c("1", "2", +"3", "4", "5", "6", "7", "8", "9", "10"), c("e", "f", "g", "e", +"f", "g", "e", "f", "g", "e", "f", "g")))) + + sig <- sigma_names(~x1 + x2 + (x1 + x2 || g), sc) + expect_equal(sig, + c("1|g", "x1|g", "x2a|g", "x2b|g")) + + nr <- get_nrandom(~x1 + x2 + (x1 + x2 || g), sc) + expect_equal(length(nr), length(sig)) + expect_equal(sum(nr), ncol(test1)) + + test2 <- as.matrix(get_Z(~0 + x2:x3 + (0 + x2:x3 | g), sc)) + sig <- sigma_names(~0 + x2:x3 + (0 + x2:x3 | g), sc) + expect_equal(sig, c("x2a:x3c|g", "x2b:x3c|g", "x2a:x3d|g", "x2b:x3d|g")) + + expect_equal(test2, structure(c(0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, +0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, +0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, +1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(10L, +12L), dimnames = list(c("1", "2", "3", "4", "5", "6", "7", "8", +"9", "10"), c("e", "f", "g", "e", "f", "g", "e", "f", "g", "e", +"f", "g")))) + + nr <- get_nrandom(~0 + x2:x3 + (0 + x2:x3 | g), sc) + expect_equal(length(nr), length(sig)) + expect_equal(sum(nr), ncol(test2)) +}) + +test_that("random factor slope estimates work", { +skip_on_ci() + +M <- 10000 +J <- 3 + +set.seed(123) +sc <- data.frame(x1 = rnorm(M), + x2 = factor(sample(letters[1:2], M, replace=TRUE)), + x3 = factor(sample(letters[3:4], M, replace=TRUE)), + g = factor(sample(letters[5:26], M, replace=TRUE))) + +X <- model.matrix(~x1 + x2, sc) + +b <- c(-0.3, 0.3, 0.5) + +b2 <- c(-0.3, 0.3, 0, 0.5) + +re <- lapply(1:length(b2), function(i){ + rnorm(length(levels(sc$g)), 0, 0.2) +}) + +Z <- as.matrix(unmarked:::get_Z(~x1 + x2 + (x1 + x2 ||g), sc)) + +Xall <- cbind(X, Z) + +ball <- c(b, unlist(re)) + +psi <- plogis(Xall %*% ball) +z <- rbinom(M, 1, psi) + +y <- matrix(NA, M, J) +for (i in 1:M){ + y[i,] <- rbinom(3, 1, 0.4) * z[i] +} + +umf <- unmarkedFrameOccu(y=y, siteCovs=sc) + +fit <- occu(~1~x1 + x2 + (x1 + x2 || g), umf) + +expect_equivalent(coef(fit), + c(b, qlogis(0.4)), tol=0.1) + +expect_equivalent(sigma(fit)$sigma, + c(0.2865,0.1859,0.3004,0.1928), tol=1e-4) + +re_est <- randomTerms(fit) +expect_equal(re_est$Level, rep(levels(sc$g), 4)) + +pr <- predict(fit, 'state', level=NULL)$Predicted + +expect_true(cor(psi, pr) > 0.95) + +# Two grouping variables +fit2 <- occu(~1~x1 + x2 + (x1 + x2 || g) + (1|x3), umf[1:100,]) + +expect_equal(nrow(sigma(fit2)), 5) + +re <- randomTerms(fit2) +expect_equal(re$Groups[89:90], c("x3", "x3")) +})