diff --git a/Code/DHARMaPackageSupport/phylolm/phylolm.R b/Code/DHARMaPackageSupport/phylolm/phylolm.R
index 50fa277e..a5be655d 100644
--- a/Code/DHARMaPackageSupport/phylolm/phylolm.R
+++ b/Code/DHARMaPackageSupport/phylolm/phylolm.R
@@ -1,6 +1,31 @@
+
+devtools::install_github("lamho86/phylolm")
+
 library(phylolm)
 
+# GLM
+
+set.seed(123456)
+tre = rtree(50)
+x = rTrait(n=1,phy=tre)
+X = cbind(rep(1,50),x)
+y = rbinTrait(n=1,phy=tre, beta=c(-1,0.5), alpha=1 ,X=X)
+dat = data.frame(trait01 = y, predictor = x)
+fit = phyloglm(trait01~predictor,phy=tre,data=dat,boot=100)
+
+summary(fit)
+coef(fit)
+vcov(fit)
+
+DHARMa:::checkModel(fit, stop = F)
+
+res = DHARMa::simulateResiduals(fit, plot = T)
+
+fitBoot = phylolm:::update(fit, boot = nsim, save = T)
 
+phylolm:::up
+
+# LM 
 
 set.seed(123456)
 tre = rcoal(60)
@@ -15,28 +40,9 @@ dat = data.frame(trait=y[taxa],pred=x[taxa])
 fit = phylolm(trait~pred,data=dat,phy=tre,model="lambda")
 summary(fit)
 
-# adding measurement errors and bootstrap
-z <- y + rnorm(60,0,1)
-dat = data.frame(trait=z[taxa],pred=x[taxa])
-fit = phylolm(trait~pred,data=dat,phy=tre,model="BM",measurement_error=TRUE,boot=100)
-summary(fit)
-
-
+res = DHARMa::simulateResiduals(fit, plot = T)
+testSpatialAutocorrelation(res, )
 
-phyloglm
+res = DHARMa::simulateResiduals(fit, plot = T, rotation = "estimated")
 
 
-
-set.seed(123456)
-tre = rtree(50)
-x = rTrait(n=1,phy=tre)
-X = cbind(rep(1,50),x)
-y = rbinTrait(n=1,phy=tre, beta=c(-1,0.5), alpha=1 ,X=X)
-dat = data.frame(trait01 = y, predictor = x)
-fit = phyloglm(trait01~predictor,phy=tre,data=dat,boot=100)
-summary(fit)
-coef(fit)
-vcov(fit)
-
-predict(fit)
-simulate(fit)
diff --git a/DHARMa/DESCRIPTION b/DHARMa/DESCRIPTION
index bbd05125..6d4bc705 100644
--- a/DHARMa/DESCRIPTION
+++ b/DHARMa/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: DHARMa
 Title: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models
-Version: 0.4.6.x
+Version: 0.4.6.1
 Date: 2022-09-08
 Authors@R: c(person("Florian", "Hartig", email = "florian.hartig@biologie.uni-regensburg.de", role = c("aut", "cre"), comment=c(ORCID="0000-0002-6255-9059")), person("Lukas", "Lohse", role = "ctb"))
 Description: The 'DHARMa' package uses a simulation-based approach to create
diff --git a/DHARMa/NAMESPACE b/DHARMa/NAMESPACE
index 2789fc85..5678a366 100644
--- a/DHARMa/NAMESPACE
+++ b/DHARMa/NAMESPACE
@@ -1,13 +1,20 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(getFamily,default)
+S3method(getFamily,phyloglm)
+S3method(getFamily,phylolm)
 S3method(getFitted,HLfit)
 S3method(getFitted,MixMod)
 S3method(getFitted,default)
 S3method(getFitted,gam)
+S3method(getFitted,phyloglm)
+S3method(getFitted,phylolm)
 S3method(getFixedEffects,MixMod)
 S3method(getFixedEffects,default)
 S3method(getObservedResponse,HLfit)
 S3method(getObservedResponse,default)
+S3method(getObservedResponse,phyloglm)
+S3method(getObservedResponse,phylolm)
 S3method(getPearsonResiduals,default)
 S3method(getPearsonResiduals,gam)
 S3method(getRefit,HLfit)
@@ -24,6 +31,8 @@ S3method(getSimulations,gam)
 S3method(getSimulations,glmmTMB)
 S3method(getSimulations,lmerMod)
 S3method(getSimulations,negbin)
+S3method(getSimulations,phyloglm)
+S3method(getSimulations,phylolm)
 S3method(hist,DHARMa)
 S3method(plot,DHARMa)
 S3method(plot,DHARMaBenchmark)
@@ -32,6 +41,7 @@ S3method(residuals,DHARMa)
 export(benchmarkRuntime)
 export(createDHARMa)
 export(createData)
+export(getFamily)
 export(getFitted)
 export(getFixedEffects)
 export(getObservedResponse)
@@ -56,6 +66,7 @@ export(testGeneric)
 export(testOutliers)
 export(testOverdispersion)
 export(testOverdispersionParametric)
+export(testPhylogeneticAutocorrelation)
 export(testQuantiles)
 export(testResiduals)
 export(testSimulatedResiduals)
diff --git a/DHARMa/R/compatibility.R b/DHARMa/R/compatibility.R
index 69ae989a..90af28f5 100644
--- a/DHARMa/R/compatibility.R
+++ b/DHARMa/R/compatibility.R
@@ -39,6 +39,7 @@ checkModel <- function(fittedModel, stop = F){
   if(!(class(fittedModel)[1] %in% getPossibleModels())){
     if(stop == FALSE) warning("DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work!")
     else stop("DHARMa: fittedModel not in class of supported models")
+    
   }
   
   # if(hasNA(fittedModel)) message("It seems there were NA values in the data used for fitting the model. This can create problems if you supply additional data to DHARMa functions. See ?checkModel for details")
@@ -58,7 +59,7 @@ checkModel <- function(fittedModel, stop = F){
 #' returns a list of supported model classes 
 #' 
 #' @keywords internal
-getPossibleModels<-function()c("lm", "glm", "negbin", "lmerMod", "lmerModLmerTest", "glmerMod", "gam", "bam", "glmmTMB", "HLfit", "MixMod")
+getPossibleModels<-function()c("lm", "glm", "negbin", "lmerMod", "lmerModLmerTest", "glmerMod", "gam", "bam", "glmmTMB", "HLfit", "MixMod", "phylolm", "phyloglm")
 
 
 weightsWarning = "Model was fit with prior weights. These will be ignored in the simulation. See ?getSimulations for details"
@@ -214,6 +215,22 @@ getPearsonResiduals <- function (object, ...) {
   UseMethod("getPearsonResiduals", object)
 }
 
+#' Get model family
+#'
+#' Wrapper to get the family of a fitted model
+#' 
+#' @param object a fitted model
+#' @param ... additional parameters to be passed on
+#'
+#' @seealso \code{\link{getObservedResponse}}, \code{\link{getSimulations}}, \code{\link{getRefit}}, \code{\link{getFixedEffects}}, \code{\link{getFitted}}
+#'
+#' @author Florian Hartig
+#' @export
+getFamily <- function (object, ...) {
+  UseMethod("getFamily", object)
+}
+
+# nObs 
 
 # default -----------------------------------------------------------------
 
@@ -238,9 +255,6 @@ getObservedResponse.default <- function (object, ...){
     
     out = out[,1]
   }
-  
-
-  
   return(out)
 }
 
@@ -343,6 +357,14 @@ getPearsonResiduals.default <- function (object, ...){
 #   if(length(x) < as.numeric(x[length(x) ])) return(TRUE)
 #   else return(FALSE)
 # }
+#
+
+
+#' @rdname getFamily
+#' @export
+getFamily.default <- function (object,...){
+  family(object)
+}
 
 
 
@@ -381,6 +403,7 @@ hasWeigths.lm <- function(object, ...){
 #' @rdname getSimulations
 #' @export
 getSimulations.negbin<- function (object, nsim = 1, type = c("normal", "refit"), ...){
+  type <- match.arg(type)
   if("(weights)" %in% colnames(model.frame(object))) warning(weightsWarning)
   getSimulations.default(object = object, nsim = nsim, type = type, ...)
 }
@@ -512,6 +535,7 @@ getRefit.glmmTMB <- function(object, newresp, ...){
 #' @export
 getSimulations.glmmTMB <- function (object, nsim = 1, type = c("normal", "refit"), ...){
   
+  type <- match.arg(type)
   if("(weights)" %in% colnames(model.frame(object)) & ! family(object)$family %in% c("binomial", "betabinomial")) warning(weightsWarning)
   
   type <- match.arg(type)
@@ -669,10 +693,98 @@ getResiduals.MixMod <- function (object,...){
   residuals(object, type = "subject_specific")
 }
 
+####### phylolm / phyloglm #########
 
 
-####### sjSDM #########
+#' @rdname getObservedResponse
+#' @export
+getObservedResponse.phylolm <- function (object, ...){
+  out = object$y
+  return(out)
+}
 
 
+#' @rdname getSimulations
+#' @export
+getSimulations.phylolm <- function(object, nsim = 1, type = c("normal", "refit"), ...){
+  type <- match.arg(type)
+  
+  fitBoot = update(object, boot = nsim, save = T)
+  out = fitBoot$bootdata
+  
+  if(type == "normal"){
+    if(!is.matrix(out)) out = data.matrix(out)
+  }else{
+    out = as.data.frame(out)
+  }
+  
+  return(out)
+}
 
+#' @rdname getFitted
+#' @export
+getFitted.phylolm  <- function (object,...){
+  return(object$fitted.values)
+}
 
+
+#' @rdname getFamily
+#' @export
+getFamily.phylolm <- function (object,...){
+  out = list()
+  out$family = "gaussian" 
+  return(out)
+}
+
+
+#' @rdname getObservedResponse
+#' @export
+getObservedResponse.phyloglm <- function (object, ...){
+  out = object$y
+  return(out)
+}
+
+
+#' @rdname getSimulations
+#' @export
+getSimulations.phyloglm <- function(object, nsim = 1, type = c("normal", "refit"), ...){
+  type <- match.arg(type)
+  
+  fitBoot = update(object, boot = nsim, save = T)
+  out = fitBoot$bootdata
+  
+  if(type == "normal"){
+    if(!is.matrix(out)) out = data.matrix(out)
+  }else{
+    out = as.data.frame(out)
+  }
+  
+  return(out)
+}
+
+#' @rdname getFitted
+#' @export
+getFitted.phyloglm  <- function (object,...){
+  return(object$fitted.values)
+}
+
+
+
+#' @rdname getFamily
+#' @export
+getFamily.phyloglm <- function (object,...){
+  out = list()
+  # out$family = object$method
+  out$family = "integer-valued" # all families of phyloglm are integer-valued
+  return(out)
+}
+
+
+#' @rdname getFamily
+#' @export
+getFamily.phyloglm <- function (object,...){
+  out = list()
+  # out$family = object$method
+  out$family = "integer-valued" # all families of phyloglm are integer-valued
+  return(out)
+}
diff --git a/DHARMa/R/simulateResiduals.R b/DHARMa/R/simulateResiduals.R
index ab3e21d8..6a05e0e5 100644
--- a/DHARMa/R/simulateResiduals.R
+++ b/DHARMa/R/simulateResiduals.R
@@ -56,7 +56,7 @@ simulateResiduals <- function(fittedModel, n = 250, refit = F, integerResponse =
   
   out = list()
   
-  family = family(fittedModel)
+  family = getFamily(fittedModel)
   out$fittedModel = fittedModel
   out$modelClass = class(fittedModel)[1]
   
@@ -71,7 +71,7 @@ simulateResiduals <- function(fittedModel, n = 250, refit = F, integerResponse =
   if (out$nObs < length(out$observedResponse)) stop("DHARMA::simulateResiduals: nobs(model) < nrow(model.frame). A possible reason is that you have observation with zero prior weights (ir binomial with n=0) in your data. Calculating residuals in this case wouldn't be sensible. Please remove zero-weight observations from your data and re-fit your model! If you believe that this is not the reason, please file an issue under https://github.com/florianhartig/DHARMa/issues")
   
   if(is.null(integerResponse)){
-    if (family$family %in% c("binomial", "poisson", "quasibinomial", "quasipoisson", "Negative Binom", "nbinom2", "nbinom1", "genpois", "compois", "truncated_poisson", "truncated_nbinom2", "truncated_nbinom1", "betabinomial", "Poisson", "Tpoisson", "COMPoisson", "negbin", "Tnegbin") | grepl("Negative Binomial",family$family) ) integerResponse = TRUE
+    if (family$family %in% c("integer-valued", "binomial", "poisson", "quasibinomial", "quasipoisson", "Negative Binom", "nbinom2", "nbinom1", "genpois", "compois", "truncated_poisson", "truncated_nbinom2", "truncated_nbinom1", "betabinomial", "Poisson", "Tpoisson", "COMPoisson", "negbin", "Tnegbin") | grepl("Negative Binomial",family$family) ) integerResponse = TRUE
     else integerResponse = FALSE
   }
   out$integerResponse = integerResponse
diff --git a/DHARMa/R/tests.R b/DHARMa/R/tests.R
index 35f95689..8ac8d0c7 100644
--- a/DHARMa/R/tests.R
+++ b/DHARMa/R/tests.R
@@ -750,6 +750,71 @@ testSpatialAutocorrelation <- function(simulationOutput, x = NULL, y  = NULL, di
 }
 
 
+#' Test for phylogenetic autocorrelation
+#'
+#' This function performs a Moran's I test for phylogenetic autocorrelation on the calculated quantile residuals
+#'
+#' @param simulationOutput an object of class DHARMa, either created via \code{\link{simulateResiduals}} for supported models or by \code{\link{createDHARMa}} for simulations created outside DHARMa, or a supported model. Providing a supported model directly is discouraged, because simulation settings cannot be changed in this case.
+#' @param tree phylogenetic tree
+#' @param alternative a character string specifying whether the test should test if observations are "greater", "less" or "two.sided" compared to the simulated null hypothesis
+#' @param plot whether to plot output
+#' @details The function performs Moran.I test from the package ape on the DHARMa residuals. If a distance matrix (distMat) is provided, calculations will be based on this distance matrix, and x,y coordinates will only used for the plotting (if provided). If distMat is not provided, the function will calculate the euclidean distances between x,y coordinates, and test Moran.I based on these distances.
+#' 
+#' If plot = T, a plot will be produced showing each residual with at its x,y position, colored according to the residual value. Residuals with 0.5 are colored white, everything below 0.5 is colored increasinly red, everything above 0.5 is colored increasingly blue. 
+#'
+#' Testing for spatial autocorrelation requires unique x,y values - if you have several observations per location, either use the recalculateResiduals function to aggregate residuals per location, or extract the residuals from the fitted object, and plot / test each of them independently for spatially repeated subgroups (a typical scenario would repeated spatial observation, in which case one could plot / test each time step separately for temporal autocorrelation). Note that the latter must be done by hand, outside testSpatialAutocorrelation.
+#'
+#' @note Standard DHARMa simulations from models with (temporal / spatial / phylogenetic) conditional autoregressive terms will still have the respective temporal / spatial / phylogenetic correlation in the DHARMa residuals, unless the package you are using is modelling the autoregressive terms as explicit REs and is able to simulate conditional on the fitted REs. This has two consequences
+#'
+#' 1. If you check the residuals for such a model, they will still show significant autocorrelation, even if the model fully accounts for this structure.
+#'
+#' 2. Because the DHARMa residuals for such a model are not statistically independent any more, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.
+#'
+#' There are three (non-exclusive) routes to address these issues when working with spatial / temporal / other autoregressive models:
+#'
+#' 1. Simulate conditional on the fitted CAR structures (see conditional simulations in the help of [simulateResiduals])
+#'
+#' 2. Rotate simulations prior to residual calculations (see parameter rotation in [simulateResiduals])
+#'
+#' 3. Use custom tests / plots that explicitly compare the correlation structure in the simulated data to the correlation structure in the observed data.
+#'
+#' @author Florian Hartig
+#' @seealso \code{\link{testResiduals}}, \code{\link{testUniformity}}, \code{\link{testOutliers}}, \code{\link{testDispersion}}, \code{\link{testZeroInflation}}, \code{\link{testGeneric}}, \code{\link{testTemporalAutocorrelation}}, \code{\link{testSpatialAutocorrelation}}, \code{\link{testQuantiles}}, \code{\link{testCategorical}}
+#' @example inst/examples/testPhylogeneticAutocorrelationHelp.R
+#' @export
+testPhylogeneticAutocorrelation <- function(simulationOutput, 
+                                            tree, 
+                                            alternative = c("two.sided", "greater", "less"), 
+                                            plot = T){
+  
+  alternative <- match.arg(alternative)
+  data.name = deparse(substitute(simulationOutput)) # needs to be before ensureDHARMa
+  simulationOutput = ensureDHARMa(simulationOutput, convert = T)
+  
+  # calculate distance matrix
+  distMat <- cophenetic(tree)
+  invDistMat <- 1/distMat
+  diag(invDistMat) <- 0
+  
+  MI = ape::Moran.I(simulationOutput$scaledResiduals, weight = invDistMat, alternative = alternative)
+  
+  out = list()
+  out$statistic = c(observed = MI$observed, expected = MI$expected, sd = MI$sd)
+  out$method = "DHARMa Moran's I test for phylogenetic autocorrelation"
+  out$alternative = "Phylogenetic autocorrelation"
+  out$p.value = MI$p.value
+  out$data.name = data.name
+  
+  class(out) = "htest"
+  
+  if(plot == T ) {
+    # todo
+  }
+  return(out)
+}
+
+
+
 getP <- function(simulated, observed, alternative, plot = F, ...){
 
   if(alternative == "greater") p = mean(simulated >= observed)
diff --git a/DHARMa/inst/examples/testPhylogeneticAutocorrelationHelp.R b/DHARMa/inst/examples/testPhylogeneticAutocorrelationHelp.R
new file mode 100644
index 00000000..6db581cd
--- /dev/null
+++ b/DHARMa/inst/examples/testPhylogeneticAutocorrelationHelp.R
@@ -0,0 +1,26 @@
+# https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00044.x
+
+library(DHARMa)
+library(phylolm)
+
+set.seed(123)
+tre = rcoal(60)
+b0=0; b1=1;
+x <- runif(length(tre$tip.label), 0,1)
+y <- b0 + b1*x + 
+  rTrait(n=1, phy=tre,model="BM",
+         parameters=list(ancestral.state=0,sigma2=10))
+dat = data.frame(trait=y,pred=x)
+
+fit = lm(trait~pred,data=dat)
+res = simulateResiduals(fit, plot = T)
+
+testPhylogeneticAutocorrelation(res, tree = tre)
+
+
+fit = phylolm(trait~pred,data=dat,phy=tre,model="BM")
+summary(fit)
+
+res = DHARMa::simulateResiduals(fit, plot = T)
+res = DHARMa::simulateResiduals(fit, plot = T, rotation = "estimated")
+
diff --git a/DHARMa/man/getFamily.Rd b/DHARMa/man/getFamily.Rd
new file mode 100644
index 00000000..4d9203ea
--- /dev/null
+++ b/DHARMa/man/getFamily.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/compatibility.R
+\name{getFamily}
+\alias{getFamily}
+\alias{getFamily.default}
+\alias{getFamily.phylolm}
+\alias{getFamily.phyloglm}
+\title{Get model family}
+\usage{
+getFamily(object, ...)
+
+\method{getFamily}{default}(object, ...)
+
+\method{getFamily}{phylolm}(object, ...)
+
+\method{getFamily}{phyloglm}(object, ...)
+
+\method{getFamily}{phyloglm}(object, ...)
+}
+\arguments{
+\item{object}{a fitted model}
+
+\item{...}{additional parameters to be passed on}
+}
+\description{
+Wrapper to get the family of a fitted model
+
+checks if the fitted model excluded NA values
+}
+\details{
+Checks if the fitted model excluded NA values
+}
+\seealso{
+\code{\link{getObservedResponse}}, \code{\link{getSimulations}}, \code{\link{getRefit}}, \code{\link{getFixedEffects}}, \code{\link{getFitted}}
+}
+\author{
+Florian Hartig
+}
diff --git a/DHARMa/man/getFitted.Rd b/DHARMa/man/getFitted.Rd
index 5c95d18b..0a5f1209 100644
--- a/DHARMa/man/getFitted.Rd
+++ b/DHARMa/man/getFitted.Rd
@@ -6,6 +6,8 @@
 \alias{getFitted.gam}
 \alias{getFitted.HLfit}
 \alias{getFitted.MixMod}
+\alias{getFitted.phylolm}
+\alias{getFitted.phyloglm}
 \title{Get fitted / predicted values}
 \usage{
 getFitted(object, ...)
@@ -17,6 +19,10 @@ getFitted(object, ...)
 \method{getFitted}{HLfit}(object, ...)
 
 \method{getFitted}{MixMod}(object, ...)
+
+\method{getFitted}{phylolm}(object, ...)
+
+\method{getFitted}{phyloglm}(object, ...)
 }
 \arguments{
 \item{object}{a fitted model}
diff --git a/DHARMa/man/getObservedResponse.Rd b/DHARMa/man/getObservedResponse.Rd
index 19f80f29..00bd172f 100644
--- a/DHARMa/man/getObservedResponse.Rd
+++ b/DHARMa/man/getObservedResponse.Rd
@@ -4,6 +4,8 @@
 \alias{getObservedResponse}
 \alias{getObservedResponse.default}
 \alias{getObservedResponse.HLfit}
+\alias{getObservedResponse.phylolm}
+\alias{getObservedResponse.phyloglm}
 \title{Get model response}
 \usage{
 getObservedResponse(object, ...)
@@ -11,6 +13,10 @@ getObservedResponse(object, ...)
 \method{getObservedResponse}{default}(object, ...)
 
 \method{getObservedResponse}{HLfit}(object, ...)
+
+\method{getObservedResponse}{phylolm}(object, ...)
+
+\method{getObservedResponse}{phyloglm}(object, ...)
 }
 \arguments{
 \item{object}{a fitted model}
diff --git a/DHARMa/man/getRefit.Rd b/DHARMa/man/getRefit.Rd
index e2a7efe1..853a5686 100644
--- a/DHARMa/man/getRefit.Rd
+++ b/DHARMa/man/getRefit.Rd
@@ -30,13 +30,9 @@ getRefit(object, newresp, ...)
 }
 \description{
 Wrapper to refit a fitted model
-
-checks if the fitted model excluded NA values
 }
 \details{
 The purpose of this wrapper is to standardize the refit of a model. The behavior of this function depends on the supplied model. When available, it uses the refit method, otherwise it will use update. For glmmTMB: since version 1.0, glmmTMB has a refit function, but this didn't work, so I switched back to this implementation, which is a hack based on the update function.
-
-Checks if the fitted model excluded NA values
 }
 \examples{
 testData = createData(sampleSize = 400, family = gaussian())
diff --git a/DHARMa/man/getSimulations.Rd b/DHARMa/man/getSimulations.Rd
index a6601a40..e8f03d8e 100644
--- a/DHARMa/man/getSimulations.Rd
+++ b/DHARMa/man/getSimulations.Rd
@@ -9,6 +9,8 @@
 \alias{getSimulations.glmmTMB}
 \alias{getSimulations.HLfit}
 \alias{getSimulations.MixMod}
+\alias{getSimulations.phylolm}
+\alias{getSimulations.phyloglm}
 \title{Get model simulations}
 \usage{
 getSimulations(object, nsim = 1, type = c("normal", "refit"), ...)
@@ -33,6 +35,12 @@ getSimulations(object, nsim = 1, type = c("normal", "refit"), ...)
 
 \method{getSimulations}{MixMod}(object, nsim = 1, type = c("normal",
   "refit"), ...)
+
+\method{getSimulations}{phylolm}(object, nsim = 1, type = c("normal",
+  "refit"), ...)
+
+\method{getSimulations}{phyloglm}(object, nsim = 1, type = c("normal",
+  "refit"), ...)
 }
 \arguments{
 \item{object}{a fitted model}
diff --git a/DHARMa/man/testPhylogeneticAutocorrelation.Rd b/DHARMa/man/testPhylogeneticAutocorrelation.Rd
new file mode 100644
index 00000000..9b0b964b
--- /dev/null
+++ b/DHARMa/man/testPhylogeneticAutocorrelation.Rd
@@ -0,0 +1,76 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tests.R
+\name{testPhylogeneticAutocorrelation}
+\alias{testPhylogeneticAutocorrelation}
+\title{Test for phylogenetic autocorrelation}
+\usage{
+testPhylogeneticAutocorrelation(simulationOutput, tree,
+  alternative = c("two.sided", "greater", "less"), plot = T)
+}
+\arguments{
+\item{simulationOutput}{an object of class DHARMa, either created via \code{\link{simulateResiduals}} for supported models or by \code{\link{createDHARMa}} for simulations created outside DHARMa, or a supported model. Providing a supported model directly is discouraged, because simulation settings cannot be changed in this case.}
+
+\item{tree}{phylogenetic tree}
+
+\item{alternative}{a character string specifying whether the test should test if observations are "greater", "less" or "two.sided" compared to the simulated null hypothesis}
+
+\item{plot}{whether to plot output}
+}
+\description{
+This function performs a Moran's I test for phylogenetic autocorrelation on the calculated quantile residuals
+}
+\details{
+The function performs Moran.I test from the package ape on the DHARMa residuals. If a distance matrix (distMat) is provided, calculations will be based on this distance matrix, and x,y coordinates will only used for the plotting (if provided). If distMat is not provided, the function will calculate the euclidean distances between x,y coordinates, and test Moran.I based on these distances.
+
+If plot = T, a plot will be produced showing each residual with at its x,y position, colored according to the residual value. Residuals with 0.5 are colored white, everything below 0.5 is colored increasinly red, everything above 0.5 is colored increasingly blue.
+
+Testing for spatial autocorrelation requires unique x,y values - if you have several observations per location, either use the recalculateResiduals function to aggregate residuals per location, or extract the residuals from the fitted object, and plot / test each of them independently for spatially repeated subgroups (a typical scenario would repeated spatial observation, in which case one could plot / test each time step separately for temporal autocorrelation). Note that the latter must be done by hand, outside testSpatialAutocorrelation.
+}
+\note{
+Standard DHARMa simulations from models with (temporal / spatial / phylogenetic) conditional autoregressive terms will still have the respective temporal / spatial / phylogenetic correlation in the DHARMa residuals, unless the package you are using is modelling the autoregressive terms as explicit REs and is able to simulate conditional on the fitted REs. This has two consequences
+\enumerate{
+\item If you check the residuals for such a model, they will still show significant autocorrelation, even if the model fully accounts for this structure.
+\item Because the DHARMa residuals for such a model are not statistically independent any more, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.
+}
+
+There are three (non-exclusive) routes to address these issues when working with spatial / temporal / other autoregressive models:
+\enumerate{
+\item Simulate conditional on the fitted CAR structures (see conditional simulations in the help of \link{simulateResiduals})
+\item Rotate simulations prior to residual calculations (see parameter rotation in \link{simulateResiduals})
+\item Use custom tests / plots that explicitly compare the correlation structure in the simulated data to the correlation structure in the observed data.
+}
+}
+\examples{
+# https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00044.x
+
+library(DHARMa)
+library(phylolm)
+
+set.seed(123)
+tre = rcoal(60)
+b0=0; b1=1;
+x <- runif(length(tre$tip.label), 0,1)
+y <- b0 + b1*x + 
+  rTrait(n=1, phy=tre,model="BM",
+         parameters=list(ancestral.state=0,sigma2=10))
+dat = data.frame(trait=y,pred=x)
+
+fit = lm(trait~pred,data=dat)
+res = simulateResiduals(fit, plot = T)
+
+testPhylogeneticAutocorrelation(res, tree = tre)
+
+
+fit = phylolm(trait~pred,data=dat,phy=tre,model="BM")
+summary(fit)
+
+res = DHARMa::simulateResiduals(fit, plot = T)
+res = DHARMa::simulateResiduals(fit, plot = T, rotation = "estimated")
+
+}
+\seealso{
+\code{\link{testResiduals}}, \code{\link{testUniformity}}, \code{\link{testOutliers}}, \code{\link{testDispersion}}, \code{\link{testZeroInflation}}, \code{\link{testGeneric}}, \code{\link{testTemporalAutocorrelation}}, \code{\link{testSpatialAutocorrelation}}, \code{\link{testQuantiles}}, \code{\link{testCategorical}}
+}
+\author{
+Florian Hartig
+}