From 684496cec413c7566ea695ef3f1304983c52238c Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 10 Dec 2023 14:56:09 -0500 Subject: [PATCH 1/6] Initial working lme4 replacement --- DESCRIPTION | 1 - R/getDesign.R | 4 +- R/mixedModelTools.R | 241 +++++++++++++++++++++++++++-------- R/predict.R | 4 +- tests/testthat/test_pcount.R | 14 +- 5 files changed, 201 insertions(+), 63 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d59b9851..ecfedd14 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,7 +22,6 @@ Depends: R (>= 4.0) Imports: graphics, lattice, - lme4, MASS, Matrix, methods, diff --git a/R/getDesign.R b/R/getDesign.R index 8597b2d6..7505bf45 100644 --- a/R/getDesign.R +++ b/R/getDesign.R @@ -10,8 +10,8 @@ setGeneric("handleNA", function(umf, ...) standardGeneric("handleNA")) setMethod("getDesign", "unmarkedFrame", function(umf, formula, na.rm=TRUE) { - detformula <- lme4::nobars(as.formula(formula[[2]])) - stateformula <- lme4::nobars(as.formula(paste("~", formula[3], sep=""))) + detformula <- remove_bars(as.formula(formula[[2]])) + stateformula <- remove_bars(as.formula(paste("~", formula[3], sep=""))) detVars <- all.vars(detformula) M <- numSites(umf) diff --git a/R/mixedModelTools.R b/R/mixedModelTools.R index 8c6a96d2..4581f1e2 100644 --- a/R/mixedModelTools.R +++ b/R/mixedModelTools.R @@ -1,50 +1,212 @@ -get_xlev <- function(data, model_frame){ - fac_col <- data[, sapply(data, is.factor), drop=FALSE] - xlevs <- lapply(fac_col, levels) - xlevs[names(xlevs) %in% names(model_frame)] -} - +# Generate required random effects info---------------------------------------- +# Sort-of drop-in replacement for lme4::mkReTrms get_reTrms <- function(formula, data, newdata=NULL){ - fb <- lme4::findbars(formula) - mf <- model.frame(lme4::subbars(formula), data, na.action=stats::na.pass) - if(is.null(newdata)) return(lme4::mkReTrms(fb, mf)) - new_mf <- model.frame(stats::terms(mf), newdata, na.action=stats::na.pass, - xlev=get_xlev(data, mf)) - lme4::mkReTrms(fb, new_mf, drop.unused.levels=FALSE) + if(!has_random(formula)){ + stop("No random effect terms in formula", call.=FALSE) + } + cnms <- get_cnms(formula) + # TODO: check these are factors + flist <- lapply(unique(names(cnms)), function(x){ + out <- data[[x]] + if(!is.factor(out)) out <- factor(out) + out + }) + names(flist) <- unique(names(cnms)) + list(Z = get_Z(formula, data), cnms = cnms, flist=flist) } +#Create Z random effects matrix from a formula, possibly using newdata--------- get_Z <- function(formula, data, newdata=NULL){ - if(is.null(lme4::findbars(formula))){ + # If no random effects in formula, return an empty Matrix + if(!has_random(formula)){ if(is.null(newdata)){ return(Matrix::Matrix(matrix(0, nrow=nrow(data), ncol=0),sparse=TRUE)) } else{ return(Matrix::Matrix(matrix(0, nrow=nrow(newdata), ncol=0),sparse=TRUE)) } } - check_formula(formula, data) - Zt <- get_reTrms(formula, data, newdata)$Zt - Z <- t(as.matrix(Zt)) + + # Get new formula + bars <- find_bars(formula) + new_form <- bars_to_formula(bars) + + # Get partial Z for each bar expression + form_list <- lapply(1:length(bars), function(x) bars_to_formula(bars[x])) + Z_parts <- lapply(form_list, get_partial_Z, data=data, newdata=newdata) + + # Create model frame + #mf <- model.frame(new_form, data, na.action=stats::na.pass) + #if(!is.null(newdata)){ + # mf <- model.frame(stats::terms(mf), newdata, na.action=stats::na.pass, + # xlev=get_xlev(data, mf)) + #} + + # Create sparse Z matrix + #Z <- model.matrix(new_form, mf) + Z <- do.call(cbind, Z_parts) + colnames(Z) <- Z_colnames(formula, data) Matrix::Matrix(Z, sparse=TRUE) } -get_group_vars <- function(formula){ - rand <- lme4::findbars(formula) - ifelse(is.null(rand), 0, length(rand)) +# Determine if formula has random effects specified---------------------------- +has_random <- function(form){ + length(find_bars(form)) > 0 } -get_nrandom <- function(formula, data){ - rand <- lme4::findbars(formula) - if(length(rand)==0) return(as.array(0)) +# Find 'bar' random effect parts of formula (i.e., the (x|y) structures)------- +# Operates recursively +find_bars <- function(form){ + out <- NULL + if(is.name(form)) return(NULL) + if(form[[1]] == as.name("(")) return(form) + if(is.call(form)){ + out <- lapply(form, find_bars) + } + unlist(out) +} + +# Convert bar components into new formula-------------------------------------- +# E.g. ~x + (1|g) becomes ~g - 1 +bars_to_formula <- function(bars){ + bar_terms <- lapply(bars, bars_to_terms) + check_duplicate_terms(bar_terms) + as.formula(str2lang(paste("~", paste(bar_terms, collapse = " + "), "- 1"))) +} + +# Translate bar components into standard formula terms------------------------- +bars_to_terms <- function(bars){ + info <- get_bar_info(bars) + new_terms <- sapply(info$LHS, function(x){ + if(x == "1") return(info$RHS) + paste0(x, ":", info$RHS) + }) + paste(new_terms, collapse=" + ") +} + +# Organize information in a bar expression into a list------------------------- +get_bar_info <- function(bar){ + out <- list( + operator = deparse(bar[[2]][[1]]), + LHS = terms_in_bar(bar, RHS=FALSE), + RHS = terms_in_bar(bar, RHS=TRUE) + ) + check_bar_info(out) +} + +terms_in_bar <- function(bars, RHS=FALSE){ + bars_sub <- bars[[2]][[2]] + if(RHS) bars_sub <- bars[[2]][[3]] + form <- formula(substitute(~X, list(X=bars_sub))) + trms <- attr(terms(form), "term.labels") + int <- attr(terms(form), "intercept") + if(int == 1 & !RHS) trms <- c("1", trms) + trms +} - out <- sapply(rand, function(x){ - col_nm <- as.character(x[[3]]) - length(unique(data[[col_nm]])) +# Check bar info to make sure the expression are allowed----------------------- +# For example unmarked doesn't support correlated random effects with | +check_bar_info <- function(info){ + if(info$operator == "|" & length(info$LHS) > 1){ + stop("Correlated random effects not supported, use || instead of |", call.=FALSE) + } + if(any(grepl(":", info$RHS)) | any(grepl("/", info$RHS))){ + stop("Nested random effects notation (: or /) not supported", call.=FALSE) + } + stopifnot(length(info$RHS) == 1) + info +} + +# Check terms to make sure there aren't any duplicates------------------------- +# E.g. as a result of a formula like ~ (1|g) + (x||) where the intercept +# is also implied in the second bar expression +check_duplicate_terms <- function(bar_terms){ + all_terms <- lapply(bar_terms, function(x){ + unlist(strsplit(x, " + ", fixed=TRUE)) + }) + all_terms <- unlist(all_terms) + dups <- duplicated(all_terms) + if(any(dups)){ + stop("Formula implies duplicate terms: ", paste0(all_terms[dups], collapse=", "), + call.=FALSE) + } + invisible() +} + +# Get partial Z for a given bar expression------------------------------------- +get_partial_Z <- function(formula, data, newdata){ + mf <- model.frame(formula, data, na.action=stats::na.pass) + if(!is.null(newdata)){ + mf <- model.frame(stats::terms(mf), newdata, na.action=stats::na.pass, + xlev=get_xlev(data, mf)) + } + model.matrix(formula, mf) +} + +# Get levels of factor--------------------------------------------------------- +# For use in specifying model frame +get_xlev <- function(data, model_frame){ + fac_col <- data[, sapply(data, is.factor), drop=FALSE] + xlevs <- lapply(fac_col, levels) + xlevs[names(xlevs) %in% names(model_frame)] +} + +# Generate colnames for Z to match what lme4::mkReTrms does-------------------- +Z_colnames <- function(formula, data){ + bars <- find_bars(formula) + info <- lapply(bars, get_bar_info) + groups <- sapply(info, function(x) x$RHS) + nreps <- sapply(info, function(x) length(x$LHS)) + lvls <- lapply(1:length(groups), function(i){ + group_dat <- data[[groups[i]]] + if(!is.factor(group_dat)) group_dat <- as.factor(group_dat) + rep(levels(group_dat), nreps[i]) + }) + unlist(lvls) +} + +# Get 'cnms' - random effect names - as with lme4::mkReTrms-------------------- +get_cnms <- function(formula){ + bars <- find_bars(formula) + info <- lapply(bars, get_bar_info) + cnms <- lapply(info, function(x){ + out <- lapply(x$LHS, function(y) ifelse(y == "1", "(Intercept)", y)) + names(out) <- rep(x$RHS, length(out)) + out + }) + do.call(c, cnms) +} + +# Get number of random effect SDs/variances------------------------------------ +get_nrandom <- function(formula, data){ + if(!has_random(formula)) return(as.array(0)) + cnms <- get_cnms(formula) + out <- sapply(names(cnms), function(x){ + length(unique(data[[x]])) }) as.array(out) } -has_random <- function(formula){ - length(lme4::findbars(formula)) > 0 +# Get number of grouping variables--------------------------------------------- +get_group_vars <- function(formula){ + if(!has_random(formula)) return(0) + cnms <- get_cnms(formula) + length(cnms) +} + +# Check if function has no support for random effects-------------------------- +check_no_support <- function(formula_list){ + has_bars <- any(sapply(formula_list, has_random)) + if(has_bars){ + stop("This function does not support random effects", call.=FALSE) + } +} + +# Remove all bar components from a formula------------------------------------- +remove_bars <- function(formula){ + s1 <- gsub("\\([^)]+\\|[^)]+\\) ?\\+?", "", deparse1(formula)) + s2 <- gsub(" \\+ {0,}$", "", s1) + if(s2 == "~") return(~1) + as.formula(str2lang(s2)) } sigma_names <- function(formula, data){ @@ -56,22 +218,6 @@ sigma_names <- function(formula, data){ nms } -check_formula <- function(formula, data){ - rand <- lme4::findbars(formula) - if(is.null(rand)) return(invisible()) - - char <- paste(formula, collapse=" ") - if(grepl(":|/", char)){ - stop("Nested random effects (using / and :) are not supported", - call.=FALSE) - } - theta <- get_reTrms(formula, data)$theta - if(0 %in% theta){ - stop("Correlated slopes and intercepts are not supported. Use || instead of |.", - call.=FALSE) - } -} - split_formula <- function(formula){ if(length(formula) != 3) stop("Double right-hand side formula required") char <- lapply(formula, function(x){ @@ -135,7 +281,7 @@ get_randvar_info <- function(tmb_report, type, formula, data){ 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(re$Z)) } get_fixed_names <- function(tmb_report){ @@ -158,13 +304,6 @@ print_randvar_info <- function(object){ #cat(paste0("Number of obs: ",object$n_obs,", groups: ",group_info,"\n")) } -check_no_support <- function(formula_list){ - has_bars <- any(sapply(formula_list, function(x) !is.null(lme4::findbars(x)))) - if(has_bars){ - stop("This function does not support random effects", call.=FALSE) - } -} - fit_TMB <- function(model, data, params, random, starts, method, ...){ diff --git a/R/predict.R b/R/predict.R index 4f3b5fb5..10bf4fcb 100644 --- a/R/predict.R +++ b/R/predict.R @@ -55,7 +55,7 @@ setMethod("predict", "unmarkedFit", # This function makes sure factor levels in newdata match, and that # any functions in the formula are handled properly (e.g. scale) make_mod_matrix <- function(formula, data, newdata, re.form=NULL){ - form_nobars <- lme4::nobars(formula) + form_nobars <- remove_bars(formula) mf <- model.frame(form_nobars, data, na.action=stats::na.pass) X.terms <- stats::terms(mf) fac_cols <- data[, sapply(data, is.factor), drop=FALSE] @@ -65,7 +65,7 @@ make_mod_matrix <- function(formula, data, newdata, re.form=NULL){ #X <- model.matrix(X.terms, newdata, xlev=xlevs) X <- model.matrix(form_nobars, nmf) offset <- model.offset(nmf) - if(is.null(re.form) & !is.null(lme4::findbars(formula))){ + if(is.null(re.form) & has_random(formula)){ Z <- get_Z(formula, data, newdata) X <- cbind(X, Z) } diff --git a/tests/testthat/test_pcount.R b/tests/testthat/test_pcount.R index 6db99667..35c9069a 100644 --- a/tests/testthat/test_pcount.R +++ b/tests/testthat/test_pcount.R @@ -93,7 +93,7 @@ test_that("pcount predict works",{ test_that("pcount can fit models with random effects",{ -set.seed(35) +set.seed(12) nSites <- 300 nVisits <- 3 x <- rnorm(nSites) # a covariate @@ -101,7 +101,7 @@ beta0 <- 0 beta1 <- 0.4 ran <- rnorm(100, 0, 1) -group <- factor(as.character(rep(1:100, each=3))) +group <- factor(as.character(rep(1:100, each=3)), levels=1:100) ran_ind <- as.numeric(group) lambda <- exp(beta0 + beta1*x + @@ -124,19 +124,19 @@ expect_is(fm, "unmarkedFitPCount") fmr <- pcount(~visit~x+(1|group), umf, K=25) -expect_equivalent(coef(fmr), c(0.05397,0.3197,-0.8760,1.3668,2.078), - tol=1e-3) +expect_equivalent(coef(fmr), c(-0.13648,0.3355,-0.9202,1.4680,2.47077), + tol=1e-4) expect_true(inherits(sigma(fmr), 'data.frame')) -expect_equal(sigma(fmr)$sigma, 1.05945, tol=1e-5) +expect_equal(sigma(fmr)$sigma, 1.038204, tol=1e-5) pr <- predict(fmr, "state") expect_equivalent(as.numeric(pr[1,]), - c(1.037050,0.58179,0.3453,3.1140), tol=1e-3) + c(0.31838, 0.20989, 0.087458,1.15905), tol=1e-3) pr2 <- predict(fmr, "state", re.form=NA) expect_equivalent(as.numeric(pr2[1,]), - c(1.48366,0.2011,1.1374,1.93255), tol=1e-3) + c(0.53086,0.090768,0.37970,0.74220), tol=1e-3) pr3 <- predict(fmr, "det") expect_true(inherits(pr3, "data.frame")) From 68cca855fdbf51906a2c255e24b5f0df67434f97 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 10 Dec 2023 19:05:43 -0500 Subject: [PATCH 2/6] Handle factor random slopes and add tests --- R/gdistremoval.R | 8 +-- R/mixedModelTools.R | 41 ++++++++------- R/power.R | 2 +- R/simulate.R | 4 +- tests/testthat/test_mixed_models.R | 80 ++++++++++++++++++++++++++++++ 5 files changed, 111 insertions(+), 24 deletions(-) create mode 100644 tests/testthat/test_mixed_models.R diff --git a/R/gdistremoval.R b/R/gdistremoval.R index b2536251..db5439b2 100644 --- a/R/gdistremoval.R +++ b/R/gdistremoval.R @@ -199,19 +199,19 @@ setMethod("getDesign", "unmarkedFrameGDR", if(return.frames) return(list(sc=sc, ysc=ysc, oc=oc)) - lam_fixed <- lme4::nobars(formula$lambdaformula) + lam_fixed <- remove_bars(formula$lambdaformula) Xlam <- model.matrix(lam_fixed, model.frame(lam_fixed, sc, na.action=NULL)) - phi_fixed <- lme4::nobars(formula$phiformula) + phi_fixed <- remove_bars(formula$phiformula) Xphi <- model.matrix(phi_fixed, model.frame(phi_fixed, ysc, na.action=NULL)) - dist_fixed <- lme4::nobars(formula$distanceformula) + dist_fixed <- remove_bars(formula$distanceformula) Xdist <- model.matrix(dist_fixed, model.frame(dist_fixed, ysc, na.action=NULL)) - rem_fixed <- lme4::nobars(formula$removalformula) + rem_fixed <- remove_bars(formula$removalformula) Xrem <- model.matrix(rem_fixed, model.frame(rem_fixed, oc, na.action=NULL)) diff --git a/R/mixedModelTools.R b/R/mixedModelTools.R index 4581f1e2..d8dad59a 100644 --- a/R/mixedModelTools.R +++ b/R/mixedModelTools.R @@ -28,11 +28,10 @@ get_Z <- function(formula, data, newdata=NULL){ # Get new formula bars <- find_bars(formula) - new_form <- bars_to_formula(bars) - + check_duplicate_terms(bars) + # Get partial Z for each bar expression - form_list <- lapply(1:length(bars), function(x) bars_to_formula(bars[x])) - Z_parts <- lapply(form_list, get_partial_Z, data=data, newdata=newdata) + Z_parts <- lapply(bars, get_partial_Z, data=data, newdata=newdata) # Create model frame #mf <- model.frame(new_form, data, na.action=stats::na.pass) @@ -42,7 +41,7 @@ get_Z <- function(formula, data, newdata=NULL){ #} # Create sparse Z matrix - #Z <- model.matrix(new_form, mf) + #Z <- model.matrix(new_form, mf, contrasts.arg=cont) Z <- do.call(cbind, Z_parts) colnames(Z) <- Z_colnames(formula, data) Matrix::Matrix(Z, sparse=TRUE) @@ -56,26 +55,25 @@ has_random <- function(form){ # Find 'bar' random effect parts of formula (i.e., the (x|y) structures)------- # Operates recursively find_bars <- function(form){ - out <- NULL + if(is.null(form)) return(NULL) if(is.name(form)) return(NULL) if(form[[1]] == as.name("(")) return(form) if(is.call(form)){ - out <- lapply(form, find_bars) + out <- return(unlist(lapply(form, find_bars))) } - unlist(out) + NULL } # Convert bar components into new formula-------------------------------------- -# E.g. ~x + (1|g) becomes ~g - 1 -bars_to_formula <- function(bars){ - bar_terms <- lapply(bars, bars_to_terms) - check_duplicate_terms(bar_terms) +# E.g. (1|g) becomes ~g - 1 +bar_to_formula <- function(bar){ + bar_terms <- bar_to_terms(bar) as.formula(str2lang(paste("~", paste(bar_terms, collapse = " + "), "- 1"))) } # Translate bar components into standard formula terms------------------------- -bars_to_terms <- function(bars){ - info <- get_bar_info(bars) +bar_to_terms <- function(bar){ + info <- get_bar_info(bar) new_terms <- sapply(info$LHS, function(x){ if(x == "1") return(info$RHS) paste0(x, ":", info$RHS) @@ -119,7 +117,8 @@ check_bar_info <- function(info){ # Check terms to make sure there aren't any duplicates------------------------- # E.g. as a result of a formula like ~ (1|g) + (x||) where the intercept # is also implied in the second bar expression -check_duplicate_terms <- function(bar_terms){ +check_duplicate_terms <- function(bars){ + bar_terms <- lapply(bars, bar_to_terms) all_terms <- lapply(bar_terms, function(x){ unlist(strsplit(x, " + ", fixed=TRUE)) }) @@ -132,14 +131,22 @@ check_duplicate_terms <- function(bar_terms){ invisible() } + # Get partial Z for a given bar expression------------------------------------- -get_partial_Z <- function(formula, data, newdata){ +get_partial_Z <- function(bar, data, newdata){ + info <- get_bar_info(bar) + formula <- bar_to_formula(bar) mf <- model.frame(formula, data, na.action=stats::na.pass) if(!is.null(newdata)){ mf <- model.frame(stats::terms(mf), newdata, na.action=stats::na.pass, xlev=get_xlev(data, mf)) } - model.matrix(formula, mf) + has_fac <- any(sapply(info$LHS[info$LHS != "1"], function(x) is.factor(data[[x]]))) + if(has_fac){ + stop("unmarked does not support random slopes for R factors.\n", + "Try converting the variable to a series of indicator variables instead.", call.=FALSE) + } + out <- model.matrix(formula, mf) } # Get levels of factor--------------------------------------------------------- diff --git a/R/power.R b/R/power.R index 3240e8e2..f0aa741d 100644 --- a/R/power.R +++ b/R/power.R @@ -163,7 +163,7 @@ check_coefs <- function(coefs, fit, template=FALSE){ # If there are random effects, adjust the expected coefficient names # to remove the b vector and add the grouping covariate name - rand <- lapply(formulas, lme4::findbars) + rand <- lapply(formulas, find_bars) if(!all(sapply(rand, is.null))){ stopifnot(all(required_subs %in% names(formulas))) rvar <- lapply(rand, function(x) unlist(lapply(x, all.vars))) diff --git a/R/simulate.R b/R/simulate.R index a8887cbc..efc372a4 100644 --- a/R/simulate.R +++ b/R/simulate.R @@ -91,7 +91,7 @@ setMethod("simulate", "character", #replace_sigma <- function(coefs, fit){ # required_subs <- names(fit@estimates@estimates) # formulas <- sapply(names(fit), function(x) get_formula(fit, x)) -# rand <- lapply(formulas, lme4::findbars) +# rand <- lapply(formulas, find_bars) # if(!all(sapply(rand, is.null))){ # rvar <- lapply(rand, function(x) unlist(lapply(x, all.vars))) # for (i in required_subs){ @@ -108,7 +108,7 @@ setMethod("simulate", "character", generate_random_effects <- function(coefs, fit){ required_subs <- names(fit@estimates@estimates) formulas <- sapply(names(fit), function(x) get_formula(fit, x)) - rand <- lapply(formulas, lme4::findbars) + rand <- lapply(formulas, find_bars) if(!all(sapply(rand, is.null))){ rvar <- lapply(rand, function(x) unlist(lapply(x, all.vars))) for (i in required_subs){ diff --git a/tests/testthat/test_mixed_models.R b/tests/testthat/test_mixed_models.R new file mode 100644 index 00000000..f0f7d10e --- /dev/null +++ b/tests/testthat/test_mixed_models.R @@ -0,0 +1,80 @@ +context("mixed model tools") + +test_that("get_reTrms matches lme4::mkReTrms", { + + skip_if(!requireNamespace("lme4", quietly=TRUE), + "lme4 package unavailable") + + dat <- data.frame(x = rnorm(20), y = rnorm(20), z = factor(sample(letters[1:3], 20, replace=T)), + group = factor(sample(letters[4:6], 20, replace=T)), + id = factor(sample(letters[7:9], 20, replace=T))) + + form1 <- ~x + (1|group) + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) + + form1 <- ~x + (x||group) + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) + + form1 <- ~x + (x||group) + (1|id) + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) + + form1 <- ~x + (x||group) + (y||id) + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) + + form1 <- ~x + (x*y||group) + (y||id) + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) + + form1 <- ~(1|group) + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) + + form1 <- ~(1|group) + x + r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r2 <- get_reTrms(form1, dat) + expect_identical(r2$Z, Matrix::t(r1$Zt)) + expect_identical(r1$cnms, r2$cnms) + attributes(r1$flist)$assign <- NULL + expect_identical(r1$flist, r2$flist) +}) + +test_that("get_reTrms errors correctly", { + form1 <- ~x + (x+z||group) + (y||id) + expect_error(get_reTrms(form1, dat)) + + form1 <- ~x + (x|group:id) + expect_error(get_reTrms(form1, dat)) + + form1 <- ~x + (x|group/id) + expect_error(get_reTrms(form1, dat)) + + expect_identical(find_bars(NULL), NULL) +}) From 9d26eb95b11731ccdf9828d94d9681694df2e003 Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 10 Dec 2023 20:26:24 -0500 Subject: [PATCH 3/6] Fix tests --- R/mixedModelTools.R | 4 +- tests/testthat/test_mixed_models.R | 107 +++++++++++++++-------------- 2 files changed, 57 insertions(+), 54 deletions(-) diff --git a/R/mixedModelTools.R b/R/mixedModelTools.R index d8dad59a..985d06e3 100644 --- a/R/mixedModelTools.R +++ b/R/mixedModelTools.R @@ -95,8 +95,8 @@ terms_in_bar <- function(bars, RHS=FALSE){ bars_sub <- bars[[2]][[2]] if(RHS) bars_sub <- bars[[2]][[3]] form <- formula(substitute(~X, list(X=bars_sub))) - trms <- attr(terms(form), "term.labels") - int <- attr(terms(form), "intercept") + trms <- attr(stats::terms(form), "term.labels") + int <- attr(stats::terms(form), "intercept") if(int == 1 & !RHS) trms <- c("1", trms) trms } diff --git a/tests/testthat/test_mixed_models.R b/tests/testthat/test_mixed_models.R index f0f7d10e..4cf3b88a 100644 --- a/tests/testthat/test_mixed_models.R +++ b/tests/testthat/test_mixed_models.R @@ -2,68 +2,71 @@ context("mixed model tools") test_that("get_reTrms matches lme4::mkReTrms", { - skip_if(!requireNamespace("lme4", quietly=TRUE), - "lme4 package unavailable") - + #skip_if(!requireNamespace("lme4", quietly=TRUE), + # "lme4 package unavailable") + set.seed(123) dat <- data.frame(x = rnorm(20), y = rnorm(20), z = factor(sample(letters[1:3], 20, replace=T)), group = factor(sample(letters[4:6], 20, replace=T)), id = factor(sample(letters[7:9], 20, replace=T))) - + + load('lme4_output.Rdata') form1 <- ~x + (1|group) - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + #l1 <- lme4::mkReTrms(lme4::findbars(form1), dat) + r1 <- get_reTrms(form1, dat) + expect_identical(r1$Z, Matrix::t(l1$Zt)) + expect_identical(r1$cnms, l1$cnms) + attributes(l1$flist)$assign <- NULL + expect_identical(r1$flist, l1$flist) - form1 <- ~x + (x||group) - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + form2 <- ~x + (x||group) + #l2 <- lme4::mkReTrms(lme4::findbars(form2), dat) + r2 <- get_reTrms(form2, dat) + expect_identical(r2$Z, Matrix::t(l2$Zt)) + expect_identical(r2$cnms, l2$cnms) + attributes(l2$flist)$assign <- NULL + expect_identical(r2$flist, l2$flist) - form1 <- ~x + (x||group) + (1|id) - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + form3 <- ~x + (x||group) + (1|id) + #l3 <- lme4::mkReTrms(lme4::findbars(form3), dat) + r3 <- get_reTrms(form3, dat) + expect_identical(r3$Z, Matrix::t(l3$Zt)) + expect_identical(r3$cnms, l3$cnms) + attributes(l3$flist)$assign <- NULL + expect_identical(r3$flist, l3$flist) - form1 <- ~x + (x||group) + (y||id) - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + form4 <- ~x + (x||group) + (y||id) + #l4 <- lme4::mkReTrms(lme4::findbars(form4), dat) + r4 <- get_reTrms(form4, dat) + expect_identical(r4$Z, Matrix::t(l4$Zt)) + expect_identical(r4$cnms, l4$cnms) + attributes(l4$flist)$assign <- NULL + expect_identical(r4$flist, l4$flist) - form1 <- ~x + (x*y||group) + (y||id) - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + form5 <- ~x + (x*y||group) + (y||id) + #l5 <- lme4::mkReTrms(lme4::findbars(form5), dat) + r5 <- get_reTrms(form5, dat) + expect_identical(r5$Z, Matrix::t(l5$Zt)) + expect_identical(r5$cnms, l5$cnms) + attributes(l5$flist)$assign <- NULL + expect_identical(r5$flist, l5$flist) - form1 <- ~(1|group) - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + form6 <- ~(1|group) + #l6 <- lme4::mkReTrms(lme4::findbars(form6), dat) + r6 <- get_reTrms(form6, dat) + expect_identical(r6$Z, Matrix::t(l6$Zt)) + expect_identical(r6$cnms, l6$cnms) + attributes(l6$flist)$assign <- NULL + expect_identical(r6$flist, l6$flist) - form1 <- ~(1|group) + x - r1 <- lme4::mkReTrms(lme4::findbars(form1), dat) - r2 <- get_reTrms(form1, dat) - expect_identical(r2$Z, Matrix::t(r1$Zt)) - expect_identical(r1$cnms, r2$cnms) - attributes(r1$flist)$assign <- NULL - expect_identical(r1$flist, r2$flist) + form7 <- ~(1|group) + x + #l7 <- lme4::mkReTrms(lme4::findbars(form7), dat) + r7 <- get_reTrms(form7, dat) + expect_identical(r7$Z, Matrix::t(l7$Zt)) + expect_identical(r7$cnms, l7$cnms) + attributes(l7$flist)$assign <- NULL + expect_identical(r7$flist, l7$flist) + + #save(l1,l2,l3,l4,l5,l6,l7, file='lme4_output.Rdata') }) test_that("get_reTrms errors correctly", { From 55061193ccf8ffdf74991336dba03233e835dfab Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Sun, 10 Dec 2023 20:43:00 -0500 Subject: [PATCH 4/6] Add missing reference objects --- tests/testthat/lme4_output.Rdata | Bin 0 -> 2701 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/testthat/lme4_output.Rdata diff --git a/tests/testthat/lme4_output.Rdata b/tests/testthat/lme4_output.Rdata new file mode 100644 index 0000000000000000000000000000000000000000..a211494cfd01ac87e26f8b3f32741adee9835aef GIT binary patch literal 2701 zcmV;83Uc)yiwFP!000001MOUYa8$(|-@JGCE_b;cG(ktu>4ZvYs$!$^s|-{Mj0&dm zLy#Xj;F(-Po5DY7hr5YkaV*N`5> zUr}($;#($6MJS-*NC^s{Fb=ik4fszb{uV(c`e-4yqzV5i$1*CPXnr&*B{W3zdkL4$ zm(zrN`muJ@MmxHUYf#HS04_-ekHx?e!cvMQoF6elWmuxLZVbx^EF-ae4od}=N^|g* z&h$A=&9n-1Ar>jcMOxy`i4<)@qs{6t-g@|wfG-XB&xh=B&BPJ)iH1abd!l|-D&5}H z(kNsH?u(}qG_9OIH0k9tnkF_Et7~m;i`S)BHMOKuiM3{KC?~W~PFo<;q0gLX&gadi ztsyM~b!gnWCK=kuD6z8NLSf*c20TZt?XY%YcCaoPPvJ~$B$pVjZ=C&&I9{ZT*chJh zHStEWOni9*3g@R0rn`ijUNW^hDyTjz4Qv4&>aZGY|7kg{K0<6akUSmo0@gfImgEN@ zA5@#B2Xm;P%5=;mJ)EQC(0iZV2-K(XvYg8dMssdcmNz?i_&Qh^uoQII>=ED*H7=bm zrwjRXx-IKyBohxNk}2e6AH(xmhi_^XvTu|iSPH&)u(7@Mfi~-dS_+g^&uvL3+UpW+ z=`WiN#>6fi+OX@$Aem;SXM<>wZMr(jL08?E3+$vpEuqaKjzl(_S}H4hk=B!Kry?<$ zgRrd~RXKb^)k9UhI>X%Q9(3PQAOfDtiZm2WuTG@naHiHawSbmrx3)o%Wb3*_JDBv^ zc=LVraaO${*_47vpsuAk1y8mj&E;ql&>?_ID{V_chtmF}en19}~NCIhHH1jKOjhmd|6U!tw;C5-+xBY~f9=8a-q(JRWo@26of$hS-o)Zv(e4ei z$1G>@{oi|Z>=(y1viOCzlNS!Z{UnRO`$*^Zy**o5`)6nDcq?B08y5e^CtXKMI%cr= z(Z}~3*!t{FR{x`?yAOSH%y<_6OQ`0l|AddR_*zG$FnPyJr)6Wva)8p0-5?l{gKUSXWg~{<@L(@=iUY;eZBqKDcSsRu1pnSk0G5*8W2UHqbEE)=K33Erl4FUH9 z+zxOzsG9-q1-KQ=od7oi+y`(QsJp=21aJ?iTfp1_>IMw@@$YDPIh0~1l-E>UTa|Hu zE}#n{d~vnmzQ!^)^@3%lqXL#FrRqdxAcP9oLi`Tp#dN+d-=&Z*79kalGj zCfrYKM87n`m(@JK??<~dT?U(^@CiZ@MnNdTDF{Va1)&J9AQWJh3~tF_mlS?MD8euZ zMK}hb2+JT8;F%PrK`6pC2nEMiFFiwPXgN`ZOxa_>z<#nGE-qYYQ1<$kd`#y#W z5YYay+s-iLaD6y79^qiv@6IFxRz(35*&zAOMBxYLlrBD$E*qqy zMeie)b#z&@rh0Djb<`(R>-C;Pcj~#=F{K0N5&mxPu^MiQefv68HbcTjutKd-Ie#sxTNeSpLSKD z#T`;k>D24vU)j8Ol9XHe?a%gp{G~Saq=u?wY=~3NuQ2pGwPB1 ziWaShq&Tsm=EMSMpg1vW+4aIu0-jt{=MgG0qFBg?p@wpbE>jw+=rTo@*)+6x@Sw7x zN<)3Fp}BOp!4lU_iQqj%Xb{ z+P)FGW-LOAtQRTX0P3CQ?UxVjYr=ODk`uleWG6ghbhja7Mf~Lm0FRXrbxRl9cPlUJ zS2j+s?ZX>*vho-Eyf7ZG^U&Av%hm;+OH8Nb9j&oJ>k*4eh-mT_A|av4+lMvCeCv>q z&$Ez8Jl}z3F_slr?!`h}PoixE5N8~$F=S^073R#%!F4#hWj{oiG*FbEDniaL?|$UTZ!2#9>z& z9-M&gb3fJRcI()84+$i6!*_WPx=!mF(^+8&&bk9OQ%rd6@^ zmD2b?-~YMuD^IN7S$G9kLhCa+KYQ)kg>7%Y;70Q^r{~bREuH7T{jeMD56;oeuNrF} zxqIWmOM4B_kN?)oo0|$`ckMxB*8)LK3{Lp}NR++)9z%uK=e>XI&PxU4aAWN!w6GkK zl>^r<8KdLjSmE*7p<<|!sJf;YHrEt}&{$utDac)fB|*sy*JDI>$zdo8 zxo9Zlp*C~=BIuNc`cFgEC5O^bpKItarI3||7K?^5-iq>3z%v1V1iTRAdnk{?_!-K( zFg}IyB#i$6UW4)#jE4Yz0eAz-2QZvZ@w>W)8}8R|L+xk!%MVo=>OTz~uBW^+U*PSc zxE!vA`NR?VlGDUwFPF*&vRx`WUM>$XtKsWw=gCNZ?X3JUf7+6vKFx-Vd{#fnt$vZQ zu Date: Tue, 30 Jan 2024 15:06:50 -0500 Subject: [PATCH 5/6] More tests --- DESCRIPTION | 4 ++-- R/mixedModelTools.R | 2 +- tests/testthat/lme4_output.Rdata | Bin 2701 -> 2818 bytes tests/testthat/test_mixed_models.R | 27 ++++++++++++++++++++++++++- 4 files changed, 29 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d42b5b99..93852c73 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.4.1.9001 -Date: 2024-01-23 +Version: 1.4.1.9002 +Date: 2024-01-30 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( diff --git a/R/mixedModelTools.R b/R/mixedModelTools.R index 985d06e3..db9d9e56 100644 --- a/R/mixedModelTools.R +++ b/R/mixedModelTools.R @@ -1,6 +1,6 @@ # Generate required random effects info---------------------------------------- # Sort-of drop-in replacement for lme4::mkReTrms -get_reTrms <- function(formula, data, newdata=NULL){ +get_reTrms <- function(formula, data){ if(!has_random(formula)){ stop("No random effect terms in formula", call.=FALSE) } diff --git a/tests/testthat/lme4_output.Rdata b/tests/testthat/lme4_output.Rdata index a211494cfd01ac87e26f8b3f32741adee9835aef..4449b57da6627051efea326c2b91292403980dbf 100644 GIT binary patch delta 2804 zcmVMKV(VuY39@3K5uV-_uFs3-)~P^I)CM@!IeRT zkczyjf>a;AM84(AZl62@A+LfX9^^#<94hgf@mDRr3L=j>PRRAd@K+Ubd@m&_We$feeQ_ zXNo?ax1BeJv<%dtapRnH(MHCKl}jxYI__`5W7OV%33DfA`%1>V z@FaDKjo|^`5NRdT#Fy8gaCR7BxINtPlA%>mPW54GU$?K!>zb_TFr5Rkj)~p58p}0U#$&k-%O|ma z)M5D)mWfy<>0@A8>l%`;buHgoFFP|(4`tK+|9bGbCwBds#b16nb@1hHr&*hq&v#82 z8;P-aZ>Vp}yz#48{O~s)pYWN9tt@_}Q7Z#~wt`(Xc0*8auWdtQ&!|Axi? z@nP=?PxovVKl$W=BRhYzm(_p&JAKE0zA}Cii~q&n@XUV#r&#>;p4p9Srd-eBzgqtM ziRq6V$UWcD#~3@pbH{=v9;kPjb(1lgJou+9G0k{>O@Ask8;>T zyoT~(I$xJ>)2h5Cq}$!MLed9W`~0(L+;xS-!azAPV#IMYmz4hSqux5StXs+{oqJ`{i`x%Qm2ykJ{mH=(zIXT? z8ON0V_taNvcdfl7<(e*z**2?o-#IDg^yXa+&jv%^ka17x<9|y&e&0*gQV#0f9dAFK zx%8rxi+Xq5&s!GFy$G7+O$+)ZZYpb>ge{kct!g5$%``tIYJVMABm*wSz)Z^a;;`!__ z%5h%uk9nI3K@N9s8z;|TR;vbHpWM^ZY=#_iU!k(~kX%lz*qoRb4HqY7e=NIM=!(~o zi)!6Nmy9SEGNRZ}&ZWy-4Rz@r83O_aq+4)AYw6JrjL_9%K~iM1Nb!16ewuesK6FVF zzEhE$@a;lo!ZS=)8$w3Je_w?F@E8eEyL7RAyYjO8%KGWFeRu;;R(`0@3;l6A4g>AK zY+m5G#B^HT(i&^D9G+Y$y*ZjxLTfWRd*_~ce@7SIbB!I%%{HC? ze)rL~Mf(`dDx0J!zF++F(A_!OSw@?~=6U4An+yNA=8JYTI@@I5Q7>MdI+}ZhQAVrN%{blP-~HoV zjqDeHQkos;^}oK*e_6k1$rE<8J8bUhLpLeyC!acQOT&W`(0A#F>Voe5wzNDr3;X-` z*FLxDUR#=;Iu+@+W`=*iBAk1*RYu!`^YQ+kyDF7BHosCDAL#o(_k8}T&3kjN;7VwH zM(3w**tn$Qji1`l{LJY&wrNMt<*z+zNBe_wblXeXhR5#Pe|mIej{yepYrU+wsX=Df z9z=F5;NwK&g#SgN?Dh9B%Dq1C{9|@rY9NQ}>mZ?pVUw&KzII6;EjPy+hsO>z#YUp; zF~z7krqG1O269Y6&LRv6YS&CLX&^@#npJllMr4N^noA*<4TW56Ggn#!ovWdxr=jj4 zhpV9_uA!rpe?oROv|Kcl@m7?N0-g!@BjAM?-$QvE#?MgRh4CqrCt>^t@EVk_U_1oy z3&0ytK7iqTir?L1xY2$LS8PACwER$4LrYIXN9!(c%@;U3D6WE|VLow8w&YAP*~y`@ zhD?XbmWRv3%xaYMvGX(}KX!KgvC_09#oo<|Mn1cLf649sBV%sex&OQEn!8d6N5*yf z|Nr0Z>8>l{uFYW%U01|)MO;?|Cebx_uDNqv5eKfwj6AN$P>r+~SsUragnJH$xr z%l!T;GQR2NP&YxV delta 2686 zcmV-^3W4>47L65uABzY80000000Zq@e{fXA9pAin_bzw295g{k(CLIqX{utQ@~aF~ z3XBS-@iZ!aLM9ZCQL;rpy5af3ZO6!wd4)>PbK~qK_&WVA-ALn|0%~ZDxYY6G%6)D zMD%+Jm(G{dgnat3cGN~Yx{Pa3%Rc}vNd}L_z!JhziY1&MF+ycnqO@)d%LptZv3w3o z1(r&4@RrVh^f^t-QLvFC}aoji>DGat(-nI>E$z;CN>wVYi(|e*QHlAwWL#twPtQ8C$vybTOiY+ z&zxw^=gp_BAuR)SXxzFc8QRDwv9jMnVc?+#JV&j6?XY%YcCaoPPvJ~$B$pVjZ=C&& zI9{ZT*chJhHStEWOni9*3g@R0rn`ijUNW^hDyTjz4Qv4&>aZGY|7kg{K0<6akUSmo z0@gfImgEN@A5@#B2Xm;P%5=;mJ)EQC(0iZV2-K(XvYg8dMssdcmNz?i_&Qh^uoQII z>=EF95j8HIFQ*Inbh<6;Xe1L4CXy-SWgo-yS%+_G6|!%XAXo~%c(Ad(^?^3)gIWrd zRnKinC)(>0ZRszY4aUSS9on$#$RL?!re}j_k!`v<%0XA%mkaEqK`o)pB925hn_4O> zdy&?YZKonJnuD;d9aTAeL)Al7ygI|&=^k`{-%=m~p390f6iu&Aq~ma=);6_(mT0%O zL6KzZx6@f=HmQr8xypwj#~tXe~)(-}_qfM!WyS_GUGGNU+sl zsDfMw`4ES4kS^xCwI`O7wT}vuZ`zG$FnPyJr)6Wva)8p0-5?l{gKUR7sbwRE z$g-g23h<3vl!g>Qi+9ua%9Gc1tFHBg60z}bR{xz&{$+* zpf6#$2FnC2HCXT*qsdslYR;+6!jN`l7AD+JY(#&*G{TqFJiqTpyEI(}o22jwLJ>wm zD8eZSMOX!)2(KU%V3rJS$zYchenBY0FbG9B2B8ScAQa%46sAEa!Zip5*d~Q<5DG9( zgmZ(ADcrd1yxQe;pA+8G;4uZyv-A5th6)hS{;}K6FywH3ICv)(dlT_ihqUczL&6^6 zVAy}}&Ljg?MFA7pAoFK}?<1CVbXl~fdT#P{)F)Ky^`1j_>bcip zC}2@6z+$MOt%ANvLlyM(uA#i9iNb59Q@W_pnx8eactFswwV}>lf*UzQ+ajfy*8q1l zob=ZN@ilQXu@~ZI0kV=HWpwR;L?Jtlv=)EtCC(EkPWfVcVt=yw0I0-#uj7OEW#e`0 zUz`sZ55#G{qcwx}qf{q|6(s}^xd^X9Ciz2N0JsjB#1An9Q;|ue&BQVX%WYU{%{8)F zvnY~Vv*^GRK}BaC7sxoG?3v5{G=J@`6*8_Ud+5ZEuR3{kvy3y!p1yy>=F{K0N5+30 zWl!E&`BGr#78!?>J&qPF|J{}M%DANLC!cmzp~W3iPU+O^<6qglcaoG_`t8s5e*B~T z@5?x*?7t_yUb$`6IVsn4cEqOXmAg(!Ij47SuX!dC{jQ9A${zk(>e0Jj9Vz9YKG^c! zleu%RNV%vFM*qBF{;XT1oYbj#^{an3cCMCkQzybp=gyk6M#@q3?Af{DjbmSzaaGyv zUB_=*GG(fav&z1)@py7`cejkY%D%aLbL1!2jUHSME8u?m=C?-_dE>GIoRE6_yFfgj z9YY1qOa8GRGf~Ll`gU>hEM~R5@AXNaj%G9Jk^71kt%syIv7zR~0%)K(F>8O>^}hH{E7QyQx1GDVl!G_-i|pt7M#Lw&BHz)S!G0Zaoh3cwsNh5(oV ziuf6_2Lw-%Iz!}uya8bYk_N;~kue}(K)QpDXdOM;z7e`+EJBK`7b)HV>Ye88mk;f0 z!gmsq6TTT_Cp=?xw;^Oj{N;ZL0FRXrbxRl9cPlUJS2j+s?ZX>*vho-Eyf7ZG^U&Av z%hm;+OH8Nb9j&oJ>k*4eh-mT_A|av4+lMvCeCv>q&$Ez8Jl}z3F_slr?!`h}PoixE z5N8~$F=ypRvGQh4yP9vOt^kqc%d5&&JMAcx1YOf=>v>LmeJyH*zFIEYgzem=vG&n z91g`dcm3j@`|c>v&NA8@4$lKe-kJBum0xqC(b*yUzJB)d^ufX_j51oC9>&?8o{s0X z)w2KilhW)!umAPs_UeE63m)VJ^I99R~jCifbMfY)#rBXcBSRPS=imP zyYl&Uce&E^)Tzk6J1zG6rLn@RtuopkoR5!o++L3z%v1V1iTRAdnk{?_!-K(Fg}IyB#i$6UW4)#jE4Yz z0eAz-2QZvZ@w>W)8}8R|L+xk!%MVo=>OTz~uBW^+U*PScxE!vA`NR?VlGDUwFPF*& zvRx`WUM>$XtKsWw=gCNZ?X3JUf7+6vKFx-Vd{#fnt$sR@v9Ryd?`|t|=ZEmpxK987 szq>s}xgu(B4sfVk5#@>~R|FPOnLB0flq=%oilEQ_2dXeqJ7aAC07}DBUjP6A diff --git a/tests/testthat/test_mixed_models.R b/tests/testthat/test_mixed_models.R index 4cf3b88a..ca7cb05a 100644 --- a/tests/testthat/test_mixed_models.R +++ b/tests/testthat/test_mixed_models.R @@ -1,4 +1,5 @@ context("mixed model tools") +skip_on_cran() test_that("get_reTrms matches lme4::mkReTrms", { @@ -65,11 +66,35 @@ test_that("get_reTrms matches lme4::mkReTrms", { expect_identical(r7$cnms, l7$cnms) attributes(l7$flist)$assign <- NULL expect_identical(r7$flist, l7$flist) + + # Check that unused factor levels aren't dropped + dat2 <- data.frame(x = c(0.1, 0.2, -0.1), + group = factor(c("a","b","c"), levels=c("a","b","c","d"))) + form8 <- ~x + (1|group) + #l8 <- lme4::mkReTrms(lme4::findbars(form8), dat2, drop.unused.levels=FALSE) + r8 <- get_reTrms(form8, dat2) + expect_identical(r8$Z, Matrix::t(l8$Zt)) + + # Check that get_Z handles newdata + form9 <- ~x + (x||group) + (1|id) + nd <- data.frame(x=c(0.5,1), group=c("e","d"), id=c("i","i")) + Z <- get_Z(form9, dat, newdata=nd) + expect_equivalent(as.matrix(Z), matrix(c(0,1,0,0,0.5,0,0,0,1, + 1,0,0,1,0,0,0,0,1), nrow=2, byrow=T)) + + # New level + nd <- data.frame(x=c(0.5,1), group=c("a","d"), id=c("i","i")) + expect_error(get_Z(form9, dat, newdata=nd)) - #save(l1,l2,l3,l4,l5,l6,l7, file='lme4_output.Rdata') + #save(l1,l2,l3,l4,l5,l6,l7,l8, file='lme4_output.Rdata') }) test_that("get_reTrms errors correctly", { + set.seed(123) + dat <- data.frame(x = rnorm(20), y = rnorm(20), z = factor(sample(letters[1:3], 20, replace=T)), + group = factor(sample(letters[4:6], 20, replace=T)), + id = factor(sample(letters[7:9], 20, replace=T))) + form1 <- ~x + (x+z||group) + (y||id) expect_error(get_reTrms(form1, dat)) From 20a7cae2d08dc9390f7110eed20288751d4266cf Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Tue, 30 Jan 2024 15:48:24 -0500 Subject: [PATCH 6/6] Fix occuCOP and IDS --- R/IDS.R | 2 ++ R/occuCOP.R | 4 ++-- tests/testthat/test_IDS.R | 3 +++ 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/R/IDS.R b/R/IDS.R index 38338b12..c0b15346 100644 --- a/R/IDS.R +++ b/R/IDS.R @@ -73,6 +73,8 @@ IDS <- function(lambdaformula = ~1, formlist <- list(lam=lambdaformula, ds=form_hds, pc=form_pc, oc=form_oc, phi=availformula) + check_no_support(formlist) + stopifnot(inherits(dataDS, "unmarkedFrameDS")) stopifnot(inherits(dataPC, c("unmarkedFramePCount", "NULL"))) stopifnot(inherits(dataOC, c("unmarkedFrameOccu", "NULL"))) diff --git a/R/occuCOP.R b/R/occuCOP.R index f22c2354..07875e3f 100644 --- a/R/occuCOP.R +++ b/R/occuCOP.R @@ -79,7 +79,7 @@ setMethod( # Occupancy submodel ------------------------------------------------------- # Retrieve the fixed-effects part of the formula - psiformula <- lme4::nobars(as.formula(formlist$psiformula)) + psiformula <- remove_bars(as.formula(formlist$psiformula)) psiVars <- all.vars(psiformula) # Retrieve the site covariates @@ -109,7 +109,7 @@ setMethod( # Detection submodel ------------------------------------------------------- # Retrieve the fixed-effects part of the formula - lambdaformula <- lme4::nobars(as.formula(formlist$lambdaformula)) + lambdaformula <- remove_bars(as.formula(formlist$lambdaformula)) lambdaVars <- all.vars(lambdaformula) # Retrieve the observation covariates diff --git a/tests/testthat/test_IDS.R b/tests/testthat/test_IDS.R index 61d348db..7565b25c 100644 --- a/tests/testthat/test_IDS.R +++ b/tests/testthat/test_IDS.R @@ -72,6 +72,9 @@ test_that("IDS can fit models with covariates", { unitsOut="kmsq") expect_equal(length(coef(mod_sim)), 4) expect_equal(length(coef(mod_sep)), 5) + + # Doesn't support random effects + expect_error(update(mod_sep, lambdaformula = ~elev + (1|group))) }) test_that("IDS can fit models with occupancy data", {