Genearlized Linear Models

A Generalized Outline

1. Why use GLMs? An Intro to Entropy
   
   
2. Logistic Regression Verus Linear Regression
   
   
3. Generalized Linear Models
   
   
4. Poisson Regression (Poisson Error, Long Link)

What is Maximum Entropy?

Maximum Entropy Principle

-McElreath 2017

Why are we thinking about MaxEnt?

• MaxEnt distributions have the widest spread - conservative
  
  
• Nature tends to favor maximum entropy distributions
  • It’s just natural probability
    
    
• The foundation of Generalized Linear ModelDistributions
  
  
• Leads to useful distributions once we impose constraints

Information Entropy

\[H(p) = - \sum p_i log \, p_i\]

• Measure of uncertainty
  
  
• If more events possible, it increases
  
  
• Nature finds the distribution with the largest entropy, given constraints of distribution

Maximum Entropy and Coin Flips

• Let’s say you are flipping a fair (p=0.5) coin twice
  
  
• What is the maximum entropy distribution of # Heads?
  
  
• Possible Outcomes: TT, HT, TH, HH
  • which leads to 0, 1, 2 heads
    
    
• Constraint is, with p=0.5, the average outcome is 1 heads

The Binomial Possibilities

Let’s compareother distributions meeting constraint using Entropy

Remember, we must average 1 Heads, so,
sum(distribution * 0,1,1,2) = 1

\[H = - \sum{p_i log p_i}\]

What about other p’s and draws?

Assume 2 draws, p=0.7, make 1000 simulated distributions

OK, what about the Gaussian?

Maximum Entropy Distributions

How Distributions are Coupled

A Generalized Outline

1. Why use GLMs? An Intro to Entropy
   
   
2. Logistic Regression Verus Linear Regression
   
   
3. Generalized Linear Models
   
   
4. Poisson Regression (Poisson Error, Long Link)

Infection by Cryptosporidium

Cryptosporidum Infection Rates

This is not linear or gaussian

Why?

The General Linear Model

\[\Large \boldsymbol{Y_i} = \boldsymbol{\beta X_i} + \boldsymbol{\epsilon} \]





\[\Large \epsilon \sim \mathcal{N}(0,\sigma^{2})\]

The General Linear Model

Likelihood:
\[\Large Y_i \sim \mathcal{N}(\hat{Y_i},\sigma^{2})\]



Data Generating Process:
\[\Large \boldsymbol{\hat{Y}_{i}} = \boldsymbol{\beta X_i} \]

The General(ized) Linear Model

Likelihood:
\[\Large Y_i \sim \mathcal{N}(\hat{Y_i},\sigma^{2})\]

Data Generating Process:
- Transformation (Identity Link):
\[\Large \hat{Y}_{i} = \eta_{i} \]
- Linear Equation:
\[\Large \boldsymbol{\eta_{i}} = \boldsymbol{\beta X_i} \]

A Generalized Linear Model with an Exponential Curve

Likelihood:
\[\Large Y_i \sim \mathcal{N}(\hat{Y_i},\sigma^{2})\]

Data Generating Process:
- Transformation (Log Link):
\[\Large Log(\hat{Y}_{i}) = \eta_{i} \]

- Linear Equation:
\[\Large \boldsymbol{\eta_{i}} = \boldsymbol{\beta X_i} \]

Log Link

Isn’t this just a transformation?

Likelihood:
\[\Large Y_i \sim \mathcal{N}(\hat{Y_i},\sigma^{2})\] Error is Normal

Data Generating Process:
- Transformation (Log Link):
\[\Large Log(\hat{Y}_{i}) = \eta_{i} \]

- Linear Equation:
\[\Large \boldsymbol{\eta_{i}} = \boldsymbol{\beta X_i} \]

But This is Not Normal

Binomial Distribution

\[ Y_i \sim B(prob, size) \]

• Discrete Distribution
• prob = probability of something happening (% Infected)
• size = # of discrete trials
• Used for frequency or probability data
• We estimate coefficients that influence prob

So, Y is a Logistic Curve

Generalized Linear Model with a Logit Link

Likelihood: \[\Large Y_i \sim \mathcal{B}(\hat{Y_i}, size)\]
Data Generating Process:
Logit Link Function \[\Large Logit(\hat{Y_i}) = \eta_{i}\]

Linear Function
\[\Large \boldsymbol{\eta_{i}} = \boldsymbol{\beta X_i} \]

Logitistic Regression

Generalized Linear Model with Logit Link

crypto_glm <- glm(Y/N ~ Dose,
                  weight=N,
                  family=binomial(link="logit"),
                  data=crypto)

OR, with Success and Failures

crypto_glm <- glm(cbind(Y, Y-N) ~ Dose,
                  family=binomial(link="logit"),
                  data=crypto)

Outputs

And logit coefficients

The Odds

The Meaning of a Logit Coefficient

But how do we assess assumptions?

• Should still be no fitted v. residual relationship
  
  
• But QQ plots lose meaning
  • Not a normal distribution
  • Mean scales with variance
    
    
• Also many types of residuals
  • Deviance, Pearson, raw, etc.

Randomized quantile residuals

• If model fits well, quantiles of residuals should be uniformly distributed
  
  
• I.E., for any point, if we had its distribution, there should be no bias in its quantile
  
  
• We do this via simulation
  
  
• Works for many models, and naturally via Bayesian simuation

Randomized quantile residuals: Steps

1. Get ~1000 (or more) simulations of model coefficients
   
   
2. For each response (y) value, create an empirical distribution from the simuations
   
   
3. For each response, determine it’s quantile from that empirical distribution
   
   
4. The quantiles of all y values should be uniformly distributed
   • QQ plot of a uniform distribution!
     

Randomized quantile residuals: Visualize

Randomized quantile residuals: Visualize

Randomized quantile residuals: Visualize

Randomized quantile residuals: Visualize

Randomized quantile residuals: Visualize

Quantile Residuals

library(DHARMa)

set.seed(2010)
#Get quantile residuals from simulations
simulationOutput_crypto <- simulateResiduals(crypto_glm, 
                                      n = 1000)
#plot!
plot(simulationOutput_crypto)

Quantile Residuals

A Generalized Outline

1. Why use GLMs? An Intro to Entropy
   
   
2. Logistic Regression Verus Linear Regression
   
   
3. Generalized Linear Models
   
   
4. Poisson Regression (Poisson Error, Long Link)

General Linear Models are a Special Case

Likelihood:
\[\Large Y_i \sim \mathcal{N}(\hat{Y_i},\sigma^{2})\]

Data Generating Process:
- Transformation (Identity Link):
\[\Large \hat{Y}_{i} = \eta_{i} \]
- Linear Equation:
\[\Large \boldsymbol{\eta_{i}} = \boldsymbol{\beta X_i} \] But what if We don’t want a Normal Distribution?

The Generalized Linear Model

Likelihood:
\[\boldsymbol{Y_i} = E(\boldsymbol{\hat{Y_i}}, \theta)\]
E is any distribution from the Exponential Family
\(\theta\) is an error parameter, and can be a function of Y

Data Generating Process:
- Link Function \[\boldsymbol{f(\hat{Y_i})} = \boldsymbol{\eta_i}\]

- Linear Predictor \[\boldsymbol{\eta_i} = \boldsymbol{\beta X}\]

Generalized Linear Models: Link Functions

Basic Premise:

1. We have a linear predictor, \(\eta_i = a+Bx_i\)
   
   

2. That predictor is linked to the fitted value of \(Y_i\), \(\hat{Y_i}\)
   
   

3. We call this a link function, such that \(g(\hat{Y_i}) = \eta_i\)
   

   • For example, for a linear function, \(\mu_i = \eta_i\)
     

   • For an exponential function, \(log(\mu_i) = \eta_i\)

Some Common Links

• Identity: \(\hat{Y_i} = \eta_i\) - e.g. \(\mu = a + bx\)

• Log: \(log(\hat{Y_i}) = \eta_i\) - e.g. \(\mu = e^{a + bx}\)

• Logit: \(logit(\hat{Y_i}) = \eta_i\) - e.g. \(\hat{Y_i} = \frac{e^{a + bx}}{1+e^{a + bx}}\)

• Inverse: \(\frac{1}{\hat{Y_i}} = \eta_i\) - e.g. \(\hat{Y_i} = (a + bx)^{-1}\)

Generalized Linear Models: Error

Basic Premise:

1. The error distribution is from the exponential family

   • e.g., Normal, Poisson, Binomial, and more.

   
   
   

2. For these distributions, the variance is a funciton of the fitted value on the curve: \(var(Y_i) = \theta V(\hat{Y_i})\)

   • For a normal distribution, \(var(Y_i) = \theta*1\) as \(V(\hat{Y_i})=1\)
     

   • For a poisson distribution, \(var(Y_i) = 1*\mu_i\) as \(V(\hat{Y_i})=\hat{Y_i}\)

Distributions, Canonical Links, and Dispersion

The Generalized Linear Model

Likelihood:
\[\boldsymbol{Y_i} = E(\boldsymbol{\hat{Y_i}}, \theta)\]
E is any distribution from the Exponential Family
\(\theta\) is an error parameter, and can be a function of Y

Data Generating Process:
- Link Function \[\boldsymbol{f(\hat{Y_i})} = \boldsymbol{\eta_i}\]

- Linear Predictor \[\boldsymbol{\eta_i} = \boldsymbol{\beta X}\]

A Generalized Outline

1. Why use GLMs? An Intro to Entropy
   
   
2. Logistic Regression Verus Linear Regression
   
   
3. Generalized Linear Models
   
   
4. Poisson Regression (Poisson Error, Long Link)

Poisson Regression with a Log Link

Likelihood:
\[\boldsymbol{Y_i} \sim \mathcal{P}(\lambda = \boldsymbol{\hat{Y_i}})\]

Data Generating Process: \[log(\boldsymbol{\hat{Y_i}}) = \boldsymbol{\eta_i}\]
\[\boldsymbol{\eta_i} = \boldsymbol{\beta X_i}\]

What is the relationship between kelp holdfast size and number of fronds?

What About Kelp Holdfasts?

How ’bout dem residuals?

kelp_lm <- lm(FRONDS ~ HLD_DIAM, data=kelp)

What is our data and error generating process?

What is our data and error generating process?

• Data generating process should be exponential - No values less than 1
  
  
• Error generating process should be Poisson - Count data

What is our data and error generating process?

kelp_glm <- glm(FRONDS ~ HLD_DIAM, data=kelp,
                family=poisson(link="log"))

Kelp GLM Results

LR Test

Coefficients:

Kelp GLM Results

kelp_plot +
  stat_smooth(method="glm", 
              method.args=list(family=poisson(link="log")))

Kelp GLM Results

Kelp GLM Quantile Residuals

library(DHARMa)
set.seed(2017)
simulationOutput <- simulateResiduals(kelp_glm, 
                                      n = 1000)
plot(simulationOutput)

Ruh roh! Overdispersion

• Sometimes, your variance changes faster than predicted by your error distribution
  
  
• This is called overdispersion
  
  
• We will deal with it more formally next week, but…
  
  
• Sometimes, a different error structure will do.

The Negative Binomial

• Related to the binomial (coin flips!)
  
  
• Number of failures before size successes are seen given p as the probability of success
  
  
• \(Y_i \sim NB(size, p)\)
  
  
• Variance = \(\hat{Y_i}^2 + \kappa\hat{Y_i}^2\)|$
  
  
• Increases with the square, not linearly

The Negative Binomial

Negative Binomials and Log Links…

• We can also write NBs using means and a dispersion parameter.
  
  
• If the mean is \(\mu_i\), then, the variance is: \[var(Y_i) = \mu_i + \mu_i^2/\theta\]
  
• We call \(\theta\) the dispersion parameters

A Negative Binomial GLM

Likelihood:
\[\boldsymbol{Y_i} \sim NB(\boldsymbol{\hat{Y_i}}, \boldsymbol{\theta})\]

Data Generating Process: \[log(\boldsymbol{\hat{Y_i}}) = \boldsymbol{\eta_i}\]
\[\boldsymbol{\eta_i} = \boldsymbol{\beta X_i}\]

A Negative Binomial GLM

library(MASS)
kelp_glm_nb <- glm.nb(FRONDS ~ HLD_DIAM, data=kelp)
Anova(kelp_glm_nb)

Analysis of Deviance Table (Type II tests)

Response: FRONDS
         LR Chisq Df Pr(>Chisq)    
HLD_DIAM   51.145  1  8.578e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kelp NB GLM Results

kelp_plot +
  stat_smooth(method="glm.nb")

Kelp NB GLM Checks

res <- simulateResiduals(kelp_glm_nb)

plot(res)

You Try: Wolf Inbreeding and Litter Size