Merge pull request #63 from florianhartig/OverdispersionTests

OverdispersionTests

Code/DHARMaPerformance/Power/TestPower.Rmd

@@ -0,0 +1,86 @@
+---
+title: "PowerTests"
+author: "Florian Hartig"
+date: "5/12/2018"
+output: 
+  html_document: 
+    keep_md: yes
+---
+
+
+```{r global_options, include=FALSE}
+knitr::opts_chunk$set(fig.width=7, fig.height=4.5, fig.align='center', warning=FALSE, message=FALSE, cache = T)
+```
+
+```{r, echo = F}
+library(AER)
+library(lme4)
+library(glmmTMB)
+library(DHARMa)
+```
+
+# Dispersion Tests
+
+## Benchmarking DHARMa dispersiontests against AER dispersiontest for Poisson GLM
+
+```{r}
+doCalculations <- function(control = 0){
+  testData = createData(sampleSize = 200, family = poisson(), overdispersion = control, randomEffectVariance = 0)
+  fittedModel <- glm(observedResponse ~ Environment1, data = testData, family = poisson())
+ 
+  out = list()
+  
+  res <- simulateResiduals(fittedModel = fittedModel, n = 250)
+  out$uniformTest = testUniformity(res)$p.value  
+  out$Dispersion = testDispersion(res, plot = F)$p.value  
+  out$DispersionAER = AER::dispersiontest(fittedModel)$p.value 
+  
+  res <- simulateResiduals(fittedModel = fittedModel, n = 250, refit = T)  
+  out$DispersionRefitted = testDispersion(res, plot = F)$p.value  
+  return(unlist(out))
+}
+
+# testing a single return
+doCalculations(control = 0.3)
+
+dispValues = seq(0,1.2, len = 20)
+
+# running benchmark
+out = runBenchmarks(doCalculations, controlValues = dispValues , nRep = 100,  parallel = F)
+
+tests = c("uniformity", "DHARMa disp basic" , "AER dispersiontest", "DHARMa disp refit")
+
+
+matplot(dispValues, t(out$summaries[,1,]), type = "l")
+legend("bottomright", tests, col = 1:4, lty = 1:4)
+
+```
+
+
+## Benchmarking DHARMa dispersiontests against pearson dispersiontest for glmms
+
+
+```{r}
+overdisp_fun <- function(model) {
+  ## number of variance parameters in 
+  ##   an n-by-n variance-covariance matrix
+  vpars <- function(m) {
+    nrow(m)*(nrow(m)+1)/2
+  }
+  model.df <- sum(sapply(VarCorr(model),vpars))+length(fixef(model))
+  rdf <- nrow(model.frame(model))-model.df
+  rp <- residuals(model,type="pearson")
+  Pearson.chisq <- sum(rp^2)
+  prat <- Pearson.chisq/rdf
+  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
+  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
+}
+```
+
+
+# Autocorrelation tests
+
+
+
+
+

Code/Experimental/parametricDispersionTestFragments.R

@@ -0,0 +1,117 @@
+
+############# DAS HIER WAR IN DER FUNKTION #####################
+
+if(method == "parametric"){
+  if(! class(model)[1] %in% c("glmerMod")){
+    return("DHARMa::testOverdispersionParametric currently only works for GLMMs in lme4. For Poisson GLMs, you can use AER::dispersiontest, otherwise use the non-parametric tests of DHARMa to test dispersion.")
+  }
+
+  ## number of variance parameters in
+  ##   an n-by-n variance-covariance matrix
+  vpars <- function(m) {
+    nrow(m)*(nrow(m)+1)/2
+  }
+
+  vcov(rd)
+
+  model.df <- sum(sapply(lme4::VarCorr(model),vpars))+length(lme4::fixef(model))
+  rdf <- nrow(model.frame(model))-model.df
+
+  # Residual df
+
+  rp <- residuals(model,type="pearson")
+  Pearson.chisq <- sum(rp^2)
+  prat <- Pearson.chisq/rdf
+
+  # Harrison 2014 seems to do more or less the same 
+
+  # Note: blmeco::dispersion_glmer defines the dispersion point estimate different
+
+  # computing  estimated scale  ( binomial model) following  D. Bates :
+  # That quantity is the square root of the penalized residual sum of
+  # squares divided by n, the number of observations, evaluated as:
+  # https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/015392.html
+
+  #n <- length(resid(modelglmer))
+  #return(  sqrt( sum(c(resid(modelglmer),modelglmer@u) ^2) / n ) ) 
+  #should be between, 0.75 and 1.4 if not under- or overdispersed, respectively
+
+  # Hypothesis test 
+
+  pval = switch(alternative, 
+                greater = pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE), 
+                two.sided = min(pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE),pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE))* 2, 
+                less = pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE)
+  )
+
+  out = list(
+    statistic=c(dispersion=prat, pearSS = Pearson.chisq, rdf=rdf), 
+    method = "Chisq test for overdispersion in GLMMs", 
+    data.name = as.character(model@call$family), 
+    p.value = pval)
+}
+
+if(method == "var"){
+  out = aerDispersion(object = model, alternative = alternative, ...)
+  out$method = "Cameron & Trivedi (1990)"
+}
+
+
+
+##################### ENDE ###########################
+
+
+
+
+## The following functions are only for validation purpose 
+
+#Function to calculate a point estimate of overdispersion from a mixed model object following Harrison 2014
+od.point<-function(modelobject){
+  x<-sum(resid(modelobject,type="pearson")^2)
+  rdf<-summary(modelobject)$AICtab[5]
+  return(x/rdf)
+}
+
+#Function to calculate point estimate of overdispersion from blmeco::dispersion_glmer
+od.point2<-function(modelglmer){
+  n <- length(resid(modelglmer))
+  return(sqrt(sum(c(resid(modelglmer), modelglmer@u)^2)/n))
+}
+
+# this is the AER dispersion function
+aerDispersion <- function (object, trafo = NULL, alternative = c("greater", "two.sided", "less")) 
+{
+  if (!inherits(object, "glm") || family(object)$family != 
+      "poisson") 
+    stop("only Poisson GLMs can be tested with method = var")
+  alternative <- match.arg(alternative)
+  otrafo <- trafo
+  if (is.numeric(otrafo)) 
+    trafo <- function(x) x^otrafo
+  y <- if (is.null(object$y)) 
+    model.response(model.frame(object))
+  else object$y
+  yhat <- fitted(object)
+  aux <- ((y - yhat)^2 - y)/yhat
+  if (is.null(trafo)) {
+    STAT <- sqrt(length(aux)) * mean(aux)/sd(aux)
+    NVAL <- c(dispersion = 1)
+    EST <- c(dispersion = mean(aux) + 1)
+  }
+  else {
+    auxreg <- lm(aux ~ 0 + I(trafo(yhat)/yhat))
+    STAT <- as.vector(summary(auxreg)$coef[1, 3])
+    NVAL <- c(alpha = 0)
+    EST <- c(alpha = as.vector(coef(auxreg)[1]))
+  }
+  rval <- list(statistic = c(z = STAT), 
+               p.value = switch(alternative, 
+                                greater = pnorm(STAT, lower.tail = FALSE), 
+                                two.sided = pnorm(abs(STAT), lower.tail = FALSE) * 2, 
+                                less = pnorm(STAT)
+               ), 
+               estimate = EST, 
+               null.value = NVAL, 
+               data.name = deparse(substitute(object)))
+  return(rval)
+}  

DHARMa/DESCRIPTION

@@ -1,6 +1,6 @@
 Package: DHARMa
 Title: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models
-Version: 0.1.6.3
+Version: 0.1.6.4
 Date: 2018-03-14
 Authors@R: c(person("Florian", "Hartig", email = "[email protected]", role = c("aut", "cre"), comment = "Theoretical Ecology, University of Regensburg, Regensburg, Germany"))
 Description: The 'DHARMa' package uses a simulation-based approach to create

DHARMa/NAMESPACE

@@ -14,8 +14,11 @@ export(plotResiduals)
 export(plotSimulatedResiduals)
 export(runBenchmarks)
 export(simulateResiduals)
+export(testDispersion)
+export(testGeneric)
 export(testOverdispersion)
 export(testOverdispersionParametric)
+export(testResiduals)
 export(testSimulatedResiduals)
 export(testSpatialAutocorrelation)
 export(testTemporalAutocorrelation)

DHARMa/NEWS

@@ -1,15 +1,20 @@
 NOTE: for more news about the package, see https://github.com/florianhartig/DHARMa/releases
 
-DHARMa 0.1.7
+DHARMa 0.2.0
 
 
 New features
 
-- support for glmmTMB https://github.com/florianhartig/DHARMa/issues/16
+- support for glmmTMB https://github.com/florianhartig/DHARMa/issues/16, implemented since https://github.com/florianhartig/DHARMa/releases/tag/v0.1.6.2
+
+Major changes
+
+- remodeled benchmarks functions in https://github.com/florianhartig/DHARMa/releases/tag/v0.1.6.3
+- remodeled dispersion tests, adresses https://github.com/florianhartig/DHARMa/issues/62 
 
 Minor changes
 
-- changed plot function names
+- changed plot function names in https://github.com/florianhartig/DHARMa/releases/tag/v0.1.6.1
 
 Bugfixes
 

DHARMa/R/plotResiduals.R

@@ -9,7 +9,7 @@
 #' @seealso \code{\link{plotResiduals}}, \code{\link{plotQQunif}}
 #' @export
 plotSimulatedResiduals <- function(simulationOutput, ...){
-  message("plotSimulatedResiduals is deprecated, switch your code to using the normal plot command")
+  message("plotSimulatedResiduals is deprecated, switch your code to using the plot function")
   plot(simulationOutput, ...)
 }
 

DHARMa/R/simulateResiduals.R

@@ -159,8 +159,7 @@ simulateResiduals <- function(fittedModel, n = 250, refit = F, integerResponse =
         out$refittedPredictedResponse[,i] = predict(refittedModel, type = "response")
         out$refittedFixedEffects[,i] = getFixedEffects(refittedModel)
         out$refittedResiduals[,i] = residuals(refittedModel, type = "response")
-        # try statement for glmmTMB
-        if(!out$modelClass == "glmmTMB") out$refittedPearsonResiduals[,i] = residuals(refittedModel, type = "pearson")
+        out$refittedPearsonResiduals[,i] = residuals(refittedModel, type = "pearson")
         #out$refittedRandomEffects[,i]  = ranef(refittedModel)
       }, silent = T)
     }
@@ -215,8 +214,7 @@ checkModel <- function(fittedModel){
   if(!(class(fittedModel)[1] %in% getPossibleModels())) warning("DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work!")
 
   if (class(fittedModel)[1] == "gam" ) if (class(fittedModel$family)[1] == "extended.family") stop("It seems you are trying to fit a model from mgcv that was fit with an extended.family. Simulation functions for these families are not yet implemented in DHARMa. See issue https://github.com/florianhartig/DHARMa/issues/11 for updates about this")
-
-  if (class(fittedModel)[1] == "glmmTMB" ) warning("Note that there are a few limitations for using glmmTMB with DHARMa. Please consult https://github.com/florianhartig/DHARMa/issues/16 for details.")
+
 }