Check neg binom vals

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

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'
     }

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)
+}

glmmTMB/src/glmmTMB.cpp

@@ -812,7 +812,9 @@ Type objective_function<Type>::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);

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
 }
 

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 :

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))