In this chapter, you will learn
gls or glmmTMB (Section 10.1);Workflow steps (Chapter 12): step 5 (residual checks: changing variance).
Prerequisites: the linear model and its residual checks (Chapter 2), GLMMs (Chapter 9) for the second part.
Modeling dispersion means that we describe how the expected variance changes as a function of the mean or predictor variables. The residual problem that we address by that is known as heteroskedasticity (see Section 2.5.1 and Section 2.4.1).
Heteroskedasticity means that the residual variance is not constant, for example because it depends on a predictor or on the mean. The easiest way to understand it is in the linear model. Let’s look at an extreme example that we create by simulating some data
set.seed(125)
data = data.frame(treatment = factor(rep(c("A", "B", "C"), each = 15)))
data$observation = c(7, 2 ,4)[as.numeric(data$treatment)] +
rnorm( length(data$treatment), sd = as.numeric(data$treatment)^2 )
boxplot(observation ~ treatment, data = data)
Especially p-values and confidence intervals of lm() and ANOVA can react quite strongly to such differences in residual variation. So, running a standard lm() / ANOVA on this data is not a good idea - in this case, we see that no regression effect is significant, and neither is the ANOVA (p = 0.27), although the true group means differ (7, 2, 4, we simulated them). The unequal residual variances distort the standard errors and p-values (here they cost us power); let’s see what happens when we model the variances.
fit = lm(observation ~ treatment, data = data)
summary(fit)
Call:
lm(formula = observation ~ treatment, data = data)
Residuals:
Min 1Q Median 3Q Max
-17.2897 -1.0514 0.3531 2.4465 19.8602
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.043 1.731 4.069 0.000204 ***
treatmentB -3.925 2.448 -1.603 0.116338
treatmentC -1.302 2.448 -0.532 0.597601
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.704 on 42 degrees of freedom
Multiple R-squared: 0.05973, Adjusted R-squared: 0.01495
F-statistic: 1.334 on 2 and 42 DF, p-value: 0.2744summary(aov(fit)) Df Sum Sq Mean Sq F value Pr(>F)
treatment 2 119.9 59.95 1.334 0.274
Residuals 42 1887.6 44.94 A subtlety about heteroskedasticity in the normal model is that it especially affects p-values when the residual variance changes with the predictor itself, as illustrated by the simulation below (no true effect, so p-values should be uniform)
par(mfrow = c(1, 2))
# error sd is variable, but uncorrelated with predictor
getP <- function(effect = 0, n = 1000){
sd = rgamma(n, shape = 2)
#hist(sd)
x = runif(n)
y = effect * x + rnorm(n, sd = sd)
fit <- lm(y ~ x)
table <- summary(fit)
return(table$coefficients[2,4])
}
out = replicate(1000, getP())
hist(out, breaks = 20, main = "variance independent of x")
# error sd is variable and correlated with predictor (identical in this case)
getP2 <- function(effect = 0, n = 1000){
sd = rgamma(n, shape = 2)
y = effect * sd + rnorm(n, sd = sd)
fit <- lm(y ~ sd)
table <- summary(fit)
return(table$coefficients[2,4])
}
out2 = replicate(1000, getP2())
hist(out2, breaks = 20, main = "variance depends on x")
Left: the p-values are uniform under H0, i.e. calibrated (about 5% false positives at alpha = 0.05). Right: strong excess of small p-values (here about 30% false positives): inflated type I error.
So, what can we do?
One option is to search for a transformation of the response that improves the problem - If heteroskedasticity correlates with the mean value, one can typically decrease it by some sqrt or log transformation, but this is often difficult, because this may also conflict with keeping the distribution normal.
The second, more general option is to model the variance, i.e. to fit a model in which the variance is not fixed. We will discuss two packages in R that allow you to model the dispersion.
The first (traditional) option is to use nlme::gls. GLS = Generalized Least Squares (not to be confused with GLMs, Chapter 9: GLS keeps the normal distribution, but allows non-constant variance and correlated errors). In the gls function, you can specify a dependency of the residual variance on a predictor or the response via the weights argument. There are different types of dependencies that you can specify, see ?varFunc. In our case, we will use the varIdent function, which allows you to specify a different variance per treatment.
library(nlme)
fit = gls(observation ~ treatment, data = data,
weights = varIdent(form = ~ 1 | treatment))
summary(fit)Generalized least squares fit by REML
Model: observation ~ treatment
Data: data
AIC BIC logLik
243.9258 254.3519 -115.9629
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | treatment
Parameter estimates:
A B C
1.000000 4.714512 11.821868
Coefficients:
Value Std.Error t-value p-value
(Intercept) 7.042667 0.2348387 29.989388 0.0000
treatmentB -3.925011 1.1317816 -3.467994 0.0012
treatmentC -1.302030 2.7861462 -0.467323 0.6427
Correlation:
(Intr) trtmnB
treatmentB -0.207
treatmentC -0.084 0.017
Standardized residuals:
Min Q1 Med Q3 Max
-2.4587934 -0.6241702 0.1687727 0.6524558 1.9480170
Residual standard error: 0.9095262
Degrees of freedom: 45 total; 42 residualIf we check the ANOVA, we see that, unlike before, we get a significant effect of the treatment
anova(fit)Denom. DF: 42
numDF F-value p-value
(Intercept) 1 899.3761 <.0001
treatment 2 6.0962 0.0047The second option for modeling variances is to use the glmmTMB package. Here, you can specify an extra regression formula for the dispersion. For the Gaussian family, this formula models the residual standard deviation on the log scale. If we fit this:
library(glmmTMB)
fit = glmmTMB(observation ~ treatment, data = data,
dispformula = ~ treatment)
summary(fit) Family: gaussian ( identity )
Formula: observation ~ treatment
Dispersion: ~treatment
Data: data
AIC BIC logLik -2*log(L) df.resid
248.7 259.5 -118.3 236.7 39
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.0427 0.2269 31.042 < 2e-16 ***
treatmentB -3.9250 1.0934 -3.590 0.000331 ***
treatmentC -1.3020 2.6917 -0.484 0.628581
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Dispersion model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1293 0.1826 -0.708 0.479
treatmentB 1.5506 0.2582 6.006 1.91e-09 ***
treatmentC 2.4700 0.2582 9.566 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We get 2 regression tables as outputs - one for the effects, and one for the dispersion. We see, as expected, that the dispersion is higher in groups B and C compared to A: for the Gaussian family the dispersion formula models the residual standard deviation on the log scale, so exp(1.55) = 4.7 means that the SD in B is about 4.7 times that in A. gls reports the same quantities as SD ratios relative to A (4.7 for B and 11.8 for C, with confidence intervals from intervals() on the gls fit), whereas glmmTMB additionally gives z-tests and p-values for these differences.
Take this plot of Ozone ~ Solar.R using the airquality data. Clearly there is heteroskedasticity in the relationship:
plot(Ozone ~ Solar.R, data = airquality)
We can also see this when we fit the regression model:
m1 = lm(Ozone ~ Solar.R, data = airquality)
par(mfrow = c(2, 2))
plot(m1)
We could of course consider other predictors, but let’s say we want to fit this model specifically
gls function (package nlme) with the untransformed response to make the variance dependent on Solar.R. Hint: Read in varClasses and decide how to model this.glmmTMB to model heteroskedasticity.GLM(M)s can be heteroskedastic as well, i.e. dispersion depends on some predictors. This is the same idea as the dispersion problems in GLMs (Section 9.3.2): the variance is not what the distribution assumes, but here we model it explicitly. In glmmTMB, you can make the dispersion of the negative Binomial dependent on a formula via the dispformula argument, analogous to the variance modeling with gls and glmmTMB for the linear model (Section 10.1.2).
Variance problems would show up when plotting residuals against predicted and predictors (see Section 9.3.1). In the section on recognizing overdispersion (Section 9.3.2.1), we saw some variance problems in the Salamander model. We could add a variable dispersion model via
m3 = glmmTMB(count ~ spp + mined + (1|site), family = nbinom1,
dispformula = ~ spp + mined , data = Salamanders)
summary(m3) Family: nbinom1 ( log )
Formula: count ~ spp + mined + (1 | site)
Dispersion: ~spp + mined
Data: Salamanders
AIC BIC logLik -2*log(L) df.resid
1654.4 1730.3 -810.2 1620.4 627
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site (Intercept) 0.2283 0.4778
Number of obs: 644, groups: site, 23
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.5288 0.2799 -5.462 4.70e-08 ***
sppPR -1.3304 0.3480 -3.822 0.000132 ***
sppDM 0.2695 0.2004 1.345 0.178561
sppEC-A -0.7525 0.2772 -2.714 0.006641 **
sppEC-L 0.6228 0.2109 2.952 0.003155 **
sppDES-L 0.7113 0.1976 3.600 0.000318 ***
sppDF 0.1470 0.2171 0.677 0.498259
minedno 2.1348 0.2825 7.557 4.14e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Dispersion model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.2834 0.6414 -0.442 0.6586
sppPR 0.3160 0.7501 0.421 0.6735
sppDM 0.1979 0.5712 0.346 0.7289
sppEC-A 0.3592 0.6477 0.554 0.5792
sppEC-L 1.0830 0.5215 2.077 0.0378 *
sppDES-L 0.7951 0.5370 1.481 0.1387
sppDF 0.3769 0.6109 0.617 0.5373
minedno 0.5583 0.4187 1.334 0.1823
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(DHARMa)
res = simulateResiduals(m3, plot = T)
par(mfrow = c(1, 2))
plotResiduals(res, Salamanders$spp)
plotResiduals(res, Salamanders$mined)
What can we do if, after accounting for the functional form (Section 2.5), response transformation (Section 10.1.1) and variance modeling (Section 10.1.2), the Q-Q plot of the residuals (Section 2.4.1) still shows non-normality, in particular strong outliers? Here simulated example data with strong outliers / deviations from normality:
set.seed(123)
n = 100
concentration = runif(n, -1, 1)
growth = 2 * concentration + rnorm(n, sd = 0.5) +
rbinom(n, 1, 0.05) * rnorm(n, mean = 6*concentration, sd = 6)
plot(growth ~ concentration)
Fitting the model, we see that the distribution is too wide:
fit = lm(growth ~ concentration)
par(mfrow = c(2, 2))
plot(fit)
What can we do to deal with such distributional problems and outliers?
Robust methods generally refer to methods that are robust to violation of assumptions, e.g. outliers. More specifically, standard robust regressions typically downweight data points that have a too high influence on the fit. See https://cran.r-project.org/web/views/Robust.html for a list of robust packages in R.
# This is the classic method.
library(MASS)
fit = rlm(growth ~ concentration)
summary(fit)
Call: rlm(formula = growth ~ concentration)
Residuals:
Min 1Q Median 3Q Max
-7.1986 -0.3724 0.0377 0.3391 7.0902
Coefficients:
Value Std. Error t value
(Intercept) -0.0978 0.0594 -1.6453
concentration 2.0724 0.1048 19.7721
Residual standard error: 0.534 on 98 degrees of freedom# No p-values and not sure if we can trust the confidence intervals.
# Would need to bootstrap by hand!
# This is another option that gives us p-values directly.
library(robustbase)
fit = lmrob(growth ~ concentration)
summary(fit)
Call:
lmrob(formula = growth ~ concentration)
\--> method = "MM"
Residuals:
Min 1Q Median 3Q Max
-7.2877 -0.4311 -0.0654 0.2788 7.0384
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.04448 0.05160 -0.862 0.391
concentration 2.00588 0.08731 22.974 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Robust residual standard error: 0.5549
Multiple R-squared: 0.8431, Adjusted R-squared: 0.8415
Convergence in 7 IRWLS iterations
Robustness weights:
9 observations c(27,40,47,52,56,76,80,91,100)
are outliers with |weight| = 0 ( < 0.001);
5 weights are ~= 1. The remaining 86 ones are summarized as
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.6673 0.9015 0.9703 0.9318 0.9914 0.9989
Algorithmic parameters:
tuning.chi bb tuning.psi refine.tol
1.548e+00 5.000e-01 4.685e+00 1.000e-07
rel.tol scale.tol solve.tol zero.tol
1.000e-07 1.000e-10 1.000e-07 1.000e-10
eps.outlier eps.x warn.limit.reject warn.limit.meanrw
1.000e-03 1.819e-12 5.000e-01 5.000e-01
nResample max.it best.r.s k.fast.s k.max
500 50 2 1 200
maxit.scale trace.lev mts compute.rd fast.s.large.n
200 0 1000 0 2000
psi subsampling cov
"bisquare" "nonsingular" ".vcov.avar1"
compute.outlier.stats
"SM"
seed : int(0) Quantile regressions don’t fit a line with an error spreading around it, but try to fit a quantile (e.g. the 0.5 quantile, the median) regardless of the distribution. Thus, they do not assume a particular residual distribution or constant variance, but functional form and independence still matter, and the median is not the mean.
library(qgam)
dat = data.frame(growth = growth, concentration = concentration)
fit = qgam(growth ~ concentration, data = dat, qu = 0.5) Estimating learning rate. Each dot corresponds to a loss evaluation.
qu = 0.5...............done summary(fit)
Family: elf
Link function: identity
Formula:
growth ~ concentration
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.08310 0.05709 -1.456 0.145
concentration 2.05033 0.09467 21.657 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.427 Deviance explained = 49.4%
-REML = 160.78 Scale est. = 1 n = 100lm() and ANOVA can be badly wrong, especially for the predictor the variance is related to; in our example, a non-significant ANOVA became significant once the variance was modeled.gls with varIdent or varPower, or glmmTMB with dispformula. For the Gaussian family, dispformula models the residual standard deviation on the log scale.dispformula, e.g. with the negative binomial). Check again with simulated residuals plotted against the predictors.