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.
Usage
testCategorical(simulationOutput, catPred, quantiles = c(0.25, 0.5, 0.75),
plot = TRUE, ...)Arguments
- 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.
- 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).
- quantiles
whether to draw the quantile lines.
- plot
if TRUE, the function will create an additional plot.
- ...
additional arguments to boxplot.
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.
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.
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.
Examples
testData = createData(sampleSize = 100, overdispersion = 0.5, randomEffectVariance = 0)
fittedModel <- glm(observedResponse ~ Environment1 , family = "poisson", data = testData)
simulationOutput <- simulateResiduals(fittedModel = fittedModel)
# the plot function shows 2 plots and runs 4 tests
# i) KS test i) Dispersion test iii) Outlier test iv) quantile test
plot(simulationOutput, quantreg = TRUE)
# testResiduals tests distribution, dispersion and outliers
testResiduals(simulationOutput)
#> $uniformity
#>
#> Asymptotic one-sample Kolmogorov-Smirnov test
#>
#> data: simulationOutput$scaledResiduals
#> D = 0.067079, p-value = 0.7591
#> alternative hypothesis: two-sided
#>
#>
#> $dispersion
#>
#> DHARMa nonparametric dispersion test via sd of residuals fitted vs.
#> simulated
#>
#> data: simulationOutput
#> dispersion = 1.5167, p-value < 2.2e-16
#> alternative hypothesis: two.sided
#>
#>
#> $outliers
#>
#> DHARMa bootstrapped outlier test
#>
#> data: simulationOutput
#> outliers at both margin(s) = 1, observations = 100, p-value = 0.6
#> alternative hypothesis: two.sided
#> percent confidence interval:
#> 0.00 0.02
#> sample estimates:
#> outlier frequency (expected: 0.0037 )
#> 0.01
#>
#>
####### Individual tests #######
# KS test for correct distribution of residuals
testUniformity(simulationOutput)
#>
#> Asymptotic one-sample Kolmogorov-Smirnov test
#>
#> data: simulationOutput$scaledResiduals
#> D = 0.067079, p-value = 0.7591
#> alternative hypothesis: two-sided
#>
# KS test for correct distribution within and between groups
# group specified as formula (recommended)
testCategorical(simulationOutput, ~group)
#> $uniformity
#> $uniformity$details
#> catPred: 1
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.25979, p-value = 0.4361
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 2
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.23865, p-value = 0.5426
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 3
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.28998, p-value = 0.3068
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 4
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.46799, p-value = 0.01578
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 5
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.19479, p-value = 0.7754
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 6
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.15653, p-value = 0.9365
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 7
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.15813, p-value = 0.9316
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 8
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.34788, p-value = 0.1384
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 9
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.132, p-value = 0.9854
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 10
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.21925, p-value = 0.6466
#> alternative hypothesis: two-sided
#>
#>
#> $uniformity$p.value
#> [1] 0.43611379 0.54259699 0.30679427 0.01577848 0.77536994 0.93645211
#> [7] 0.93164870 0.13835501 0.98537937 0.64660322
#>
#> $uniformity$p.value.cor
#> [1] 1.0000000 1.0000000 1.0000000 0.1577848 1.0000000 1.0000000 1.0000000
#> [8] 1.0000000 1.0000000 1.0000000
#>
#>
#> $homogeneity
#> Levene's Test for Homogeneity of Variance (center = median)
#> Df F value Pr(>F)
#> group 9 0.2544 0.9847
#> 90
#>
# group specified as variable in your environment
testCategorical(simulationOutput, testData$group)
#> $uniformity
#> $uniformity$details
#> catPred: 1
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.25979, p-value = 0.4361
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 2
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.23865, p-value = 0.5426
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 3
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.28998, p-value = 0.3068
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 4
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.46799, p-value = 0.01578
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 5
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.19479, p-value = 0.7754
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 6
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.15653, p-value = 0.9365
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 7
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.15813, p-value = 0.9316
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 8
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.34788, p-value = 0.1384
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 9
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.132, p-value = 0.9854
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 10
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.21925, p-value = 0.6466
#> alternative hypothesis: two-sided
#>
#>
#> $uniformity$p.value
#> [1] 0.43611379 0.54259699 0.30679427 0.01577848 0.77536994 0.93645211
#> [7] 0.93164870 0.13835501 0.98537937 0.64660322
#>
#> $uniformity$p.value.cor
#> [1] 1.0000000 1.0000000 1.0000000 0.1577848 1.0000000 1.0000000 1.0000000
#> [8] 1.0000000 1.0000000 1.0000000
#>
#>
#> $homogeneity
#> Levene's Test for Homogeneity of Variance (center = median)
#> Df F value Pr(>F)
#> group 9 0.2544 0.9847
#> 90
#>
# Dispersion test - for details see ?testDispersion
testDispersion(simulationOutput) # tests under and overdispersion
#>
#> DHARMa nonparametric dispersion test via sd of residuals fitted vs.
#> simulated
#>
#> data: simulationOutput
#> dispersion = 1.5167, p-value < 2.2e-16
#> alternative hypothesis: two.sided
#>
# Outlier test (number of observations outside simulation envelope)
# Use type = "boostrap" for exact values, see ?testOutliers
testOutliers(simulationOutput, type = "binomial")
#>
#> DHARMa outlier test based on exact binomial test with approximate
#> expectations
#>
#> data: simulationOutput
#> outliers at both margin(s) = 2, observations = 100, p-value = 0.1898
#> alternative hypothesis: true probability of success is not equal to 0.007968127
#> 95 percent confidence interval:
#> 0.002431337 0.070383932
#> sample estimates:
#> frequency of outliers (expected: 0.00796812749003984 )
#> 0.02
#>
# testing zero inflation
testZeroInflation(simulationOutput)
#>
#> DHARMa zero-inflation test via comparison to expected zeros with
#> simulation under H0 = fitted model
#>
#> data: simulationOutput
#> ratioObsSim = 1.1814, p-value = 0.328
#> alternative hypothesis: two.sided
#>
# testing generic summaries
countOnes <- function(x) sum(x == 1) # testing for number of 1s
testGeneric(simulationOutput, summary = countOnes) # 1-inflation
#>
#> DHARMa generic simulation test
#>
#> data: simulationOutput
#> ratioObsSim = 0.913, p-value = 0.576
#> alternative hypothesis: two.sided
#>
testGeneric(simulationOutput, summary = countOnes, alternative = "less") # 1-deficit
#>
#> DHARMa generic simulation test
#>
#> data: simulationOutput
#> ratioObsSim = 0.913, p-value = 0.288
#> alternative hypothesis: less
#>
means <- function(x) mean(x) # testing if mean prediction fits
testGeneric(simulationOutput, summary = means)
#>
#> DHARMa generic simulation test
#>
#> data: simulationOutput
#> ratioObsSim = 1.0045, p-value = 0.96
#> alternative hypothesis: two.sided
#>
spread <- function(x) sd(x) # testing if mean sd fits
testGeneric(simulationOutput, summary = spread)
#>
#> DHARMa generic simulation test
#>
#> data: simulationOutput
#> ratioObsSim = 1.185, p-value = 0.024
#> alternative hypothesis: two.sided
#>