18.4: Generalized Linear Squares
- Page ID
- 45261
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)Introduction
Draft
With access to powerful computers and better algorithms, we can move past the classical ANOVA and ordinary least squares approaches to linear models. We have discussed general linear models, but here we introduce generalized linear models, GLM. What follows is just a brief foray; for more — and better! discussion, see Zuur et al (2009).
Model variances
Data from Corn and Hiesey (1973) ohia.RData
> head(ohia) Site Height Width 1 M-1 12.5567 19.1264 2 M-1 13.2019 13.1547 3 M-1 8.0699 16.0320 4 M-1 6.0952 22.8586 5 M-1 11.3879 11.0105 6 M-1 12.2242 21.8102
# ignore the variance issue
> AnovaModel.1 <- aov(Height ~ Site, data = ohia); summary(AnovaModel.1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 2 4070 2034.8 22.63 0.000000131 ***
Residuals 47 4227 89.9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Alternatively, use gls(). Default fits by restricted maximum likelihood, REML. That is, it’s the likelihood of linear combinations of the original data.
>model.aov.1 <- gls(Height ~ Site, data = ohia)
Generalized least squares fit by REML
Model: Height ~ Site
Data: ohia
AIC BIC logLik
361.1312 368.5318 -176.5656
Coefficients:
Value Std.Error t-value p-value
(Intercept) 15.313745 2.120550 7.221591 0.0000
Site[T.M-2] 19.261000 2.998911 6.422666 0.0000
Site[T.M-3] 2.924215 3.672900 0.796160 0.4299
Correlation:
(Intr) S[T.M-2
Site[T.M-2] -0.707
Site[T.M-3] -0.577 0.408
Standardized residuals:
Min Q1 Med Q3 Max
-1.9832938 -0.5020880 -0.1850871 0.5017636 3.0850635
Residual standard error: 9.483388
Degrees of freedom: 50 total; 47 residual
Code for the plot in Figure \(\PageIndex{1}\):
par(mfrow = c(1, 2)) plot(residuals(model.aov.1) ~ Site, pch=19, cex=1.5,col=”blue”, data = ohia) plot(residuals(model.aov.1) ~ fitted(model.aov.1), pch=19, cex=1.5, col=”blue”, ylab=””, data=ohia) mtext(“ANOVA Ohi`a Height, 3 Maui sites “, side = 3, line = -3, outer = TRUE)
# test equal variances, Height
> leveneTest(Height ~ Site, data=ohia, center="median")
Levene's Test for Homogeneity of Variance (center = "median")
Df F value Pr(>F)
group 2 2.1663 0.1259
47
We would conclude no significant departures from equal variances.
> bartlett.test(Height ~ Site, data=ohia) Bartlett test of homogeneity of variances data: Height by Site Bartlett's K-squared = 10.373, df = 2, p-value = 0.005592
Bartlett’s test is sensitive to deviations from normality.
Include variances as part of model
> model.aov.3 <- gls(Height ~ Site, data = ohia, weights = varIdent(form = ~1|Site)); summary(model.aov.3)
Generalized least squares fit by REML
Model: Height ~ Site
Data: ohia
AIC BIC logLik
354.421 365.5219 -171.2105
varIdent permits variances for each group to vary. Results from R continue below.
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Site
Parameter estimates:
M-1 M-2 M-3
1.0000000 1.0396880 0.3471771
We see here that comparisons were carried out versus the M-1 site.
Coefficients:
Value Std.Error t-value p-value
(Intercept) 15.313745 2.280931 6.713812 0.0000
Site[T.M-2] 19.261000 3.290358 5.853770 0.0000
Site[T.M-3] 2.924215 2.541027 1.150800 0.2556
Marginal differences between M-1 and M-2 for height were significantly different, but not between the M-1 and M-3 site.
Correlation:
(Intr) S[T.M-2
Site[T.M-2] -0.693
Site[T.M-3] -0.898 0.622
Standardized residuals:
Min Q1 Med Q3 Max
-1.7734556 -0.6909962 -0.2108834 0.5801370 2.7586550
Residual standard error: 10.20064
Degrees of freedom: 50 total; 47 residual
# Test the models
> anova(model.aov.1, model.aov.3)
Model df AIC BIC logLik Test L.Ratio p-value
model.aov.1 1 4 361.1312 368.5318 -176.5656
model.aov.3 2 6 354.4210 365.5219 -171.2105 1 vs 2 10.7102 0.0047
Although additional degrees of freedom are required, note that this model (model.aov.3) has higher (better!) log likelihood (-171.21) than model.aov.1, the gls model lacking a fit for different variances (-176.57). Introduce a test of the hypothesis that the two models are equal by comparing the log (natural) likelihoods, the log likelihood ratio test, LRT.
\[LRT = -2 \cdot \ln \left(\frac{LL \ model_{aov.1}}{LL \ model_{aov.3}}\right) = -2 \cdot \ln [(LL \ model_{aov.1}) - (LL \ model_{aov.3})] \nonumber\]
The LRT follows a chi-square distribution (per Wilk’s theorem). If there was no advantage to fitting for unequal variances, then the model fit would not be improved and p-value of the LRT would not be less than 5%.
Conclusion
You can see why this approach, modeling versus separate test of assumptions would be the preferred way to go. We get a better fitting model, cf discussion in
Another example, same data set.
# ignore variances, Width
model.aov.2 <- gls(Width ~ Site, data = ohia); summary(model.aov.2)
Figure 2.
par(mfrow = c(1, 2)) plot(residuals(model.aov.2) ~ Site, pch=19, cex=1.5,col=”blue”, data = ohia) plot(residuals(model.aov.2) ~ fitted(model.aov.2), pch=19, cex=1.5, col=”blue”, ylab=””, data=ohia) mtext(“ANOVA Ohi`a Width 3 Maui sites “, side = 3, line = -3, outer = TRUE)
# test equal variances, Width
Tapply(Width ~ Site, var, na.action=na.omit, data=ohia) # variances by group leveneTest(Width ~ Site, data=ohia, center=”median”) Tapply(Width ~ Site, var, na.action=na.omit, data=ohia) # variances by group bartlett.test(Width ~ Site, data=ohia)
# model the variances, Height
library(nlme) model.aov.3 <- gls(Height ~ Site, data = ohia, weights = varIdent(form = ~1|Site)); summary(model.aov.3) par(mfrow = c(1, 2)) plot(residuals(model.aov.3) ~ Site, pch=19, cex=1.5,col=”red”, data = ohia) plot(residuals(model.aov.3) ~ fitted(model.aov.3), pch=19, cex=1.5, col=”red”, ylab=””, data=ohia) mtext(“GLS Ohi`a Height 3 Maui sites “, side = 3, line = -3, outer = TRUE)
# Test the models
anova(model.aov.1, model.aov.3)
# model the variances, Width
model.aov.4 <- gls(Width ~ Site, data = ohia, weights = varIdent(form = ~1|Site)); summary(model.aov.4) par(mfrow = c(1, 2)) plot(residuals(model.aov.4) ~ Site, pch=19, cex=1.5,col=”red”, data = ohia) plot(residuals(model.aov.4) ~ fitted(model.aov.4), pch=19, cex=1.5, col=”red”, ylab=””, data=ohia) mtext(“GLS Ohi`a Width 3 Maui sites “, side = 3, line = -3, outer = TRUE)
# Test the models
anova(model.aov.2, model.aov.4)
Model correlated residuals
[pending]
Questions
[pending]