Binomial Logistic Regression

Y ~ Binomial(m, π)

\(P(Y=y)={m\choose y} \pi^y (1-\pi)^{m-y}\) v.s. Bernoulli => \(\pi\) remains the same in binomial, while in Bernoulli case it is independent for each \(X_i\)

\(E(Y)=m\pi\), \(var(Y)=m\pi(1-\pi)\)

Consider modelling \(Y.m\) be the proportion of success out of \(m\) independent Bernoulli trials \(E(Y/m)=\pi, var(Y/m)=\pi(1-\pi)/m\)

library(Sleuth3)
library(ggplot2)
krunnit = case2101

Extinct = krunnit$Extinct # number of "success"
AtRisk = krunnit$AtRisk
Area = krunnit$Area

pis = Extinct/AtRisk
NExtinct = AtRisk - Extinct # number of "failure"
logitpi = log(pis/(1-(pis)))
logarea = log(Area)
case2101

Let \(\pi_i\) be the probability of extinction for each island, assume that this is the same for each species for bird on a particular island

Assume species survival is independent.
Then \(Y_i\sim Binomial(m_i,\pi_i)\)

Unlike binary logistic model for Bernoulli distribution, we can estimate \(\pi_i\) from the data.

Observed response proportion \(\bar{\pi_i} = y_i/m_i\)

Observed or empirical logits: (S-"saturated)

ggplot(krunnit, aes(x=Area, y=logitpi)) + geom_point()

​

ggplot(krunnit, aes(x=log(Area), y=logitpi)) + geom_point()

​

Then the proposed model is based on plot

# cbind(Extinct, NExtinct ~ log(Area) is the way to represent it in R
fitbl <- glm(cbind(Extinct, NExtinct)~log(Area), family=binomial, data=krunnit)
summary(fitbl)

Call:
glm(formula = cbind(Extinct, NExtinct) ~ log(Area), family = binomial, 
    data = krunnit)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.71726  -0.67722   0.09726   0.48365   1.49545

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.19620    0.11845 -10.099  < 2e-16 ***
log(Area)   -0.29710    0.05485  -5.416 6.08e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 45.338  on 17  degrees of freedom
Residual deviance: 12.062  on 16  degrees of freedom
AIC: 75.394

Number of Fisher Scoring iterations: 4

anova(fitbl, test="Chisq")

vcov(fitbl) # variance co-variance matrix

# 95% CI for beta
CL = cbind(bhat=coef(fitbl), confint.default(fitbl))
CL

Model summary

# observations: 18
# coefficients: 2
fitted model \(logit(\hat\pi) = -1.196-0.297\log(Area)\)

Wald procedures still work the same

\(H_0:\beta_1=0\)
\(z=\frac{\hat\beta_1}{se(\hat\beta_1)}\sim N(0,1)\)
CI: \(\hat\beta_1\pm t_{\alpha,1}se(\hat\beta_1)\)

Interpretation of beta1

Model:

CHanging \(x\) by a factor of \(h\) will change the odds by a multiplicative factor of \(h^{\beta_1}\)

Example: Halving island area changes odds by a factor of \(0.5^{-0.2971} = 1.23\)

Therefore, the odds of extinction on a smaller island are 123% of the odds of extinction on an island double its size.

Halving of area is associated with an increase in the odds of extinction by an estimated 23%. An approximate 95% confidence interval for the percentage change in odds is 14%-32%

# CI of beta's for halfing area
0.5^(CL)

phats = predict.glm(fitbl, type="response")
options(digits=4)
rbind(Extinct, NExtinct, pis, phats)

Question: estimate the probability of extinction for a species on the Ulkokrunni island (Area \(= 185.5 km^2\))

\(logit(\hat\pi_{M,1}=)-1.196 - 0.297\log(185.5) = -2.75\)

\(\hat\pi_{M,1} = 0.06\)

logit_result = -1.196 - 0.297 * log(185.5)
pi_m1 = exp(logit_result) / (1+ exp(logit_result))
print(logit_result)
print(pi_m1)

[1] -2.747
[1] 0.06024

Diagnostics

Model assumptions

• Underlying probability model is Binomial variance non constant, is a function of the mean
• Observation independent
• The form of the model is correct
• Liner relationship between logits and explanatory variables
• All relevant variables are included, irrelevant ones excluded
• Sample size is large enough for valid inferences
• outliers

Saturated Model

Model that fits exactly with the data
Most general model possible for the data

Consider one explanatory variable, X with \(n\) unique levels for the outcome, \(Y\sim Binomial(m,\pi)\)

Saturated Model as many parameter coefficients as \(n\)

Fitted Model nested within a FULL model, has (p+1) parameters

Null ModelL intercept only

Checking model adequacy: Form of the model - Deviance Goodness of Fit Test

Form of Hypothesis: $H_0: $ reduced model, $H_a: $ Full model

Deviance GOF test compares the fitted model \(M\) to the saturated model \(S\)

Test Statistic:

This is an asymptotic approximation, so it works better if each \(m_i>5\)

summary(fitbl)

Call:
glm(formula = cbind(Extinct, NExtinct) ~ log(Area), family = binomial, 
    data = krunnit)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7173  -0.6772   0.0973   0.4837   1.4954

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.1962     0.1184  -10.10  < 2e-16 ***
log(Area)    -0.2971     0.0549   -5.42  6.1e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 45.338  on 17  degrees of freedom
Residual deviance: 12.062  on 16  degrees of freedom
AIC: 75.39

Number of Fisher Scoring iterations: 4

Example Deviance GOF with R-output \(H_0:\) Fitted model \(logit(\pi)~\log(Area)\)
$H_a: $ Saturated model

Test statistic: \(Deviance = 12.062\) (Residual deviance)

Distribution: \(Deviance\sim \chi^2_{18-2}\)

p-value: \(P(\chi^2_{16}\leq 12.062) = 0.74\)

Conclusion: the data are consistent with \(H_0\); the simpler model with linear function of log(Area) is adequate (fit as well as the saturated model)

# calculate p-value for example
1 - pchisq(12.062, 16)

0.739700862578766

Small deviance leads to larger p-value and vice versa.

Large p-value means: - fitted model is adequate - Test is not powerful enough to detect inadequacies

Small p-value means: - fitted model is not adequate, consider a more complex model with more explanatory variables or higher order terms and so on OR - response distribution is not adequately modeled by the Binomial distribution, OR - There are severe outliers

Outliers

Response (raw) residuals

Standardized residuals
- Pearson Residuals: uses estimate of s.d. of \(Y\)

• Deviance Residuals: defined so that the sum of the square of the residuals is the deviance

# response residuals
rres = residuals(fitbl, type="response")
# pearson residuals
pres = residuals(fitbl, type="pearson")
# Deviance residuals
dres = residuals(fitbl, type="deviance")
cbind(pis, phats, rres, pres, dres)

par(mfrow=c(1,3))
plot(log(Area), rres)
plot(log(Area), pres)
plot(log(Area), dres)

​

Pearson vs Deviance Residuals

• Both asymptotically to \(N(0,1)\)
• Possible outlier if \(|res|>2\)
• Definite outlier if \(|res|>3\)
• Under small \(n\) \(D_{res}\) closer to \(N(0,1)\) then \(P_{res}\)
• When \(\hat\pi\) close to extreme (0, 1), \(P_{res}\) are unstable

Problems and Solutions

Issues Related to General Linear Models

extrapolation

• Don't make inferences/predictions outside range of observed data
• Model may no longer be appropriate

Multicollinearity
Consequences includes

• unstable fitted equation
• Coefficient that should be statistically significant while not
• Coefficient may have wrong sign
• Sometimes large \(se(\beta)\)
• Sometimes numerical procedure to find MLEs does not converge

Influential Points

Model: overfit issue, should build model on training data (cross validation)

Problems specific to Logistic Regression

Extra-binomial variation (over dispersion) - variance of \(Y_i\) greater than \(m_i\pi_i(1-\pi_i)\) - Does not bias \(\hat\beta\)'s but se of \(\hat\beta\)'s will be too small

Solution: add one more parameter to the model, \(\phi\)-dispersion parameter. Then \(var(Y_i)=\phi m_i\pi_i(1-\pi_i)\)

Complete separation

• One or a linear combination of explanatory variables perfectly predict whether \(Y=1\lor 0\)
• In Binary response, when \(y_i=1,\hat y_i=1\), then \(\sum_{i=1}^n\{y_i\log(\hat y_i + (1-y_i)\log(1-\hat y_i))=0\}\)
• MLE cannot be computed

Quasi-complete separation:

• explanatory variables predict \(Y=1\lor 0\) almost perfectly
• MLE are numerically unstable

Solution: simplify the model. Try other methods like penalized maximum likelihood, exact logistic regression, Bayesian methods

Conclusions about Logistic Regression

False Logistic regression describes population proportion of probability as a linear function of explanatory variables.

Non-linear \(\hat\pi = e^{\hat\mu}/(1+e^{\hat\mu})\)

Hence Logistic regression is a nonlinear regression model

Extra-binomial variation

Suppose \(X_1,...,X_m\) are not independent but identically distributed \(\sim Bernoulli(\pi)\). Further suppose \(\forall X_i,X_j. cov(X_1,X_j)=\rho, \rho>0\).

Let \(Y_1=\sum_i^{m_1}X_i\)

Then

Therefore, the model for variance:

estimate of \(\hat\psi\) by scaled Pearson chi-square statistic

\(\hat\psi>>1\) indicates evidence of over-dispersion

\(\psi\) does not effect \(E(Y_i)\), hence using over-dispersion does not change \(\hat\beta\)

The standard errors under this assumption is \(se_{\psi}(\hat\beta)=\sqrt{\hat\psi}se(\hat\beta)\)

The following two ways are equivalent

# manually conpute the psi hat
psihat = sum(residuals(fitbl, type="pearson")^2/fitbl$df.residual)
psihat

0.73257293737338

# apply the dispersion
summary(fitbl, dispersion=psihat)

Call:
glm(formula = cbind(Extinct, NExtinct) ~ log(Area), family = binomial, 
    data = krunnit)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7173  -0.6772   0.0973   0.4837   1.4954

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.1962     0.1014  -11.80  < 2e-16 ***
log(Area)    -0.2971     0.0469   -6.33  2.5e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 0.7326)

    Null deviance: 45.338  on 17  degrees of freedom
Residual deviance: 12.062  on 16  degrees of freedom
AIC: 75.39

Number of Fisher Scoring iterations: 4

# use quasi binomial family
fitbl2 <- glm(cbind(Extinct, NExtinct)~log(Area), family=quasibinomial, data=krunnit)
fitbl2

Call:  glm(formula = cbind(Extinct, NExtinct) ~ log(Area), family = quasibinomial, 
    data = krunnit)

Coefficients:
(Intercept)    log(Area)  
     -1.196       -0.297

Degrees of Freedom: 17 Total (i.e. Null);  16 Residual
Null Deviance:      45.3 
Residual Deviance: 12.1     AIC: NA

Using Logistic Regression for Classification

Want: \(y^*\mid (x_1^*, ...,x_p^*)= \mathbb{I}\)

Do: calculate \(\hat\pi_M^*\) be the estimated probability that \(y^*=1\) based on the fitted model given \(X_i=x_i^*\), then predict that \(y^* = \mathbb{I}(\hat\pi_M^* \text{is large enough})\)

Need: a good cutoff \(\hat\pi_M^*\)

Classification

Try different cutoffs and see which gives fewest incorrect classifications

• Useful if proportion of 1's and 0's in date reflect their relative proportions in the population
• Likely to overestimate. To overcome,cross-validation (training group vs. validation group)

Confusion Matrix

https://en.wikipedia.org/wiki/Confusion_matrix

Choose a cutoff probability based on one of the 5 criteria for success of classification that is most important to you

Examples

• High sensitivity makes good screening test
• High specificity makes a good confirmatory test
• A screening test followed by a confirmatory test is good (but expensive) diagnostic procedure