general tests for phylolm support

florianhartig · Oct 1, 2024 · 0d4b615 · 0d4b615
1 parent 1e65ada
commit 0d4b615
Show file tree

Hide file tree

Showing 5 changed files with 785 additions and 213 deletions.
diff --git a/DHARMa/R/compatibility.R b/DHARMa/R/compatibility.R
@@ -709,10 +709,10 @@ getObservedResponse.phylolm <- function (object, ...){
 #' @rdname getSimulations
 #' @export
 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"){
@@ -729,7 +729,7 @@ getSimulations.phylolm <- function(object, nsim = 1, type = c("normal", "refit")
 getRefit.phylolm <- function(object, newresp, ...){
   newData <- model.frame(object)
   newData[,1] = newresp
-  refittedModel = update(object, data = newData)
+  refittedModel = update(object, data = newData, ...)
 }
 
 
@@ -761,10 +761,11 @@ 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"){
@@ -786,7 +787,7 @@ getRefit.phyloglm <- function(object, newresp, ...){
   names(newData) <- terms
   newData[,1] = newresp
 
-  refittedModel = update(object, data = newData)
+  refittedModel = update(object, data = newData, ...)
   return(refittedModel)
 }
 

diff --git a/DHARMa/R/simulateResiduals.R b/DHARMa/R/simulateResiduals.R
@@ -41,7 +41,10 @@
 #' @example inst/examples/simulateResidualsHelp.R
 #' @import stats
 #' @export
-simulateResiduals <- function(fittedModel, n = 250, refit = FALSE, integerResponse = NULL, plot = FALSE, seed = 123, method = c("PIT", "traditional"), rotation = NULL, ...){
+simulateResiduals <- function(fittedModel, n = 250, refit = FALSE,
+                              integerResponse = NULL, plot = FALSE, seed = 123,
+                              method = c("PIT", "traditional"), rotation = NULL,
+                              ...){
 
   ######## general assertions and startup calculations ##########
 
@@ -88,11 +91,15 @@ simulateResiduals <- function(fittedModel, n = 250, refit = FALSE, integerRespon
 
     out$method = method
 
-    out$simulatedResponse = getSimulations(fittedModel, nsim = n, type = "normal", ...)
+    out$simulatedResponse <- getSimulations(fittedModel, nsim = n,
+                                           type = "normal", ...)
 
     checkSimulations(out$simulatedResponse, out$nObs, out$nSim)
 
-    out$scaledResiduals = getQuantile(simulations = out$simulatedResponse , observed = out$observedResponse , integerResponse = integerResponse, method = method, rotation = rotation)
+    out$scaledResiduals <- getQuantile(simulations = out$simulatedResponse ,
+                                      observed = out$observedResponse ,
+                                      integerResponse = integerResponse,
+                                      method = method, rotation = rotation)
 
     ######## refit = T ##################
   } else {
@@ -102,30 +109,33 @@ simulateResiduals <- function(fittedModel, n = 250, refit = FALSE, integerRespon
     # Adding new outputs
 
     out$refittedPredictedResponse <- matrix(nrow = out$nObs, ncol = n )
-    out$refittedFixedEffects <- matrix(nrow = length(out$fittedFixedEffects), ncol = n )
+    out$refittedFixedEffects <- matrix(nrow = length(out$fittedFixedEffects),
+                                       ncol = n )
     #out$refittedRandomEffects <- matrix(nrow = length(out$fittedRandomEffects), ncol = n )
-    out$refittedResiduals = matrix(nrow = out$nObs, ncol = n)
-    out$refittedPearsonResiduals = matrix(nrow = out$nObs, ncol = n)
+    out$refittedResiduals <- matrix(nrow = out$nObs, ncol = n)
+    out$refittedPearsonResiduals <- matrix(nrow = out$nObs, ncol = n)
 
-    out$simulatedResponse = getSimulations(fittedModel, nsim = n, type = "refit", ...)
+    out$simulatedResponse <- getSimulations(fittedModel, nsim = n,
+                                           type = "refit", ...)
 
     for (i in 1:n){
 
-      simObserved = out$simulatedResponse[[i]]
+      simObserved <- out$simulatedResponse[[i]]
 
       try({
 
         # for testing
         # if (i==3) stop("x")
         # Note: also set silent = T for production
 
-        refittedModel = getRefit(fittedModel, simObserved)
+        refittedModel <- getRefit(fittedModel, simObserved, ...)
 
-        out$refittedPredictedResponse[,i] = getFitted(refittedModel)
-        out$refittedFixedEffects[,i] = getFixedEffects(refittedModel)
-        out$refittedResiduals[,i] = getResiduals(refittedModel)
-        out$refittedPearsonResiduals[,i] = residuals(refittedModel, type = "pearson")
-        #out$refittedRandomEffects[,i]  = ranef(refittedModel)
+        out$refittedPredictedResponse[,i] <- getFitted(refittedModel)
+        out$refittedFixedEffects[,i] <- getFixedEffects(refittedModel)
+        out$refittedResiduals[,i] <- getResiduals(refittedModel)
+        out$refittedPearsonResiduals[,i] <- residuals(refittedModel,
+                                                     type = "pearson")
+        #out$refittedRandomEffects[,i]  <- ranef(refittedModel)
       }, silent = TRUE)
     }
 
@@ -147,7 +157,11 @@ simulateResiduals <- function(fittedModel, n = 250, refit = FALSE, integerRespon
 
     ######### residual calculations ###########
 
-    out$scaledResiduals = getQuantile(simulations = out$refittedResiduals, observed = out$fittedResiduals, integerResponse = integerResponse, method = "traditional", rotation = rotation)
+    out$scaledResiduals = getQuantile(simulations = out$refittedResiduals,
+                                      observed = out$fittedResiduals,
+                                      integerResponse = integerResponse,
+                                      method = "traditional",
+                                      rotation = rotation)
   }
 
   ########### Wrapup ############