Merge pull request #445 from florianhartig/406-errorMessageQgam

improving error message for qgam + plot fixes #406
florianhartig · Oct 15, 2024 · 74c39e5 · 74c39e5
2 parents 38a313b + 9465f37
commit 74c39e5
Show file tree

Hide file tree

Showing 3 changed files with 61 additions and 12 deletions.
diff --git a/Code/DHARMaIssues/406.R b/Code/DHARMaIssues/406.R
@@ -0,0 +1,36 @@
+## improving the information about errors in qgam package to calculate quantile regression of the residuals
+# from issue #406
+
+library(DHARMa)
+library(glmmTMB)
+
+# reproducible example from the issue
+habitat <- c(rep("Town", 6), rep("Country", 6), rep("City", 6))
+location <- c(rep(c("a","b","c"), 6))
+beetles <- c(275, 89, 776, 325, 3302, 3335, 288, 566, 1546, 80, 282, 2529, 550, 479, 983, 372, 450, 902)
+testData <- data.frame(habitat, location, beetles)
+
+# it doesn't throw an error with 18 data points
+M2_Total <- glmmTMB(beetles ~ habitat + (1 | location),
+                    family = "nbinom2", data=testData)
+simulationOutput_M2_Total <- simulateResiduals(fittedModel = M2_Total,  plot = F)
+plot(simulationOutput_M2_Total)
+
+# with 14 data points, we can see the warning messages, but quantiles are plotted:
+M3_Total <- glmmTMB(beetles ~ habitat + (1 | location),
+                    family = "nbinom2", data=testData[1:14,])
+simulationOutput_M3_Total <- simulateResiduals(fittedModel = M3_Total,  plot = F)
+plot(simulationOutput_M3_Total)
+
+# with 13 data points, we can see the error messages + warnings, just ONE quantile is plotted:
+M4_Total <- glmmTMB(beetles ~ habitat + (1 | location),
+                    family = "nbinom2", data=testData[1:13,])
+simulationOutput_M4_Total <- simulateResiduals(fittedModel = M4_Total,  plot = F)
+plot(simulationOutput_M4_Total)
+
+# with 10 data points, we can see the error messages + warnings, no quantiles plotted:
+M5_Total <- glmmTMB(beetles ~ habitat + (1 | location),
+                    family = "nbinom2", data=testData[1:10,])
+simulationOutput_M5_Total <- simulateResiduals(fittedModel = M5_Total,  plot = F)
+plot(simulationOutput_M5_Total)
+
diff --git a/DHARMa/R/plots.R b/DHARMa/R/plots.R
@@ -181,7 +181,10 @@ plotQQunif <- function(simulationOutput, testUniformity = TRUE, testOutliers = T
 #' @seealso [plotQQunif], [testQuantiles], [testOutliers]
 #' @example inst/examples/plotsHelp.R
 #' @export
-plotResiduals <- function(simulationOutput, form = NULL, quantreg = NULL, rank = TRUE, asFactor = NULL, smoothScatter = NULL, quantiles = c(0.25, 0.5, 0.75), absoluteDeviation = FALSE, ...){
+plotResiduals <- function(simulationOutput, form = NULL, quantreg = NULL,
+                          rank = TRUE, asFactor = NULL, smoothScatter = NULL,
+                          quantiles = c(0.25, 0.5, 0.75),
+                          absoluteDeviation = FALSE, ...){
 
 
   ##### Checks #####
@@ -233,13 +236,17 @@ plotResiduals <- function(simulationOutput, form = NULL, quantreg = NULL, rank =
 
   # categorical plot
   if(is.factor(pred)){
-    testCategorical(simulationOutput = simulationOutput, catPred = pred, quantiles = quantiles)
+    testCategorical(simulationOutput = simulationOutput, catPred = pred,
+                    quantiles = quantiles)
   }
   # smooth scatter
   else if (smoothScatter == TRUE) {
     defaultCol = ifelse(res == 0 | res == 1, 2,blackcol)
-    do.call(graphics::smoothScatter, append(list(x = pred, y = res , ylim = c(0,1), axes = FALSE, colramp = colorRampPalette(c("white", "darkgrey"))),a))
-    points(pred[defaultCol == 2], res[defaultCol == 2], col = .Options$DHARMaSignalColor, cex = 0.5)
+    do.call(graphics::smoothScatter, append(list(x = pred, y = res ,
+                                                 ylim = c(0,1), axes = FALSE,
+                    colramp = colorRampPalette(c("white", "darkgrey"))),a))
+    points(pred[defaultCol == 2], res[defaultCol == 2],
+            col = .Options$DHARMaSignalColor, cex = 0.5)
 
     axis(1)
     axis(2, at=c(0, quantiles, 1))
@@ -266,14 +273,18 @@ plotResiduals <- function(simulationOutput, form = NULL, quantreg = NULL, rank =
       title(main = main, cex.main = 1)
       abline(h = quantiles, col = "black", lwd = 0.5, lty = 2)
       try({
-        lines(smooth.spline(pred, res, df = 10), lty = 2, lwd = 2, col = .Options$DHARMaSignalColor)
+        lines(smooth.spline(pred, res, df = 10), lty = 2, lwd = 2,
+              col = .Options$DHARMaSignalColor)
         abline(h = 0.5, col = .Options$DHARMaSignalColor, lwd = 2)
       }, silent = TRUE)
     }else{
 
       out = testQuantiles(res, pred, quantiles = quantiles, plot = FALSE)
 
-
+      if(is.na(out$p.value)){
+        main = paste(main, "Some quantile regressions failed", sep = "\n")
+        maincol = .Options$DHARMaSignalColor
+      } else{
       if(any(out$pvals < 0.05, na.rm = TRUE)){
         main = paste(main, "Quantile deviations detected (red curves)", sep ="\n")
         if(out$p.value <= 0.05){
@@ -286,7 +297,7 @@ plotResiduals <- function(simulationOutput, form = NULL, quantreg = NULL, rank =
         main = paste(main, "No significant problems detected", sep ="\n")
         maincol = "black"
       }
-
+    }
 
       title(main = main, cex.main = 0.8,
             col.main = maincol)
@@ -303,8 +314,6 @@ plotResiduals <- function(simulationOutput, form = NULL, quantreg = NULL, rank =
         lines(out$predictions$pred, out$predictions[,2*i], col = lineCol, lwd = 2)
       }
 
-      # legend("bottomright", c(paste("Quantile test: p=", round(out$p.value, digits = 5)), paste("Deviation ", ifelse(out$p.value < 0.05, "significant", "n.s."))), text.col = ifelse(out$p.value < 0.05, .Options$DHARMaSignalColor, "black" ), bty="n")
-
     }
   }
   invisible(out)

diff --git a/DHARMa/R/tests.R b/DHARMa/R/tests.R
@@ -119,7 +119,8 @@ testQuantiles <- function(simulationOutput, predictor = NULL, quantiles = c(0.25
                                                     data = datTemp,
                                                     qu = quantiles[i])), silent = T)
       if(inherits(quantResult, "try-error")){
-        message("Unable to calculate quantile regression for quantile ", quantiles[i], ". Possibly to few (unique) data points / predictions. Will be ommited in plots and significance calculations.")
+        message("\n DHARMa: qgam was unable to calculate quantile regression for quantile ",
+                quantiles[i], ". Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations. \n")
       } else {
         x = summary(quantileFits[[i]])
         pval[i] = min(p.adjust(c(x$p.table[1,4], x$s.table[1,4]), method = "BH")) # correction for test on slope and intercept
@@ -335,14 +336,17 @@ testOutliers <- function(simulationOutput, alternative = c("two.sided", "greater
 #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
 #' @example inst/examples/testsHelp.R
 #' @export
-testCategorical <- function(simulationOutput, catPred, quantiles = c(0.25, 0.5, 0.75), plot = TRUE){
+testCategorical <- function(simulationOutput, catPred,
+                            quantiles = c(0.25, 0.5, 0.75), plot = TRUE){
 
   simulationOutput = ensureDHARMa(simulationOutput, convert = T)
 
   catPred = as.factor(catPred)
   out = list()
 
-  out$uniformity$details = suppressWarnings(by(simulationOutput$scaledResiduals, catPred, ks.test, 'punif', simplify = TRUE))
+  out$uniformity$details = suppressWarnings(by(simulationOutput$scaledResiduals,
+                                               catPred, ks.test, 'punif',
+                                               simplify = TRUE))
   out$uniformity$p.value = rep(NA, nlevels(catPred))
   for(i in 1:nlevels(catPred)) out$uniformity$p.value[i] = out$uniformity$details[[i]]$p.value
   out$uniformity$p.value.cor = p.adjust(out$uniformity$p.value)