OSCR

Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § STAR★METHODS › QUANTIFICATION AND STATISTICAL ANALYSIS ↔ DHARMa/R/tests.R, lines 402–520 · score 0.53 · zero inflation, glmmTMB, outliers, Pearson, nonparametric, ratios
  2. [2] § STAR★METHODS › QUANTIFICATION AND STATISTICAL ANALYSIS ↔ glmmTMB/R/diagnose.R, lines 1–35 · score 0.51 · zero inflation, Diagnostic, glmmTMB, logit, coefficients, ratios

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 903 lines · 56 KB · gnu · 1 match

  1. #' DHARMa general residual test
  2. #'
  3. #' Calls uniformity, dispersion and outliers tests.
  4. #'
  5. #' This function is a wrapper for the various test functions implemented in DHARMa. Currently, this function calls the functions [testUniformity], [testDispersion], and [testOutliers]. All other tests (see list below) have to be called by hand.
  6. #'
  7. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  8. #' @param plot if TRUE, plot functions of the tests are called.
  9. #' @author Florian Hartig
  10. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  11. #' @example inst/examples/testsHelp.R
  12. #' @export
  13. testResiduals <- function(simulationOutput, plot = TRUE){
  15. opar = par(mfrow = c(1,3))
  16. on.exit(par(opar))
  17. out = list()
  18. out$uniformity = testUniformity(simulationOutput, plot = plot)
  19. out$dispersion = testDispersion(simulationOutput, plot = plot)
  20. out$outliers = testOutliers(simulationOutput, plot = plot)
  22. #print(out) # do we need it?
  23. return(out)
  24. }
  26. #' Residual tests
  27. #'
  28. #' @details Deprecated, switch your code to using the [testResiduals] function
  29. #'
  30. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  31. #' @author Florian Hartig
  32. #' @export
  33. testSimulatedResiduals <- function(simulationOutput){
  34. message("testSimulatedResiduals is deprecated, switch your code to using the testResiduals function")
  35. testResiduals(simulationOutput)
  36. }
  39. #' Test for overall uniformity
  40. #'
  41. #' This function tests the overall uniformity of the simulated residuals in a DHARMa object.
  42. #'
  43. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  44. #' @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. See [stats::ks.test] for details.
  45. #' @param plot if TRUE, plots calls [plotQQunif] as well.
  46. #' @details The function applies a [stats::ks.test] for uniformity on the simulated residuals.
  47. #' @author Florian Hartig
  48. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  49. #' @example inst/examples/testsHelp.R
  50. #' @export
  51. testUniformity <- function(simulationOutput, alternative = c("two.sided", "less", "greater"), plot = TRUE){
  53. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  55. out <- suppressWarnings(ks.test(simulationOutput$scaledResiduals, 'punif', alternative = alternative))
  56. if(plot == T) plotQQunif(simulationOutput = simulationOutput)
  57. return(out)
  58. }
  61. # Experimental
  62. testBivariateUniformity <- function(simulationOutput, alternative = c("two.sided", "less", "greater"), plot = TRUE){
  64. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  66. #out <- suppressWarnings(ks.test(simulationOutput$scaledResiduals, 'punif', alternative = alternative))
  67. #if(plot == T) plotQQunif(simulationOutput = simulationOutput)
  68. out = NULL
  69. return(out)
  70. }
  74. #' Test for quantiles
  75. #'
  76. #' This function fits quantile regressions on the residuals, and compares their location to the expected location.
  77. #'
  78. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  79. #' @param predictor an optional predictor variable to be used, instead of the predicted response (default). See details.
  80. #' @param rank if TRUE, the values provided in predictor will be rank transformed. This will usually make patterns easier to spot visually, especially if the distribution of the predictor is skewed. If form is a factor, this has no effect.
  81. #' @param quantiles the quantiles to be tested.
  82. #' @param plot if TRUE, the function will create an additional plot.
  83. #' @details The function fits quantile regressions (via package qgam) on the residuals, and compares their location to the expected location (because of the uniform distribution, the expected location is 0.5 for the 0.5 quantile).
  84. #'
  85. #' A significant p-value for the splines means the fitted spline deviates from a flat line at the expected location.
  86. #'
  87. #' The p-values of the intercept and splines are combined into a total p-value via Benjamini & Hochberg adjustment to control the FDR.
  88. #'
  89. #' Predictor can be a formula (e.g. predictor = ~predictor), in which case NAs are handled automatically (recommended). Predictor can also be a variable in your environment (e.g. predictor = data$predictor), but then you need to remove rows that were excluded by the model due to NAs by hand. For more details and a more flexible syntax for predictors, see the help of the argument `form` in [plotResiduals].
  90. #'
  91. #' When plotting (plot = TRUE), the shaded gray areas indicate 95% confidence intervals of the quantile estimates (1.96 * standard error).
  92. #'
  93. #' @author Florian Hartig
  94. #' @example inst/examples/testQuantilesHelp.R
  95. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  96. #' @export
  97. testQuantiles <- function(simulationOutput, predictor = NULL, rank = TRUE,
  98. quantiles = c(0.25,0.5,0.75), plot = TRUE){
  100. if(plot == F){
  102. out = list()
  103. out$data.name = deparse(substitute(simulationOutput))
  105. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  106. res = simulationOutput$scaledResiduals
  108. if(inherits(predictor, "formula")) predictor = getFormulaPredictors(simulationOutput, predictor)[[1]]
  109. pred = ensurePredictor(simulationOutput, predictor)
  110. if (rank == TRUE) pred = rankTransform(pred)
  112. dat = data.frame(res = simulationOutput$scaledResiduals, pred = pred)
  114. quantileFits <- list()
  115. pval = rep(NA, length(quantiles))
  116. predictions = data.frame(pred = sort(dat$pred))
  117. predictions = cbind(predictions, matrix(ncol = 2 * length(quantiles),
  118. nrow = nrow(dat)))
  119. for(i in 1:length(quantiles)){
  120. datTemp = dat
  121. datTemp$res = datTemp$res - quantiles[i]
  123. # settings for k = the dimension of the basis used to represent the smooth term.
  124. # see https://github.com/mfasiolo/qgam/issues/37
  125. dimSmooth = min(length(unique(datTemp$pred)), 10)
  126. quantResult = try(capture.output(quantileFits[[i]] <-
  127. qgam::qgam(res ~ s(pred, k = dimSmooth),
  128. data = datTemp,
  129. qu = quantiles[i])), silent = T)
  130. if(inherits(quantResult, "try-error")){
  131. message("\n DHARMa: qgam was unable to calculate quantile regression for quantile ",
  132. quantiles[i], ". Possibly too few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations. \n")
  133. } else {
  134. x = summary(quantileFits[[i]])
  135. pval[i] = min(p.adjust(c(x$p.table[1,4], x$s.table[1,4]), method = "BH")) # correction for test on slope and intercept
  136. quantPre = predict(quantileFits[[i]], newdata = predictions, se = T)
  137. predictions[, 2*i] = quantPre$fit + quantiles[i]
  138. predictions[, 2*i + 1] = 1.96*quantPre$se.fit
  139. }
  140. }
  142. out$method = "Test for location of quantiles via qgam"
  143. out$alternative = "both"
  144. out$pvals = pval
  145. out$p.value = min(p.adjust(pval, method = "BH")) # correction for multiple quantile tests
  146. out$predictions = predictions
  147. out$qgamFits = quantileFits
  149. class(out) = "htest"
  151. } else if(plot == T) {
  152. if(is.null(predictor)) {
  153. out <- plotResiduals(simulationOutput = simulationOutput, form = NULL,
  154. rank = rank, quantiles = quantiles, quantreg = TRUE)
  155. } else {
  156. out <- plotResiduals(simulationOutput = simulationOutput, form = predictor,
  157. rank = rank, quantiles = quantiles, quantreg = TRUE)
  158. }
  160. }
  161. return(out)
  162. }
  165. #unif.2017YMi(X, type = c("Q1", "Q2", "Q3"), lower = rep(0, ncol(X)),upper = rep(1, ncol(X)))
  167. #' Test for outliers
  168. #'
  169. #' This function tests if the number of observations outside the simulation envelope are larger or smaller than expected
  170. #'
  171. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  172. #' @param alternative a character string specifying whether the test should test if observations are "greater", "less" or "two.sided" (default) compared to the simulated null hypothesis.
  173. #' @param margin whether to test for outliers only at the lower, only at the upper, or both sides (default) of the simulated data distribution.
  174. #' @param type either default, bootstrap or binomial. See details.
  175. #' @param nBoot number of bootstrap replicates. Only used at type = "bootstrap".
  176. #' @param plot if TRUE, the function will create an additional plot.
  177. #' @param plotBoostrap if plot should be produced of outlier frequencies calculated under the bootstrap.
  178. #' @details DHARMa residuals are created by simulating from the fitted model, and comparing the simulated values to the observed data. It can occur that all simulated values are higher or smaller than the observed data, in which case they get the residual value of 0 and 1, respectively. I refer to these values as simulation outliers, or simply outliers.
  179. #'
  180. #' Because no data was simulated in the range of the observed values, we don't know "how strongly" these values deviate from the model expectation, so the term "outlier" should be used with a grain of salt. It is not a judgment about the magnitude of the residual deviation, but simply a dichotomous sign that we are outside the simulated range. Moreover, the number of outliers will decrease as we increase the number of simulations.
  181. #'
  182. #' To test if the outliers are a concern, testOutliers implements 2 options (bootstrap, binomial), which can be chosen via the parameter "type". The third option (default) chooses bootstrap for integer-valued distributions with nObs < 500, and else binomial.
  183. #'
  184. #' The binomial test considers that under the null hypothesis that the model is correct, and for continuous distributions (i.e. data and the model distribution are identical and continuous), the probability that a given observation is higher than all simulations is 1/(nSim +1), and binomial distributed. The testOutlier function can test this null hypothesis via type = "binomial". In principle, it would be nice if we could extend this idea to integer-valued distributions, which are randomized via the PIT procedure (see [simulateResiduals]); the rate of "true" outliers is more difficult to calculate, and in general not 1/(nSim +1). The testOutlier function implements a small tweak that calculates the rate of residuals that are closer than 1/(nSim+1) to the 0/1 border, which roughly occur at a rate of nData /(nSim +1). This approximate value, however, is generally not exact, and may be particularly off non-bounded integer-valued distributions (such as Poisson or Negative Binomial).
  185. #'
  186. #' For this reason, the testOutlier function implements an alternative procedure that uses the bootstrap to generate a simulation-based expectation for the outliers. It is recommended to use the bootstrap for integer-valued distributions (and integer-valued only, because it has no advantage for continuous distributions, ideally with reasonably high values of nSim and nBoot (I recommend at least 1000 for both)). Because of the high runtime, however, this option is switched off for type = default when nObs > 500.
  187. #'
  188. #' Both binomial and bootstrap generate a null expectation, and then test for an excess or lack of outliers. Per default, testOutliers() looks for both, so if you get a significant p-value, you have to check if you have too many or too few outliers. An excess of outliers is to be interpreted as too many values outside the simulation envelope. This could be caused by overdispersion, or by what we classically call outliers. A lack of outliers would be caused, for example, by underdispersion.
  189. #'
  190. #'
  191. #' @author Florian Hartig
  192. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  193. #' @export
  194. testOutliers <- function(simulationOutput, alternative = c("two.sided", "greater", "less"), margin = c("both", "upper", "lower"), type = c("default","bootstrap", "binomial"), nBoot = 100, plot = TRUE, plotBoostrap = FALSE){
  196. # check inputs
  197. alternative = match.arg(alternative)
  198. margin = match.arg(margin)
  199. type = match.arg(type)
  200. data.name = deparse(substitute(simulationOutput)) # remember: needs to be called before ensureDHARMa
  201. simulationOutput = ensureDHARMa(simulationOutput, convert = "Model")
  203. if(type == "default"){
  204. if(simulationOutput$integerResponse == FALSE) type = "binomial"
  205. else{
  206. if(simulationOutput$nObs > 500) type = "binomial"
  207. else type = "bootstrap"
  208. }
  209. }
  211. # using the binomial test, not exact
  212. if(type == "binomial"){
  214. # calculation of outliers
  215. if(margin == "both") outliers = sum(simulationOutput$scaledResiduals < (1/(simulationOutput$nSim+1))) + sum(simulationOutput$scaledResiduals > (1-1/(simulationOutput$nSim+1)))
  216. if(margin == "upper") outliers = sum(simulationOutput$scaledResiduals > (1-1/(simulationOutput$nSim+1)))
  217. if(margin == "lower") outliers = sum(simulationOutput$scaledResiduals < (1/(simulationOutput$nSim+1)))
  219. # calculations of trials and H0
  220. outFreqH0 = 1/(simulationOutput$nSim +1) * ifelse(margin == "both", 2, 1)
  221. trials = simulationOutput$nObs
  223. out = binom.test(outliers, trials, p = outFreqH0, alternative = alternative)
  225. # overwrite information in binom.test
  226. out$method = "DHARMa outlier test based on exact binomial test with approximate expectations"
  227. out$data.name = data.name
  228. out$margin = margin
  229. names(out$statistic) = paste("outliers at", margin, "margin(s)")
  230. names(out$parameter) = "observations"
  231. names(out$estimate) = paste("frequency of outliers (expected:", out$null.value,")")
  233. if (simulationOutput$integerResponse == T & out$p.value < 0.05) message("DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details")
  235. if(plot == T) {
  237. hist(simulationOutput, main = "")
  239. main = ifelse(out$p.value <= 0.05,
  240. "Outlier test significant",
  241. "Outlier test n.s.")
  243. title(main = main, cex.main = 1,
  244. col.main = ifelse(out$p.value <= 0.05, "red", "black"))
  246. }
  248. } else {
  250. if(margin == "both") outliers = mean(simulationOutput$scaledResiduals == 0) +
  251. mean(simulationOutput$scaledResiduals == 1)
  252. if(margin == "upper") outliers = mean(simulationOutput$scaledResiduals == 1)
  253. if(margin == "lower") outliers = mean(simulationOutput$scaledResiduals == 0)
  256. # Bootstrapping to compare to expected
  258. simIndices = 1:simulationOutput$nSim
  259. nSim = simulationOutput$nSim
  260. if(simulationOutput$refit == T){
  261. simResp = simulationOutput$refittedResiduals
  262. } else {
  263. simResp = simulationOutput$simulatedResponse
  264. }
  265. resMethod = simulationOutput$method
  266. resInteger = simulationOutput$integerResponse
  269. if (nBoot > nSim){
  270. message("DHARMa::testOutliers: nBoot > nSim does not make much sense, thus changed to nBoot = nSim. If you want to increase nBoot, increase nSim in DHARMa::simulateResiduals as well.")
  271. nBoot = nSim
  272. }
  274. frequBoot <- rep(NA, nBoot)
  276. for (i in 1:nBoot){
  278. #sel = -i
  279. sel = sample(simIndices[-i], size = nSim, replace = T)
  281. residuals <- getQuantile(simulations = simResp[,sel],
  282. observed = simResp[,i],
  283. integerResponse = resInteger,
  284. method = resMethod)
  286. if(margin == "both") frequBoot[i] = mean(residuals == 1) + mean(residuals == 0)
  287. else if(margin == "upper") frequBoot[i] = mean(residuals == 1)
  288. else if(margin == "lower") frequBoot[i] = mean(residuals == 0)
  289. }
  291. out = list()
  292. class(out) = "htest"
  293. out$alternative = alternative
  294. out$p.value = getP(frequBoot, outliers, alternative = alternative)
  296. out$conf.int = quantile(frequBoot, c(0.025, 0.975))
  298. out$data.name = data.name
  299. out$margin = margin
  301. out$method = "DHARMa bootstrapped outlier test"
  303. out$statistic = outliers * simulationOutput$nObs
  304. names(out$statistic) = paste("outliers at", margin, "margin(s)")
  305. out$parameter = simulationOutput$nObs
  306. names(out$parameter) = "observations"
  307. out$estimate = outliers
  308. names(out$estimate) = paste("outlier frequency (expected:", mean(frequBoot),")")
  310. if(plotBoostrap == T){
  311. hist(frequBoot, xlim = range(frequBoot, outliers), col = "lightgrey", main = "Bootstrapped outlier frequency")
  312. abline(v = mean(frequBoot), col = 1, lwd = 2)
  313. abline(v = outliers, col = "red", lwd = 2)
  315. # legend("center", c(paste("p=", round(out$p.value, digits = 5)), paste("Deviation ", ifelse(out$p.value < 0.05, "significant", "n.s."))), text.col = ifelse(out$p.value < 0.05, "red", "black" ))
  316. }
  318. if(plot == T) {
  320. hist(simulationOutput, main = "")
  322. main = ifelse(out$p.value <= 0.05,
  323. "Outlier test significant",
  324. "Outlier test n.s.")
  326. title(main = main, cex.main = 1,
  327. col.main = ifelse(out$p.value <= 0.05, "red", "black"))
  329. }
  330. }
  331. return(out)
  332. }
  336. #' Test for categorical dependencies
  337. #'
  338. #' This function tests if there are problems in a res ~ group structure. It performs two tests: test for within-group uniformity, and test for between-group homogeneity of variances.
  339. #'
  340. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  341. #' @param catPred a categorical predictor with the same dimensions as the residuals in simulationOutput. If specified as formula, e.g. catPred = ~group, NAs are handled automatically (recommended).
  342. #' @param quantiles whether to draw the quantile lines.
  343. #' @param plot if TRUE, the function will create an additional plot.
  344. #' @param ... additional arguments to boxplot.
  345. #' @details The function tests for two common problems: are residuals within each group distributed according to model assumptions, and is the variance between groups heterogeneous.
  346. #'
  347. #' The test for within-group uniformity is performed via multiple KS-tests, with adjustment of p-values for multiple testing. If the plot is drawn, problematic groups are highlighted in red, and a corresponding message is displayed in the plot.
  348. #'
  349. #' The test for homogeneity of variances is done with a Levene test. A significant p-value means that group variances are not constant. In this case, you should consider modelling variances, e.g. via ~dispformula in glmmTMB.
  350. #'
  351. #' @author Florian Hartig
  352. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  353. #' @example inst/examples/testsHelp.R
  354. #' @export
  355. testCategorical <- function(simulationOutput, catPred,
  356. quantiles = c(0.25, 0.5, 0.75), plot = TRUE, ...){
  358. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  360. a = list(...)
  361. a$xlab = paste(checkDots("xlab", deparse(substitute(catPred)), ...), "(catPred)")
  363. if(!inherits(catPred, "formula")) {
  364. catPred = catPred # old default
  365. } else {
  366. catPred = getFormulaPredictors(simulationOutput, formula = catPred)[[1]] # formula syntax
  367. }
  369. catPred = as.factor(catPred)
  370. out = list()
  372. out$uniformity$details = suppressWarnings(by(simulationOutput$scaledResiduals,
  373. catPred, ks.test, 'punif',
  374. simplify = TRUE))
  375. out$uniformity$p.value = rep(NA, nlevels(catPred))
  376. for(i in 1:nlevels(catPred)) out$uniformity$p.value[i] = out$uniformity$details[[i]]$p.value
  377. out$uniformity$p.value.cor = p.adjust(out$uniformity$p.value)
  379. if(nlevels(catPred) > 1) out$homogeneity = leveneTest_formula(simulationOutput$scaledResiduals ~ catPred)
  381. if(plot == T){
  382. boxplot(simulationOutput$scaledResiduals ~ catPred, ylim = c(0,1), xlab = a$xlab, axes = FALSE, col = ifelse(out$uniformity$p.value.cor < 0.05, "red", "lightgrey"))
  383. axis(1, at = 1:nlevels(catPred), levels(catPred))
  384. axis(2, at=c(0, quantiles, 1))
  385. abline(h = quantiles, lty = 2)
  386. }
  388. title(ifelse(any(out$uniformity$p.value.cor < 0.05), "Within-group deviations from uniformity significant (red)", "Within-group deviation from uniformity n.s."),
  389. col.main = ifelse(any(out$uniformity$p.value.cor < 0.05), "red", "black"),
  390. line = 1, cex.main = 0.8)
  392. if(length(out) > 1) {
  393. title(ifelse(out$homogeneity$`Pr(>F)`[1] < 0.05, "Levene Test for homogeneity of variance significant", "Levene Test for homogeneity of variance n.s."),
  394. col.main = ifelse(out$homogeneity$`Pr(>F)`[1] < 0.05, "red", "black"), cex.main = 0.8)
  395. }
  397. return(out)
  399. }
  402. #' DHARMa dispersion tests
  403. #'
  404. #' This function performs simulation-based tests for over / underdispersion. If type = "DHARMa" (default and recommended), simulation-based dispersion tests are performed. Their behavior differs depending on whether simulations are done with refit = F, or refit = T, and whether data is simulated conditional (see below). If type = "PearsonChisq", a chi2 test on Pearson residuals is performed.
  405. #'
  406. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  407. #' @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. Greater corresponds to testing only for overdispersion. It is recommended to keep the default setting (testing for both over and underdispersion).
  408. #' @param plot whether to provide a plot for the results.
  409. #' @param type which test to run. Default is DHARMa, other options are PearsonChisq (see details).
  410. #' @param ... arguments to pass on to [testGeneric].
  411. #'
  412. #' @details Over / underdispersion means that the observed data is more / less dispersed than expected under the fitted model. There is no unique way to test for dispersion problems, and there are a number of different dispersion tests implemented in various R packages. For detailed recommendations see also Leite M.S., Rettelbach, D. and Hartig, F. (2025), preprint available at https://doi.org/10.32942/X23M14.
  413. #'
  414. #' The testDispersion function implements several dispersion tests:
  415. #'
  416. #' **Simulation-based dispersion tests (type == "DHARMa")**
  417. #'
  418. #' If type = "DHARMa" (default and recommended), simulation-based dispersion tests are performed. Their behavior differs depending on whether simulations are done with refit = F, or refit = T.
  419. #'
  420. #' **Important:** for either refit = T or F, the results of type = "DHARMa" dispersion test will differ depending on whether simulations are done conditional (= conditional on fitted random effects) or unconditional (= REs are re-simulated). The default as of DHARMa 0.5.0 is conditional simulations, which substantially increase the power and sensitivity of the dispersion test (for details see [simulateResiduals] and the DHARMa vignette).
  421. #'
  422. #' If refit = F, the function uses [testGeneric] to compare the variance of the observed raw residuals (i.e. var(observed - predicted), displayed as a red line) against the variance of the simulated residuals (i.e. var(simulated - predicted), histogram). The variances are scaled to the mean simulated variance. A significant ratio > 1 indicates overdispersion, a significant ratio < 1 underdispersion.
  423. #'
  424. #' If refit = T, the function compares the approximate deviance (via squared Pearson residuals) with the same quantity from the models refitted with simulated data. Applying this is much slower than the previous alternative. Given the computational cost, I would suggest that most users will be satisfied with the standard dispersion test.
  425. #'
  426. #' **Analytical dispersion tests (type == "PearsonChisq")**
  427. #'
  428. #' This is the test described in https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion, identical to performance::check_overdispersion. Works only if the fitted model provides df.residual and Pearson residuals.
  429. #'
  430. #' The test statistic is biased to lower values under quite general conditions, and will therefore tend to test significant for underdispersion. It is recommended to use this test only for overdispersion, i.e. use alternative == "greater". Also, obviously, it requires that Pearson residuals are available for the chosen model, which will not be the case for all models / packages.
  431. #'
  432. #' @note For particular model classes / situations, there may be more powerful and thus preferable options over the DHARMa test. The advantage of the DHARMa test is that it directly targets the spread of the data (unlike other tests such as dispersion/df, which essentially measure fit and may thus be triggered by problems other than dispersion as well), and it makes practically no assumptions about the fitted model, other than the availability of simulations.
  433. #'
  434. #' @author Florian Hartig
  435. #' @references Leite, MS, Rettelbach, D. & Hartig, F. 2025. Dispersion tests in generalised linear mixed-effects models - a methods comparison and practical guide. EcoEvoRxiv Preprint: https://doi.org/10.32942/X23M14.
  436. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  437. #' @example inst/examples/testDispersionHelp.R
  438. #' @export
  439. testDispersion <- function(simulationOutput, alternative = c("two.sided", "greater", "less"), plot = T, type = c("DHARMa", "PearsonChisq"), ...){
  441. alternative <- match.arg(alternative)
  442. type <- match.arg(type)
  444. out = list()
  445. out$data.name = deparse(substitute(simulationOutput))
  447. if(type == "DHARMa"){
  449. simulationOutput = ensureDHARMa(simulationOutput, convert = "Model")
  451. if(simulationOutput$refit == F){
  453. expectedVar = sd(simulationOutput$simulatedResponse)^2
  454. spread <- function(x) var(x - simulationOutput$fittedPredictedResponse) / expectedVar
  455. out = testGeneric(simulationOutput, summary = spread, alternative = alternative, methodName = "DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated", plot = plot, ...)
  456. names(out$statistic) = "dispersion"
  457. } else {
  459. rp <- getPearsonResiduals(simulationOutput$fittedModel)
  460. observed = sum(rp^2)
  461. expected = apply(simulationOutput$refittedPearsonResiduals^2 , 2, sum)
  462. out$statistic = c(dispersion = observed / mean(expected))
  463. names(out$statistic) = "dispersion"
  464. out$method = "DHARMa nonparametric dispersion test via mean deviance residual fitted vs. simulated-refitted"
  466. p = getP(simulated = expected, observed = observed, alternative = alternative)
  468. out$alternative = alternative
  469. out$p.value = p
  470. class(out) = "htest"
  472. if(plot == T) {
  473. #plotTitle = gsub('(.{1,50})(\\s|$)', '\\1\n', out$method)
  474. xLabel = paste("Simulated values, red line = fitted model. p-value (",out$alternative, ") = ", out$p.value, sep ="")
  476. hist(expected, xlim = range(expected, observed, na.rm=T ), col = "lightgrey", main = "", xlab = xLabel, breaks = 20, cex.main = 1)
  477. abline(v = observed, lwd= 2, col = "red")
  479. main = ifelse(out$p.value <= 0.05,
  480. "Dispersion test significant",
  481. "Dispersion test n.s.")
  483. title(main = main, cex.main = 1,
  484. col.main = ifelse(out$p.value <= 0.05, "red", "black"))
  485. }
  486. }
  488. } else if(type == "PearsonChisq"){
  490. if("DHARMa" %in% class(simulationOutput)){
  491. model = simulationOutput$fittedModel
  492. }
  493. else model = simulationOutput
  495. if(!alternative == "greater" & class(model)[1] %in% c("lmerMod", "lmerModLmerTest", "glmerMod", "bam", "glmmTMB", "HLfit", "MixMod")) {
  496. message("Note that the Chi2 test on Pearson residuals is biased for MIXED models towards underdispersion. Tests with alternative = two.sided or less are therefore not reliable. If you have random effects in your model, we recommend to test only with alternative = 'greater', i.e. test for overdispersion, or else use the DHARMa default tests which are unbiased. See help for details.")}
  498. rp <- getPearsonResiduals(model)
  500. rdf <- df.residual(model)
  501. Pearson.chisq <- sum(rp^2)
  502. prat <- Pearson.chisq/rdf
  503. if(alternative == "greater") pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  504. else if (alternative == "less") pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE)
  505. else if (alternative == "two.sided") pval <- min(min(pchisq(Pearson.chisq, df=rdf, lower.tail=TRUE), pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)) * 2,1)
  507. out$statistic = prat
  508. names(out$statistic) = "dispersion"
  509. out$parameter = rdf
  510. names(out$parameter) = "df"
  511. out$method = "Parametric dispersion test via mean Pearson-chisq statistic"
  512. out$alternative = alternative
  513. out$p.value = pval
  514. class(out) = "htest"
  515. # c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
  516. return(out)
  517. }
  519. return(out)
  520. }
  522. #' Simulated overdisperstion tests (deprecated)
  523. #'
  524. #' @details Deprecated, switch your code to using the [testDispersion] function
  525. #'
  526. #' @param simulationOutput an object of class DHARMa with simulated quantile residuals, either created via [simulateResiduals] or by [createDHARMa] for simulations created outside DHARMa
  527. #' @param ... additional arguments to [testDispersion]
  528. #' @export
  529. testOverdispersion <- function(simulationOutput, ...){
  530. message("testOverdispersion is deprecated, switch your code to using the testDispersion function")
  531. testDispersion(simulationOutput, ...)
  532. }
  534. #' Parametric overdisperstion tests (deprecated)
  535. #'
  536. #' @details Deprecated, switch your code to using the [testDispersion] function.
  537. #'
  538. #' @param ... arguments will be ignored, the parametric tests is no longer recommend
  539. #' @export
  540. testOverdispersionParametric <- function(...){
  541. message("testOverdispersionParametric is deprecated - switch your code to using the testDispersion function")
  542. return(0)
  543. }
  548. #' Tests for zero-inflation
  549. #'
  550. #' This function compares the observed number of zeros with the zeros expected from simulations.
  551. #'
  552. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  553. #' @param ... further arguments to [testGeneric].
  554. #' @details Zero-inflation means that the observed data contains more zeros than would be expected under the fitted model. Zero-inflation must always be accessed with respect to a particular model, so the mere fact that there are many zeros in the observed data is not an indication of zero-inflation, see Warton, D. I. (2005). Many zeros does not mean zero inflation: comparing the goodness-of-fit of parametric models to multivariate abundance data. Environmetrics 16(3), 275-289.
  555. #'
  556. #' The testZeroInflation function simulates new datasets from the fitted model and compares this null distribution (gray histogram in the plot) with the observed values (red line in the plot). Technically, it is a wrapper for [testGeneric], with the summary argument set to function(x) sum(x == 0). The test statistic is the ratio of observed to simulated zeros. A value < 1 means that the observed data have fewer zeros than expected, a value > 1 means that they have more zeros than expected (aka zero inflation). By default, the function tests both sides, so it would also test for fewer zeros than expected.
  557. #'
  558. #' @note Zero-inflation can occur for a number of reasons other than an underlying data generating process corresponding to a ZIP model. Vice versa, it is very well possible that no zero-inflation will be observed when fitting models to data derived from a ZIP process. The latter is due to the fact that excess zeros can often be explained by other model parameters, such as the theta parameter in the negative binomial.
  559. #'
  560. #' For this reason, results of the zero-inflation test should be interpreted as a residual pattern that can have many reasons, not as a decision criterion for whether or not to fit a ZIP model. To decide whether to add a ZIP term, I would advise relying on appropriate model selection techniques such as AIC, BIC, WAIC, Bayes factor, or LRT. Note that these tests are often not reliable in GLMMs because it is difficult to determine the df spent by the different models. The [simulateLRT] function in DHARMa provides a nonparametric alternative to obtain p-values for LRTs in nested models with unknown df.
  561. #'
  562. #' @author Florian Hartig
  563. #' @example inst/examples/testsHelp.R
  564. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  565. #' @export
  566. testZeroInflation <- function(simulationOutput, ...){
  567. countZeros <- function(x) sum( x == 0)
  568. testGeneric(simulationOutput = simulationOutput, summary = countZeros, methodName = "DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model", ... )
  569. }
  572. #' Test for a generic summary statistic based on simulated data
  573. #'
  574. #' This function tests if a user-defined summary differs when applied to simulated / observed data.
  575. #'
  576. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  577. #' @param summary a function that can be applied to simulated / observed data. See examples below.
  578. #' @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.
  579. #' @param plot whether to plot the simulated summary.
  580. #' @param methodName name of the test (will be used in plot).
  581. #'
  582. #' @details This function applies a user-defined summary to the simulated / observed data of a DHARMa object and then performs a hypothesis test using the ratio Obs / Sim as the test statistic.
  583. #'
  584. #' The summary is applied directly to the data and not to the residuals, but it can easily be remodeled to apply summaries to the residuals by simply defining something like f = function(x) summary (x - predictions), as done in [testDispersion].
  585. #'
  586. #' @note The summary function you specify will be applied to the data as it appears in your fitted model, which may not always be what you want.
  587. #'
  588. #' As an example, consider the case where we want to test for n-inflation in k/n data. If you provide your data via cbind (k, n-k), you have to test for n-inflation, but if you provide your data via k/n and weights = n, you should test for 1-inflation. When in doubt, check how the data is represented internally in model.frame(model) or via simulate(model).
  589. #'
  590. #' @export
  591. #' @author Florian Hartig
  592. #' @example inst/examples/testsHelp.R
  593. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  594. testGeneric <- function(simulationOutput, summary, alternative = c("two.sided", "greater", "less"), plot = T, methodName = "DHARMa generic simulation test"){
  596. out = list()
  597. out$data.name = deparse(substitute(simulationOutput))
  599. simulationOutput = ensureDHARMa(simulationOutput, convert = "Model")
  601. alternative <- match.arg(alternative)
  603. observed = summary(simulationOutput$observedResponse)
  605. simulated = apply(simulationOutput$simulatedResponse, 2, summary)
  607. p = getP(simulated = simulated, observed = observed, alternative = alternative)
  609. out$statistic = c(ratioObsSim = observed / mean(simulated))
  610. out$method = methodName
  611. out$alternative = alternative
  612. out$p.value = p
  615. class(out) = "htest"
  617. if(plot == T) {
  618. plotTitle = gsub('(.{1,50})(\\s|$)', '\\1\n', methodName)
  619. xLabel = paste("Simulated values, red line = fitted model. p-value (",out$alternative, ") = ", out$p.value, sep ="")
  620. hist(simulated, xlim = range(simulated, observed, na.rm=T ), col = "lightgrey", main = plotTitle, xlab = xLabel, breaks = max(round(simulationOutput$nSim / 5), 20), cex.main = 0.8)
  621. abline(v = observed, lwd= 2, col = "red")
  622. }
  623. return(out)
  624. }
  627. #' Test for temporal autocorrelation
  628. #'
  629. #' This function performs a standard test for temporal autocorrelation on the simulated residuals.
  630. #'
  631. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or by [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.
  632. #' @param time time variable in the same order as the datapoints and with the same length as the residuals in simulationOutput. If specified as formula, e.g. time = ~time, NAs are handled automatically (recommended).
  633. #' @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.
  634. #' @param plot whether to plot output.
  635. #' @details The function performs a Durbin-Watson test on the uniformly scaled residuals, and plots the residuals against time. The DB test was originally designed for normal residuals. In simulations, I didn't see a problem with this setting though. The alternative is to transform the uniform residuals to normal residuals and perform the DB test on those.
  636. #'
  637. #' Testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use [recalculateResiduals] function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testTemporalAutocorrelation.
  638. #'
  639. #' @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:
  640. #'
  641. #' 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.
  642. #'
  643. #' 2. Because the DHARMa residuals for such a model are not statistically independent anymore, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.
  644. #'
  645. #' There are three (non-exclusive) routes to address these issues when working with spatial / temporal / other autoregressive models:
  646. #'
  647. #' 1. Simulate conditional on the fitted CAR structures (see conditional simulations in the help of [simulateResiduals]).
  648. #'
  649. #' 2. Rotate simulations prior to residual calculations (see parameter rotation in [simulateResiduals]).
  650. #'
  651. #' 3. Use custom tests / plots that explicitly compare the correlation structure in the simulated data to the correlation structure in the observed data.
  652. #'
  653. #' @author Florian Hartig
  654. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  655. #' @example inst/examples/testTemporalAutocorrelationHelp.R
  656. #' @export
  657. testTemporalAutocorrelation <- function(simulationOutput, time, alternative = c("two.sided", "greater", "less"), plot = TRUE){
  659. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  661. alternative <- match.arg(alternative)
  663. if(is.null(time)){
  664. time = sample.int(simulationOutput$nObs, simulationOutput$nObs)
  665. message("DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point")
  666. }
  668. if(!inherits(time, "formula")){
  669. # old default
  670. # To avoid Issue #190
  671. if (length(time) != length(residuals(simulationOutput))) stop("Dimensions of time don't match the dimension of the residuals. Remove rows with NAs or use formula syntax for group to handle NAs automatically.")
  672. } else {
  673. time = getFormulaPredictors(simulationOutput, formula = time)[[1]] # formula syntax
  674. }
  676. # actually not sure if this is neccessary for dwtest, but seems better to aggregate
  677. if(any(duplicated(time))) stop("testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use the recalculateResiduals function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testTemporalAutocorrelation.")
  679. out = lmtest::dwtest(simulationOutput$scaledResiduals ~ 1, order.by = time, alternative = alternative)
  681. if(plot == T) {
  682. oldpar <- par(mfrow = c(1,2))
  683. on.exit(par(oldpar))
  685. plot(simulationOutput$scaledResiduals[order(time)] ~ time[order(time)],
  686. type = "l", ylab = "Scaled residuals", xlab = "Time", main = "Residuals vs. time", ylim = c(0,1))
  688. abline(h=c(0.5))
  689. abline(h=c(0,0.25,0.75,1), lty = 2 )
  691. acf(simulationOutput$scaledResiduals[order(time)], main = "Autocorrelation", ylim = c(-1,1))
  692. legend("topright",
  693. c(paste(out$method, " p=", round(out$p.value, digits = 5)),
  694. paste("Deviation ", ifelse(out$p.value < 0.05, "significant", "n.s."))),
  695. text.col = ifelse(out$p.value < 0.05, "red", "black" ), bty="n")
  697. }
  699. return(out)
  700. }
  703. #' Test for distance-based spatial (or similar type) autocorrelation
  704. #'
  705. #' This function performs a Moran's I test for distance-based spatial (or similar type) autocorrelation on the calculated quantile residuals.
  706. #'
  707. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or via [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.
  708. #' @param x the x coordinate, in the same order as the data points. Must be specified unless distMat is provided.
  709. #' @param y the y coordinate, in the same order as the data points. Must be specified unless distMat is provided.
  710. #' @param distMat optional distance matrix. Must be provided as class "matrix/array" (class "dist" is not supported). If distMat is not provided, euclidean distances based on x and y will be calculated. See details for explanation.
  711. #' @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.
  712. #' @param plot if T, and if x and y is provided, plot the output (see details).
  713. #'
  714. #' @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 be 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.
  715. #'
  716. #' If plot = T, a plot will be produced showing each residual at its x,y position, colored according to the residual value. Residuals with 0.5 are colored white, everything below 0.5 is colored increasingly red, everything above 0.5 is colored increasingly blue.
  717. #'
  718. #' 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 be 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].
  719. #'
  720. #' @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:
  721. #'
  722. #' 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.
  723. #'
  724. #' 2. Because the DHARMa residuals for such a model are not statistically independent anymore, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.
  725. #'
  726. #' There are three (non-exclusive) routes to address these issues when working with spatial / temporal / phylogenetic / other autoregressive models:
  727. #'
  728. #' 1. Simulate conditional on the fitted CAR structures (see conditional simulations in the help of [simulateResiduals]).
  729. #'
  730. #' 2. Rotate simulations prior to residual calculations (see parameter rotation in [simulateResiduals]).
  731. #'
  732. #' 3. Use custom tests / plots that explicitly compare the correlation structure in the simulated data to the correlation structure in the observed data.
  733. #'
  734. #' @author Florian Hartig
  735. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  736. #' @import grDevices
  737. #' @example inst/examples/testSpatialAutocorrelationHelp.R
  738. #' @export
  739. testSpatialAutocorrelation <- function(simulationOutput, x = NULL, y = NULL, distMat = NULL, alternative = c("two.sided", "greater", "less"), plot = TRUE){
  741. alternative <- match.arg(alternative)
  742. data.name = deparse(substitute(simulationOutput)) # needs to be before ensureDHARMa
  743. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  745. # Assertions
  747. if( (!is.null(x) | !is.null(y)) & !is.null(distMat) ){
  748. if(isTRUE(plot)) {message("Both coordinates and distMat provided, calculations will be done based on the distance matrix, coordinates will only be used for plotting.")}
  750. else { # set x,y to NULL so they don't cause errors if wrong dimensions
  751. x = NULL
  752. y = NULL
  753. message("Both coordinates and distMat provided, calculations will be done based on the distance matrix, coordinates will be ignored.")
  754. }
  755. }
  757. if( (is.null(x) | is.null(y)) & is.null(distMat) ) stop("You need to provide either x,y coordinates or a distMatrix.")
  759. # check if x,y are formulas
  760. if(!inherits(x, "formula") | !inherits(y, "formula")) {
  761. # old default
  762. # To avoid Issue #190
  763. if (is.null(distMat) & !is.null(x) & length(x) != length(residuals(simulationOutput)) | is.null(distMat) & !is.null(y) & length(y) != length(residuals(simulationOutput))) stop("Dimensions of x / y coordinates don't match the dimension of the residuals. Remove rows with NAs or specify x,y as formula to handle NAs automatically.")
  765. if(inherits(x, "formula") != inherits(y, "formula")) stop("Please specify x and y the same way, either both as formula (e.g.x = ~x) or both as vector (e.g. x = data$x).")
  766. } else {
  767. # formula syntax
  768. x = getFormulaPredictors(simulationOutput, formula = x)[[1]]
  769. y = getFormulaPredictors(simulationOutput, formula = y)[[1]]
  770. }
  772. # if not provided, create distance matrix based on x and y
  773. if(is.null(distMat)) distMat <- as.matrix(dist(cbind(x, y)))
  775. # check for duplicates in x,y and distMat (off-diagonal zero values) (#73)
  776. testDistMat = distMat
  777. diag(testDistMat) = NA # set diagonal to NA
  778. if(any(duplicated(cbind(x,y))) | any(testDistMat == 0, na.rm = TRUE)) stop("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 be 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.")
  780. invDistMat <- 1/distMat
  781. diag(invDistMat) <- 0
  783. MI = ape::Moran.I(simulationOutput$scaledResiduals, weight = invDistMat, alternative = alternative)
  785. out = list()
  786. out$statistic = c(observed = MI$observed, expected = MI$expected, sd = MI$sd)
  787. out$method = "DHARMa Moran's I test for distance-based autocorrelation"
  788. out$alternative = "Distance-based autocorrelation"
  789. out$p.value = MI$p.value
  790. out$data.name = data.name
  792. class(out) = "htest"
  794. if(plot == T & !is.null(x) & !is.null(y)) {
  795. opar = par(no.readonly = TRUE)
  796. on.exit(par(opar))
  798. col = colorRamp(c("red", "white", "blue"))(simulationOutput$scaledResiduals)
  799. layout(matrix(c(1, 2), ncol = 2), widths = c(6, 1))
  801. # scatterplot
  802. par(mar = c(5, 4, 4, 1))
  803. plot(x, y,
  804. col = rgb(col, maxColorValue = 255),
  805. main = "Standardized residuals in space",
  806. cex.main = 0.8)
  808. mtext(paste0(out$method, "\n",
  809. " p = ", round(out$p.value, digits = 5), ", Deviation ", ifelse(out$p.value < 0.05, "significant", "n.s.")),
  810. col = ifelse(out$p.value < 0.05, "red", "black"),
  811. cex = 0.8,
  812. side = 3, line = 0)
  814. # legend
  815. par(mar = c(1, 0.5, 2, 2))
  816. plot.new()
  817. plot.window(xlim = c(0, 1), ylim = c(-1, 1))
  819. legend_image = as.raster(matrix(colorRampPalette(c("blue", "white", "red"))(200), ncol = 1))
  820. rasterImage(legend_image,
  821. xleft = -0.5, xright = 1,
  822. ybottom = -0.1, ytop = 0.5)
  824. axis(4, at = c(-0.1, 0.5), labels = c("0", "1"), las = 1)
  825. mtext("Standardized \n residuals", side = 3, line = -5,cex = 0.8,)
  827. # TODO implement correlogram
  828. }
  830. return(out)
  831. }
  834. #' Test for phylogenetic autocorrelation
  835. #'
  836. #' This function performs a Moran's I test for phylogenetic autocorrelation on the calculated quantile residuals.
  837. #'
  838. #' @param simulationOutput an object of class DHARMa, either created via [simulateResiduals] for supported models or via [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.
  839. #' @param tree a phylogenetic tree object.
  840. #' @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 of no phylogenetic correlation.
  841. #'
  842. #' @details The function performs Moran.I test from the package ape on the DHARMa residuals, based on the phylogenetic distance matrix internally created from the provided tree. For custom distance matrices, you can use [testSpatialAutocorrelation].
  843. #'
  844. #' @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:
  845. #'
  846. #' 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.
  847. #'
  848. #' 2. Because the DHARMa residuals for such a model are not statistically independent anymore, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.
  849. #'
  850. #' There are three (non-exclusive) routes to address these issues when working with spatial / temporal / phylogenetic autoregressive models:
  851. #'
  852. #' 1. Simulate conditional on the fitted CAR structures (see conditional simulations in the help of [simulateResiduals]).
  853. #'
  854. #' 2. Rotate simulations prior to residual calculations (see parameter rotation in [simulateResiduals]).
  855. #'
  856. #' 3. Use custom tests / plots that explicitly compare the correlation structure in the simulated data to the correlation structure in the observed data.
  857. #'
  858. #' @author Florian Hartig
  859. #' @seealso [testResiduals], [testUniformity], [testOutliers], [testDispersion], [testZeroInflation], [testGeneric], [testTemporalAutocorrelation], [testSpatialAutocorrelation], [testQuantiles], [testCategorical]
  860. #' @example inst/examples/testPhylogeneticAutocorrelationHelp.R
  861. #' @export
  862. testPhylogeneticAutocorrelation <- function(simulationOutput,
  863. tree,
  864. alternative = c("two.sided", "greater", "less")){
  866. alternative <- match.arg(alternative)
  867. data.name = deparse(substitute(simulationOutput)) # needs to be before ensureDHARMa
  868. simulationOutput = ensureDHARMa(simulationOutput, convert = T)
  870. # calculate distance matrix
  871. distMat <- cophenetic(tree)
  872. invDistMat <- 1/distMat
  873. diag(invDistMat) <- 0
  875. MI = ape::Moran.I(simulationOutput$scaledResiduals, weight = invDistMat, alternative = alternative)
  877. out = list()
  878. out$statistic = c(observed = MI$observed, expected = MI$expected, sd = MI$sd)
  879. out$method = "DHARMa Moran's I test for phylogenetic autocorrelation"
  880. out$alternative = "Phylogenetic autocorrelation"
  881. out$p.value = MI$p.value
  882. out$data.name = data.name
  884. class(out) = "htest"
  885. return(out)
  886. }
  890. getP <- function(simulated, observed, alternative, plot = FALSE, ...){
  892. if(alternative == "greater") p = mean(simulated >= observed)
  893. if(alternative == "less") p = mean(simulated <= observed)
  894. if(alternative == "two.sided") p = min(min(mean(simulated <= observed), mean(simulated >= observed) ) * 2,1)
  896. if(plot == T){
  897. hist(simulated, xlim = range(simulated, observed), col = "lightgrey", main = "Distribution of test statistic \n grey = simulated, red = observed", ...)
  898. abline(v = mean(simulated), col = 1, lwd = 2)
  899. abline(v = observed, col = "red", lwd = 2)
  900. }
  902. return(p)
  903. }

tests.R at commit c9dbcee, under gnu · at the source

Overview

Authors: Gustavo A Rodriguez1, Andrew Aoun1,2,3, Eva F Rothenberg1, C Oliver Shetler1, Lorenzo Posani4,5,6, Srujan V Vajram1, Thomas Tedesco1, Anurag Sharma2, Stefano Fusi4,5,6,7, S Abid Hussaini1,2,8,9
ORCID iDs: S Abid Hussaini
  1. Taub Institute for Research on Alzheimer’s Disease and the Aging Brain, Columbia University Irving Medical Center, New York, NY 10032, USA
  2. Burke Neurological Institute, White Plains, NY 10605, USA
  3. Center for Neural Science, New York University, New York, NY 10003, USA
  4. Department of Neuroscience, Columbia University Irving Medical Center, New York, NY 10027, USA
  5. Center for Theoretical Neuroscience, Columbia University, New York, NY 10027, USA
  6. Zuckerman Mind Brain Behavior Institute, Columbia University, New York, NY 10027, USA
  7. Kavli Institute for Brain Science, Columbia University, New York, NY 10027, USA
  8. Department of Pathology and Cell Biology, Columbia University Irving Medical Center, New York, NY 10032, USA
  9. Lead contact
Institutions: Columbia University Irving Medical Center (United States); Columbia University (United States); New York University (United States); Mortimer B. Zuckerman Mind Brain Behavior Institute (United States)
Journal: Cell reports, volume 45, issue 6, article 117505
Dates: published online 6 June 2026; in print 23 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.celrep.2026.117505 · PMID 42250220 · PMCID PMC13560757 · OpenAlex W7163717440
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: extracellular electrophysiology (units, LFP) (modality), mouse (organism), Alzheimer's / dementia (population)
Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions
Keywords: Amyloid beta, Medial Entorhinal Cortex, Aβ Plaque, In Vivo Electrophysiology, Neuronal Hyperactivity, Spatial Remapping, App Knock-in, Cp: Neuroscience, Spatial Decoding, App Nl-g-f, Ratemap Stability
MeSH: Aging*, Amyloid beta-Protein Precursor*, Entorhinal Cortex*, Neurons*, Alzheimer Disease, Amyloid beta-Peptides, Animals, Gene Knock-In Techniques, Interneurons, Male, Mice, Mice, Inbred C57BL, Mice, Transgenic (* major topic)
Topic: Memory and Neural Mechanisms (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIA NIH HHS (RF1 AG080818, R01 AG050425, K01 AG068598, R01 AG064066); BrightFocus Foundation (A2019382F); National Institute on Aging (R01AG064066, R01AG050425, K01AG068598, RF1AG080818)
Citations: cited by 1 paper (Europe PMC); 142 references in the paper
Research resources: RRID:AB_141930, RRID:AB_143165, RRID:AB_2619988, RRID:AB_2797642, RRID:IMSR_JAX:000664, MATLAB RRID:SCR_001622, R v4.4.1 RRID:SCR_001905, Fiji (ImageJ) RRID:SCR_002285, GraphPad Prism v10.4.0 RRID:SCR_002798, Python RRID:SCR_008394

Abstract

Advanced amyloid beta (Aβ) pathology is associated with aberrant neuronal network activity and cognitive impairment in preclinical Alzheimer’s disease (AD) models. Here, we assess Aβ pathology’s impact on spatial information processing in the medial entorhinal cortex (MEC) of 18-month AppNL-G-F/NL-G-F knock-in (APP KI) mice during exploration of open field arenas. Spatial information scores are decreased in APP KI MEC neurons versus age-matched controls. Border cell firing preferences are unstable across sessions and grid cell spatial periodicity is disrupted. Ratemap stability analysis using the Earth Mover’s Distance indicates increased instability in spatially tuned APP KI neurons. Spatial decoding analysis indicates deficits in position and speed coding in APP KI mice across all comparisons. Additionally, APP KI mice display a mild hyperactive phenotype driven by narrow-spiking putative interneurons. These findings tie Aβ-associated dysregulation in neuronal firing to disruptions in spatial information processing that may underlie cognitive deficits associated with AD.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

Zenodo 19040105

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the resources table
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
1 file
At the source: zenodo.org/records/19040105

klusta-team/klustakwik

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d8b750107fef53996957d774dfa3ec71f2442462, 21 April 2015
Languages: C/C++ (8), C++ (8), Python (2)
Size: 26 files, 18 scripts
Software Heritage: archived
Found in: the resources table
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
19 files
At the source: github.com/klusta-team/klustakwik

glmmTMB/glmmTMB

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: db725d143d529c93c55cf7baa7402dc5135fe6ff, 23 September 2026
Languages: R (170), JavaScript (20), C/C++ (3), C++ (2), Shell (2), Quarto (2), Jupyter (1)
Size: 487 files, 200 scripts
Software Heritage: archived
Found in: the resources table
Holds: README, environment (glmmTMB/DESCRIPTION, misc/dockerfile), tests, continuous integration, documentation, 18 notebooks
Not found: license file, CITATION.cff
Tools: glmmTMB (126 files), lme4 (32 files), tidyverse (15 files), ggplot2 (10 files), car (7 files), emmeans (7 files), nlme (6 files), mgcv (5 files), broom (4 files), easystats (3 files), reshape2 (3 files), brms (2 files), lmerTest (2 files), metafor (2 files), patchwork (2 files), Stan (2 files), NumPy (1 file), SymPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
201 files
At the source: github.com/glmmTMB/glmmTMB

florianhartig/DHARMa

License: gnu
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c9dbcee2eec12904e505dd8c7c4f1f8d61ec1f79, 25 September 2026
Languages: R (142)
Size: 305 files, 142 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, environment (DHARMa/DESCRIPTION), tests, continuous integration, documentation, 18 notebooks
Not found: license file, CITATION.cff
Tools: lme4 (55 files), glmmTMB (27 files), mgcv (15 files), JAGS (5 files), tidyverse (5 files), brms (4 files), easystats (4 files), nlme (4 files), ggplot2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
143 files
At the source: github.com/florianhartig/DHARMa

hussainilab/rodriguez-et-al-cell-reports-2026

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 7772cfa1cce93cad101419229805151e27f5beea, 14 May 2026
Size: 6 files, 0 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
1 file
At the source: github.com/HussainiLab/Rodriguez-et-al-Cell-Reports-2026/tree/Data

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 5 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 360 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data and code availability

Source data for all figures will be made freely available on the Hussaini Lab GitHub page (https://doi.org/10.5281/zenodo.19040105).

All original code used for data analysis will be made freely available on the Hussaini Lab GitHub page.

Any additional information regarding data is available from the lead contact upon request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Cell Press
  • Authors: added S Abid Hussaini (0000-0002-3921-9021); removed S Abid Hussaini

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 11 keywords, 13 MeSH terms, 3 funders, 140 references, 10 RRIDs.

Cite

This paper

Rodriguez, G. A., Aoun, A., Rothenberg, E. F., Shetler, C. O., Posani, L., Vajram, S. V., Tedesco, T., Sharma, A., Fusi, S., & Hussaini, S. A. (2026). Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice. Cell reports, 45(6), 117505. https://doi.org/10.1016/j.celrep.2026.117505

BibTeX

@article{rodriguez2026impaired,
author = {Rodriguez, Gustavo A and Aoun, Andrew and Rothenberg, Eva F and Shetler, C Oliver and Posani, Lorenzo and Vajram, Srujan V and Tedesco, Thomas and Sharma, Anurag and Fusi, Stefano and Hussaini, S Abid},
title = {{Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice}},
journal = {Cell reports},
year = {2026},
month = jun,
volume = {45},
number = {6},
pages = {117505},
publisher = {Cell Press},
issn = {2211-1247},
doi = {10.1016/j.celrep.2026.117505},
url = {https://doi.org/10.1016/j.celrep.2026.117505},
pmid = {42250220},
pmcid = {PMC13560757}
}

RIS

TY - JOUR
AU - Rodriguez, Gustavo A
AU - Aoun, Andrew
AU - Rothenberg, Eva F
AU - Shetler, C Oliver
AU - Posani, Lorenzo
AU - Vajram, Srujan V
AU - Tedesco, Thomas
AU - Sharma, Anurag
AU - Fusi, Stefano
AU - Hussaini, S Abid
TI - Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice
T2 - Cell reports
J2 - Cell Rep
PY - 2026
DA - 2026/06/06
VL - 45
IS - 6
SP - 117505
SN - 2211-1247
PB - Cell Press
DO - 10.1016/j.celrep.2026.117505
UR - https://doi.org/10.1016/j.celrep.2026.117505
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.celrep.2026.117505",
"type": "article-journal",
"title": "Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice",
"container-title": "Cell reports",
"author": [
{
"family": "Rodriguez",
"given": "Gustavo A"
},
{
"family": "Aoun",
"given": "Andrew"
},
{
"family": "Rothenberg",
"given": "Eva F"
},
{
"family": "Shetler",
"given": "C Oliver"
},
{
"family": "Posani",
"given": "Lorenzo"
},
{
"family": "Vajram",
"given": "Srujan V"
},
{
"family": "Tedesco",
"given": "Thomas"
},
{
"family": "Sharma",
"given": "Anurag"
},
{
"family": "Fusi",
"given": "Stefano"
},
{
"family": "Hussaini",
"given": "S Abid"
}
],
"container-title-short": "Cell Rep",
"volume": "45",
"issue": "6",
"page": "117505",
"DOI": "10.1016/j.celrep.2026.117505",
"PMID": "42250220",
"PMCID": "PMC13560757",
"ISSN": "2211-1247",
"publisher": "Cell Press",
"URL": "https://doi.org/10.1016/j.celrep.2026.117505",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
6
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.celrep.2026.117646 [code]
Medial entorhinal-hippocampal desynchronization parallels the emergence of memory impairment in a mouse model of Alzheimer's disease pathology.
Journal: Cell reports
In common: car, broom, emmeans, 7 other tools, Alzheimer's / dementia, mouse, 19 references
[2] doi:10.1038/s41467-026-69866-3 [code]
Selective weakening of population-coupled synaptic activity in vivo in a mouse model of amyloid-beta pathology.
Journal: Nature communications
In common: patchwork, ggplot2, tidyverse, Alzheimer's / dementia, mouse, 17 references
[3] doi:10.1016/j.isci.2026.116747 [code]
Age and loneliness relate to reduced trust learning and alterations in amygdala function.
Journal: iScience
In common: Stan, brms, nlme, 9 other tools
[4] doi:10.1111/ejn.70480 [code]
Astrocyte Proximity Protects Synapses From Human Amyloid-Beta Induced Degeneration in a Mouse Ex Vivo Model of Early Alzheimer's Disease.
Journal: The European journal of neuroscience
In common: emmeans, lmerTest, lme4, 6 other tools, Alzheimer's / dementia, mouse, 6 references
[5] doi:10.1038/s43587-026-01096-0 [code]
Neuronal APOE4-induced early hippocampal network hyperexcitability in Alzheimer's disease pathogenesis.
Journal: Nature aging
In common: pandas, Matplotlib, NumPy, Alzheimer's / dementia, mouse, 10 references
[6] doi:10.1038/s41467-026-73865-9 [code]
Histamine shapes the neurocomputational dynamics of human learning.
Journal: Nature communications
In common: metafor, brms, easystats, 10 other tools
[7] doi:10.1002/hbm.70605 [code]
BrainEnrich: Revealing Biological Insights for Imaging-Derived Phenotypes Through Transcriptomic Enrichment.
Journal: Human brain mapping
In common: JAGS, nlme, easystats, 10 other tools
[8] doi:10.64898/2026.03.02.709173 [code]
Corpus Callosum Dysgenesis impairs metacognition: evidence from multi-modality and multi-cohort replications
Journal: bioRxiv (preprint)
In common: JAGS, Stan, brms, 7 other tools
[9] doi:10.1038/s41467-025-62798-4 [code]
Distinct manifestations of excitatory-inhibitory imbalance associated with amyloid-β and tau in patients with Alzheimer’s disease
Journal: n/a
In common: Alzheimer's / dementia, 11 references
[10] doi:10.1126/sciadv.aec9291 [code]
Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks.
Journal: Science advances
In common: Stan, nlme, easystats, 8 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.