diff --git a/glmmTMB/R/family.R b/glmmTMB/R/family.R index 984ac56d..46a93a64 100644 --- a/glmmTMB/R/family.R +++ b/glmmTMB/R/family.R @@ -21,16 +21,33 @@ in_glm_fit <- function() { identical(up_two[[1]], quote(glm.fit)) } -make_family <- function(x, link) { +make_family <- function(x, link, needs_nonneg = FALSE, needs_int = FALSE) { x <- c(x, list(link=link), make.link(link)) ## stubs for Effect.default/glm.fit if (is.null(x$aic)) { x <- c(x,list(aic=function(...) NA_real_)) } if (is.null(x$initialize)) { - ## should handle log-links adequately - x <- c(x,list(initialize=expression({mustart <- y+0.1}))) + x <- c(x,list(initialize= + substitute(env = list(FAMILY=x$family), + expr = expression({ + ## should handle log-links adequately + mustart <- y+0.1 + if (needs_int) { + if (any(abs(y - round(y)) > 0.001)) { + warning(gettextf("non-integer counts in a %s response variable", + FAMILY), domain = NA) + } + } + if (needs_nonneg) { + if (any(y < 0)) { + warning(gettextf("negative values in a %s response variable", + FAMILY), domain = NA) + } + } + })))) } + if (is.null(x$dev.resids)) { x <- c(x,list(dev.resids=function(y,mu,wt) { if (in_glm_fit()) { @@ -182,7 +199,7 @@ compois <- function(link="log") { if (length(phi)==1) phi <- rep(phi, length=length(mu)) .Call("compois_calc_var", mu, 1/phi, PACKAGE="glmmTMB") }) - return(make_family(r,link)) + return(make_family(r, link, needs_nonneg = TRUE, needs_int = TRUE)) } #' @rdname nbinom2 @@ -192,7 +209,7 @@ truncated_compois <- function(link="log") { variance=function(mu,phi) { stop("variance for truncated compois family not yet implemented") }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE, needs_int = TRUE)) } #' @rdname nbinom2 @@ -202,7 +219,7 @@ genpois <- function(link="log") { variance=function(mu,phi) { mu*phi }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE, needs_int = TRUE)) } #' @rdname nbinom2 @@ -212,7 +229,7 @@ truncated_genpois <- function(link="log") { variance=function(mu,phi) { stop("variance for truncated genpois family not yet implemented") }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE, needs_int = TRUE)) } #' @rdname nbinom2 @@ -222,7 +239,7 @@ truncated_poisson <- function(link="log") { variance=function(lambda) { (lambda+lambda^2)/(1-exp(-lambda)) - lambda^2/((1-exp(-lambda))^2) }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE, needs_int = TRUE)) } #' @rdname nbinom2 @@ -248,7 +265,7 @@ truncated_nbinom2 <- function(link="log") { (a*(1-pnbinom(c,mu=mu,size=theta))) return(mu_star + c*(mu_star-mu) +mu_star*mu*(1+1/a)-mu_star^2) }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE, needs_int = TRUE)) } #' @rdname nbinom2 @@ -258,7 +275,7 @@ truncated_nbinom1 <- function(link="log") { variance=function(mu,alpha) { stop("variance for truncated nbinom1 family not yet implemented") }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE, needs_int = TRUE)) } ## similar to mgcv::betar(), but simplified. @@ -287,8 +304,6 @@ beta_family <- function(link="logit") { return(make_family(r,link)) } -## fixme: better name? - #' @rdname nbinom2 #' @export ## variance= (Wikipedia) @@ -304,7 +319,7 @@ betabinomial <- function(link="logit") { variance = function(mu, phi) { mu*(1-mu) }, - initialize = binomial()$initialize) + initialize = our_binom_initialize(binomial()$initialize)) return(make_family(r,link)) } @@ -315,7 +330,7 @@ tweedie <- function(link="log") { variance = function(mu, phi, power) { phi * mu ^ power }) - return(make_family(r,link)) + return(make_family(r,link, needs_nonneg = TRUE)) } #' @rdname nbinom2 diff --git a/glmmTMB/R/glmmTMB.R b/glmmTMB/R/glmmTMB.R index a6bd577c..3a85bbc5 100644 --- a/glmmTMB/R/glmmTMB.R +++ b/glmmTMB/R/glmmTMB.R @@ -1241,6 +1241,11 @@ glmmTMB <- function( ## and then only to check whether it's NULL or not ... etastart <- mustart <- NULL + ## hack initialization method to flag negative values + if (family$family == "binomial") { + family$initialize <- our_binom_initialize(family$family) + } + if (!is.null(family$initialize)) { local(eval(family$initialize)) ## 'local' so it checks but doesn't modify 'y' and 'weights' } diff --git a/glmmTMB/R/utils.R b/glmmTMB/R/utils.R index 958dab7a..c300f846 100644 --- a/glmmTMB/R/utils.R +++ b/glmmTMB/R/utils.R @@ -682,3 +682,17 @@ get_family <- function(family) { } return(family) } + + +## add negative-value check to binomial initialization method +our_binom_initialize <- function(family) { + newtest <- substitute( + ## added test for glmmTMB + if (any(y<0)) { + stop(sprintf('negative values not allowed for the %s family', FAMILY)) + } + , list(FAMILY=family)) + b0 <- binomial()$initialize + b0[[length(b0)+1]] <- newtest + return(b0) +} diff --git a/glmmTMB/src/glmmTMB.cpp b/glmmTMB/src/glmmTMB.cpp index a3b3fc6f..2adda1e3 100644 --- a/glmmTMB/src/glmmTMB.cpp +++ b/glmmTMB/src/glmmTMB.cpp @@ -812,7 +812,9 @@ Type objective_function::operator() () s3 = sqrt(s1); //log-scale sd tmp_loglik = zt_lik_zero(yobs(i), dnorm(log(yobs(i)), s2, s3, true) - log(yobs(i))); - SIMULATE{yobs(i) = exp(rnorm(s2, s3));} // untested + SIMULATE{ + yobs(i) = exp(rnorm(s2, s3)); + } // untested break; case t_family: s1 = (yobs(i) - mu(i))/phi(i); diff --git a/glmmTMB/tests/testthat/test-basics.R b/glmmTMB/tests/testthat/test-basics.R index a4efb4b3..77729555 100644 --- a/glmmTMB/tests/testthat/test-basics.R +++ b/glmmTMB/tests/testthat/test-basics.R @@ -2,8 +2,10 @@ stopifnot(require("testthat"), require("glmmTMB")) +## drop (unimportant) info that may not match across versions drop_version <- function(obj) { obj$modelInfo$packageVersion <- NULL + obj$modelInfo$family$initialize <- NULL ## updated initialization expressions obj } diff --git a/glmmTMB/tests/testthat/test-families.R b/glmmTMB/tests/testthat/test-families.R index 1b06cec2..d66424b3 100644 --- a/glmmTMB/tests/testthat/test-families.R +++ b/glmmTMB/tests/testthat/test-families.R @@ -15,7 +15,6 @@ simfun0 <- function(beta=c(2,1), return(data.frame(x,f,mu)) } -context("alternative binomial specifications") test_that("binomial", { load(system.file("testdata","radinger_dat.RData",package="lme4")) radinger_dat <<- radinger_dat ## global assignment for testthat @@ -62,29 +61,44 @@ test_that("binomial", { expect_warning( glmmTMB(prop~1, family=binomial()) ) ## Warning as glm x <- c(1, 2, 3) ## weights=1 => x > weights ! expect_error ( glmmTMB(x~1, family=binomial(), - data = data.frame(x))) ## Error as glm + data = data.frame(x))) ## Error as glm + }) -context("non-integer count warnings") -test_that("count distributions", { - dd <- data.frame(y=c(0.5,1,1,1)) - for (f in c("binomial","betabinomial","poisson", - "genpois", - ## "compois", ## fails anyway ... +## check for negative values +test_that("detect negative values in two-column binomial response", { + x <- matrix(c(-1, 1, 2, 2, 3, 4), nrow = 3) + expect_error(glmmTMB(x~1, family=binomial(), data = NULL), + "negative values not allowed") +}) + +count_dists <-c("poisson", "genpois", "compois", "truncated_genpois", - # "truncated_compois", - "nbinom1", - "nbinom2" - # why do these truncated cases fail? - ##, "truncated_nbinom1", - ##"truncated_nbinom2" - )) { - expect_warning(m <- glmmTMB(y~1,data=dd,family=f), - "non-integer") + "nbinom1", "nbinom2", + "truncated_nbinom1", + "truncated_nbinom2" + ) + +binom_dists <- c("binomial", "betabinomial") + +test_that("count distributions", { + dd <- data.frame(y=c(0.5, rep(1:4, c(9, 2, 2, 2)))) + for (f in count_dists) { + expect_warning(glmmTMB(y~1, data=dd, family=f), "non-integer") + } +}) + +test_that("binom-type distributions", { + dd <- data.frame(y=c(0.5, rep(1:4, c(9, 2, 2, 2)))/10) + for (f in binom_dists) { + expect_warning(glmmTMB(y~1, + weights = rep(10, nrow(dd)), + data=dd, family=f), "non-integer") } }) -context("fitting exotic families") + + test_that("beta", { skip_on_cran() set.seed(101) @@ -333,7 +347,6 @@ test_that("truncated_genpois",{ }) -context("trunc compois") ##Compois test_that("truncated_compois",{ skip_on_cran() @@ -345,7 +358,6 @@ test_that("truncated_compois",{ expect_equal(predict(tcmp1,type="response")[1:2], c(19.4, 7.3), tolerance = 1e-6) }) -context("compois") test_that("compois", { skip_on_cran() # cmpdat <<- data.frame(f=factor(rep(c('a','b'), 10)), @@ -356,7 +368,6 @@ test_that("compois", { expect_equal(predict(cmp1,type="response")[1:2], c(19.4, 7.3), tolerance = 1e-6) }) -context("genpois") test_that("genpois", { skip_on_cran() gendat <<- data.frame(y=c(11,10,9,10,9,8,11,7,9,9,9,8,11,10,11,9,10,7,13,9)) @@ -365,7 +376,6 @@ test_that("genpois", { expect_equal(sigma(gen1), 0.235309, tolerance = 1e-6) }) -context("tweedie") test_that("tweedie", { skip_on_cran() ## Boiled down tweedie:::rtweedie : diff --git a/misc/neg_binom_ex.R b/misc/neg_binom_ex.R new file mode 100644 index 00000000..8a61a6ab --- /dev/null +++ b/misc/neg_binom_ex.R @@ -0,0 +1,31 @@ +## from https://github.com/glmmTMB/glmmTMB/issues/988 + +library(glmmTMB) +set.seed(101); dd <- data.frame(x = 1:10) +m1 <- glmmTMB(cbind(x, 10-x) ~ 1, family = betabinomial, data = dd) +sessionInfo() + +dd2 <- rbind(dd, data.frame(x=-1)) +update(m1, data = dd2) ## NOT detected +update(m1, data = dd2, family = binomial) ## gives weird/infinite answers + +try(glm(cbind(x, 10-x) ~ 1, family = binomial, data = dd2)) + + + +our_binomInitialize <- function(family) { + newtest <- substitute( + ## added test for glmmTMB + if (any(y<0)) { + stop(sprintf('negative values not allowed in %s responses', FAMILY)) + } + , list(FAMILY=family)) + b0 <- binomial()$initialize + b0[[3]] <- newtest + return(b0) +} + +b_new <- betabinomial() +b_new$initialize <- our_binomInitialize("betabinomial") + +try(update(m1, data = dd2, family = b_new))