Merge pull request #441 from florianhartig/MelinaCheckPackageSupport

Melina check package support

.github/workflows/R-CMD-check.yaml

@@ -49,13 +49,15 @@ jobs:
           needs: check
           working-directory: DHARMa
       # second check with different rcmdcheck arguments.
+
       - uses: r-lib/actions/check-r-package@v2
         with:
           upload-snapshots: true
           working-directory: DHARMa
           args: 'c("--no-multiarch", "--no-manual")'
+          build_args: c('--compact-vignettes=gs+qpdf')
 
       - uses: r-lib/actions/check-r-package@v2
         with:
-          upload-snapshots: true
+          upload-snapshots: false
           working-directory: DHARMa

DHARMa/DESCRIPTION

@@ -6,8 +6,9 @@ Authors@R: c(person("Florian", "Hartig", email = "[email protected]
 Description: The 'DHARMa' package uses a simulation-based approach to create
     readily interpretable scaled (quantile) residuals for fitted (generalized) linear mixed
     models. Currently supported are linear and generalized linear (mixed) models from 'lme4'
-    (classes 'lmerMod', 'glmerMod'), 'glmmTMB' 'GLMMadaptive' and 'spaMM', generalized additive models
-    ('gam' from 'mgcv'), 'glm' (including 'negbin' from 'MASS', but excluding quasi-distributions) and
+    (classes 'lmerMod', 'glmerMod'), 'glmmTMB', 'GLMMadaptive' and 'spaMM'; phylogenetic 
+    linear models from 'phylolm'(classes 'phylolm' and 'phyloglm'); generalized additive 
+    models ('gam' from 'mgcv'); 'glm' (including 'negbin' from 'MASS', but excluding quasi-distributions) and
     'lm' model classes. Moreover, externally created simulations, e.g. posterior predictive simulations
     from Bayesian software such as 'JAGS', 'STAN', or 'BUGS' can be processed as well.
     The resulting residuals are standardized to values between 0 and 1 and can be interpreted
@@ -30,16 +31,17 @@ Imports:
     lme4
 Suggests:
     knitr,
-    testthat,
+    testthat (>= 3.0.0),
     rmarkdown,
     KernSmooth,
     sfsmisc,
     MASS,
     mgcv,
-    mgcViz (>= 	0.1.9),
+    mgcViz (>= 0.1.9),
     spaMM (>= 3.2.0),
     GLMMadaptive,
-    glmmTMB (>= 1.1.2.3)
+    glmmTMB (>= 1.1.2.3),
+    phylolm (>= 2.6.5)
 Enhances:
     phyr,
     rstan,
@@ -49,7 +51,8 @@ License: GPL (>= 3)
 URL: http://florianhartig.github.io/DHARMa/
 LazyData: TRUE
 BugReports: https://github.com/florianhartig/DHARMa/issues
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Roxygen: list(old_usage = TRUE, markdown = TRUE)
 VignetteBuilder: knitr
 Encoding: UTF-8
+Config/testthat/edition: 3

DHARMa/NAMESPACE

@@ -22,6 +22,8 @@ S3method(getRefit,MixMod)
 S3method(getRefit,default)
 S3method(getRefit,glmmTMB)
 S3method(getRefit,lm)
+S3method(getRefit,phyloglm)
+S3method(getRefit,phylolm)
 S3method(getResiduals,MixMod)
 S3method(getResiduals,default)
 S3method(getSimulations,HLfit)

DHARMa/NEWS.md

@@ -1,6 +1,23 @@
 NOTE: for more news about the package, see https://github.com/florianhartig/DHARMa/releases
 
 
+# DHARMa 0.4.7
+
+## New Features
+
+* Includes support for the package phylolm
+* new function for testing residual phylogenetic autocorrelation testPylogeneticAutocorrelation()
+
+## Major change
+
+* new optional argument 'rotation' in simulateResiduals() for the rotation of the residual space prior to calculating the quantile residuals to account for residual covariance as created by temporal, spatial or phylogenetic autocorrelation. 
+
+## Minor changes
+
+* color bind friendly plots as default
+* documentation: improved help files and vignette
+
+
 # DHARMa 0.4.6
 
 ## New Features
@@ -38,7 +55,7 @@ NOTE: for more news about the package, see https://github.com/florianhartig/DHAR
 * re-introduced glmmTMB to suggests
 * phyr moved to enhances
 * re-modelled package unit tests
-* added RStan, CmdStanR, rjags, BayesianTools to enhances, as they coudld be used with DHARMa
+* added RStan, CmdStanR, rjags, BayesianTools to enhances, as they could be used with DHARMa
 * moved parallel calculations in runBenchmark to R native parallel functions
 
 ## New features
@@ -169,7 +186,6 @@ This is actually a bugfix release for 0.3.4, but on reflection I decided that 0.
 ## New features
 
 * added smooth scatter in plotResiduals https://github.com/florianhartig/DHARMa/commit/da01d8c7a9a74558817e4a73fe826084164cf05d
-
 * glmmTMB now fully supported through new compulsory version 1.0 of glmmTMB, which includes the re.form argument in the simulations required by DHARMa https://github.com/florianhartig/DHARMa/pull/140
 
 
@@ -306,7 +322,6 @@ This is actually a bugfix release for 0.3.4, but on reflection I decided that 0.
 ## Minor changes
 
 * plotResiduals includes support for factors
-
 * updates to the help
 
 
@@ -318,7 +333,6 @@ This is actually a bugfix release for 0.3.4, but on reflection I decided that 0.
 # DHARMa 0.1.1
 
 * including now the negative binomial models from MASS and lme4, as well as the possibility to create synthetic data from the negative binomial family
-
 * includes a createDHARMa function that allows using the plot functions of DHARMa also with externally created simualtions, for example for Bayesian predictive simulations
 
 

DHARMa/R/compatibility.R

@@ -247,13 +247,13 @@ getObservedResponse.default <- function (object, ...){
   }
 
   if(is.matrix(out)){
-    # case scaled variables or something like that 
+    # case scaled variables or something like that
     if(ncol(out) == 1){
       out = as.vector(out)
     } else if(ncol(out) == 2) {
       # case k/n binomial
       if(!(family(object)$family %in% c("binomial", "betabinomial"))) securityAssertion("nKcase - wrong family")
-      out = out[,1]      
+      out = out[,1]
     }
     else securityAssertion("Response in the model is a matrix with > 2 dim")
   }
@@ -708,10 +708,11 @@ getObservedResponse.phylolm <- function (object, ...){
 
 #' @rdname getSimulations
 #' @export
-getSimulations.phylolm <- function(object, nsim = 1, type = c("normal", "refit"), ...){
+getSimulations.phylolm <- function(object, nsim = 1, type = c("normal", "refit"),
+                                    ...){
   type <- match.arg(type)
 
-  fitBoot = update(object, boot = nsim, save = T)
+  fitBoot = update(object, boot = nsim, save = T, ...)
   out = fitBoot$bootdata
 
   if(type == "normal"){
@@ -723,6 +724,16 @@ getSimulations.phylolm <- function(object, nsim = 1, type = c("normal", "refit")
   return(out)
 }
 
+#' @rdname getRefit
+#' @export
+getRefit.phylolm <- function(object, newresp, ...){
+  newData <- model.frame(object)
+  newData[,1] = newresp
+  refittedModel = update(object, data = newData, ...)
+}
+
+
+
 #' @rdname getFitted
 #' @export
 getFitted.phylolm  <- function (object,...){
@@ -747,12 +758,14 @@ getObservedResponse.phyloglm <- function (object, ...){
 }
 
 
+
 #' @rdname getSimulations
 #' @export
-getSimulations.phyloglm <- function(object, nsim = 1, type = c("normal", "refit"), ...){
+getSimulations.phyloglm <- function(object, nsim = 1,
+                                    type = c("normal", "refit"), ...){
   type <- match.arg(type)
 
-  fitBoot = update(object, boot = nsim, save = T)
+  fitBoot = update(object, boot = nsim, save = T, ...)
   out = fitBoot$bootdata
 
   if(type == "normal"){
@@ -764,24 +777,30 @@ getSimulations.phyloglm <- function(object, nsim = 1, type = c("normal", "refit"
   return(out)
 }
 
-#' @rdname getFitted
+
+#' @rdname getRefit
 #' @export
-getFitted.phyloglm  <- function (object,...){
-  return(object$fitted.values)
+getRefit.phyloglm <- function(object, newresp, ...){
+  #object phyloglm doesn't have a model.frame object
+  terms <-  as.character(formula(object))[-1]
+  newData <- model.frame(object$y ~ object$X[,-1])
+  names(newData) <- terms
+  newData[,1] = newresp
+
+  refittedModel = update(object, data = newData, ...)
+  return(refittedModel)
 }
 
 
 
-#' @rdname getFamily
+#' @rdname getFitted
 #' @export
-getFamily.phyloglm <- function (object,...){
-  out = list()
-  # out$family = object$method
-  out$family = "integer-valued" # all families of phyloglm are integer-valued
-  return(out)
+getFitted.phyloglm  <- function (object,...){
+  return(object$fitted.values)
 }
 
 
+
 #' @rdname getFamily
 #' @export
 getFamily.phyloglm <- function (object,...){
@@ -790,3 +809,4 @@ getFamily.phyloglm <- function (object,...){
   out$family = "integer-valued" # all families of phyloglm are integer-valued
   return(out)
 }
+

DHARMa/R/helper.R

@@ -1,6 +1,6 @@
 #' Modified ECDF function.
 #'
-#' @details ensures symmetric ECDF (standard ECDF is <), and that 0 / 1 values are only produced if the data is strictly < > than the observed data.
+#' @details Ensures symmetric ECDF (standard ECDF is <), and that 0 / 1 values are only produced if the data is strictly < > than the observed data.
 #'
 #' @keywords internal
 DHARMa.ecdf <- function (x)
@@ -20,7 +20,7 @@ DHARMa.ecdf <- function (x)
 
 
 
-#' calculate quantiles
+#' Calculate Residual Quantiles
 #'
 #' Calculates residual quantiles from a given simulation.
 #'
@@ -43,12 +43,14 @@ DHARMa.ecdf <- function (x)
 #'
 #' **Rotation (optional)**
 #'
-#' The getQuantile function includes an additional option to rotate residuals. This option should ONLY be used when the fitted model includes a particular residuals covariance structure, such as an AR1 or a spatial CAR model.
+#' The getQuantile function includes an additional option to rotate residuals prior to calculating the quantile residuals. This option should ONLY be used when the fitted model includes a particular residuals covariance structure, such as an AR1 or a spatial or phylogenetic CAR model.
 #'
 #' For these models, residuals calculated from unconditional simulations will include the specified covariance structure, which will trigger e.g. temporal autocorrelation tests and can inflate type I errors of other tests. The idea of the rotation is to rotate the residual space according to the covariance structure of the fitted model, such that the rotated residuals are conditional independent (provided the fitted model is correct).
 #'
 #' If the residual covariance of the fitted model at the response scale can be extracted (e.g. when fitting gls type models), it would be best to extract it and provide this covariance matrix to the rotation option. If that is not the case, providing the argument "estimated" to rotation will estimate the covariance from the data simulated by the model. This is probably without alternative for GLMMs, where the covariance at the response scale is likely not known / provided, but note, that this approximation will tend to have considerable error and may be slow to compute for high-dimensional data. If you try to estimate the rotation from simulations, you should set n as high as possible! See [testTemporalAutocorrelation] for a practical example.
 #'
+#' The rotation of residuals implemented here is similar to the Variogram.lme() and Variongram.gls() functions in nlme package using the argument resType = "normalized".
+#'
 #' @references
 #'
 #' Smith, J. Q. "Diagnostic checks of non-standard time series models." Journal of Forecasting 4.3 (1985): 283-291.
@@ -58,12 +60,13 @@ DHARMa.ecdf <- function (x)
 #' Warton, David I., Loïc Thibaut, and Yi Alice Wang. "The PIT-trap—A “model-free” bootstrap procedure for inference about regression models with discrete, multivariate responses." PloS one 12.7 (2017).
 #'
 #' @export
-getQuantile <- function(simulations, observed, integerResponse, method = c("PIT", "traditional"), rotation = NULL){
+getQuantile <- function(simulations, observed, integerResponse,
+                        method = c("PIT", "traditional"), rotation = NULL){
 
   method = match.arg(method)
 
   n = length(observed)
-  if (nrow(simulations) != n) stop("DHARMa::getquantile: wrong dimension of simulations")
+  if(nrow(simulations) != n) stop("DHARMa::getquantile: wrong dimension of simulations")
   nSim = ncol(simulations)
 
 
@@ -77,7 +80,7 @@ getQuantile <- function(simulations, observed, integerResponse, method = c("PIT"
 
       values = as.vector(simulations)[duplicated(as.vector(simulations))]
       if(length(values) > 0){
-        if (all(values%%1==0)){
+        if(all(values%%1 == 0)){
           integerResponse = T
           message("Model family was recognized or set as continuous, but duplicate values were detected in the simulation - changing to integer residuals (see ?simulateResiduals for details).")
         } else {