Skip to content

Commit

Permalink
Fix methods and add tests
Browse files Browse the repository at this point in the history
  • Loading branch information
kenkellner committed Oct 25, 2024
1 parent 7270224 commit 0fefc91
Show file tree
Hide file tree
Showing 3 changed files with 140 additions and 14 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -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(
Expand Down
30 changes: 18 additions & 12 deletions R/mixedModelTools.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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){
Expand All @@ -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)
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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)

Expand Down
120 changes: 120 additions & 0 deletions tests/testthat/test_random_effects.R
Original file line number Diff line number Diff line change
@@ -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"))
})

0 comments on commit 0fefc91

Please sign in to comment.