diff --git a/DHARMa/R/compatibility.R b/DHARMa/R/compatibility.R
index 650ed594..42a62dac 100644
--- 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
index e9fc57fa..b2368b0f 100644
--- 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,16 +109,18 @@ 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({
 
@@ -119,13 +128,14 @@ simulateResiduals <- function(fittedModel, n = 250, refit = FALSE, integerRespon
         # 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 ############
diff --git a/DHARMa/tests/testthat/_snaps/testTests.new.md b/DHARMa/tests/testthat/_snaps/testTests.new.md
new file mode 100644
index 00000000..d26ddaf3
--- /dev/null
+++ b/DHARMa/tests/testthat/_snaps/testTests.new.md
@@ -0,0 +1,654 @@
+# tests work
+
+    Code
+      testUniformity(simulationOutput, plot = FALSE)
+    Output
+      
+      	Asymptotic one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D = 0.071007, p-value = 0.2655
+      alternative hypothesis: two-sided
+      
+
+---
+
+    Code
+      testUniformity(simulationOutput, plot = FALSE, alternative = "less")
+    Output
+      
+      	Asymptotic one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D^- = 0.071007, p-value = 0.1331
+      alternative hypothesis: the CDF of x lies below the null hypothesis
+      
+
+---
+
+    Code
+      testUniformity(simulationOutput, plot = FALSE, alternative = "greater")
+    Output
+      
+      	Asymptotic one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D^+ = 0.017586, p-value = 0.8836
+      alternative hypothesis: the CDF of x lies above the null hypothesis
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE)
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.1461, p-value = 0.248
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE, alternative = "less")
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.1461, p-value = 0.876
+      alternative hypothesis: less
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE, alternative = "greater")
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.1461, p-value = 0.124
+      alternative hypothesis: greater
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE, alternative = "two.sided")
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.1461, p-value = 0.248
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testResiduals(simulationOutput, plot = FALSE)
+    Output
+      $uniformity
+      
+      	Asymptotic one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D = 0.071007, p-value = 0.2655
+      alternative hypothesis: two-sided
+      
+      
+      $dispersion
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.1461, p-value = 0.248
+      alternative hypothesis: two.sided
+      
+      
+      $outliers
+      
+      	DHARMa bootstrapped outlier test
+      
+      data:  simulationOutput
+      outliers at both margin(s) = 0, observations = 200, p-value = 1
+      alternative hypothesis: two.sided
+       percent confidence interval:
+       0.00 0.01
+      sample estimates:
+      outlier frequency (expected: 0.0029 ) 
+                                          0 
+      
+      
+
+---
+
+    Code
+      testZeroInflation(simulationOutput, plot = FALSE)
+    Output
+      
+      	DHARMa zero-inflation test via comparison to expected zeros with
+      	simulation under H0 = fitted model
+      
+      data:  simulationOutput
+      ratioObsSim = 1.0501, p-value = 0.616
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testZeroInflation(simulationOutput, plot = FALSE, alternative = "less")
+    Output
+      
+      	DHARMa zero-inflation test via comparison to expected zeros with
+      	simulation under H0 = fitted model
+      
+      data:  simulationOutput
+      ratioObsSim = 1.0501, p-value = 0.744
+      alternative hypothesis: less
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = countOnes, plot = FALSE)
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = 1.0037, p-value = 1
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = countOnes, plot = FALSE, alternative = "less")
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = 1.0037, p-value = 0.532
+      alternative hypothesis: less
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = means, plot = FALSE)
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = 1.0098, p-value = 0.904
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = spread, plot = FALSE)
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = 1.068, p-value = 0.304
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testUniformity(simulationOutput, plot = FALSE)
+    Output
+      
+      	Exact one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D = 0.12286, p-value = 0.9932
+      alternative hypothesis: two-sided
+      
+
+---
+
+    Code
+      testUniformity(simulationOutput, plot = FALSE, alternative = "less")
+    Output
+      
+      	Exact one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D^- = 0.10586, p-value = 0.7447
+      alternative hypothesis: the CDF of x lies below the null hypothesis
+      
+
+---
+
+    Code
+      testUniformity(simulationOutput, plot = FALSE, alternative = "greater")
+    Output
+      
+      	Exact one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D^+ = 0.12286, p-value = 0.685
+      alternative hypothesis: the CDF of x lies above the null hypothesis
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE)
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.0862, p-value = 0.8
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE, alternative = "less")
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.0862, p-value = 0.6
+      alternative hypothesis: less
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE, alternative = "greater")
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.0862, p-value = 0.4
+      alternative hypothesis: greater
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE, alternative = "two.sided")
+    Output
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.0862, p-value = 0.8
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testResiduals(simulationOutput, plot = FALSE)
+    Output
+      $uniformity
+      
+      	Exact one-sample Kolmogorov-Smirnov test
+      
+      data:  simulationOutput$scaledResiduals
+      D = 0.12286, p-value = 0.9932
+      alternative hypothesis: two-sided
+      
+      
+      $dispersion
+      
+      	DHARMa nonparametric dispersion test via sd of residuals fitted vs.
+      	simulated
+      
+      data:  simulationOutput
+      dispersion = 1.0862, p-value = 0.8
+      alternative hypothesis: two.sided
+      
+      
+      $outliers
+      
+      	DHARMa bootstrapped outlier test
+      
+      data:  simulationOutput
+      outliers at both margin(s) = 0, observations = 200, p-value = 1
+      alternative hypothesis: two.sided
+       percent confidence interval:
+       0.0 0.1
+      sample estimates:
+      outlier frequency (expected: 0.004 ) 
+                                         0 
+      
+      
+
+---
+
+    Code
+      testZeroInflation(simulationOutput, plot = FALSE)
+    Output
+      
+      	DHARMa zero-inflation test via comparison to expected zeros with
+      	simulation under H0 = fitted model
+      
+      data:  simulationOutput
+      ratioObsSim = NaN, p-value = 1
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testZeroInflation(simulationOutput, plot = FALSE, alternative = "less")
+    Output
+      
+      	DHARMa zero-inflation test via comparison to expected zeros with
+      	simulation under H0 = fitted model
+      
+      data:  simulationOutput
+      ratioObsSim = NaN, p-value = 1
+      alternative hypothesis: less
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = countOnes, plot = FALSE)
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = NaN, p-value = 1
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = countOnes, plot = FALSE, alternative = "less")
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = NaN, p-value = 1
+      alternative hypothesis: less
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = means, plot = FALSE)
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = 1.0098, p-value = 0.904
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testGeneric(simulationOutput, summary = spread, plot = FALSE)
+    Output
+      
+      	DHARMa generic simulation test
+      
+      data:  simulationOutput
+      ratioObsSim = 1.1306, p-value = 0.52
+      alternative hypothesis: two.sided
+      
+
+---
+
+    Code
+      testDispersion(simulationOutput, plot = FALSE)
+    Output
+      
+      	DHARMa nonparametric dispersion test via mean deviance residual fitted
+      	vs. simulated-refitted
+      
+      data:  simulationOutput
+      dispersion = 1.0869, p-value = 0.344
+      alternative hypothesis: two.sided
+      
+
+# correlation tests work
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, x = testData$x, y = testData$y,
+      plot = FALSE)
+    Output
+      
+      	DHARMa Moran's I test for distance-based autocorrelation
+      
+      data:  simulationOutput
+      observed = 0.0082334, expected = -0.0050251, sd = 0.0119231, p-value =
+      0.2661
+      alternative hypothesis: Distance-based autocorrelation
+      
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, x = testData$x, y = testData$y,
+      plot = FALSE, alternative = "two.sided")
+    Output
+      
+      	DHARMa Moran's I test for distance-based autocorrelation
+      
+      data:  simulationOutput
+      observed = 0.0082334, expected = -0.0050251, sd = 0.0119231, p-value =
+      0.2661
+      alternative hypothesis: Distance-based autocorrelation
+      
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, distMat = dM, plot = FALSE)
+    Output
+      
+      	DHARMa Moran's I test for distance-based autocorrelation
+      
+      data:  simulationOutput
+      observed = 0.0082334, expected = -0.0050251, sd = 0.0119231, p-value =
+      0.2661
+      alternative hypothesis: Distance-based autocorrelation
+      
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, distMat = dM, plot = FALSE,
+        alternative = "two.sided")
+    Output
+      
+      	DHARMa Moran's I test for distance-based autocorrelation
+      
+      data:  simulationOutput
+      observed = 0.0082334, expected = -0.0050251, sd = 0.0119231, p-value =
+      0.2661
+      alternative hypothesis: Distance-based autocorrelation
+      
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, plot = FALSE, x = testData$x[1:10],
+      y = testData$y[1:9])
+    Condition
+      Warning in `cbind()`:
+      number of rows of result is not a multiple of vector length (arg 2)
+      Error in `testSpatialAutocorrelation()`:
+      ! Dimensions of x / y coordinates don't match the dimension of the residuals
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput[1:10], plot = FALSE, x = testData$x[
+        1:10], y = testData$y[1:10])
+    Condition
+      Error in `ensureDHARMa()`:
+      ! wrong argument to function, simulationOutput must be a DHARMa object or a numeric vector of quantile residuals!
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, distMat = dM, plot = FALSE, x = testData$
+        x)
+    Message
+      both coordinates and distMat provided, calculations will be done based on the distance matrix, coordinates will only be used for plotting
+    Output
+      
+      	DHARMa Moran's I test for distance-based autocorrelation
+      
+      data:  simulationOutput
+      observed = 0.0082334, expected = -0.0050251, sd = 0.0119231, p-value =
+      0.2661
+      alternative hypothesis: Distance-based autocorrelation
+      
+
+---
+
+    Code
+      testSpatialAutocorrelation(simulationOutput, distMat = dM, plot = FALSE, y = testData$
+        y)
+    Message
+      both coordinates and distMat provided, calculations will be done based on the distance matrix, coordinates will only be used for plotting
+    Output
+      
+      	DHARMa Moran's I test for distance-based autocorrelation
+      
+      data:  simulationOutput
+      observed = 0.0082334, expected = -0.0050251, sd = 0.0119231, p-value =
+      0.2661
+      alternative hypothesis: Distance-based autocorrelation
+      
+
+---
+
+    Code
+      testTemporalAutocorrelation(simulationOutput, plot = FALSE, time = testData$
+      time)
+    Output
+      
+      	Durbin-Watson test
+      
+      data:  simulationOutput$scaledResiduals ~ 1
+      DW = 2.0905, p-value = 0.5202
+      alternative hypothesis: true autocorrelation is not 0
+      
+
+---
+
+    Code
+      testTemporalAutocorrelation(simulationOutput, plot = FALSE, time = testData$
+      time, alternative = "greater")
+    Output
+      
+      	Durbin-Watson test
+      
+      data:  simulationOutput$scaledResiduals ~ 1
+      DW = 2.0905, p-value = 0.7399
+      alternative hypothesis: true autocorrelation is greater than 0
+      
+
+# test phylogenetic autocorrelation
+
+    Code
+      restest
+    Output
+      
+      	DHARMa Moran's I test for phylogenetic autocorrelation
+      
+      data:  res
+      observed = 0.851667, expected = -0.016949, sd = 0.088733, p-value <
+      2.2e-16
+      alternative hypothesis: Phylogenetic autocorrelation
+      
+
+# testOutliers
+
+    Code
+      testOutliers(simulationOutput, plot = F, margin = "lower")
+    Output
+      
+      	DHARMa outlier test based on exact binomial test with approximate
+      	expectations
+      
+      data:  simulationOutput
+      outliers at lower margin(s) = 4, observations = 1000, p-value = 0.8037
+      alternative hypothesis: true probability of success is not equal to 0.003984064
+      95 percent confidence interval:
+       0.001090908 0.010209665
+      sample estimates:
+      frequency of outliers (expected: 0.00398406374501992 ) 
+                                                       0.004 
+      
+
+---
+
+    Code
+      testOutliers(simulationOutput, plot = F, alternative = "two.sided", margin = "lower")
+    Output
+      
+      	DHARMa outlier test based on exact binomial test with approximate
+      	expectations
+      
+      data:  simulationOutput
+      outliers at lower margin(s) = 4, observations = 1000, p-value = 0.8037
+      alternative hypothesis: true probability of success is not equal to 0.003984064
+      95 percent confidence interval:
+       0.001090908 0.010209665
+      sample estimates:
+      frequency of outliers (expected: 0.00398406374501992 ) 
+                                                       0.004 
+      
+
+---
+
+    Code
+      testOutliers(simulationOutput, plot = F, margin = "upper")
+    Output
+      
+      	DHARMa outlier test based on exact binomial test with approximate
+      	expectations
+      
+      data:  simulationOutput
+      outliers at upper margin(s) = 4, observations = 1000, p-value = 0.8037
+      alternative hypothesis: true probability of success is not equal to 0.003984064
+      95 percent confidence interval:
+       0.001090908 0.010209665
+      sample estimates:
+      frequency of outliers (expected: 0.00398406374501992 ) 
+                                                       0.004 
+      
+
diff --git a/DHARMa/tests/testthat/testModelPhylolm.R b/DHARMa/tests/testthat/testModelPhylolm.R
deleted file mode 100644
index 99105354..00000000
--- a/DHARMa/tests/testthat/testModelPhylolm.R
+++ /dev/null
@@ -1,187 +0,0 @@
-
-skip_on_cran()
-skip_on_ci()
-
-set.seed(1234)
-
-doPlots = F
-
-# Test functions --------------------------------------------------------------
-checkOutput <- function(simulationOutput){
-
-  # print(simulationOutput)
-
-  if(any(simulationOutput$scaledResiduals < 0)) stop()
-  if(any(simulationOutput$scaledResiduals > 1)) stop()
-  if(any(is.na(simulationOutput$scaledResiduals))) stop()
-
-  if(length(simulationOutput$scaledResiduals) !=
-     length(simulationOutput$observedResponse)) stop()
-  if(length(simulationOutput$fittedPredictedResponse) != length(simulationOutput$observedResponse)) stop()
-
-}
-
-expectDispersion <- function(x, answer = T){
-  res <- simulateResiduals(x)
-  if (answer) expect_lt(testDispersion(res, plot = doPlots)$p.value, 0.05)
-  else expect_gt(testDispersion(res, plot = doPlots)$p.value, 0.05)
-}
-
-# pylolm/phyloglms works -------------------------------------------------------
-#
-# test_that("phylolm works",
-#           {
-#             # LM
-#             set.seed(123456)
-#             tre = ape::rcoal(60)
-#             taxa = sort(tre$tip.label)
-#             b0 = 0; b1 = 1;
-#             x <- phylolm::rTrait(n = 1, phy = tre, model = "BM",
-#                                  parameters = list(ancestral.state = 0,
-#                                                    sigma2 = 10))
-#             y <- b0 + b1*x +
-#               phylolm::rTrait(n = 1, phy = tre, model = "lambda",
-#                               parameters = list(ancestral.state = 0, sigma2 = 1,
-#                                                 lambda = 0.5))
-#             testData = data.frame(trait = y[taxa], predictor = x[taxa],
-#                                   x = runif(length(y)), y= runif(length(y)))
-#             fittedModel = phylolm::phylolm(trait ~ predictor, data = testData,
-#                                            phy = tre, model = "lambda")
-#
-#             t = getObservedResponse(fittedModel)
-#             expect_true(is.vector(t))
-#             expect_true(is.numeric(t))
-#
-#             x = getSimulations(fittedModel, 1)
-#             expect_true(is.matrix(x))
-#             expect_true(ncol(x) == 1)
-#
-#             x = getSimulations(fittedModel, 2)
-#             expect_true(is.numeric(x))
-#             expect_true(is.matrix(x))
-#             expect_true(ncol(x) == 2)
-#
-#             x = getSimulations(fittedModel, 1, type = "refit")
-#             expect_true(is.data.frame(x))
-#
-#             x = getSimulations(fittedModel, 2, type = "refit")
-#             expect_true(is.data.frame(x))
-#
-#             fittedModel2 = getRefit(fittedModel,x[[1]])
-#             expect_false(any(getFixedEffects(fittedModel) -
-#                              getFixedEffects(fittedModel2) > 0.5)) # doesn't work for some models
-#
-#             simulationOutput <- simulateResiduals(fittedModel = fittedModel,  n = 200)
-#
-#             checkOutput(simulationOutput)
-#
-#             if(doPlots) plot(simulationOutput, quantreg = F)
-#
-#             expect_gt(testOutliers(simulationOutput, plot = doPlots)$p.value, 0.001)
-#             expect_gt(testDispersion(simulationOutput, plot = doPlots)$p.value, 0.001)
-#             expect_gt(testUniformity(simulationOutput = simulationOutput,
-#                                      plot = doPlots)$p.value, 0.001)
-#             expect_gt(testZeroInflation(simulationOutput = simulationOutput,
-#                                         plot = doPlots)$p.value, 0.001)
-#             expect_gt(testTemporalAutocorrelation(simulationOutput = simulationOutput,
-#                                                   time = testData$time,
-#                                                   plot = doPlots)$p.value, 0.001)
-#             expect_gt(testSpatialAutocorrelation(simulationOutput = simulationOutput,
-#                                                  x = testData$x, y = testData$y,
-#                                                  plot = F)$p.value, 0.001)
-#
-#             simulationOutput <- recalculateResiduals(simulationOutput, group = testData$group)
-#             expect_gt(testDispersion(simulationOutput, plot = doPlots)$p.value, 0.001)
-#
-#             simulationOutput2 <- simulateResiduals(fittedModel = fittedModel,
-#                                                    refit = T, n = 100)
-#
-#             checkOutput(simulationOutput2)
-#             if(doPlots) plot(simulationOutput2, quantreg = F)
-#
-#             # note that the pearson test is biased, therefore have to test greater
-#             #expect_gt(testDispersion(simulationOutput2, plot = doPlots, alternative = "greater")$p.value, 0.001)
-#             x = testDispersion(simulationOutput2, plot = doPlots)
-#
-#             simulationOutput3 <- recalculateResiduals(simulationOutput2, group = testData$group)
-#             #expect_gt(testDispersion(simulationOutput3, plot = doPlots, alternative = "greater")$p.value, 0.001)
-#             x = testDispersion(simulationOutput3, plot = doPlots)
-#           })
-
-test_that("phyloglm works",
-          {
-            #GLM
-            set.seed(123456)
-            tre = ape::rtree(50)
-            x = phylolm::rTrait(n = 1, phy = tre)
-            X = cbind(rep(1, 50), x)
-            y = phylolm::rbinTrait(n = 1, phy = tre, beta = c(-1,0.5), alpha = 1,
-                                   X = X)
-            testData = data.frame(trait = y, predictor = x,
-                                  x = runif(length(y)), y = runif(length(y)))
-            fittedModel = phylolm::phyloglm(trait ~ predictor, btol=20,
-                                            phy = tre, data = testData)
-
-            t = getObservedResponse(fittedModel)
-            expect_true(is.vector(t))
-            expect_true(is.numeric(t))
-
-            x = getSimulations(fittedModel, 1)
-            expect_true(is.matrix(x))
-            expect_true(ncol(x) == 1)
-
-            x = getSimulations(fittedModel, 2)
-            expect_true(is.numeric(x))
-            expect_true(is.matrix(x))
-            expect_true(ncol(x) == 2)
-
-            x = getSimulations(fittedModel, 1, type = "refit")
-            expect_true(is.data.frame(x))
-
-            x = getSimulations(fittedModel, 2, type = "refit")
-            expect_true(is.data.frame(x))
-
-            fittedModel2 = getRefit(fittedModel,x[[1]])
-            expect_false(any(getFixedEffects(fittedModel) -
-                               getFixedEffects(fittedModel2) > 0.5)) # doesn't work for some models
-
-            simulationOutput <- simulateResiduals(fittedModel = fittedModel,  n = 200)
-
-            checkOutput(simulationOutput)
-
-            if(doPlots) plot(simulationOutput, quantreg = F)
-
-            expect_gt(testOutliers(simulationOutput, plot = doPlots)$p.value, 0.001)
-            expect_gt(testDispersion(simulationOutput, plot = doPlots)$p.value, 0.001)
-            expect_gt(testUniformity(simulationOutput = simulationOutput,
-                                     plot = doPlots)$p.value, 0.001)
-            expect_gt(testZeroInflation(simulationOutput = simulationOutput,
-                                        plot = doPlots)$p.value, 0.001)
-            expect_gt(testTemporalAutocorrelation(simulationOutput = simulationOutput,
-                                                  time = testData$time,
-                                                  plot = doPlots)$p.value, 0.001)
-            expect_gt(testSpatialAutocorrelation(simulationOutput = simulationOutput,
-                                                 x = testData$x, y = testData$y,
-                                                 plot = F)$p.value, 0.001)
-
-            simulationOutput <- recalculateResiduals(simulationOutput, group = testData$group)
-            expect_gt(testDispersion(simulationOutput, plot = doPlots)$p.value, 0.001)
-
-            simulationOutput2 <- simulateResiduals(fittedModel = fittedModel,
-                                                   refit = T, n = 100)
-
-            checkOutput(simulationOutput2)
-            if(doPlots) plot(simulationOutput2, quantreg = F)
-
-            # note that the pearson test is biased, therefore have to test greater
-            #expect_gt(testDispersion(simulationOutput2, plot = doPlots, alternative = "greater")$p.value, 0.001)
-            x = testDispersion(simulationOutput2, plot = doPlots)
-
-            simulationOutput3 <- recalculateResiduals(simulationOutput2, group = testData$group)
-            #expect_gt(testDispersion(simulationOutput3, plot = doPlots, alternative = "greater")$p.value, 0.001)
-            x = testDispersion(simulationOutput3, plot = doPlots)
-
-
-
-
-          })
diff --git a/DHARMa/tests/testthat/testModelTypes.R b/DHARMa/tests/testthat/testModelTypes.R
index b04d532a..bac7cfd6 100644
--- a/DHARMa/tests/testthat/testModelTypes.R
+++ b/DHARMa/tests/testthat/testModelTypes.R
@@ -1,8 +1,8 @@
 
 skip_on_cran()
-skip_on_ci()
+#skip_on_ci()
 
-set.seed(123)
+set.seed(1234)
 
 doPlots = F
 
@@ -27,7 +27,7 @@ expectDispersion <- function(x, answer = T){
   else expect_gt(testDispersion(res, plot = doPlots)$p.value, 0.05)
 }
 
-runEverything = function(fittedModel, testData, DHARMaData = T,
+runEverything = function(fittedModel, testData, DHARMaData = T, phy = NULL,
                          expectOverdispersion = F){
 
   t = getObservedResponse(fittedModel)
@@ -122,6 +122,48 @@ testData$binomial_nk_weights2$prop = testData$binomial_nk_weights2$observedRespo
 
 
 
+testData$poisson1 = createData(sampleSize = 500, overdispersion = 0,
+                               randomEffectVariance = 0.000, family = poisson())
+testData$poisson2 = createData(sampleSize = 200, overdispersion = 2,
+                               randomEffectVariance =0.000, family = poisson())
+testData$poisson3 = createData(sampleSize = 500, overdispersion = 0.5,
+                               randomEffectVariance = 0.000, family = poisson())
+
+testData$poisson_weights = createData(sampleSize = 200, overdispersion = 0.5,
+                                      randomEffectVariance = 0.5, family = poisson())
+testData$weights = rep(c(1,1.1), each = 100)
+testData$poisson_weights$weights = testData$weights
+
+# testData --------------------------------------------------------------------
+testData = list()
+testData$lm = createData(sampleSize = 200, fixedEffects = c(1,0),
+                         overdispersion = 0, randomEffectVariance = 0,
+                         family = gaussian())
+testData$lmm = createData(sampleSize = 200,
+                          overdispersion = 0, randomEffectVariance = 0.5,
+                          family = gaussian())
+
+
+testData$binomial_10 = createData(sampleSize = 200, randomEffectVariance = 0,
+                                  family = binomial())
+testData$binomial_yn = createData(sampleSize = 200, fixedEffects = c(1,0),
+                                  overdispersion = 0, randomEffectVariance = 0,
+                                  family = binomial(), factorResponse = T)
+testData$binomial_nk_matrix = createData(sampleSize = 200, overdispersion = 0,
+                                         randomEffectVariance = 0, family = binomial(),
+                                         binomialTrials = 20)
+testData$binomial_nk_weights = createData(sampleSize = 200, overdispersion = 0,
+                                          randomEffectVariance = 0, family = binomial(),
+                                          binomialTrials = 20)
+testData$binomial_nk_weights$prop = testData$binomial_nk_weights$observedResponse1 / 20
+
+testData$binomial_nk_weights2 = createData(sampleSize = 200, overdispersion = 1,
+                                           randomEffectVariance = 0, family = binomial(),
+                                           binomialTrials = 20)
+testData$binomial_nk_weights2$prop = testData$binomial_nk_weights2$observedResponse1 / 20
+
+
+
 testData$poisson1 = createData(sampleSize = 500, overdispersion = 0,
                                randomEffectVariance = 0.000, family = poisson())
 testData$poisson2 = createData(sampleSize = 200, overdispersion = 2,
@@ -399,7 +441,7 @@ test_that("glmmTMB works",
 )
 
 
-# spaMM::HLfit works  --------------------------------------------------------------
+# spaMM::HLfit  --------------------------------------------------------------
 
 test_that("spaMM::HLfit works",
           {
@@ -452,7 +494,7 @@ test_that("spaMM::HLfit works",
 )
 
 
-# GLMMadaptive works  --------------------------------------------------------------
+# GLMMadaptive  --------------------------------------------------------------
 
 test_that("GLMMadaptive works",
           {
@@ -511,6 +553,54 @@ test_that("GLMMadaptive works",
           }
 )
 
+rm(testData)
+
+
+# phylolm::pylolm/phyloglm -----------------------------------------------
+
+# OBS: somehow, the base function update() in getSimulations.phylolm() and
+# getRefit.phylolm() searches for the objects in the global env, then the
+# test_that wasn't working. A workaround was to put the objects of the tree and
+# the dataset in the global env.
+
+test_that("phylolm works",
+          {
+            # LM
+            set.seed(123456)
+            tre1 <<- ape::rcoal(60) # global env
+            taxa = sort(tre1$tip.label)
+            b0 = 0; b1 = 1;
+            x <- phylolm::rTrait(n = 1, phy = tre1, model = "BM",
+                                 parameters = list(ancestral.state = 0,
+                                                   sigma2 = 10))
+            y <- b0 + b1*x +
+              phylolm::rTrait(n = 1, phy = tre1, model = "lambda",
+                              parameters = list(ancestral.state = 0, sigma2 = 1,
+                                                lambda = 0.5))
+            testData <<- data.frame(trait = y[taxa], predictor = x[taxa],
+                                    x = runif(length(y)), y= runif(length(y)))
+            fittedModel = phylolm::phylolm(trait ~ predictor, data = testData,
+                                           phy = tre1, model = "lambda")
+
+            runEverything(fittedModel, testData = testData, phy = tre1)
+          })
+
+test_that("phyloglm works",
+          {
+            #GLM
+            set.seed(123456)
+            tre <<- ape::rtree(50) # global env
+            x = phylolm::rTrait(n = 1, phy = tre)
+            X = cbind(rep(1, 50), x)
+            y = phylolm::rbinTrait(n = 1, phy = tre, beta = c(-1,0.5), alpha = 1,
+                                   X = X)
+            testData <<- data.frame(trait = y, predictor = x,
+                                    x = runif(length(y)), y = runif(length(y)))
+            fittedModel = phylolm::phyloglm(trait ~ predictor,
+                                            phy = tre, data = testData)
+
+            runEverything(fittedModel,  testData = testData, phy = tre)
+          })
 
 
 # isNA works? ------------------------------------------------------------------