129 phylolm

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)

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 = "[email protected]", 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

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)

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)
+}

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

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)