This function compares the observed number of zeros with the zeros expected from simulations.
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.
- ...
further arguments to testGeneric.
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.
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.
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.
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.
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.11641, p-value = 0.133
#> alternative hypothesis: two-sided
#>
#>
#> $dispersion
#>
#> DHARMa nonparametric dispersion test via sd of residuals fitted vs.
#> simulated
#>
#> data: simulationOutput
#> dispersion = 1.2702, p-value = 0.144
#> alternative hypothesis: two.sided
#>
#>
#> $outliers
#>
#> DHARMa bootstrapped outlier test
#>
#> data: simulationOutput
#> outliers at both margin(s) = 0, observations = 100, p-value = 1
#> alternative hypothesis: two.sided
#> percent confidence interval:
#> 0.00 0.02
#> sample estimates:
#> outlier frequency (expected: 0.0029 )
#> 0
#>
#>
####### Individual tests #######
# KS test for correct distribution of residuals
testUniformity(simulationOutput)
#>
#> Asymptotic one-sample Kolmogorov-Smirnov test
#>
#> data: simulationOutput$scaledResiduals
#> D = 0.11641, p-value = 0.133
#> 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.25837, p-value = 0.4429
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 2
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.2083, p-value = 0.7053
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 3
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.20257, p-value = 0.7355
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 4
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.17073, p-value = 0.887
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 5
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.19125, p-value = 0.7932
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 6
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.22906, p-value = 0.5936
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 7
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.14971, p-value = 0.9547
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 8
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.21376, p-value = 0.6762
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 9
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.13682, p-value = 0.9793
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 10
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.32628, p-value = 0.1899
#> alternative hypothesis: two-sided
#>
#>
#> $uniformity$p.value
#> [1] 0.4429026 0.7053317 0.7354573 0.8870295 0.7932444 0.5936246 0.9546701
#> [8] 0.6761893 0.9792570 0.1898504
#>
#> $uniformity$p.value.cor
#> [1] 1 1 1 1 1 1 1 1 1 1
#>
#>
#> $homogeneity
#> Levene's Test for Homogeneity of Variance (center = median)
#> Df F value Pr(>F)
#> group 9 0.0377 1
#> 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.25837, p-value = 0.4429
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 2
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.2083, p-value = 0.7053
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 3
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.20257, p-value = 0.7355
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 4
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.17073, p-value = 0.887
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 5
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.19125, p-value = 0.7932
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 6
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.22906, p-value = 0.5936
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 7
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.14971, p-value = 0.9547
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 8
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.21376, p-value = 0.6762
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 9
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.13682, p-value = 0.9793
#> alternative hypothesis: two-sided
#>
#> ------------------------------------------------------------
#> catPred: 10
#>
#> Exact one-sample Kolmogorov-Smirnov test
#>
#> data: dd[x, ]
#> D = 0.32628, p-value = 0.1899
#> alternative hypothesis: two-sided
#>
#>
#> $uniformity$p.value
#> [1] 0.4429026 0.7053317 0.7354573 0.8870295 0.7932444 0.5936246 0.9546701
#> [8] 0.6761893 0.9792570 0.1898504
#>
#> $uniformity$p.value.cor
#> [1] 1 1 1 1 1 1 1 1 1 1
#>
#>
#> $homogeneity
#> Levene's Test for Homogeneity of Variance (center = median)
#> Df F value Pr(>F)
#> group 9 0.0377 1
#> 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.2702, p-value = 0.144
#> 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) = 1, observations = 100, p-value = 0.5507
#> alternative hypothesis: true probability of success is not equal to 0.007968127
#> 95 percent confidence interval:
#> 0.000253146 0.054459385
#> sample estimates:
#> frequency of outliers (expected: 0.00796812749003984 )
#> 0.01
#>
# testing zero inflation
testZeroInflation(simulationOutput)
#>
#> DHARMa zero-inflation test via comparison to expected zeros with
#> simulation under H0 = fitted model
#>
#> data: simulationOutput
#> ratioObsSim = 1.1058, p-value = 0.488
#> 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.79533, p-value = 0.192
#> alternative hypothesis: two.sided
#>
testGeneric(simulationOutput, summary = countOnes, alternative = "less") # 1-deficit
#>
#> DHARMa generic simulation test
#>
#> data: simulationOutput
#> ratioObsSim = 0.79533, p-value = 0.096
#> 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.005, p-value = 0.968
#> 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.1021, p-value = 0.264
#> alternative hypothesis: two.sided
#>