Generalized linear models
A screen-reader-friendly transcript of the lecture slides: all slide text, tables, figure descriptions, and speaker notes.
About these notes
Each section is one slide. Figures carry text descriptions where available. For the concepts with worked examples, also read the plain-language summary and lab for this session.Slide 1: Generalized Linear Models
Some material from Olivier Gimenez 2020
Slide 2: General Linear Model
Often written as:
But even clearer:
From this we can see assumption of constant variance
y ~ N(β0 + β1x1 + β2x2 + β3x3 + …, σ)
- y = β0 + β1x1 + β2x2 + β3x3 + … + ε
- ε ~ N(0, σ)
Slide 3
Slide from Federico Vegetti
Slide 4
Slide from Federico Vegetti
Slide 5: Regression is ONE KIND of General Linear Model
Slide 6: What does it mean to be linear?
- #1 is a linear model.
- The COEFFICIENTS need to be linear. Not the predictors!
Slide 7: All models have stochastic part for errors and deterministic “scientific” model defined by both algebra and geometry
What is the associated R code for this?
Slide 8: General Linear Model
y ~ N(β0 + β1x1 + β2x2 + β3x3 + …, σ)
Slide 9: What influences partying in college?
- partyHr
- student’s party hours per week
- drinks
- average number of drinks per week
- Is there also a difference among male and female students?
From Yuk Tung Liu’s STAT 390EF
Speaker notes for this slide
What kind of model is this? (regression) How would bad scientists interpret? Let’s interpret correctly
Slide 10: What influences partying in college?
- What is the average number of party hours for males and females who don’t drink?
- What is the average number of party hours for males and females who have an average of 5 drinks per week?
- What does the algebra and geometry (graph) of this model look like?
Group Work
Speaker notes for this slide
What is this model called in stats classes? (ANCOVA)
Slide 11
- Can change the reference level to female
- survey$gender2 <- relevel(survey$gender,"female")
Slide 12: Interaction term
- What is the interpretation?
- What is the average number of party hours for males and females who have an average of 5 drinks per week?
- What does the algebra and geometry (graph) of this model look like?
Group Work
Slide 13: More than two levels – Drinking by religion
Start with reference level “Christian”
Slide 14: More than two levels – Drinking by religion
Now with reference level “Atheist”
- Why might this model not by ideal?
- Is there a problem with releveling to assess all pairwise interactions?
Speaker notes for this slide
Remember to think about sample size
Slide 15: Sometimes we can’t be Normal
(No text on this slide.)
Slide 16: Titanic Data – Try a General Linear Model
- Survival (1/0) for passengers with known age, sex, and class
- We could try fitting the linear model for class (crew, first, second, third)
- titanic.lm <- lm(survived ~ class, data = titanic)
- But the residuals, of course, are not remotely normally distributed.
- And model is grossly misspecified because “survival” can only be 0 or 1
The actual probability model under the hood
Slide 17: Titanic Data – General Linear Model
Let’s interpret this output
Slide 18: Titanic Data – Generalized Linear Model
Slide 19: From probability to odds to log of odds
- Imagine the probability of success of some event is .8. Then the probability of failure is 1 – .8 = .2. The odds of success is the ratio of the probability of success over the probability of failure. 0.8/0.2 = 4. That is to say that the odds of success are 4 to 1. If the probability of success is .5, i.e., 50-50 percent chance, then the odds of success is 1 to 1.
- The transformation from probability to odds is a monotonic transformation, meaning the odds increase as the probability increases or vice versa. Probability ranges from 0 and 1. Odds range from 0 and positive infinity. Log odds from negative Inf to Inf
- MUCH easier to interpret probability
Slide 20: Titanic Data – Generalized Linear Model
- The logit link function. If:
- Then what is p? Find it.
Group Work
Slide 21: Titanic Data – Generalized Linear Model
- The logit link function. If:
- Then what is p?
- Note: We have transformed p to the whole real line, and can now fit a linear combination of our predictors, which in this case are discrete factors
Slide 22: Titanic Data – Let’s look at the data
Slide 23: Binomial GLM (i.e. Logistic Regression)
How do we interpret this?
Slide 24: Interpreting the effect size with odds ratio
Group Work
Show that this is also true for any continuous variables, x, by using the odds ratio of x + 1 to x
Slide 25: Questions
- Based on the model output, what is the probability that a member of the crew, first, second, third class survived? How much higher is the probability of survival for first class relative to crew?
- How should these probabilities relate to the data in this model?
- What are the p-values telling you?
Group Work
Slide 26
Slide 27
So this model was trivial for estimation but was still needed for inference!
Slide 28
Slide 29: Did males have higher survival than females?
(No text on this slide.)
Slide 30
Slide 31
Slide 32: Did survival vary by class and sex?
(No text on this slide.)
Slide 33
Slide 34: How do we interpret this?
- What are the log(odds) for first class males and females?
- What are the corresponding survival probabilities?
- What are the log(odds) for male and female crew?
- What are the corresponding survival probabilities?
Group Work
Slide 35
Slide 36
Slide 37: Instead of looking at individual passengers, can’t we look at proportion by group?
(No text on this slide.)
Slide 38
Not 1 anymore
Bernoulli becomes a Binomial
- Response ~ Distribution(Mean Response)
- Yi ~ Binomial(Ni, pi)
- logit(pi) = a + b xi
- pi = logit−1(a + b xi) = exp(a + b xi) / (1 + exp(a + b xi))
Slide 39
Check to make sure results are identical. I’ve made it not 100% easy to check
Slide 40: Models for count data
- Poisson
- Quasipoisson
- Negative Binomial
- Binomial
Unbounded
Bounded
Slide 41
Oncorhynchus nerka (sockeye salmon)
Does environmental DNA (eDNA) contain enough information to estimate salmon escapement and smolt abundance?
=
One person, one sample /day
Several people, all day
Speaker notes for this slide
A minor revolution in ecology is using eDNA from water to sample aquatic species. I have started a project with He Daming of the Asian International Rivers Center at Yunnan University to use this technology to map fish populations.
Slide 42
Levi. T., Allen, J. M., Bell, D., Joyce, J. Russel, J. R., Tallmon, D. A., Vulstek, S. C., Yahan, Y., Yu, D. W. 2019. Environmental DNA for the enumeration and management of Pacific salmon. Molecular Ecology Resources
Slide 43
- Validation data from Auke Creek Station in Juneau, Alaska
- Daily counts of
- returning Coho and Sockeye adults
- out-migrating Sockeye smolts
Slide 44
Sockeye returning adults
Slide 45
Sockeye returning adults
Slide 46
Sockeye returning adults
Slide 47
Sockeye returning adults
Water flux (cubic_ft/sec) X DNA concentration (ng/ul) = ng/sec X cubic_ft/ul
Slide 48
Sockeye returning adults
eDNA rate: ng/sec
Slide 49
Sockeye out-migrating smolts
Slide 50: Can we predict salmon counts using eDNA
Let’s start with Poisson model using flow-corrected eDNA concentration as a predictor.
Slide 51: What is happening here?
- Data are Poisson distributed with varying rate parameter
- This parameter is the mean (and variance)
- We model the mean as a linear function of predictors, but we need a link function
- Why do we use log link here?
Slide 52: Model output
- glm(formula = Count ~ Qcorr_qPCR, family = "poisson", data = sockeye.adult[sockeye.adult$Year ==
- "2015", ])
- Deviance Residuals:
- Min 1Q Median 3Q Max
- -13.5661 -11.7512 -3.2418 3.5480 19.4317
- Coefficients:
- Estimate Std. Error z value Pr(>|z|)
- (Intercept) 4.29048 0.03329 128.90 <2e-16 ***
- Qcorr_qPCR 3.82647 0.10393 36.82 <2e-16 ***
- ---
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- (Dispersion parameter for poisson family taken to be 1)
- Null deviance: 2666.3 on 15 degrees of freedom
- Residual deviance: 1470.1 on 14 degrees of freedom
- AIC: 1555.9
- Number of Fisher Scoring iterations: 5
Slide 53: Model output
- glm(formula = Count ~ Qcorr_qPCR, family = "poisson", data = sockeye.adult[sockeye.adult$Year ==
- "2016", ])
- Deviance Residuals:
- Min 1Q Median 3Q Max
- -14.2212 -6.4894 -6.4878 -0.8424 36.7864
- Coefficients:
- Estimate Std. Error z value Pr(>|z|)
- (Intercept) 3.04662 0.03573 85.27 <2e-16 ***
- Qcorr_qPCR 4.33630 0.06636 65.34 <2e-16 ***
- ---
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- (Dispersion parameter for poisson family taken to be 1)
- Null deviance: 6890.4 on 38 degrees of freedom
- Residual deviance: 2918.3 on 37 degrees of freedom
- AIC: 3025.1
- Number of Fisher Scoring iterations: 6
Slide 54: How does this look? Concerns? What now? Are we done?
(No text on this slide.)
Slide 55: How does this look? Concerns? What now? Are we done?
- Poisson assumes variance = mean
- Overdispersion evident if residual deviance > degrees of freedom
- We clearly have overdispersion here because deviance is 100 times larger than degrees of freedom
- Let’s deal with the overdispersion later.
- Is the model otherwise ok?
Slide 56: 1. X vs. Y
- Why isn’t this a line if we fit
- glm(Count ~ Qcorr_qPCR , family = "poisson")
- CI is far too narrow because of overdispersion (variance > mean)
- Nothing else obvious
Slide 57: 2. Good models make good predictions
- Ah, look at that. We are consistently overpredicting the low/no counts!
- Why? Let’s think more carefully
- Is this model mispecified somehow?
Slide 58: What is being fit here exactly?
- What would eDNA rate have to be to match zero counts? Do you see it?
- Yeah, that’s a problem! How do we fix it?
Slide 59: eDNA rate should be transformed from positive to the whole real line! How?
- Log(eDNA) > 0 when eDNA > 1
- Log(eDNA) = 0 when eDNA = 1
- Log(eDNA) < 0 when eDNA = <1
- Log(eDNA) = -Inf when eDNA = 0
- Log transforming a positive variable stretches it to the real line.
- But what functional form are we really fitting here?
Slide 60: log(Y) = a + b*log(X)
(No text on this slide.)
Slide 61: New model output
- glm(formula = Count ~ log(Qcorr_qPCR), family = "poisson", data = sockeye.adult[sockeye.adult$Year ==
- "2015", ])
- Deviance Residuals:
- Min 1Q Median 3Q Max
- -14.5709 -6.1344 -1.5773 3.4050 17.3989
- Coefficients:
- Estimate Std. Error z value Pr(>|z|)
- (Intercept) 6.42652 0.04458 144.16 <2e-16 ***
- log(Qcorr_qPCR) 0.62802 0.02284 27.49 <2e-16 ***
- ---
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- (Dispersion parameter for poisson family taken to be 1)
- Null deviance: 2666.3 on 15 degrees of freedom
- Residual deviance: 1088.8 on 14 degrees of freedom
- AIC: 1174.6
- Number of Fisher Scoring iterations: 6
Slide 62: New model output
- glm(formula = Count ~ log(Qcorr_qPCR), family = "poisson", data = sockeye.adult[sockeye.adult$Year ==
- "2016", ])
- Deviance Residuals:
- Min 1Q Median 3Q Max
- -12.2184 -2.6595 -0.9382 -0.4217 21.8418
- Coefficients:
- Estimate Std. Error z value Pr(>|z|)
- (Intercept) 6.25952 0.02922 214.24 <2e-16 ***
- log(Qcorr_qPCR) 0.78915 0.01716 45.99 <2e-16 ***
- ---
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- (Dispersion parameter for poisson family taken to be 1)
- Null deviance: 6890.4 on 38 degrees of freedom
- Residual deviance: 1174.6 on 37 degrees of freedom
- AIC: 1281.4
- Number of Fisher Scoring iterations: 5
Slide 63: Good models make good predictions
- Much better at predicting zeros now
- But still very narrow confidence intervals that don’t make sense
- Clear overdispersion because deviance >> degrees of freedom
- Now let’s fix the overdispersion problem
Slide 64: Quasipoisson GLM
- The Normal distribution has the mean and variance as separate parameters, but this is unusual
- Poisson distribution has variance equal to the mean
- Quasipoisson allows the variance to be proportional to the mean
- This ”dispersion parameter” is the constant of proportionality that defines how much larger or smaller the variance is than the mean
Slide 65: Quasipoisson GLM
- glm(formula = Count ~ log(Qcorr_qPCR), family = "quasipoisson",
- data = sockeye.adult[sockeye.adult$Year == "2015", ])
- Deviance Residuals:
- Min 1Q Median 3Q Max
- -14.5709 -6.1344 -1.5773 3.4050 17.3989
- Coefficients:
- Estimate Std. Error t value Pr(>|t|)
- (Intercept) 6.4265 0.3940 16.310 1.67e-10 ***
- log(Qcorr_qPCR) 0.6280 0.2019 3.111 0.00767 **
- ---
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- (Dispersion parameter for quasipoisson family taken to be 78.1201)
- Null deviance: 2666.3 on 15 degrees of freedom
- Residual deviance: 1088.8 on 14 degrees of freedom
- AIC: NA
- Number of Fisher Scoring iterations: 6
Slide 66: Quasipoisson figures
Slide 67: Quasipoisson figures
Slide 68: This is great, but….
- No AIC.
- We’ll talk about AIC later, but if you need AIC then we need another distribution
- Which one?
Slide 69: Negative binomial
Can’t use glm(). Need to use glm.nb() in the MASS package
- sockeye.adult.2015.reg.NB <- glm.nb(Count ~ log(Qcorr_qPCR) ,
- data = sockeye.adult[sockeye.adult$Year == "2015" , ])
- sockeye.adult.2016.reg.NB <- glm.nb(Count ~ log(Qcorr_qPCR) ,
- data = sockeye.adult[sockeye.adult$Year == "2016" , ])
- Default log link
- Still use logarithm here
Slide 70: Negative binomial
We now have AIC
- glm.nb(formula = Count ~ log(Qcorr_qPCR), data = sockeye.adult[sockeye.adult$Year ==
- "2015", ], init.theta = 0.7864342616, link = log)
- Deviance Residuals:
- Min 1Q Median 3Q Max
- -2.7537 -0.6590 -0.4000 0.2669 1.3384
- Coefficients:
- Estimate Std. Error z value Pr(>|z|)
- (Intercept) 6.8527 0.5032 13.619 < 2e-16 ***
- log(Qcorr_qPCR) 0.8128 0.1518 5.355 8.55e-08 ***
- ---
- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- (Dispersion parameter for Negative Binomial(0.7864) family taken to be 1)
- Null deviance: 40.432 on 15 degrees of freedom
- Residual deviance: 17.362 on 14 degrees of freedom
- AIC: 169.4
Slide 71: Negative binomial
Slide 72: Negative binomial
Slide 73: Gamma
- Another common GLM.
- Useful for positive and continuous data
- What is mean variance relationship?
Slide 74: Gamma
- Another common GLM.
- Useful for positive and continuous data
- What is mean variance relationship?
- Parameters, shape-scale: k > 0 (shape), θ > 0 (scale)
- Parameters, shape-rate: α > 0 (shape), β > 0 (rate)
- Support: x ∈ (0, ∞)
- PDF, shape-scale: (1 / (Γ(k) θk)) · xk−1 · exp(−x / θ)
- PDF, shape-rate: (βα / Γ(α)) · xα−1 · exp(−β x)
- Mean: E[X] = k θ (shape-scale) = α / β (shape-rate)
- Variance: Var(X) = k θ2 (shape-scale) = α / β2 (shape-rate)
Slide 75: Interaction Terms
“Context dependence” – The effect of one predictor depends on the value of another predictor
Slide 76: Interpretation and Symmetry of Interactions
- How does the effect of x1 on y change depending on the value of x2
- How does the effect of x2 on y change depending on the value of x1
Let’s focus on deterministic part
Factor out x1
Factor out x2
Slide 77: Interpretation and Symmetry of Interactions
Vs.
Slide 78: Interpretation and Symmetry of Interactions
Vs.
Slide 79: Dependence on scaling and transformation I
- Simulate data with no interaction
- Y = 200+3*x1 - 4*x2
lm(y~x1*x2)
Slide 80: Dependence on scaling and transformation II
Log transform dependent variable
lm(log(y)~x1*x2)
What is happening and why?!
Slide 81: Additivity not preserved across nonlinear transformations!
If log(y) is a function of additive x1 and x2, will y generate an interaction?
Equations for reference, work through on paper
Slide 82: Confusion in action with major policy implications!
Q: Deer are expanding, driving decline of caribou by supporting wolves. Is this due to logging or climate or both?
(1)
(2)
(3)
Slide 83: What was the point made by Barnas et al.?
Slide 84: What was the point made by Barnas et al.?
Weird. Why? Who is right?
Slide 85: Models, especially with interactions, require algebra and especially geometry. Interpret the predictions!
- As they SHOULD, the two approaches produce IDENTICAL predictions.
- -What are the predictions?
- -Was there a transformation?
- What lesson should we take from this?
- You cannot just interpret main effects. Let’s think about why they were different