Lab Materials: 208429 Generalised Linear Models
Kamonrat Suphawan |
Department of Statistics
Faculty of Science, Chiang Mai University
Faculty of Science, Chiang Mai University
Data files
One-time setup: point R to your data folder
Every lab below reads a data file (SalesPrice.csv,
StudyHours.csv, etc.) from a DataFiles_208429
folder. This document is set up as an RStudio Project:
keep this .Rmd file and the DataFiles_208429/
folder (containing all the CSVs) side by side in the same folder, and
the relative path below will just work - no path to edit, no
setwd() needed, and it will knit identically on any
computer.
# If your .Rmd and DataFiles_208429/ folder are in the same folder, leave this as is.
data_dir <- "DataFiles_208429"
# --- Do not edit below this line ---
# Fail fast with a clear message if the data folder can't be found, instead of a
# confusing low-level error several chunks later.
if (!dir.exists(data_dir)) {
stop(
"\n\nSTOP: 'data_dir' does not point to a real folder.\n",
"Current value: ", data_dir, "\n",
"Fix this: place the DataFiles_208429 folder in the same folder as this .Rmd,\n",
"or edit data_dir in the 'setup-data-dir' chunk to the full path of that folder\n",
"on your computer (e.g. \"C:/Users/yourname/Documents/DataFiles_208429\").",
call. = FALSE
)
}Packages used in this document. Install these once before starting:
install.packages(c("MASS","aod","rms","pscl","car","lmtest","AER","nnet","readr",
"ggplot2","gridExtra","Metrics","ResourceSelection","brant","pROC"))Lab 1: Multiple Linear Regression Analysis
Objectives: Students are able to use R language to
analyse data using multiple linear regression:
1. transform qualitative independent variable into dummy variables
2. select independent variables
3. perform linear regression analysis and inference on regression
parameters
4. interpret the results
- Import data into R. ( must be in .csv form).
data = read.csv(file.path(data_dir, "SalesPrice.csv"),header=TRUE)
data# name variables for convenience
y = data$SalesPrice
x1 = data$NumberOfHousehold #quantitative independent variable
x2 = data$Location #qualitative independent variable with k=3 groups- Explore data using descriptive statistic
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 93.28 136.01 195.29 177.19 215.71 242.16hist(y,col = "darkorange",
border = "dodgerblue")
shapiro.test(y)##
## Shapiro-Wilk normality test
##
## data: y
## W = 0.92999, p-value = 0.2727summary(x1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 99.0 127.5 179.0 173.0 210.0 248.0plot(x1,y,xlab="Number Of Household",ylab="Sale price")
cor(x1,y) #Compute the r between two quantitative variables## [1] 0.9610064boxplot(y~x2,xlab="Location",ylab="Sale price")
NOTE: It is easier to work with graphs when using ggplot
data_df = data.frame(data) # set data as dataframe type
library(ggplot2) # you might need to install packgage "ggplot2" first!
ggplot(data_df,aes(x=NumberOfHousehold,y=SalesPrice))+
geom_point(aes(color=Location))
ggplot(data_df,aes(y=SalesPrice,x=Location))+
geom_boxplot(aes(color=Location))
- Since Location is qualitative independent variable with k=3
groups.
We need to transform Location into dummy variables
We create k-1=3-1=2 dummy variables, set as. x2.dummy, x3.dummy.
#create x2.dummy
x2.dummy= c() #create null vector
for(i in 1:length(data$Location)){
if(data$Location[i]=="Street") x2.dummy[i] = 1
else x2.dummy[i] = 0
}
x2.dummy## [1] 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0#create x3.dummy
x3.dummy= c()
for(i in 1:length(data$Location)){
if(data$Location[i]=="Mall") x3.dummy[i] = 1
else x3.dummy[i] = 0
}
x3.dummy## [1] 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0- Select independent variables using variable selection methods
null=lm(y~1)
full=lm(y~x1+x2.dummy+x3.dummy)
#The R function vif() [car package] can be used to detect multicollinearity in a regression model:
car::vif(full)## x1 x2.dummy x3.dummy
## 1.446801 1.881636 1.366321# Forward selection
fw.fit = step(null, scope=list(lower=null,upper=full),direction="forward")## Start: AIC=117.83
## y ~ 1
##
## Df Sum of Sq RSS AIC
## + x1 1 31278.1 2590 81.269
## + x2.dummy 1 14690.3 19178 111.302
## + x3.dummy 1 5230.6 28637 117.316
## <none> 33868 117.833
##
## Step: AIC=81.27
## y ~ x1
##
## Df Sum of Sq RSS AIC
## + x3.dummy 1 2038.47 551.29 60.063
## + x2.dummy 1 931.21 1658.56 76.585
## <none> 2589.77 81.269
##
## Step: AIC=60.06
## y ~ x1 + x3.dummy
##
## Df Sum of Sq RSS AIC
## + x2.dummy 1 84.366 466.92 59.572
## <none> 551.29 60.063
##
## Step: AIC=59.57
## y ~ x1 + x3.dummy + x2.dummyfw.fit##
## Call:
## lm(formula = y ~ x1 + x3.dummy + x2.dummy)
##
## Coefficients:
## (Intercept) x1 x3.dummy x2.dummy
## 21.958 0.868 22.101 -6.901# Backward elimination
be.fit = step(full,direction="backward")## Start: AIC=59.57
## y ~ x1 + x2.dummy + x3.dummy
##
## Df Sum of Sq RSS AIC
## <none> 466.9 59.572
## - x2.dummy 1 84.4 551.3 60.063
## - x3.dummy 1 1191.6 1658.6 76.585
## - x1 1 18527.4 18994.3 113.158be.fit##
## Call:
## lm(formula = y ~ x1 + x2.dummy + x3.dummy)
##
## Coefficients:
## (Intercept) x1 x2.dummy x3.dummy
## 21.958 0.868 -6.901 22.101# Stepwise regression
sw.fit = step(full, direction="both")## Start: AIC=59.57
## y ~ x1 + x2.dummy + x3.dummy
##
## Df Sum of Sq RSS AIC
## <none> 466.9 59.572
## - x2.dummy 1 84.4 551.3 60.063
## - x3.dummy 1 1191.6 1658.6 76.585
## - x1 1 18527.4 18994.3 113.158sw.fit##
## Call:
## lm(formula = y ~ x1 + x2.dummy + x3.dummy)
##
## Coefficients:
## (Intercept) x1 x2.dummy x3.dummy
## 21.958 0.868 -6.901 22.101- Obtain fitted equation and CI of beta
summary(sw.fit)##
## Call:
## lm(formula = y ~ x1 + x2.dummy + x3.dummy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.422 -2.989 2.243 4.572 5.852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.95824 8.78193 2.500 0.029486 *
## x1 0.86800 0.04155 20.892 3.34e-10 ***
## x2.dummy -6.90102 4.89503 -1.410 0.186240
## x3.dummy 22.10084 4.17123 5.298 0.000253 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.515 on 11 degrees of freedom
## Multiple R-squared: 0.9862, Adjusted R-squared: 0.9825
## F-statistic: 262.3 on 3 and 11 DF, p-value: 1.641e-10anova(sw.fit)confint.lm(sw.fit)## 2.5 % 97.5 %
## (Intercept) 2.6293467 41.2871246
## x1 0.7765583 0.9594473
## x2.dummy -17.6749031 3.8728631
## x3.dummy 12.9200368 31.2816515What the model is and how to read this output
Multiple linear regression models a continuous response (\(y\) = sale price) as a linear function of
one or more predictors: \[y = \beta_0 +
\beta_1 x_1 + \beta_2 x_{2.dummy} + \beta_3 x_{3.dummy} + \varepsilon,
\quad \varepsilon \sim N(0,\sigma^2)\] Each \(\beta\) is estimated by least squares, and
summary() reports the estimate, its standard error, a
t-statistic, and a p-value testing \(H_0:
\beta = 0\).
Running this exact model on SalesPrice.csv gives
approximately:
| Term | Estimate | p-value | Interpretation |
|---|---|---|---|
| (Intercept) | 21.96 | 0.029 | Expected sale price when NumberOfHousehold = 0 and the
location is the reference category (Mall was dropped/kept depending on
selection - check your own run) |
x1 (NumberOfHousehold) |
0.87 | <0.001 | Holding location fixed, each additional household in the area is associated with a 0.87 increase in sale price |
x2.dummy (Street) |
-6.90 | 0.19 | Not statistically significant - being on a Street location does not
differ detectably from the reference category once x1 is
accounted for |
x3.dummy (Mall) |
22.10 | <0.001 | Being in a Mall location is associated with a 22.10
increase in sale price relative to the reference category,
holding x1 fixed |
The model explains about 98.6% of the variance in
sale price (\(R^2 \approx 0.986\)) -
very high for this small, illustrative dataset. In general:
- A significant, positive coefficient means that
predictor is associated with higher sale price, holding the
others constant.
- A non-significant coefficient (like
x2.dummy here) means we cannot rule out that its true
effect is zero - don’t over-interpret its sign or magnitude.
- \(R^2\) tells you how much of the
variation in \(y\) the model
captures, not whether individual predictors are significant - a model
can have high \(R^2\) with some
non-significant terms, and vice versa.
- Plot the fitted:
new = data.frame(x1=x1,x2.dummy=x2.dummy,x3.dummy=x3.dummy)
yhat = predict(sw.fit,newdata=new) #compute fitted y
yhat## 1 2 3 4 5 6 7 8
## 154.8057 100.9895 132.2376 119.2176 157.4097 235.8877 199.4316 221.1317
## 9 10 11 12 13 14 15
## 229.8117 131.7274 222.4669 200.7668 237.2229 114.8345 199.8988plot(x1,y,pch=1,xlab="Number Of Household",ylab="Sale price")
points(x1,yhat,type="p",pch=20)
legend("topleft",c("observed y","predicted y"),pch=c(1,20))
- Analyse the residuals:
par(mfrow=c(2,2)) #set plot layout as 2 row 2 column
plot(sw.fit)
#get the residuals
res = resid(sw.fit)
par(mfrow=c(1,2)) #set plot layout as 1 row 2 column
plot(res)
hist(res)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.87123, p-value = 0.03518# Performs the Durbin-Watson test for autocorrelation of disturbances.
lmtest::dwtest(sw.fit)##
## Durbin-Watson test
##
## data: sw.fit
## DW = 2.611, p-value = 0.7911
## alternative hypothesis: true autocorrelation is greater than 0# Test the constant variance assumption: Breusch-Pagan Test
lmtest::bptest(sw.fit)##
## studentized Breusch-Pagan test
##
## data: sw.fit
## BP = 2.2962, df = 3, p-value = 0.5132Lab 1 (Multiple Linear Regression) Assignment
Using the SalesPrice.csv dataset:
- Create the two dummy variables for
Locationand confirm they matchx2.dummyandx3.dummyabove.
- Fit the full model and perform forward, backward, and stepwise
selection. Do all three methods agree on the same final model? If not,
which would you report and why?
- Write out the fitted regression equation from your final model,
being explicit about what each dummy variable represents.
- Check the four assumptions of linear regression (linearity,
independence, homoscedasticity, normality of residuals) using the
diagnostic plots, the Durbin-Watson test, the Breusch-Pagan test, and
the Shapiro-Wilk test. Report your conclusions for each.
- Check for multicollinearity using
car::vif(). Is multicollinearity a concern for this model?
- Identify any influential observations using Cook’s distance (hint:
plot(sw.fit, which = 4)orcooks.distance(sw.fit)), and comment on whether they should be investigated further.
Lab 2: Binary logistic regression
Objectives: Students are able to use R to analyse a
binary response variable using logistic regression:
1. fit a binary logistic regression model and interpret the coefficients
on the logit and odds scale
2. perform inference on the parameters (Wald test, likelihood ratio
test)
3. assess overall model fit (pseudo-\(R^2\), Hosmer-Lemeshow test)
4. use the fitted model to predict probabilities and classify
observations
5. check model assumptions and influential observations (linearity in
the logit, multicollinearity, Cook’s distance, leverage)
Task 1: Example I: Study hours and Exam
What the model is
Binary logistic regression models the probability of a binary outcome (\(y\) = 1 if passed, 0 if failed) as a function of predictors, using the logit (log-odds) link: \[\ln\left(\frac{P(y=1)}{1-P(y=1)}\right) = \beta_0 + \beta_1 x\] Unlike linear regression, we can’t interpret \(\beta_1\) directly as “a change in probability” - it’s a change in the log-odds. Exponentiating \(\beta_1\) converts it to an odds ratio (OR): the multiplicative change in the odds of the outcome for each one-unit increase in \(x\).
# Load the packages
library(aod) #for 'Wald test'
library(MASS)
library(rms) #for 'lrm'
library(pscl) #to get pseudo R2
data = read.csv(file.path(data_dir, "StudyHours.csv"),header=TRUE)
head(data)y = data$Pass
x = data$Hours
# Generating the frequency table
table(y)## y
## 0 1
## 10 10# 1.Building the Model
fit.logist = glm(y ~x,family=binomial(link='logit'))
fit.logist##
## Call: glm(formula = y ~ x, family = binomial(link = "logit"))
##
## Coefficients:
## (Intercept) x
## -4.078 1.505
##
## Degrees of Freedom: 19 Total (i.e. Null); 18 Residual
## Null Deviance: 27.73
## Residual Deviance: 16.06 AIC: 20.06anova(fit.logist)summary(fit.logist)##
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0777 1.7610 -2.316 0.0206 *
## x 1.5046 0.6287 2.393 0.0167 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 16.060 on 18 degrees of freedom
## AIC: 20.06
##
## Number of Fisher Scoring iterations: 5# 2. Perform overall fit test
mod <- lrm(y ~x, data = data)
mod## Logistic Regression Model
##
## lrm(formula = y ~ x, data = data)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 20 LR chi2 11.67 R2 0.589 C 0.895
## 0 10 d.f. 1 R2(1,20) 0.413 Dxy 0.790
## 1 10 Pr(> chi2) 0.0006 R2(1,15) 0.509 gamma 0.798
## max |deriv| 1e-07 Brier 0.137 tau-a 0.416
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -4.0777 1.7610 -2.32 0.0206
## x 1.5046 0.6287 2.39 0.0167# 3. Perform Wald test:
#wald.test(b = coef(fit.logist), Sigma = vcov(fit.logist), Terms = 1)
wald.test(b = coef(fit.logist), Sigma = vcov(fit.logist), Terms = 2)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 5.7, df = 1, P(> X2) = 0.017# 4. Perform lr test
null <- glm(y ~ 1, data = data,family="binomial")
summary(null)##
## Call:
## glm(formula = y ~ 1, family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.930e-17 4.472e-01 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 27.726 on 19 degrees of freedom
## AIC: 29.726
##
## Number of Fisher Scoring iterations: 2full <- glm(y ~ x, data = data,family="binomial")
summary(full)##
## Call:
## glm(formula = y ~ x, family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0777 1.7610 -2.316 0.0206 *
## x 1.5046 0.6287 2.393 0.0167 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.726 on 19 degrees of freedom
## Residual deviance: 16.060 on 18 degrees of freedom
## AIC: 20.06
##
## Number of Fisher Scoring iterations: 5lrtest(null, full)##
## Model 1: y ~ 1
## Model 2: y ~ x
##
## L.R. Chisq d.f. P
## 1.166613e+01 1.000000e+00 6.364826e-04# 5. get LL and pseudo R2
logLik(full)## 'log Lik.' -8.029878 (df=2)#To easily get a McFadden's pseudo R2 for a fitted model in R, use the "pscl" package by Simon Jackman and use the pR2 command.
pR2(full)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -8.0298785 -13.8629436 11.6661303 0.4207667 0.4419499 0.5892665mod1b <- lrm(y ~ x, data = data)
mod1b## Logistic Regression Model
##
## lrm(formula = y ~ x, data = data)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 20 LR chi2 11.67 R2 0.589 C 0.895
## 0 10 d.f. 1 R2(1,20) 0.413 Dxy 0.790
## 1 10 Pr(> chi2) 0.0006 R2(1,15) 0.509 gamma 0.798
## max |deriv| 1e-07 Brier 0.137 tau-a 0.416
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -4.0777 1.7610 -2.32 0.0206
## x 1.5046 0.6287 2.39 0.0167## CIs using profiled log-likelihood
confint(fit.logist)## 2.5 % 97.5 %
## (Intercept) -8.5853315 -1.268378
## x 0.5293557 3.145138## CIs using standard errors
confint.default(fit.logist)## 2.5 % 97.5 %
## (Intercept) -7.5291793 -0.6262476
## x 0.2723838 2.7369071logLik(fit.logist)## 'log Lik.' -8.029878 (df=2)## odds ratios only
exp(coef(fit.logist))## (Intercept) x
## 0.01694617 4.50255687## odds ratios and 95% CI
exp(cbind(OR = coef(fit.logist), confint(fit.logist)))## OR 2.5 % 97.5 %
## (Intercept) 0.01694617 0.0001868263 0.2812875
## x 4.50255687 1.6978380343 23.2228735Reading these numbers: fitting this exact model on
StudyHours.csv gives a coefficient on Hours of
about 1.50 (p = 0.017, so significant at the 5% level),
which exponentiates to an odds ratio of about 4.5. In
words: each additional hour of study multiplies the odds of passing
the exam by roughly 4.5, holding nothing else constant (there are no
other predictors here). Since the 95% CI for the OR excludes 1,
we’re confident study hours genuinely help. If the OR had been close to
1 (and its CI straddled 1), that would indicate no meaningful
association. The intercept (on the log-odds scale) represents the
log-odds of passing when Hours = 0 - exponentiating and
converting to a probability with 1/(1+exp(-intercept))
gives the baseline pass probability for zero hours of study.
pred_logit = predict(fit.logist, newdata = list(x=x), type = "link",se = TRUE)
pred_logit## $fit
## 1 2 3 4 5 6 7
## -3.3253907 -2.9492294 -2.5730680 -2.1969066 -1.8207453 -1.4445839 -1.4445839
## 8 9 10 11 12 13 14
## -1.0684226 -0.6922612 -0.3160999 0.0600615 0.4362229 0.8123842 1.1885456
## 15 16 17 18 19 20
## 1.9408683 2.3170296 2.6931910 3.0693524 3.4455137 4.1978364
##
## $se.fit
## 1 2 3 4 5 6 7 8
## 1.4715288 1.3310441 1.1947260 1.0641770 0.9417992 0.8312096 0.8312096 0.7377287
## 9 10 11 12 13 14 15 16
## 0.6685718 0.6317781 0.6330161 0.6720757 0.7430159 0.8377769 1.0722391 1.2032148
## 17 18 19 20
## 1.3398379 1.4805458 1.6242772 1.9180797
##
## $residual.scale
## [1] 1pred = predict(fit.logist, newdata = list(x=x), type = "response",se = TRUE)
pred## $fit
## 1 2 3 4 5 6 7
## 0.03471034 0.04977295 0.07089196 0.10002862 0.13934447 0.19083650 0.19083650
## 8 9 10 11 12 13 14
## 0.25570318 0.33353024 0.42162653 0.51501086 0.60735865 0.69261733 0.76648084
## 15 16 17 18 19 20
## 0.87444750 0.91027764 0.93662366 0.95561071 0.96909707 0.98519444
##
## $se.fit
## 1 2 3 4 5 6 7
## 0.04930435 0.06295253 0.07869217 0.09580030 0.11294771 0.12835367 0.12835367
## 8 9 10 11 12 13 14
## 0.14040383 0.14861538 0.15406389 0.15811139 0.16027266 0.15818702 0.14995198
## 15 16 17 18 19 20
## 0.11772013 0.09826927 0.07953248 0.06280310 0.04864376 0.02797779
##
## $residual.scale
## [1] 1pred_prob = 1/(1+exp(-pred_logit$fit))
pred_prob## 1 2 3 4 5 6 7
## 0.03471034 0.04977295 0.07089196 0.10002862 0.13934447 0.19083650 0.19083650
## 8 9 10 11 12 13 14
## 0.25570318 0.33353024 0.42162653 0.51501086 0.60735865 0.69261733 0.76648084
## 15 16 17 18 19 20
## 0.87444750 0.91027764 0.93662366 0.95561071 0.96909707 0.98519444LCI_pred = 1/(1+exp(-(pred_logit$fit-1.96*pred_logit$se.fit)))
UCI_pred = 1/(1+exp(-(pred_logit$fit+1.96*pred_logit$se.fit)))
plot(x,pred_prob,type="o",ylim=c(-0.1,1.1),xlab="Number of Study Hours",ylab="Prob of passing exam",col="blue")
points(x,y)
lines(x,LCI_pred,col="red",lty=2)
lines(x,UCI_pred,col="red",lty=2)
abline(h=0.5,col="gray")
#create data.frame
dtf = data.frame( y =y ,
x =x ,
pred_logit = pred_logit,
pred_prob = pred_prob,
LCI_pred = LCI_pred,
UCI_pred = UCI_pred
)
dtf#create another variable called "pred_pass"
# if pred_prob <= 0.5 , pred_pass = 0
# if pred_prob > 0.5 , pred_pass = 1
#create x2.dummy
pred_pass= c() #create null vector
for(i in 1:length(y)){
if(pred_prob[i] <= 0.5) pred_pass[i] = 0
else pred_pass[i] = 1
}
pred_pass## [1] 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1#add pred_pass in dtf
dtf$pred_pass = pred_pass
dtf#Classification table
xtabs(~ y + pred_pass, data = dtf)## pred_pass
## y 0 1
## 0 8 2
## 1 2 8Model diagnostics for Task 1
Fitting a model and reading summary() is not the end of
the analysis: we should also check whether the model’s assumptions are
reasonable and whether any single observation is unduly influencing the
fit.
# Hosmer-Lemeshow goodness-of-fit test
# H0: the model fits the data well (large p-value = no evidence of lack of fit)
library(ResourceSelection)
hoslem.test(fit.logist$y, fitted(fit.logist), g = 10)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: fit.logist$y, fitted(fit.logist)
## X-squared = 2.2455, df = 8, p-value = 0.9725# ROC curve and AUC: how well does the model discriminate Pass vs Fail?
library(pROC)
roc_obj <- roc(y, fitted(fit.logist))
plot(roc_obj, main = "ROC curve: Study Hours model")
auc(roc_obj)## Area under the curve: 0.895# Influence diagnostics: Cook's distance, leverage (hat values), standardized residuals
infl <- influence.measures(fit.logist)
summary(infl)## Potentially influential observations of
## glm(formula = y ~ x, family = binomial(link = "logit")) :
##
## dfb.1_ dfb.x dffit cov.r cook.d hat
## 7 0.87 -0.65 0.99 0.57_* 0.28 0.11cooksd <- cooks.distance(fit.logist)
plot(cooksd, type = "h", main = "Cook's distance", ylab = "Cook's D")
abline(h = 4/length(cooksd), col = "red", lty = 2) # common rule-of-thumb cutoff
hatvalues_ <- hatvalues(fit.logist)
plot(hatvalues_, type = "h", main = "Leverage (hat values)", ylab = "Leverage")
abline(h = 2*mean(hatvalues_), col = "red", lty = 2)
# Standardized (deviance) residuals vs fitted values - look for patterns/outliers
plot(fitted(fit.logist), rstandard(fit.logist, type = "deviance"),
xlab = "Fitted probability", ylab = "Standardized deviance residual")
abline(h = c(-2, 2), col = "red", lty = 2)
Interpretation notes:
- A Hosmer-Lemeshow p-value < 0.05 suggests the model does not fit
the data well and the functional form (e.g. linearity in the logit)
should be reconsidered.
- Points with Cook’s distance well above the 4/n cutoff, or leverage
well above 2 times the mean leverage, deserve a closer look - refit the
model without them and see whether conclusions change.
- An AUC close to 0.5 indicates the model discriminates no better than
chance; AUC close to 1 indicates excellent discrimination.
Task 2: Admission data
researcher is interested in how variables, such as GRE (Graduate
Record Exam scores), GPA (grade point average) and prestige of the
undergraduate institution, effect admission into graduate school. The
response variable, admit/don’t admit, is a binary variable.
# Load the packages
library(aod)
library(MASS)
library(rms)
data2 = read.csv(file.path(data_dir, "Admission.csv"),header=TRUE)
head(data2)summary(data2)## admit gre gpa rank
## Min. :0.0000 Min. :220.0 Min. :2.260 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:520.0 1st Qu.:3.130 1st Qu.:2.000
## Median :0.0000 Median :580.0 Median :3.395 Median :2.000
## Mean :0.3175 Mean :587.7 Mean :3.390 Mean :2.485
## 3rd Qu.:1.0000 3rd Qu.:660.0 3rd Qu.:3.670 3rd Qu.:3.000
## Max. :1.0000 Max. :800.0 Max. :4.000 Max. :4.000sapply(data2, sd)## admit gre gpa rank
## 0.4660867 115.5165364 0.3805668 0.9444602xtabs(~admit + rank, data = data2)## rank
## admit 1 2 3 4
## 0 28 97 93 55
## 1 33 54 28 12#Using the logit model
data2$rank <- factor(data2$rank)
#Perform variable selection
null <- glm(admit ~ 1, data = data2,family="binomial")
full <- glm(admit ~ gre + gpa + rank, data = data2,family="binomial")
forwards = step(null,scope=list(lower=null,upper=full), direction="forward")## Start: AIC=501.98
## admit ~ 1
##
## Df Deviance AIC
## + rank 3 474.97 482.97
## + gre 1 486.06 490.06
## + gpa 1 486.97 490.97
## <none> 499.98 501.98
##
## Step: AIC=482.97
## admit ~ rank
##
## Df Deviance AIC
## + gpa 1 462.88 472.88
## + gre 1 464.53 474.53
## <none> 474.97 482.97
##
## Step: AIC=472.88
## admit ~ rank + gpa
##
## Df Deviance AIC
## + gre 1 458.52 470.52
## <none> 462.88 472.88
##
## Step: AIC=470.52
## admit ~ rank + gpa + grebackwards = step(full)## Start: AIC=470.52
## admit ~ gre + gpa + rank
##
## Df Deviance AIC
## <none> 458.52 470.52
## - gre 1 462.88 472.88
## - gpa 1 464.53 474.53
## - rank 3 480.34 486.34bothways = step(null, list(lower=null,upper=full),direction="both")## Start: AIC=501.98
## admit ~ 1
##
## Df Deviance AIC
## + rank 3 474.97 482.97
## + gre 1 486.06 490.06
## + gpa 1 486.97 490.97
## <none> 499.98 501.98
##
## Step: AIC=482.97
## admit ~ rank
##
## Df Deviance AIC
## + gpa 1 462.88 472.88
## + gre 1 464.53 474.53
## <none> 474.97 482.97
## - rank 3 499.98 501.98
##
## Step: AIC=472.88
## admit ~ rank + gpa
##
## Df Deviance AIC
## + gre 1 458.52 470.52
## <none> 462.88 472.88
## - gpa 1 474.97 482.97
## - rank 3 486.97 490.97
##
## Step: AIC=470.52
## admit ~ rank + gpa + gre
##
## Df Deviance AIC
## <none> 458.52 470.52
## - gre 1 462.88 472.88
## - gpa 1 464.53 474.53
## - rank 3 480.34 486.34mod1b <- lrm(admit ~ gre + gpa + rank, data = data2)
mod1b## Logistic Regression Model
##
## lrm(formula = admit ~ gre + gpa + rank, data = data2)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 400 LR chi2 41.46 R2 0.138 C 0.693
## 0 273 d.f. 5 R2(5,400)0.087 Dxy 0.386
## 1 127 Pr(> chi2) <0.0001 R2(5,260)0.131 gamma 0.386
## max |deriv| 2e-06 Brier 0.195 tau-a 0.168
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -3.9900 1.1400 -3.50 0.0005
## gre 0.0023 0.0011 2.07 0.0385
## gpa 0.8040 0.3318 2.42 0.0154
## rank=2 -0.6754 0.3165 -2.13 0.0328
## rank=3 -1.3402 0.3453 -3.88 0.0001
## rank=4 -1.5515 0.4178 -3.71 0.0002#lrtest
fit <- glm(admit ~ gre + gpa + rank, data = data2, family = "binomial")
summary(fit)##
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "binomial",
## data = data2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.989979 1.139951 -3.500 0.000465 ***
## gre 0.002264 0.001094 2.070 0.038465 *
## gpa 0.804038 0.331819 2.423 0.015388 *
## rank2 -0.675443 0.316490 -2.134 0.032829 *
## rank3 -1.340204 0.345306 -3.881 0.000104 ***
## rank4 -1.551464 0.417832 -3.713 0.000205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.52
##
## Number of Fisher Scoring iterations: 4confint(fit)## 2.5 % 97.5 %
## (Intercept) -6.2716202334 -1.792547080
## gre 0.0001375921 0.004435874
## gpa 0.1602959439 1.464142727
## rank2 -1.3008888002 -0.056745722
## rank3 -2.0276713127 -0.670372346
## rank4 -2.4000265384 -0.753542605confint.default(fit)## 2.5 % 97.5 %
## (Intercept) -6.2242418514 -1.755716295
## gre 0.0001202298 0.004408622
## gpa 0.1536836760 1.454391423
## rank2 -1.2957512650 -0.055134591
## rank3 -2.0169920597 -0.663415773
## rank4 -2.3703986294 -0.732528724wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 4:6)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 20.9, df = 3, P(> X2) = 0.00011pR2(fit)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML
## -229.25874624 -249.98825878 41.45902508 0.08292194 0.09845702
## r2CU
## 0.13799580lrm(admit ~ gre + gpa + rank, data = data2)## Logistic Regression Model
##
## lrm(formula = admit ~ gre + gpa + rank, data = data2)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 400 LR chi2 41.46 R2 0.138 C 0.693
## 0 273 d.f. 5 R2(5,400)0.087 Dxy 0.386
## 1 127 Pr(> chi2) <0.0001 R2(5,260)0.131 gamma 0.386
## max |deriv| 2e-06 Brier 0.195 tau-a 0.168
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -3.9900 1.1400 -3.50 0.0005
## gre 0.0023 0.0011 2.07 0.0385
## gpa 0.8040 0.3318 2.42 0.0154
## rank=2 -0.6754 0.3165 -2.13 0.0328
## rank=3 -1.3402 0.3453 -3.88 0.0001
## rank=4 -1.5515 0.4178 -3.71 0.0002newdata1 <- with(data2, data.frame(gre = gre, gpa = gpa, rank = factor(rank)))
pred.response <- predict(fit, newdata = newdata1, type = "response")
head(pred.response)## 1 2 3 4 5 6
## 0.1726265 0.2921750 0.7384082 0.1783846 0.1183539 0.3699699newdata1$rankP <- predict(fit, newdata = newdata1, type = "response")
head(newdata1)Reading these numbers: with gre,
gpa, and rank all in the model, a
representative fit gives odds ratios of about gre =
1.00 (per 1-point increase - GRE scores range in the
hundreds, so this tiny-looking per-point effect adds up),
gpa = 2.24 (each 1-point increase in GPA
more than doubles the odds of admission), and for rank
(relative to the most prestigious rank 1 institutions): rank 2 ≈
0.51, rank 3 ≈ 0.26, rank 4 ≈
0.21. Since rank is categorical, each of
these compares that rank directly to the reference (rank 1) holding
gre and gpa fixed - e.g. students from
rank-4 institutions have roughly 21% of the admission odds of
otherwise-identical rank-1 students. Notice the ORs decrease
monotonically from rank 2 to rank 4 - exactly the pattern you’d expect
if institutional prestige matters for admission, which is a good
face-validity check on the model.
Visualizing predicted probabilities across rank (Task
2)
A useful way to communicate a fitted logistic model is to plot
predicted probabilities, with confidence bands, across a range of one
predictor - here gre - separately for each level of a
categorical predictor - here rank.
# predict across a grid of gre values, holding gpa at its mean, for every rank
newdata2 <- with(data2, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 4),
gpa = mean(gpa),
rank = factor(rep(1:4, each = 100))))
PredictedProb <- predict(fit, newdata = newdata2, type = "response", se = TRUE)
PredictedProb_LL <- PredictedProb$fit - 1.96*PredictedProb$se.fit
PredictedProb_UL <- PredictedProb$fit + 1.96*PredictedProb$se.fit
dtf <- cbind.data.frame(newdata2, data.frame(PredictedProb = PredictedProb$fit,
PredictedProb_LL = PredictedProb_LL,
PredictedProb_UL = PredictedProb_UL))
head(dtf)ggplot(dtf, aes(x = gre, y = PredictedProb)) +
geom_ribbon(aes(ymin = PredictedProb_LL, ymax = PredictedProb_UL, fill = rank), alpha = 0.2) +
geom_line(aes(colour = rank), linewidth = 1) +
labs(title = "Predicted probability of admission by GRE score and institution rank",
x = "GRE score", y = "Predicted P(admit)")
# same idea, varying gpa instead, holding gre at its mean
newdata3 <- with(data2, data.frame(gre = mean(gre),
gpa = rep(seq(from = 0, to = 4, length.out = 100), 4),
rank = factor(rep(1:4, each = 100))))
PredictedProb <- predict(fit, newdata = newdata3, type = "response", se = TRUE)
PredictedProb_LL <- PredictedProb$fit - 1.96*PredictedProb$se.fit
PredictedProb_UL <- PredictedProb$fit + 1.96*PredictedProb$se.fit
dtf2 <- cbind.data.frame(newdata3, data.frame(PredictedProb = PredictedProb$fit,
PredictedProb_LL = PredictedProb_LL,
PredictedProb_UL = PredictedProb_UL))
ggplot(dtf2, aes(x = gpa, y = PredictedProb)) +
geom_ribbon(aes(ymin = PredictedProb_LL, ymax = PredictedProb_UL, fill = rank), alpha = 0.2) +
geom_line(aes(colour = rank), linewidth = 1) +
labs(title = "Predicted probability of admission by GPA and institution rank",
x = "GPA", y = "Predicted P(admit)")
Checking multicollinearity and influential points (Task 2)
fit has three predictors (gre,
gpa, rank), so it is worth checking whether
they are highly correlated with each other (multicollinearity), which
would inflate standard errors and make coefficients unstable.
# Variance Inflation Factor: VIF > 5 (some use >10) suggests problematic multicollinearity
library(car)
vif(fit)## GVIF Df GVIF^(1/(2*Df))
## gre 1.134377 1 1.065071
## gpa 1.155902 1 1.075129
## rank 1.025759 3 1.004248# Influential observations
cooksd <- cooks.distance(fit)
plot(cooksd, type = "h", main = "Cook's distance: Admission model")
abline(h = 4/nrow(data2), col = "red", lty = 2)
# Rows flagged as influential
influential <- which(cooksd > 4/nrow(data2))
data2[influential, ]Task 3: Buying Cosmetics
# Load the packages
library(aod) #for 'Wald test'
library(MASS)
library(rms) #for 'lrm'
library(pscl) #to get pseudo R2
data3 = read.csv(file.path(data_dir, "BuyingCosmetics.csv"),header=TRUE)
head(data3)summary(data3)## Obs Buying Age Income
## Min. : 1.00 Min. :0.0000 Min. :15.00 Min. : 3.00
## 1st Qu.:15.75 1st Qu.:0.0000 1st Qu.:22.00 1st Qu.:12.00
## Median :30.50 Median :1.0000 Median :30.00 Median :20.00
## Mean :30.50 Mean :0.5333 Mean :32.35 Mean :21.63
## 3rd Qu.:45.25 3rd Qu.:1.0000 3rd Qu.:40.00 3rd Qu.:30.00
## Max. :60.00 Max. :1.0000 Max. :60.00 Max. :51.00#Using the logit model
#Perform variable selection
null <- glm(Buying ~ 1, data = data3,family="binomial")
full <- glm(Buying ~ Age + Income, data = data3,family="binomial")
summary(full)##
## Call:
## glm(formula = Buying ~ Age + Income, family = "binomial", data = data3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.56798 1.22985 -3.714 0.000204 ***
## Age 0.16265 0.05200 3.128 0.001761 **
## Income -0.01801 0.04232 -0.426 0.670439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 82.911 on 59 degrees of freedom
## Residual deviance: 59.067 on 57 degrees of freedom
## AIC: 65.067
##
## Number of Fisher Scoring iterations: 5lrtest(null, full)##
## Model 1: Buying ~ 1
## Model 2: Buying ~ Age + Income
##
## L.R. Chisq d.f. P
## 2.384415e+01 2.000000e+00 6.642161e-06wald.test(b = coef(full), Sigma = vcov(full), Terms = 1)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 13.8, df = 1, P(> X2) = 2e-04wald.test(b = coef(full), Sigma = vcov(full), Terms = 2)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 9.8, df = 1, P(> X2) = 0.0018wald.test(b = coef(full), Sigma = vcov(full), Terms = 3)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 0.18, df = 1, P(> X2) = 0.67mod = lrm(Buying ~ Age + Income, data = data3)
mod## Logistic Regression Model
##
## lrm(formula = Buying ~ Age + Income, data = data3)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 60 LR chi2 23.84 R2 0.438 C 0.835
## 0 28 d.f. 2 R2(2,60) 0.305 Dxy 0.671
## 1 32 Pr(> chi2) <0.0001 R2(2,44.8)0.386 gamma 0.672
## max |deriv| 1e-08 Brier 0.167 tau-a 0.340
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -4.5680 1.2299 -3.71 0.0002
## Age 0.1626 0.0520 3.13 0.0018
## Income -0.0180 0.0423 -0.43 0.6704-2*logLik(full)## 'log Lik.' 59.06665 (df=3)#To easily get a McFadden's pseudo R2 for a fitted model in R, use the "pscl" package by Simon Jackman and use the pR2 command.
pR2(full)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -29.5333253 -41.4553986 23.8441464 0.2875880 0.3279365 0.4378993#R2 McFadden
1-as.numeric(logLik(full))/as.numeric(logLik(null))## [1] 0.287588#r2 Cox and snell
n=dim(data3)[1]
r2_cs = 1-(as.numeric(exp(logLik(null))/exp(logLik(full))))^(2/n)
r2_cs## [1] 0.3279365#r2 Nagelkerki
r2_N = r2_cs/(1-as.numeric(exp(logLik(null)))^(2/n))
r2_N## [1] 0.4378993bothways = step(null, list(lower=null,upper=full),direction="both")## Start: AIC=84.91
## Buying ~ 1
##
## Df Deviance AIC
## + Age 1 59.250 63.250
## + Income 1 73.593 77.593
## <none> 82.911 84.911
##
## Step: AIC=63.25
## Buying ~ Age
##
## Df Deviance AIC
## <none> 59.250 63.250
## + Income 1 59.067 65.067
## - Age 1 82.911 84.911summary(bothways)##
## Call:
## glm(formula = Buying ~ Age, family = "binomial", data = data3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.51528 1.22025 -3.700 0.000215 ***
## Age 0.14847 0.03878 3.828 0.000129 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 82.911 on 59 degrees of freedom
## Residual deviance: 59.250 on 58 degrees of freedom
## AIC: 63.25
##
## Number of Fisher Scoring iterations: 5confint(bothways)## 2.5 % 97.5 %
## (Intercept) -7.2076248 -2.3490179
## Age 0.0805232 0.2349225confint.default(bothways)## 2.5 % 97.5 %
## (Intercept) -6.9069200 -2.123635
## Age 0.0724545 0.224479## odds ratios and 95% CI
exp(cbind(OR = coef(bothways), confint(bothways)))## OR 2.5 % 97.5 %
## (Intercept) 0.01094057 0.0007409149 0.09546287
## Age 1.16005425 1.0838539957 1.26481079newdata1 <- with(data3, data.frame(Age=Age, Income=Income))
pred_response <- predict(bothways, newdata = newdata1, type = "response")
head(pred_response)## 1 2 3 4 5 6
## 0.80586986 0.89712932 0.17567316 0.48468624 0.97467705 0.09209862mod1b <- lrm(Buying ~ Age, data = data3)
mod1b## Logistic Regression Model
##
## lrm(formula = Buying ~ Age, data = data3)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 60 LR chi2 23.66 R2 0.435 C 0.827
## 0 28 d.f. 1 R2(1,60) 0.315 Dxy 0.654
## 1 32 Pr(> chi2) <0.0001 R2(1,44.8)0.397 gamma 0.677
## max |deriv| 3e-09 Brier 0.167 tau-a 0.331
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -4.5153 1.2202 -3.70 0.0002
## Age 0.1485 0.0388 3.83 0.0001Reading these numbers: the full model gives odds
ratios of Age ≈ 1.18 (p = 0.002,
significant) and Income ≈ 0.98 (p = 0.67,
not significant). Stepwise selection agrees, dropping
Income and keeping Age alone with OR ≈
1.16: each additional year of age multiplies the
odds of buying the cosmetics product by about 16%, holding nothing
else in the final model. Income looked like it might matter
on its own, but once Age is in the model its effect is no
longer distinguishable from zero - a good reminder that a predictor’s
marginal association with the outcome (e.g. from a simple
boxplot) doesn’t always survive once you control for other correlated
variables.
Task 4: Newborn birth weight
A classic public-health question: which maternal factors are
associated with delivering a low birth weight baby
(LOW = 1 if birth weight < 2500g)? This example uses
age (maternal age), lwt (mother’s weight in
pounds at last menstrual period), and smoke (whether the
mother smoked during pregnancy).
library(readr)
dat1 <- read_csv(file.path(data_dir, "NewbornWeight.csv"))
#change data type
dat1$LOW <- as.factor(dat1$LOW)
dat1$smoke <- as.factor(dat1$smoke)
#descriptive stat
summary(dat1)## Obs LOW age lwt smoke
## Min. : 1.00 0:14 Min. :24.00 Min. :100.0 0:14
## 1st Qu.: 8.25 1:16 1st Qu.:28.00 1st Qu.:120.2 1:16
## Median :15.50 Median :33.50 Median :142.0
## Mean :15.50 Mean :32.33 Mean :144.0
## 3rd Qu.:22.75 3rd Qu.:36.00 3rd Qu.:165.5
## Max. :30.00 Max. :39.00 Max. :237.0boxplot(dat1$age, main = "Maternal age")
boxplot(dat1$age~dat1$LOW, xlab = "Low birth weight", ylab = "Maternal age")
boxplot(lwt~smoke, data=dat1, xlab = "Smoker", ylab = "Mother's weight (lwt)")
plot(dat1$age, dat1$lwt, xlab = "Maternal age", ylab = "Mother's weight (lwt)")
table(dat1$LOW, dat1$smoke)##
## 0 1
## 0 11 3
## 1 3 13aggregate(age~smoke, data=dat1, FUN=summary)#1. binary logistic regression
null <- glm(LOW ~ 1, data=dat1, family = "binomial")
fit <- glm(LOW ~ age+lwt+smoke, data=dat1, family = "binomial")
summary(fit)##
## Call:
## glm(formula = LOW ~ age + lwt + smoke, family = "binomial", data = dat1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.35323 4.50831 -2.075 0.0380 *
## age 0.21000 0.11519 1.823 0.0683 .
## lwt 0.01103 0.01565 0.705 0.4811
## smoke1 2.18617 1.00810 2.169 0.0301 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.455 on 29 degrees of freedom
## Residual deviance: 25.865 on 26 degrees of freedom
## AIC: 33.865
##
## Number of Fisher Scoring iterations: 4#2. likelihood ratio test
library(rms)
lrm(LOW ~ age + lwt + smoke, data = dat1)## Logistic Regression Model
##
## lrm(formula = LOW ~ age + lwt + smoke, data = dat1)
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 30 LR chi2 15.59 R2 0.541 C 0.875
## 0 14 d.f. 3 R2(3,30) 0.343 Dxy 0.750
## 1 16 Pr(> chi2) 0.0014 R2(3,22.4)0.430 gamma 0.750
## max |deriv| 3e-09 Brier 0.130 tau-a 0.386
##
## Coef S.E. Wald Z Pr(>|Z|)
## Intercept -9.3532 4.5086 -2.07 0.0380
## age 0.2100 0.1152 1.82 0.0683
## lwt 0.0110 0.0157 0.70 0.4811
## smoke=1 2.1862 1.0081 2.17 0.0301anova(null,fit)#3. get LL and pseudo R2
library(pscl)
pR2(fit)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -12.9322642 -20.7276993 15.5908701 0.3760878 0.4052985 0.5412022#4. odds ratios
exp(fit$coefficients)## (Intercept) age lwt smoke1
## 8.668468e-05 1.233676e+00 1.011090e+00 8.901032e+00exp(confint(fit))## 2.5 % 97.5 %
## (Intercept) 2.360282e-09 0.2631688
## age 9.939895e-01 1.5922943
## lwt 9.813108e-01 1.0471535
## smoke1 1.317240e+00 76.6756271#5. find p_hat(Y=1) when age=27, lwt=120, for a smoker vs non-smoker
newdata1 <- data.frame(age = 27, lwt = 120, smoke = as.factor(1))
phat_smoker <- predict(fit, newdata = newdata1, type = "response")
phat_smoker## 1
## 0.4566911newdata1b <- data.frame(age = 27, lwt = 120, smoke = as.factor(0))
phat_nonsmoker <- predict(fit, newdata = newdata1b, type = "response")
phat_nonsmoker## 1
## 0.08628697#6. predicted-probability curve across maternal age, by smoking status
newdata2 <- with(dat1,
data.frame(
age = rep(seq(from=min(age),to=max(age),length.out=100),2),
lwt = mean(lwt),
smoke = as.factor(rep(0:1,each=100)))
)
PredictedProb <- predict(fit, newdata = newdata2, type = "response")
dtf_low <- cbind.data.frame(newdata2, data.frame(PredictedProb = PredictedProb))
ggplot(dtf_low, aes(x = age, y = PredictedProb)) +
geom_line(aes(colour = smoke), linewidth = 1) +
labs(title = "Predicted probability of low birth weight by maternal age and smoking",
x = "Maternal age", y = "Predicted P(low birth weight)")
#7. classification table
phat <- predict(fit, newdata=dat1, type = "response")
yhat <- ifelse(phat > 0.5, 1, 0)
xtabs(~ dat1$LOW + yhat)## yhat
## dat1$LOW 0 1
## 0 11 3
## 1 3 13Reading these numbers: with age,
lwt, and smoke all in the model, the odds
ratios are age ≈ 1.23 (p = 0.068,
borderline), lwt ≈ 1.01 (p = 0.48, not
significant), and smoke ≈ 8.90 (p = 0.030,
significant). The standout finding is smoking: holding maternal age
and weight fixed, smoking during pregnancy multiplies the odds of a low
birth weight baby by roughly 9-fold - a striking effect, consistent
with the well-established public-health link between smoking and fetal
growth restriction. This shows up clearly in the predicted probabilities
too: for a 27-year-old, 120lb mother, the model predicts about a
46% chance of low birth weight if she smokes, versus
only about 9% if she doesn’t - holding everything else
about her identical. lwt (mother’s pre-pregnancy weight)
turns out not to add much once age and smoking are accounted for, at
least in this small sample.
Lab 2 Assignment
Using the BuyingCosmetics.csv (Task 3),
Admission.csv (Task 2), or NewbornWeight.csv
(Task 4) dataset:
- Fit the full logistic model and select variables using stepwise
selection. Report the final model equation on the logit scale.
- Interpret the odds ratio of each significant predictor in a sentence
(e.g. “holding X constant, a one-unit increase in … multiplies the odds
of … by …”).
- Report the 95% CI for each odds ratio and state whether each
predictor is significant at the 5% level.
- Perform the Hosmer-Lemeshow test and state your conclusion about
goodness of fit.
- Construct the classification table at a cutoff of 0.5 and compute
the overall accuracy, sensitivity and specificity.
- Check for influential observations using Cook’s distance. Does
removing the most influential point change any conclusion?
- Check for multicollinearity among predictors using VIF. Are there
any predictors you would consider dropping on this basis alone?
Lab 3: Multinomial logistic regression
Objectives: Students are able to use R to analyse a
nominal (unordered) categorical response variable with more than two
categories:
1. set an appropriate reference (base) category for the response
2. fit a multinomial logistic regression model and obtain z-tests for
coefficients
3. perform variable selection using likelihood ratio tests and stepwise
AIC
4. interpret fitted coefficients as (relative-risk) odds ratios
comparing each category to the base category
5. evaluate model fit and classification accuracy
Example: Product Buying
What the model is
Multinomial logistic regression extends binary logistic regression to
a response with more than two unordered categories
(here, products A/B/C). It works by fitting \(k-1\) simultaneous logit equations, each
comparing one category to a chosen reference category: \[\ln\left(\frac{P(y = j)}{P(y =
\text{ref})}\right) = \beta_{0j} + \beta_{1j} x_1 + \dots, \quad j \neq
\text{ref}\] With product releveled so “A” is the
reference, summary(fit2) therefore reports two
sets of coefficients: one for “B vs. A” and one for “C vs. A”.
Each exponentiated coefficient is a relative risk ratio
(RRR) - the multiplicative change in the odds of choosing that
category over the reference category for a one-unit increase in
the predictor.
1.Import data
library(readr)
dt <- read_csv(file.path(data_dir, "product.csv"))
str(dt)## spc_tbl_ [400 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ obs : num [1:400] 1 2 3 4 5 6 7 8 9 10 ...
## $ product : chr [1:400] "C" "A" "C" "A" ...
## $ age : num [1:400] 57 21 66 36 23 31 37 37 55 66 ...
## $ household : num [1:400] 2 7 7 4 0 5 3 0 3 2 ...
## $ position_level: num [1:400] 2 2 2 2 2 1 3 3 3 4 ...
## $ gender : chr [1:400] "Male" "Male" "Male" "Female" ...
## $ absent : num [1:400] 10 7 1 6 11 14 12 25 3 18 ...
## - attr(*, "spec")=
## .. cols(
## .. obs = col_double(),
## .. product = col_character(),
## .. age = col_double(),
## .. household = col_double(),
## .. position_level = col_double(),
## .. gender = col_character(),
## .. absent = col_double()
## .. )
## - attr(*, "problems")=<pointer: 0x000002121dfff450>- set data structure
dt$product <- as.factor(dt$product)
dt$position_level <- as.factor(dt$position_level)
dt$gender <- as.factor(dt$gender)
str(dt)## spc_tbl_ [400 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ obs : num [1:400] 1 2 3 4 5 6 7 8 9 10 ...
## $ product : Factor w/ 3 levels "A","B","C": 3 1 3 1 1 1 1 2 3 2 ...
## $ age : num [1:400] 57 21 66 36 23 31 37 37 55 66 ...
## $ household : num [1:400] 2 7 7 4 0 5 3 0 3 2 ...
## $ position_level: Factor w/ 5 levels "1","2","3","4",..: 2 2 2 2 2 1 3 3 3 4 ...
## $ gender : Factor w/ 2 levels "Female","Male": 2 2 2 1 2 2 2 1 1 1 ...
## $ absent : num [1:400] 10 7 1 6 11 14 12 25 3 18 ...
## - attr(*, "spec")=
## .. cols(
## .. obs = col_double(),
## .. product = col_character(),
## .. age = col_double(),
## .. household = col_double(),
## .. position_level = col_double(),
## .. gender = col_character(),
## .. absent = col_double()
## .. )
## - attr(*, "problems")=<pointer: 0x000002121dfff450>- descriptive statistics
summary(dt)## obs product age household position_level
## Min. : 1.0 A:142 Min. :21.00 Min. :0.000 1: 47
## 1st Qu.:100.8 B:117 1st Qu.:30.00 1st Qu.:2.000 2: 99
## Median :200.5 C:141 Median :37.00 Median :3.000 3:135
## Mean :200.5 Mean :41.05 Mean :3.212 4: 64
## 3rd Qu.:300.2 3rd Qu.:53.00 3rd Qu.:5.000 5: 55
## Max. :400.0 Max. :67.00 Max. :7.000
## gender absent
## Female:232 Min. : 0.00
## Male :168 1st Qu.: 7.00
## Median :14.00
## Mean :14.17
## 3rd Qu.:21.00
## Max. :31.00boxplot(age~product,data=dt)
table(dt$product,dt$position_level)##
## 1 2 3 4 5
## A 14 36 44 28 20
## B 20 29 39 13 16
## C 13 34 52 23 19plot(dt$age,dt$absent)
4.set base group for y variable
dt$product <- relevel(dt$product, ref = "A")- fit multinomial logistic regression
#Important packages
library(MASS)# for multinom
library(rms) # for lrt
library(pscl) #to get pseudo R2
#package "nnet"
library(nnet)
null <- multinom(product~1,data=dt)## # weights: 6 (2 variable)
## initial value 439.444915
## final value 437.908882
## convergedfull <- multinom(product~age+household+position_level+gender+absent,data=dt)## # weights: 30 (18 variable)
## initial value 439.444915
## iter 10 value 268.455865
## iter 20 value 208.024202
## iter 30 value 207.509074
## iter 30 value 207.509072
## iter 30 value 207.509072
## final value 207.509072
## convergedsummary(full)## Call:
## multinom(formula = product ~ age + household + position_level +
## gender + absent, data = dt)
##
## Coefficients:
## (Intercept) age household position_level2 position_level3
## B -4.727535 0.2368700 -0.9628612 -0.34360355 -1.3305304
## C -10.114127 0.2647383 0.1552594 0.01216742 -0.1824063
## position_level4 position_level5 genderMale absent
## B -1.9177254 -1.9734897 -2.39588762 0.01384059
## C -0.2653796 -0.7251654 0.09274777 -0.01208771
##
## Std. Errors:
## (Intercept) age household position_level2 position_level3
## B 1.027988 0.02864214 0.13731517 0.5999365 0.5954285
## C 1.213453 0.02960232 0.09266699 0.6027884 0.5777555
## position_level4 position_level5 genderMale absent
## B 0.7839738 0.7471631 0.4633714 0.02489143
## C 0.6440092 0.7116233 0.3694506 0.02293563
##
## Residual Deviance: 415.0181
## AIC: 451.0181- Calculate z value=b/se(b) , then p-value
b = summary(full)$coefficients
se.b = summary(full)$standard.errors
z <- b/se.b
t(z)#transpose of Z## B C
## (Intercept) -4.5988214 -8.33500061
## age 8.2699833 8.94315942
## household -7.0120526 1.67545552
## position_level2 -0.5727332 0.02018522
## position_level3 -2.2345763 -0.31571534
## position_level4 -2.4461600 -0.41207432
## position_level5 -2.6413106 -1.01902993
## genderMale -5.1705554 0.25104244
## absent 0.5560382 -0.52702743# p-value for 2-tailed z test
p <- (1 - pnorm(abs(z),mean=0,sd=1))*2
t(p)#transpose of p## B C
## (Intercept) 4.248879e-06 0.00000000
## age 2.220446e-16 0.00000000
## household 2.348566e-12 0.09384489
## position_level2 5.668254e-01 0.98389562
## position_level3 2.544518e-02 0.75221859
## position_level4 1.443869e-02 0.68028495
## position_level5 8.258598e-03 0.30818876
## genderMale 2.333993e-07 0.80178130
## absent 5.781847e-01 0.59817454- test if position_level should be included when the rest are in!
#likelihood ratio test
full <- multinom(product~age+household+position_level+gender+absent,data=dt)## # weights: 30 (18 variable)
## initial value 439.444915
## iter 10 value 268.455865
## iter 20 value 208.024202
## iter 30 value 207.509074
## iter 30 value 207.509072
## iter 30 value 207.509072
## final value 207.509072
## convergedfit1 <- multinom(product~age+household+gender+absent,data=dt)## # weights: 18 (10 variable)
## initial value 439.444915
## iter 10 value 256.830887
## iter 20 value 214.842594
## iter 30 value 214.837090
## final value 214.837084
## convergedanova(fit1,full)- from
fit1, test if “absent” should be included when the rest are in
fit2 <- multinom(product~age+household+gender,data=dt)## # weights: 15 (8 variable)
## initial value 439.444915
## iter 10 value 220.300742
## iter 20 value 215.255121
## final value 215.248449
## convergedanova(fit2,fit1)- stepAIC for variable selection in
glm
stepAIC(full,direction="both")## Start: AIC=451.02
## product ~ age + household + position_level + gender + absent
##
## # weights: 27 (16 variable)
## initial value 439.444915
## iter 10 value 342.197670
## iter 20 value 334.868396
## final value 334.851268
## converged
## # weights: 27 (16 variable)
## initial value 439.444915
## iter 10 value 301.317069
## iter 20 value 276.620201
## final value 276.458307
## converged
## # weights: 18 (10 variable)
## initial value 439.444915
## iter 10 value 256.830887
## iter 20 value 214.842594
## iter 30 value 214.837090
## final value 214.837084
## converged
## # weights: 27 (16 variable)
## initial value 439.444915
## iter 10 value 260.954489
## iter 20 value 230.815041
## final value 230.581071
## converged
## # weights: 27 (16 variable)
## initial value 439.444915
## iter 10 value 229.652679
## iter 20 value 208.237597
## final value 208.149341
## converged
## Df AIC
## - absent 2 448.30
## - position_level 8 449.67
## <none> 451.02
## - gender 2 493.16
## - household 2 584.92
## - age 2 701.70
## # weights: 27 (16 variable)
## initial value 439.444915
## iter 10 value 229.652679
## iter 20 value 208.237597
## final value 208.149341
## converged
##
## Step: AIC=448.3
## product ~ age + household + position_level + gender
##
## # weights: 24 (14 variable)
## initial value 439.444915
## iter 10 value 345.161008
## iter 20 value 335.869815
## final value 335.869805
## converged
## # weights: 24 (14 variable)
## initial value 439.444915
## iter 10 value 286.695124
## iter 20 value 277.917014
## final value 277.916928
## converged
## # weights: 15 (8 variable)
## initial value 439.444915
## iter 10 value 220.300742
## iter 20 value 215.255121
## final value 215.248449
## converged
## # weights: 24 (14 variable)
## initial value 439.444915
## iter 10 value 252.565848
## iter 20 value 231.441304
## final value 231.409442
## converged
## # weights: 30 (18 variable)
## initial value 439.444915
## iter 10 value 268.455865
## iter 20 value 208.024202
## iter 30 value 207.509074
## iter 30 value 207.509072
## iter 30 value 207.509072
## final value 207.509072
## converged
## Df AIC
## - position_level 8 446.50
## <none> 448.30
## + absent 2 451.02
## - gender 2 490.82
## - household 2 583.83
## - age 2 699.74
## # weights: 15 (8 variable)
## initial value 439.444915
## iter 10 value 220.300742
## iter 20 value 215.255121
## final value 215.248449
## converged
##
## Step: AIC=446.5
## product ~ age + household + gender
##
## # weights: 12 (6 variable)
## initial value 439.444915
## iter 10 value 341.114844
## final value 340.866118
## converged
## # weights: 12 (6 variable)
## initial value 439.444915
## iter 10 value 284.455247
## iter 20 value 283.954925
## final value 283.954920
## converged
## # weights: 12 (6 variable)
## initial value 439.444915
## iter 10 value 239.284726
## iter 20 value 238.014588
## final value 238.014521
## converged
## # weights: 27 (16 variable)
## initial value 439.444915
## iter 10 value 229.652679
## iter 20 value 208.237597
## final value 208.149341
## converged
## # weights: 18 (10 variable)
## initial value 439.444915
## iter 10 value 256.830887
## iter 20 value 214.842594
## iter 30 value 214.837090
## final value 214.837084
## converged
## Df AIC
## <none> 446.50
## + position_level 8 448.30
## + absent 2 449.67
## - gender 2 488.03
## - household 2 579.91
## - age 2 693.73## Call:
## multinom(formula = product ~ age + household + gender, data = dt)
##
## Coefficients:
## (Intercept) age household genderMale
## B -5.229876 0.2247682 -0.9143186 -2.2522987
## C -10.339791 0.2608011 0.1518043 0.1087794
##
## Residual Deviance: 430.4969
## AIC: 446.4969- overall fit test
anova(null,fit2)- obtain the fitted equations
summary(fit2)## Call:
## multinom(formula = product ~ age + household + gender, data = dt)
##
## Coefficients:
## (Intercept) age household genderMale
## B -5.229876 0.2247682 -0.9143186 -2.2522987
## C -10.339791 0.2608011 0.1518043 0.1087794
##
## Std. Errors:
## (Intercept) age household genderMale
## B 0.8699672 0.02746372 0.12695214 0.4347157
## C 1.0956022 0.02854696 0.09074394 0.3624084
##
## Residual Deviance: 430.4969
## AIC: 446.4969- get odd ratios
exp(summary(fit2)$coefficients)## (Intercept) age household genderMale
## B 5.354189e-03 1.252032 0.4007896 0.1051572
## C 3.232109e-05 1.297969 1.1639324 1.1149163Reading these numbers: fitting
product ~ age + household + gender gives relative risk
ratios of roughly - for product B vs. A: age ≈
1.25 and household ≈
0.40; for product C vs. A: age ≈
1.30 and household ≈
1.16. In words: each additional year of age
multiplies the odds of choosing product B (over A) by about 1.25, and
the odds of choosing product C (over A) by about 1.30 - age pushes
customers away from A toward either alternative. household
size, however, tells a different story for each comparison: more
household members reduce the relative odds of B vs. A (RRR <
1) but increase the relative odds of C vs. A (RRR > 1) -
household size distinguishes B-buyers from C-buyers even though both
differ from A in the same direction on age. This is the key advantage
(and complexity) of multinomial models: each predictor can have
a different relationship with each non-reference category, so
always interpret coefficients “vs. the reference” rather than as one
global effect.
- Get the pseudo-R2
library(pscl)
pR2(fit2)## fitting null model for pseudo-r2
## # weights: 6 (2 variable)
## initial value 439.444915
## final value 437.908882
## converged## llh llhNull G2 McFadden r2ML r2CU
## -215.2484490 -437.9088817 445.3208655 0.5084629 0.6715275 0.7561972- Get \(p_{hat}\) and \(y_{hat}\)
phat <- predict(fit2,newdata=dt,type="probs")
head(phat)#show the first 6 rows## A B C
## 1 0.0057578835 0.1908607399 0.80338138
## 2 0.9755735475 0.0001023612 0.02432409
## 3 0.0003204954 0.0008306117 0.99884889
## 4 0.4628612448 0.2089157798 0.32822298
## 5 0.8980469293 0.0889180955 0.01303498
## 6 0.7961988323 0.0049228842 0.19887828yhat <- predict(fit2,newdata=dt)
head(yhat)## [1] C A C A A A
## Levels: A B C15.Classification table and %accuracy
tab_class <- table(dt$product,yhat)
tab_class## yhat
## A B C
## A 118 15 9
## B 28 79 10
## C 16 20 105# overall accuracy = sum of diagonal / total
accuracy <- sum(diag(tab_class))/sum(tab_class)
accuracy## [1] 0.755Lab 3 Assignment
Using the product.csv dataset:
- Refit
fit2and write out the two fitted logit equations (each non-reference category vs. the reference category “A”).
- For one significant predictor, interpret its coefficient as a
relative-risk ratio (odds ratio) for each non-reference category.
- Using
anova(), test whethergendercan be dropped fromfit2. State your conclusion.
- Report the pseudo-\(R^2\) and the
overall classification accuracy. Comment on how well the model
performs.
- Which product category is hardest for the model to classify
correctly? Look at the off-diagonal cells of the classification table to
justify your answer.
Lab 4: Ordinal logistic regression
Objectives: Students are able to use R to analyse an
ordinal (ordered categorical) response variable:
1. declare the response as an ordered factor with the correct level
ordering
2. fit a proportional-odds (ordinal logistic) model using
polr()
3. perform inference on parameters and select variables using likelihood
ratio tests and stepwise AIC
4. check the proportional-odds (parallel lines) assumption
5. interpret odds ratios and use the model for
prediction/classification
Example: Attitude towards welfare spending (World Values Survey)
The poverty.csv dataset (drawn from the World Values
Survey, also available as carData::WVS) records each
respondent’s opinion on whether government is doing “Too Little”, “About
Right”, or “Too Much” to reduce poverty (an ordered response), along
with religion, degree (has a college degree),
country, age, and gender.
What the model is
The proportional-odds (ordinal logistic) model treats the three ordered categories as arising from a single continuous latent attitude, cut into groups by two thresholds. It estimates one set of predictor coefficients shared across both thresholds, plus threshold-specific intercepts: \[\ln\left(\frac{P(y \le j)}{P(y > j)}\right) = \alpha_j - (\beta_1 x_1 + \beta_2 x_2 + \dots), \quad j = 1, 2\] The same \(\beta\) applies at both cut-points (that’s the “proportional odds” assumption, checked later with the Brant test) - only the threshold \(\alpha_j\) differs. Exponentiated coefficients are cumulative odds ratios: the multiplicative change in the odds of being in a higher category (e.g. “Too Much” rather than “About Right” or “Too Little”) for a one-unit increase in the predictor.
import data
library(readr)
dp <- read_csv(file.path(data_dir, "poverty.csv"))- set data structure
str(dp)## spc_tbl_ [5,381 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ...1 : num [1:5381] 1 2 3 4 5 6 7 8 9 10 ...
## $ poverty : chr [1:5381] "Too Little" "About Right" "Too Little" "Too Much" ...
## $ religion: chr [1:5381] "yes" "yes" "yes" "yes" ...
## $ degree : chr [1:5381] "no" "no" "no" "yes" ...
## $ country : chr [1:5381] "USA" "USA" "USA" "USA" ...
## $ age : num [1:5381] 44 40 36 25 39 80 48 32 74 30 ...
## $ gender : chr [1:5381] "male" "female" "female" "female" ...
## - attr(*, "spec")=
## .. cols(
## .. ...1 = col_double(),
## .. poverty = col_character(),
## .. religion = col_character(),
## .. degree = col_character(),
## .. country = col_character(),
## .. age = col_double(),
## .. gender = col_character()
## .. )
## - attr(*, "problems")=<pointer: 0x000002121dfffdb0>dp$poverty <- as.factor(dp$poverty)
dp$religion <- as.factor(dp$religion)
dp$degree <- as.factor(dp$degree)
dp$country <- as.factor(dp$country)
dp$gender <- as.factor(dp$gender)- descriptive statistics
dp$poverty <- ordered(dp$poverty,
levels =c("Too Little","About Right","Too Much"))
levels(dp$poverty)## [1] "Too Little" "About Right" "Too Much"ggplot(dp, aes(x = poverty, y = age, fill = poverty)) +
geom_boxplot(size = .75) +
facet_grid(country ~ gender, margins = FALSE) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
labs(title = "Age distribution by poverty-attitude category, country and gender")
#3.m0: null model
library(MASS)
m0 <- polr(poverty~1,data=dp)
#4. mf: full model
mf <- polr(poverty~religion+degree+country+age+gender,data=dp)
summary(mf)## Call:
## polr(formula = poverty ~ religion + degree + country + age +
## gender, data = dp)
##
## Coefficients:
## Value Std. Error t value
## religionyes 0.17973 0.077346 2.324
## degreeyes 0.14092 0.066193 2.129
## countryNorway -0.32235 0.073766 -4.370
## countrySweden -0.60330 0.079494 -7.589
## countryUSA 0.61777 0.070665 8.742
## age 0.01114 0.001561 7.139
## gendermale 0.17637 0.052972 3.329
##
## Intercepts:
## Value Std. Error t value
## Too Little|About Right 0.7298 0.1041 7.0128
## About Right|Too Much 2.5325 0.1103 22.9496
##
## Residual Deviance: 10402.59
## AIC: 10420.59# 5.calculate p-value
## store table
ctable <- coef(summary(mf))
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE)*2
## combined table
ctable <- cbind(ctable, "p value" = p)
ctable## Value Std. Error t value p value
## religionyes 0.17973194 0.077346025 2.323738 2.013951e-02
## degreeyes 0.14091745 0.066193093 2.128885 3.326381e-02
## countryNorway -0.32235359 0.073766012 -4.369947 1.242765e-05
## countrySweden -0.60329785 0.079493879 -7.589237 3.217957e-14
## countryUSA 0.61777260 0.070664754 8.742302 2.284084e-18
## age 0.01114091 0.001560585 7.138935 9.405696e-13
## gendermale 0.17636863 0.052972239 3.329454 8.701647e-04
## Too Little|About Right 0.72976353 0.104061619 7.012802 2.335919e-12
## About Right|Too Much 2.53247870 0.110349763 22.949562 1.488394e-116# 6. answer is age
# 7.Test if "age" should be included in the model when others are in
m1 <- polr(poverty~religion+degree+country+gender,data=dp)
#likelihood ratio test, use anova function
anova(mf,m1)## store table
ctable <- coef(summary(m1))
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE)*2
## combined table
ctable <- cbind(ctable, "p value" = p)
ctable## Value Std. Error t value p value
## religionyes 0.23479820 0.07679364 3.057521 2.231758e-03
## degreeyes 0.09304962 0.06569864 1.416310 1.566848e-01
## countryNorway -0.32957463 0.07354477 -4.481279 7.419704e-06
## countrySweden -0.60346856 0.07923965 -7.615740 2.621848e-14
## countryUSA 0.66319466 0.07033120 9.429593 4.116813e-21
## gendermale 0.18646926 0.05284173 3.528826 4.174078e-04
## Too Little|About Right 0.28016816 0.08232821 3.403064 6.663467e-04
## About Right|Too Much 2.06918584 0.08824151 23.449121 1.349496e-121#8. From m1, test the significance of "degree"
m2 <- polr(poverty~religion+country+gender,data=dp)
anova(m1,m2)#Variable selection for logistic
library(MASS)
stepAIC(mf,direction="both")#use stepwise## Start: AIC=10420.59
## poverty ~ religion + degree + country + age + gender
##
## Df AIC
## <none> 10421
## - degree 1 10423
## - religion 1 10424
## - gender 1 10430
## - age 1 10470
## - country 3 10666## Call:
## polr(formula = poverty ~ religion + degree + country + age +
## gender, data = dp)
##
## Coefficients:
## religionyes degreeyes countryNorway countrySweden countryUSA
## 0.17973194 0.14091745 -0.32235359 -0.60329785 0.61777260
## age gendermale
## 0.01114091 0.17636863
##
## Intercepts:
## Too Little|About Right About Right|Too Much
## 0.7297635 2.5324787
##
## Residual Deviance: 10402.59
## AIC: 10420.59#9.
anova(m0,mf)#11. OR and its confidence interval (CI)
ci <- confint(mf) # 95% CI of log odd
#odd
exp(coef(mf))#odds## religionyes degreeyes countryNorway countrySweden countryUSA
## 1.1968965 1.1513296 0.7244420 0.5470047 1.8547921
## age gendermale
## 1.0112032 1.1928777exp(cbind(OR = coef(mf), ci))#95% CI of odds## OR 2.5 % 97.5 %
## religionyes 1.1968965 1.0289511 1.3934593
## degreeyes 1.1513296 1.0110549 1.3106194
## countryNorway 0.7244420 0.6267512 0.8369294
## countrySweden 0.5470047 0.4678450 0.6389299
## countryUSA 1.8547921 1.6151018 2.1306390
## age 1.0112032 1.0081175 1.0143029
## gendermale 1.1928777 1.0752794 1.3234479Reading these numbers: fitting the full model on
poverty.csv gives odds ratios of roughly
religion (yes) ≈ 1.20,
degree (yes) ≈ 1.15, age ≈
1.01, gender (male) ≈
1.19, and for country (relative to
Australia): Norway ≈ 0.72, Sweden ≈
0.55, USA ≈ 1.86. In words:
holding the other variables fixed, being male multiplies the odds of
holding a “government is doing too much” (vs. “too little”/“about
right”) view by about 1.19, and each additional year of age multiplies
those odds by 1.01 (small per-year, but compounds over a lifetime).
The country effects are the most striking: Americans have nearly
1.9 times the odds of Australians of saying government
does “too much” on poverty, while Swedes have only about 0.55
times (i.e. Swedes lean much more toward “too little”) - a
pattern consistent with the two countries’ different welfare-state
traditions. Because this is a cumulative odds ratio, “higher
category” means further toward “Too Much” on the ordered scale in every
one of these comparisons simultaneously (assuming the proportional-odds
assumption below actually holds).
#12.pseudo R2
library(pscl) #to get pseudo R2
pR2(mf)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML
## -5.201296e+03 -5.370188e+03 3.377841e+02 3.144993e-02 6.084382e-02
## r2CU
## 7.041132e-02#13. getting phat
phat <- predict(mf,newdata=dp,type="probs")
head(phat)## Too Little About Right Too Much
## 1 0.3242497 0.4200439 0.2557065
## 2 0.3744021 0.4096330 0.2159649
## 3 0.3848970 0.4065882 0.2085148
## 4 0.3805578 0.4078799 0.2115623
## 5 0.3058650 0.4218762 0.2722588
## 6 0.2770756 0.4221685 0.3007559#getting yhat
yhat <- predict(mf,newdata=dp)
head(yhat)## [1] About Right About Right About Right About Right About Right About Right
## Levels: Too Little About Right Too Much#14. classification table
tab_class <- table(dp$poverty,yhat)
tab_class## yhat
## Too Little About Right Too Much
## Too Little 2190 518 0
## About Right 1483 379 0
## Too Much 377 434 0accuracy <- sum(diag(tab_class))/sum(tab_class)
accuracy## [1] 0.4774206Visualizing predicted probabilities across age (by
degree)
n <- nrow(dp)
newdata <- data.frame(religion = rep(dp$religion, 3),
degree = rep(dp$degree, 3),
country = rep(dp$country, 3),
age = rep(dp$age, 3),
gender = rep(dp$gender, 3),
group = c(rep("Too Little", n),
rep("About Right", n),
rep("Too Much", n)),
phat = c(phat[, "Too Little"],
phat[, "About Right"],
phat[, "Too Much"])
)
ggplot(newdata, aes(x = age, y = phat, color = group)) +
geom_smooth(se = FALSE) +
facet_grid(degree ~ ., scales = "free") +
labs(title = "Predicted probability of each poverty-attitude category by age and degree",
x = "Age", y = "Predicted probability")
Checking the proportional-odds assumption
The polr() model assumes the effect of each predictor is
the same across all cut-points of the ordinal response (the “parallel
lines” or “proportional odds” assumption). We can check this with the
Brant test.
library(brant)
# NOTE: brant() requires the response coded as an ordered factor, which dp$poverty already is
brant(mf)## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 261.22 7 0
## religionyes 6.07 1 0.01
## degreeyes 2.08 1 0.15
## countryNorway 83.64 1 0
## countrySweden 59.58 1 0
## countryUSA 43.85 1 0
## age 0.7 1 0.4
## gendermale 1.69 1 0.19
## --------------------------------------------
##
## H0: Parallel Regression Assumption holdsInterpretation: a significant (p < 0.05)
“Omnibus” test, or a significant test for an individual predictor,
indicates that predictor’s effect differs across cut-points and the
proportional-odds assumption is violated for it. If this happens,
consider a partial proportional-odds model (e.g. vglm() in
package VGAM) or a multinomial logistic regression
instead.
Lab 4 Assignment
Using the poverty.csv dataset:
- Confirm that
povertyis coded as an ordered factor with the levels in the correct order. Why does the order matter here (but not in Lab 3’s multinomial model)?
- Fit the full model and use
stepAIC()to select a reduced model. Report the final model.
- Interpret the odds ratio for one significant predictor in terms of
the cumulative odds of being in a “higher” poverty-attitude
category.
- Run the Brant test on your final model. Does the proportional-odds
assumption hold? If not, what would you do next?
- Report the classification accuracy and comment on whether the model
is more successful for some categories of
povertythan others.
Lab 5: Poisson regression
Objectives: Students are able to use R to analyse a
count response variable using Poisson regression:
1. fit a Poisson regression model with a log link and interpret
coefficients as rate ratios
2. perform inference and variable selection using likelihood ratio
tests
3. assess model fit using pseudo-\(R^2\) and residual analysis (response and
deviance residuals)
4. test for over-dispersion and understand its consequences for standard
errors
5. use an offset term when observations have different exposure (time,
area, population at risk)
What the model is
Poisson regression models a count response (number of events) using a log link, which guarantees predicted counts stay non-negative: \[\ln(\mu) = \beta_0 + \beta_1 x, \quad y \sim \text{Poisson}(\mu)\] Exponentiating a coefficient gives a rate ratio (RR): the multiplicative change in the expected count for a one-unit increase in the predictor. A key structural assumption of the Poisson distribution is that the mean equals the variance - when the data are more spread out than that (over-dispersion), standard errors from a plain Poisson model are too small, making predictors look more significant than they really are (this is exactly what Lab 6’s Negative Binomial model fixes).
#———————————————————-
Example1: Number of Deaths Due to AIDs
#———————————————————-
ขั้นตอนการวิเคราะห์
1. ดาวน์โหลดข้อมูล โดยให้ข้อมูลอยู่ในรูป dataframe ชื่อ dat1
library("readr")
dat1 <- read_csv(file.path(data_dir, "NumberDeathsAIDs.csv"))
dat1 <- data.frame(dat1)
head(dat1)
str(dat1)2.วิเคราะห์ข้อมูลเบื้องต้นด้วย สถิติเชิงพรรณา
summary(dat1)## y x
## Min. : 0.00 Min. :0.000
## 1st Qu.: 2.25 1st Qu.:1.442
## Median :13.50 Median :2.013
## Mean :15.64 Mean :1.799
## 3rd Qu.:24.50 3rd Qu.:2.374
## Max. :45.00 Max. :2.639library("ggplot2")
g1<-ggplot(dat1,aes(x=x,y=y))+
geom_point()
g2<-ggplot(dat1,aes(x=y))+
geom_histogram(binwidth=5)
library("gridExtra")
grid.arrange(g1,g2,nrow=1)
3.สมมติ ทำการพิจารณาความสัมพันธ์ระหว่าง x และ y ด้วยการวิเคราะห์การถดถอยเชิงเส้นอย่างง่าย (Simple linear regression analysis)
fit.lm <- lm(y~x,data=dat1)
summary(fit.lm)##
## Call:
## lm(formula = y ~ x, data = dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.659 -5.877 -1.726 6.273 16.168
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -12.621 6.362 -1.984 0.070624 .
## x 15.708 3.266 4.809 0.000427 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.116 on 12 degrees of freedom
## Multiple R-squared: 0.6584, Adjusted R-squared: 0.6299
## F-statistic: 23.13 on 1 and 12 DF, p-value: 0.0004267pred.lm <- predict(fit.lm)
dat1$pred.lm <- pred.lm
ggplot(dat1,aes(x=x,y=y))+
geom_point()+
geom_line(aes(x=x,y=pred.lm),color="blue")
4. ทำการวิเคราะห์ความสัมพันธ์ระหว่าง x และ y ด้วยการวิเคราะห์การถดถอยปัวซง
(Poisson regression analysis)
#สร้างตัวแบบ null และตัวแบบที่มีตัวแปรอิสระ
library("MASS") #package for polr function
null <- glm(y~1, data=dat1, family=poisson(link="log"))
fit.pr <- glm(y~x, data=dat1, family=poisson(link="log"))
summary(fit.pr)##
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"), data = dat1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.9442 0.5116 -3.80 0.000145 ***
## x 2.1748 0.2150 10.11 < 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: 207.272 on 13 degrees of freedom
## Residual deviance: 17.092 on 12 degrees of freedom
## AIC: 74.019
##
## Number of Fisher Scoring iterations: 4#ทดสอบว่าตัวแปรอิสระสมควรอยู่ในตัวแบบหรือไม่ ด้วยการใช้ Lr test หรือ Deviance test
anova(null,fit.pr)library("lmtest")
lrtest(null,fit.pr)#หา ช่วงความเชื่อมั่นของค่า beta และค่า exp(beta)
confint(fit.pr)## 2.5 % 97.5 %
## (Intercept) -2.999299 -0.9925442
## x 1.771631 2.6151271exp(confint(fit.pr))## 2.5 % 97.5 %
## (Intercept) 0.04982198 0.3706325
## x 5.88043832 13.6689531Reading these numbers: fitting y ~ x on
the AIDS deaths data gives a coefficient on x of about
2.17 (highly significant), which exponentiates to a
rate ratio of about 8.8. Since x here is
already on a transformed (log-time) scale, this says each one-unit
increase in x multiplies the expected number of deaths by roughly
8.8 - a dramatic rate of increase, consistent with an epidemic in
its early exponential-growth phase. Compare this fit to the simple
linear model (fit.lm) from the previous step: Poisson
naturally enforces non-negative, curved (multiplicative) growth in the
predicted counts, which is usually a much better match for
epidemic-style count data than a straight line.
#get pseudo R2
library("pscl")
pR2(fit.pr)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -35.0094212 -130.0997464 190.1806504 0.7309032 0.9999987 0.9999987#ตรวจสอบ over-dispersion
library("AER")
dispersiontest(fit.pr,trafo=1)##
## Overdispersion test
##
## data: fit.pr
## z = 0.53497, p-value = 0.2963
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 0.1956112#คำนวณค่า พยากรณ์ของ mu
pred.pr <- predict(fit.pr, newdata=dat1, type = "response")
dat1$pred.pr <- pred.pr5.วิเคราะห์ค่าส่วนเหลือของตัวแบบการถดถอยปัวซง
#residuals analysis
dat1$res.pr <- residuals(fit.pr,type="response")
ggplot(dat1,aes(x=pred.pr,y=res.pr))+
geom_point()
#compute deviance residuals
poisson.dev <- function (y, mu) {
# unit deviance
2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
}
expected <- dat1$y
estimates <- dat1$pred.pr
dev.res <- sqrt(poisson.dev(expected, estimates)) *
ifelse(expected > estimates, 1, -1)
dat1$dev.res <- dev.res
ggplot(dat1,aes(y=dev.res,x=pred.pr))+
geom_point()
6.เปรียบเทียบค่า mean square error ของทั้งสองตัวแบบ
library("Metrics")
with(dat1,rmse(y,pred.lm))
with(dat1,rmse(y,pred.pr))#—————————————
Example 2: Number of awards
#—————————————
1. Download dataset และให้อยู่ในรูปแบบ data.frame ชื่อ dat2
dataset <- read.csv(file.path(data_dir, "numberawards.csv"))
dat2 <- dataset
str(dat2)## 'data.frame': 200 obs. of 4 variables:
## $ id : int 45 108 15 67 153 51 164 133 2 53 ...
## $ num_awards: int 0 0 0 0 0 0 0 0 0 0 ...
## $ prog : int 3 1 3 3 3 1 3 3 3 3 ...
## $ math : int 41 41 44 42 40 42 46 40 33 46 ...dat2$prog <- factor(dat2$prog,
levels = 1:3,
labels = c("General", "Academic", "Vocational")
)
str(dat2)## 'data.frame': 200 obs. of 4 variables:
## $ id : int 45 108 15 67 153 51 164 133 2 53 ...
## $ num_awards: int 0 0 0 0 0 0 0 0 0 0 ...
## $ prog : Factor w/ 3 levels "General","Academic",..: 3 1 3 3 3 1 3 3 3 3 ...
## $ math : int 41 41 44 42 40 42 46 40 33 46 ...แสดงการวิเคราะห์ข้อมูลเบื้องต้นด้วยสถิติเชิงพรรณา พร้อมทั้งอธิบายแผนภาพหรือตารางที่ได้
ทำการสร้างตัวแบบด้วยการถดถอยปัวซง ทีี่มีตัวแปรตามคือ num_awards และชุดของตัวแปรอิสระ ดังนี้
- ตัวแบบที่ไม่มีตัวแปรอิสระ
- ตัวแบบที่มีเฉพาะตัวแปรอิสระ prog
- ตัวแบบที่มีเฉพาะตัวแปรอิสระ math
- ตัวแบบที่มีตัวแปรอิสระ prog และ math
- ตัวแบบที่ไม่มีตัวแปรอิสระ
จากข้อ 3 ตัวแบบใดมีความเหมาะสมที่สุดเพราะเหตุใด ให้แสดงการทดสอบเพื่อสนับสนุนคำตอบ
จากตัวแบบที่ได้ในข้อ 4 จงเขียนสมการพยากรณ์ค่า mu พร้อมทั้งอธิบายความหมายของค่าที่ได้
ให้ทำการตรวจสอบความเหมาะสมของตัวแบบ
แสดงการวิเคราะห์ค่าส่วนเหลือของตัวแบบที่ได้
สร้างแผนภาพเพื่อแสดงค่าพยากรณ์ของ mu จากคะแนนสอบ math จำแนกตาม prog
m0 <- glm(num_awards ~ 1, family=poisson(link="log"), data = dat2)
m1 <- glm(num_awards ~ prog, family=poisson(link="log"), data = dat2)
m2 <- glm(num_awards ~ math, family=poisson(link="log"), data = dat2)
m3 <- glm(num_awards ~ math+prog, family=poisson(link="log"), data = dat2)
lrtest(m0,m1)lrtest(m0,m2)lrtest(m1,m3)lrtest(m2,m3)fit <- m3
summary(fit)##
## Call:
## glm(formula = num_awards ~ math + prog, family = poisson(link = "log"),
## data = dat2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.24712 0.65845 -7.969 1.60e-15 ***
## math 0.07015 0.01060 6.619 3.63e-11 ***
## progAcademic 1.08386 0.35825 3.025 0.00248 **
## progVocational 0.36981 0.44107 0.838 0.40179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 287.67 on 199 degrees of freedom
## Residual deviance: 189.45 on 196 degrees of freedom
## AIC: 373.5
##
## Number of Fisher Scoring iterations: 6#หา ช่วงความเชื่อมั่นของค่า beta และค่า exp(beta)
confint(fit)## 2.5 % 97.5 %
## (Intercept) -6.57106799 -3.98529805
## math 0.04949674 0.09107698
## progAcademic 0.43504394 1.85586906
## progVocational -0.49020399 1.26602301exp(confint(fit))## 2.5 % 97.5 %
## (Intercept) 0.001400301 0.0185869
## math 1.050742170 1.0953533
## progAcademic 1.545030949 6.3972554
## progVocational 0.612501440 3.5467192#get pseudo R2
library("pscl")
pR2(fit)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -182.7522516 -231.8635588 98.2226145 0.2118112 0.3880551 0.4304109#check over-dispersion
library(AER)
dispersiontest(fit,trafo=1)##
## Overdispersion test
##
## data: fit
## z = 0.53224, p-value = 0.2973
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 0.04725442#prediction
# NOTE: predict.glm() has no "data" argument - use "newdata" instead.
# (Writing data = ... here is silently ignored, so double-check argument names
# whenever a predict() call looks like it "worked" without errors.)
pred <- predict(m3, newdata = dat2, type = "response", se.fit = TRUE)
dataset$pred <- pred$fit
ggplot(dataset,aes(x=math,y=pred,color=prog))+
geom_line()+
geom_point(aes(x=math,y=num_awards,color=prog))
#compute deviance residuals
poisson.dev <- function (y, mu) {
# unit deviance
2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
}
expected <- dataset$num_awards
estimates <- dataset$pred
dev.res <- sqrt(poisson.dev(expected, estimates)) *
ifelse(expected > estimates, 1, -1)
dataset$dev.res <- dev.res
ggplot(dataset,aes(y=dev.res,x=pred))+
geom_point()
Reading these numbers: the final model
(num_awards ~ math + prog) gives rate ratios of
math ≈ 1.07 and, relative to the “General”
program, prog = Academic ≈ 2.96 and
prog = Vocational ≈ 1.45 (not significant,
p = 0.40). In words: each additional point on the math exam
multiplies the expected number of awards won by about 7%, and
students in the Academic program are expected to win roughly 3 times
as many awards as otherwise-identical students in the General
program, holding math score fixed. The lrtest()
results above confirm both math and prog
improve the model significantly over simpler alternatives (comparing
m1 vs m3 and m2 vs
m3), which is why the combined model m3 is
preferred. Because the math coefficient is on the log scale,
its effect compounds - a student scoring 20 points higher on the math
exam is predicted to win \(e^{0.070 \times 20}
\approx 4\) times as many awards, not just 20 times the per-point
effect.
Influence diagnostics for the Poisson model (Example 2)
As with logistic regression, individual observations can have a disproportionate effect on a Poisson fit. Cook’s distance and leverage generalize directly to GLMs.
cooksd <- cooks.distance(fit)
plot(cooksd, type = "h", main = "Cook's distance: awards model")
abline(h = 4/nrow(dat2), col = "red", lty = 2)
hatvalues_ <- hatvalues(fit)
plot(hatvalues_, type = "h", main = "Leverage: awards model")
abline(h = 2*mean(hatvalues_), col = "red", lty = 2)
# Rows flagged as influential
which(cooksd > 4/nrow(dat2))## 27 28 36 54 58 87 154 157 181 190 194 199
## 27 28 36 54 58 87 154 157 181 190 194 199Using an offset term when exposure varies
So far each row/unit has been observed for the same amount of
“exposure” (the same time period, the same population, the same area).
When units differ in exposure - e.g. number of claims per policy-year,
number of accidents per 10,000 vehicle-km, number of infections per
person-years at risk - we model the rate rather than
the raw count by adding offset(log(exposure)) to the linear
predictor. This constrains the coefficient on log(exposure)
to be exactly 1, which is the correct way to account for exposure (as
opposed to simply adding exposure as an ordinary
predictor).
# Example: dat_claims has columns "claims" (count), "policy_years" (exposure),
# and predictor "age_group"
fit.offset <- glm(claims ~ age_group + offset(log(policy_years)),
family = poisson(link = "log"), data = dat_claims)
summary(fit.offset)
# exp(coefficient) is now interpreted as a RATE RATIO (claims per policy-year),
# not a simple count ratio
exp(coef(fit.offset))Why this matters: if you fit
glm(claims ~ age_group, ...) without the offset, you are
implicitly assuming every policyholder was observed for the same length
of time - which is rarely true in insurance, epidemiology, or
reliability data. Omitting the offset when exposure varies biases the
comparison between groups.
A note on sample size for count models
Because Poisson (and later, Negative Binomial and zero-inflated)
models estimate more than just a mean and a variance, small samples can
produce unstable coefficient estimates and unreliable standard errors,
especially when:
- the total number of events (not just the number of rows) is small -
as a rule of thumb, aim for at least 10 events per predictor in the
model (similar in spirit to the “10 events per variable” guideline for
logistic regression);
- the data contain many zeros relative to the mean count (low average
count), which increases the risk of needing a zero-inflated model rather
than a plain Poisson/NB model - see Labs 7-8;
- you are fitting many predictors (or many
factorlevels) relative to the number of observations, in which case stepwise selection results become unstable and should be validated (e.g. via cross-validation or a held-out test set).
When in doubt, examine summary() for implausibly large
standard errors or coefficients, and consider simplifying the model or
collecting more data before trusting the inference.
Lab 5 Assignment
Using the NumberDeathsAIDs.csv (Example 1) or
numberawards.csv (Example 2) dataset:
- Fit the null and full Poisson models. Use the likelihood ratio test
to decide which predictors belong in the final model.
- Interpret
exp(coefficient)for one significant predictor as a rate ratio.
- Test for over-dispersion using
dispersiontest(). What would you do differently if the test were significant? (Preview: see Lab 6.)
- Plot response and deviance residuals against fitted values. Is there
any pattern suggesting lack of fit?
- Identify any influential observations using Cook’s distance, and
comment on whether removing them changes your conclusions.
- Suppose each observational unit had actually been observed for a
different length of time (an
exposurevariable). Rewrite your model formula to correctly account for this using an offset term.
Lab 6: Negative Binomial regression
Objectives: Students are able to use R to analyse
over-dispersed count data using Negative Binomial regression:
1. recognize over-dispersion in a Poisson model and understand its
consequences
2. fit a Negative Binomial regression model using
glm.nb()
3. compare Poisson and Negative Binomial models formally (likelihood
ratio test on the dispersion parameter)
4. perform variable selection and interpret coefficients as rate
ratios
5. check residuals and influential observations for the Negative
Binomial fit
Example: Number of articles published by the scientist
What the model is
The Negative Binomial distribution can be thought of as a Poisson
distribution where the rate \(\mu\)
itself varies randomly across individuals (a Poisson-Gamma mixture),
which lets the variance exceed the mean: \[y
\sim \text{NB}(\mu, \theta), \quad \text{Var}(y) = \mu +
\frac{\mu^2}{\theta}, \quad \ln(\mu) = \beta_0 + \beta_1 x_1 +
\dots\] The extra parameter \(\theta\) (sometimes written \(1/\alpha\)) controls how much more
spread-out the counts are than a Poisson would predict - as \(\theta \to \infty\), the model collapses
back to ordinary Poisson (Var = Mean). A small \(\theta\) means strong over-dispersion;
glm.nb() estimates \(\theta\) from the data alongside the usual
regression coefficients. Interpretation of exp(coefficient)
as a rate ratio works exactly the same way as in
Poisson regression - the only thing that changes is that the standard
errors (and therefore the p-values and confidence intervals) are now
correctly widened to reflect the true, greater variability in the
data.
The articles.csv dataset records the number of articles
(Articles) published by PhD biochemists in their final 3
years of study (Long, 1990), along with Female (gender),
Married (marital status), kid5 (number of
children under 6), Phdprestige (prestige of PhD program),
and Mentor (articles published by mentor).
- Import and inspect the data
dt <- read.csv(file.path(data_dir, "articles.csv"))
# the source file has a trailing comma producing an extra empty column - drop it
dt <- dt[, c("Articles","Female","Married","kid5","Phdprestige","Mentor")]
str(dt)## 'data.frame': 915 obs. of 6 variables:
## $ Articles : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Female : chr "Men" "Women" "Women" "Men" ...
## $ Married : chr "Married" "Single" "Single" "Married" ...
## $ kid5 : int 0 0 0 1 0 2 0 2 0 0 ...
## $ Phdprestige: num 2.52 2.05 3.75 1.18 3.75 3.59 3.19 2.96 4.62 1.25 ...
## $ Mentor : int 7 6 6 3 26 2 3 4 6 0 ...#check variable type
dt$Female <- as.factor(dt$Female)
dt$Married <- as.factor(dt$Married)
str(dt)## 'data.frame': 915 obs. of 6 variables:
## $ Articles : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Female : Factor w/ 2 levels "Men","Women": 1 2 2 1 2 2 2 1 1 2 ...
## $ Married : Factor w/ 2 levels "Married","Single": 1 2 2 1 2 1 2 1 2 1 ...
## $ kid5 : int 0 0 0 1 0 2 0 2 0 0 ...
## $ Phdprestige: num 2.52 2.05 3.75 1.18 3.75 3.59 3.19 2.96 4.62 1.25 ...
## $ Mentor : int 7 6 6 3 26 2 3 4 6 0 ...- Descriptive statistics
head(dt)summary(dt)## Articles Female Married kid5 Phdprestige
## Min. : 0.000 Men :494 Married:606 Min. :0.0000 Min. :0.750
## 1st Qu.: 0.000 Women:421 Single :309 1st Qu.:0.0000 1st Qu.:2.260
## Median : 1.000 Median :0.0000 Median :3.150
## Mean : 1.693 Mean :0.4951 Mean :3.103
## 3rd Qu.: 2.000 3rd Qu.:1.0000 3rd Qu.:3.920
## Max. :19.000 Max. :3.0000 Max. :4.620
## Mentor
## Min. : 0.000
## 1st Qu.: 3.000
## Median : 6.000
## Mean : 8.767
## 3rd Qu.:12.000
## Max. :77.000hist(dt$Articles, xlab = "Number of Articles", main = "")
plot(dt$Phdprestige, dt$Articles, xlab = "PhD program prestige", ylab = "Articles")
boxplot(dt$Articles ~ dt$Female, xlab = "Female", ylab = "Articles")
boxplot(dt$Articles ~ dt$Married, xlab = "Married", ylab = "Articles")
aggregate(Articles~Female, data=dt, FUN=summary)aggregate(Articles~Married, data=dt, FUN=summary)# a quick look at the mean-variance relationship: for Poisson we expect mean ~= variance
mean(dt$Articles)## [1] 1.692896var(dt$Articles)## [1] 3.709742Notice the variance is far larger than the mean - a first warning sign that a Poisson model will not fit well.
- Fit the Poisson model first (as a baseline) and perform variable selection
null.pr <- glm(Articles~1, data=dt, family='poisson')
full.pr <- glm(Articles~Female+Married+kid5+Phdprestige+Mentor,
data=dt, family='poisson')
library(MASS)
stepAIC(full.pr, direction='both') #stepwise## Start: AIC=3314.11
## Articles ~ Female + Married + kid5 + Phdprestige + Mentor
##
## Df Deviance AIC
## - Phdprestige 1 1634.6 3312.3
## <none> 1634.4 3314.1
## - Married 1 1640.8 3318.5
## - Female 1 1651.5 3329.2
## - kid5 1 1656.5 3334.2
## - Mentor 1 1766.2 3444.0
##
## Step: AIC=3312.35
## Articles ~ Female + Married + kid5 + Mentor
##
## Df Deviance AIC
## <none> 1634.6 3312.3
## + Phdprestige 1 1634.4 3314.1
## - Married 1 1640.8 3316.6
## - Female 1 1651.8 3327.5
## - kid5 1 1656.7 3332.4
## - Mentor 1 1776.7 3452.5##
## Call: glm(formula = Articles ~ Female + Married + kid5 + Mentor, family = "poisson",
## data = dt)
##
## Coefficients:
## (Intercept) FemaleWomen MarriedSingle kid5 Mentor
## 0.49735 -0.22530 -0.15218 -0.18499 0.02576
##
## Degrees of Freedom: 914 Total (i.e. Null); 910 Residual
## Null Deviance: 1817
## Residual Deviance: 1635 AIC: 3312# stepwise selection points to dropping Phdprestige
fit.pr <- glm(Articles~Female+Married+kid5+Mentor,
data=dt, family='poisson')
summary(fit.pr)##
## Call:
## glm(formula = Articles ~ Female + Married + kid5 + Mentor, family = "poisson",
## data = dt)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.49735 0.05224 9.520 < 2e-16 ***
## FemaleWomen -0.22530 0.05461 -4.125 3.70e-05 ***
## MarriedSingle -0.15218 0.06107 -2.492 0.0127 *
## kid5 -0.18499 0.04014 -4.609 4.05e-06 ***
## Mentor 0.02576 0.00195 13.212 < 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: 1817.4 on 914 degrees of freedom
## Residual deviance: 1634.6 on 910 degrees of freedom
## AIC: 3312.3
##
## Number of Fisher Scoring iterations: 5- Test for over-dispersion
library(AER)
dispersiontest(fit.pr, trafo = 1)##
## Overdispersion test
##
## data: fit.pr
## z = 5.8377, p-value = 2.647e-09
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 0.8218008# a quick rule-of-thumb check: residual deviance / df should be close to 1 for Poisson
fit.pr$deviance / fit.pr$df.residual## [1] 1.796272A dispersion test p-value below 0.05 (and a deviance/df ratio well above 1) confirms over-dispersion: the Poisson model understates the standard errors, which can make predictors look more significant than they really are. This is a strong signal to move to a Negative Binomial model.
- Fit the Negative Binomial model
null.nb <- glm.nb(Articles~1, data=dt)
full.nb <- glm.nb(Articles~Female+Married+kid5+Phdprestige+Mentor, data=dt)
# variable selection using stepAIC()
stepAIC(full.nb, direction = "both")## Start: AIC=3133.92
## Articles ~ Female + Married + kid5 + Phdprestige + Mentor
##
## Df AIC
## - Phdprestige 1 3132.1
## <none> 3133.9
## - Married 1 3135.3
## - Female 1 3140.7
## - kid5 1 3143.0
## - Mentor 1 3203.1
##
## Step: AIC=3132.1
## Articles ~ Female + Married + kid5 + Mentor
##
## Df AIC
## <none> 3132.1
## - Married 1 3133.3
## + Phdprestige 1 3133.9
## - Female 1 3138.9
## - kid5 1 3141.2
## - Mentor 1 3206.9##
## Call: glm.nb(formula = Articles ~ Female + Married + kid5 + Mentor,
## data = dt, init.theta = 2.264120079, link = log)
##
## Coefficients:
## (Intercept) FemaleWomen MarriedSingle kid5 Mentor
## 0.45027 -0.21667 -0.14694 -0.17680 0.02943
##
## Degrees of Freedom: 914 Total (i.e. Null); 910 Residual
## Null Deviance: 1109
## Residual Deviance: 1004 AIC: 3134fit.nb <- glm.nb(Articles~Female+Married+kid5+Mentor, data=dt)
summary(fit.nb)##
## Call:
## glm.nb(formula = Articles ~ Female + Married + kid5 + Mentor,
## data = dt, init.theta = 2.26411693, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.450272 0.071833 6.268 3.65e-10 ***
## FemaleWomen -0.216673 0.072624 -2.983 0.002850 **
## MarriedSingle -0.146944 0.081765 -1.797 0.072312 .
## kid5 -0.176797 0.052826 -3.347 0.000818 ***
## Mentor 0.029430 0.003108 9.470 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.2641) family taken to be 1)
##
## Null deviance: 1108.9 on 914 degrees of freedom
## Residual deviance: 1004.4 on 910 degrees of freedom
## AIC: 3134.1
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.264
## Std. Err.: 0.271
##
## 2 x log-likelihood: -3122.096# the estimated dispersion parameter theta
fit.nb$theta## [1] 2.264117- Compare Poisson vs Negative Binomial formally
Because the Poisson model is nested inside the Negative Binomial
model as theta goes to infinity, we can use a likelihood ratio test (the
test statistic follows a 50:50 mixture of chi-sq_0 and chi-sq_1, so
pchisq(..., df=1)/2 gives the correct one-sided
p-value).
LR <- 2 * (logLik(fit.nb) - logLik(fit.pr))
LR## 'log Lik.' 180.2526 (df=6)pchisq(as.numeric(LR), df = 1, lower.tail = FALSE) / 2## [1] 2.134203e-41# Equivalently and more conveniently:
library(pscl)
odTest(fit.nb) # tests H0: Poisson is adequate (theta -> infinity) vs Negative Binomial## Likelihood ratio test of H0: Poisson, as restricted NB model:
## n.b., the distribution of the test-statistic under H0 is non-standard
## e.g., see help(odTest) for details/references
##
## Critical value of test statistic at the alpha= 0.05 level: 2.7055
## Chi-Square Test Statistic = 180.2526 p-value = < 2.2e-16A significant result confirms that the Negative Binomial model fits significantly better than the Poisson model, i.e. the over-dispersion is real and should be modeled.
- Interpret coefficients as rate ratios and obtain confidence intervals
exp(coef(fit.nb))## (Intercept) FemaleWomen MarriedSingle kid5 Mentor
## 1.5687391 0.8051929 0.8633420 0.8379497 1.0298674exp(confint(fit.nb))## 2.5 % 97.5 %
## (Intercept) 1.3572663 1.8118634
## FemaleWomen 0.6982060 0.9283549
## MarriedSingle 0.7355413 1.0131409
## kid5 0.7548603 0.9294277
## Mentor 1.0231431 1.0367616#get pseudo R2
pR2(fit.nb)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML
## -1561.0481131 -1609.9367432 97.7772602 0.0303668 0.1013489
## r2CU
## 0.1044435Interpretation example:
exp(coefficient) for Mentor tells you the
multiplicative change in the expected number of articles for
each additional article published by the mentor, holding gender, marital
status, and number of young children constant.
Reading these numbers: the fitted rate ratios are
Female (Women) ≈ 0.81,
Married (Single) ≈ 0.86, kid5
≈ 0.84, and Mentor ≈
1.03. In words: holding the other variables fixed,
women are expected to publish about 19% fewer articles than men, and
each additional young child at home is associated with a further 16%
reduction in expected articles - while each additional article
published by the mentor multiplies the student’s own expected output by
about 3% (small per-unit, but mentors range widely - a mentor with
20 vs. 0 articles corresponds to roughly \(1.03^{20} \approx 1.8\times\) the expected
output). The estimated dispersion is \(\hat\theta \approx 2.26\), confirming
genuine over-dispersion (a Poisson model would require \(\theta = \infty\)). Comparing standard
errors side by side, the Negative Binomial gives noticeably
larger SEs than the Poisson fit for every coefficient
(e.g. Mentor: 0.0031 vs. 0.0019) - this is the
over-dispersion correction in action, and it’s the reason the Poisson
model’s p-values were overly optimistic.
- Residual and influence diagnostics
dt$pred.nb <- predict(fit.nb, newdata = dt, type = "response")
dt$res.nb <- residuals(fit.nb, type = "deviance")
plot(dt$pred.nb, dt$res.nb, xlab = "Fitted values", ylab = "Deviance residuals",
main = "Residuals vs Fitted: Negative Binomial model")
abline(h = 0, col = "red", lty = 2)
cooksd <- cooks.distance(fit.nb)
plot(cooksd, type = "h", main = "Cook's distance: Negative Binomial model")
abline(h = 4/nrow(dt), col = "red", lty = 2)
Lab 6 Assignment
Using the articles.csv dataset:
- Fit both a Poisson and a Negative Binomial model with the same
predictors (
Female,Married,kid5,Phdprestige,Mentor). Report the dispersion test result and theodTest()result. What do you conclude?
- Perform variable selection on the Negative Binomial model using
stepAIC(). Report your final model.
- Interpret the rate ratio (
exp(coefficient)) for one significant predictor.
- Compare the standard errors of the Poisson and Negative Binomial
models for the same predictors. Which model gives more trustworthy
inference here, and why?
- Check the residual plot and Cook’s distance for the final Negative
Binomial model. Are there any observations you would investigate
further?
Lab 7: Zero-Inflated Poisson (ZIP) regression
Objectives: Students are able to recognize and model
count data with excess zeros:
1. distinguish “structural” zeros (can never occur) from “sampling”
zeros (happened not to occur) conceptually
2. fit a Zero-Inflated Poisson model using zeroinfl(), with
separate predictors for the count part and the zero-inflation part
3. interpret the two parts of a ZIP model (logit part vs. Poisson part)
separately
4. formally compare a ZIP model against an ordinary Poisson model (Vuong
test)
Example: Scientific productivity (articles data)
What the model is
A Zero-Inflated Poisson model assumes the population is a mixture of
two latent groups: “always-zero” individuals (who structurally cannot
have the event - e.g. a biochemist who has left research and will never
publish again) and “count” individuals whose counts follow an ordinary
Poisson distribution (which can, of course, still include sampling zeros
- someone who simply hasn’t published yet). The overall probability of
observing \(y=0\) is therefore higher
than a plain Poisson would predict: \[P(y=0)
= \pi + (1-\pi)e^{-\mu}, \qquad P(y=k) = (1-\pi)\frac{e^{-\mu}\mu^k}{k!}
\text{ for } k \ge 1\] where \(\pi\) (the “zero-inflation probability”)
and \(\mu\) (the Poisson mean for the
count part) are each modeled with their own logit and log link
respectively, each with its own set of predictors. This is why
zeroinfl() reports two full sets of
coefficients: exp(count coefficient) is a rate
ratio (as in ordinary Poisson) for those in the “could publish” group,
while exp(zero coefficient) is an odds ratio for the
probability of being in the “always-zero” group.
The articles.csv dataset records the number of articles
(Articles) published by PhD biochemists in their final 3
years of study (Long, 1990), along with Female (gender),
Married (marital status), kid5 (number of
children under 6), Phdprestige (prestige of PhD program),
and Mentor (articles published by mentor). Many biochemists
published zero articles, which is a classic setting for zero-inflated
models.
- Import and inspect the data
dat_art <- read.csv(file.path(data_dir, "articles.csv"))
# the source file has a trailing comma producing an extra empty column - drop it
dat_art <- dat_art[, c("Articles","Female","Married","kid5","Phdprestige","Mentor")]
dat_art$Female <- factor(dat_art$Female)
dat_art$Married <- factor(dat_art$Married)
str(dat_art)## 'data.frame': 915 obs. of 6 variables:
## $ Articles : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Female : Factor w/ 2 levels "Men","Women": 1 2 2 1 2 2 2 1 1 2 ...
## $ Married : Factor w/ 2 levels "Married","Single": 1 2 2 1 2 1 2 1 2 1 ...
## $ kid5 : int 0 0 0 1 0 2 0 2 0 0 ...
## $ Phdprestige: num 2.52 2.05 3.75 1.18 3.75 3.59 3.19 2.96 4.62 1.25 ...
## $ Mentor : int 7 6 6 3 26 2 3 4 6 0 ...summary(dat_art$Articles)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.693 2.000 19.000table(dat_art$Articles)##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 16 19
## 275 246 178 84 67 27 17 12 1 2 1 1 2 1 1# visualize the excess of zeros
hist(dat_art$Articles, breaks = 0:20,
main = "Number of articles published", xlab = "Number of articles")
- Fit an ordinary Poisson model as a baseline
fit.pois <- glm(Articles ~ Female + Married + kid5 + Phdprestige + Mentor,
family = poisson(link = "log"), data = dat_art)
summary(fit.pois)##
## Call:
## glm(formula = Articles ~ Female + Married + kid5 + Phdprestige +
## Mentor, family = poisson(link = "log"), data = dat_art)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.459809 0.093331 4.927 8.36e-07 ***
## FemaleWomen -0.224593 0.054613 -4.112 3.92e-05 ***
## MarriedSingle -0.155247 0.061374 -2.530 0.0114 *
## kid5 -0.184882 0.040127 -4.607 4.08e-06 ***
## Phdprestige 0.012840 0.026395 0.486 0.6266
## Mentor 0.025542 0.002006 12.732 < 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: 1817.4 on 914 degrees of freedom
## Residual deviance: 1634.4 on 909 degrees of freedom
## AIC: 3314.1
##
## Number of Fisher Scoring iterations: 5library(AER)
dispersiontest(fit.pois, trafo = 1)##
## Overdispersion test
##
## data: fit.pois
## z = 5.7824, p-value = 3.682e-09
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 0.8245435- Fit the Zero-Inflated Poisson model
zeroinfl() fits two linked sub-models at once: the count
model (Poisson, for biochemists who could publish), and the
zero-inflation model (a logistic model for the probability of being a
“structural” zero, i.e. someone who will never publish in this window).
The formula y ~ x1 + x2 | z1 + z2 separates the two:
predictors before | go in the count model, predictors after
| go in the zero-inflation model.
library(pscl)
fit.zip <- zeroinfl(Articles ~ Female + Married + kid5 + Phdprestige + Mentor
| Female + Married + kid5 + Phdprestige + Mentor,
data = dat_art, dist = "poisson")
summary(fit.zip)##
## Call:
## zeroinfl(formula = Articles ~ Female + Married + kid5 + Phdprestige +
## Mentor | Female + Married + kid5 + Phdprestige + Mentor, data = dat_art,
## dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -2.3253 -0.8652 -0.2826 0.5404 7.2976
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.744572 0.110282 6.752 1.46e-11 ***
## FemaleWomen -0.209143 0.063405 -3.299 0.000972 ***
## MarriedSingle -0.103753 0.071111 -1.459 0.144560
## kid5 -0.143321 0.047429 -3.022 0.002513 **
## Phdprestige -0.006161 0.031009 -0.199 0.842522
## Mentor 0.018098 0.002294 7.888 3.07e-15 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.930844 0.469691 -1.982 0.04750 *
## FemaleWomen 0.109754 0.280084 0.392 0.69516
## MarriedSingle 0.354026 0.317612 1.115 0.26500
## kid5 0.217096 0.196482 1.105 0.26920
## Phdprestige 0.001176 0.145271 0.008 0.99354
## Mentor -0.134103 0.045243 -2.964 0.00304 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 21
## Log-likelihood: -1605 on 12 Df- Interpret the two parts separately
# Count model: exp(coefficient) is a rate ratio for those NOT in the "always-zero" group
exp(coef(fit.zip, model = "count"))## (Intercept) FemaleWomen MarriedSingle kid5 Phdprestige
## 2.1055406 0.8112792 0.9014483 0.8664759 0.9938584
## Mentor
## 1.0182625# Zero-inflation model: exp(coefficient) is an odds ratio for being a "structural" zero
exp(coef(fit.zip, model = "zero"))## (Intercept) FemaleWomen MarriedSingle kid5 Phdprestige
## 0.3942207 1.1160033 1.4247923 1.2424635 1.0011763
## Mentor
## 0.8745003Reading these numbers: in the count
part, rate ratios are close to those from the plain NB model in
Lab 6 - e.g. Mentor ≈ 1.02,
kid5 ≈ 0.87 - describing the publication
rate among biochemists who are still active researchers. In the
zero-inflation part, the odds ratios describe who is
more likely to be a “structural zero”: Married (Single) ≈
1.43 and kid5 ≈ 1.24 both
increase the odds of never publishing, while
Mentor ≈ 0.88 decreases it -
i.e. having a more productive mentor makes a student less likely to
fall into the “will never publish” group, on top of raising the
publication rate for those who do publish. This is the key insight a
plain Poisson or NB model cannot give you: it separates “what drives
whether someone publishes at all” from “what drives how
many articles they publish, given that they publish.”
- Compare the ZIP model to the plain Poisson model
Because ZIP and Poisson are not nested in the usual sense, we use the Vuong test rather than a likelihood ratio test.
vuong(fit.pois, fit.zip)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -4.180506 model2 > model1 1.4543e-05
## AIC-corrected -3.638557 model2 > model1 0.00013709
## BIC-corrected -2.332751 model2 > model1 0.00983060A significant Vuong test favouring fit.zip indicates the
zero-inflated model fits significantly better than the plain Poisson
model, i.e. there is evidence of excess zeros beyond what Poisson alone
can explain. Running this comparison here gives an AIC-corrected
z-statistic favouring the ZIP model (p < 0.001) - confirming that
modelling the excess zeros explicitly is worthwhile, on top of the
over-dispersion issue Lab 6 already flagged.
- Prediction
# predicted mean count (combines both parts)
pred.zip <- predict(fit.zip, type = "response")
head(pred.zip)## 1 2 3 4 5 6
## 2.0379328 1.3230829 1.3087242 1.4399136 2.3632355 0.8548015# predicted probability of being a structural zero
pred.zero <- predict(fit.zip, type = "zero")
head(pred.zero)## 1 2 3 4 5 6
## 0.13393555 0.21938213 0.21972458 0.24700144 0.01890346 0.34279049Lab 7 Assignment
Using the articles.csv dataset:
- Fit a ZIP model using
kid5andMentoras predictors in both the count part and the zero-inflation part. Report and interpret the zero-inflation coefficients.
- Compare this ZIP model to the corresponding plain Poisson model
using the Vuong test. State your conclusion.
- Which variables are most associated with being a “structural zero”
(never publishing)? Which are most associated with the rate of
publishing among those who do publish?
- Obtain the predicted mean number of articles for a hypothetical
unmarried, childless, female biochemist with an average
PhdprestigeandMentorvalue. Compare this to a married male biochemist with two young children and the samePhdprestigeandMentor.
Lab 8: Zero-Inflated Negative Binomial (ZINB) regression
Objectives: Students are able to model count data
that shows both excess zeros AND over-dispersion beyond what Poisson
allows:
1. recognize when a ZIP model’s count part is still over-dispersed
2. fit a Zero-Inflated Negative Binomial model using
zeroinfl(..., dist = "negbin")
3. compare ZIP vs. ZINB, and ZINB vs. plain Negative Binomial, to choose
the best-fitting model
4. interpret and report the final chosen model
Example: Scientific productivity (continued)
What the model is
The Zero-Inflated Negative Binomial model combines both extensions at
once: it has the same two-part “always-zero vs. count” mixture structure
as ZIP, but replaces the Poisson count part with a Negative Binomial, so
it can absorb over-dispersion and excess zeros
simultaneously: \[P(y=0) = \pi +
(1-\pi)\left(\frac{\theta}{\theta+\mu}\right)^{\theta}, \qquad P(y=k) =
(1-\pi)\,\text{NB}(k;\mu,\theta) \text{ for } k \ge 1\] As with
ZIP, exp(count coefficient) is a rate ratio and
exp(zero coefficient) is an odds ratio for structural-zero
membership; the extra theta in the count part behaves
exactly as in Lab 6’s plain NB model, absorbing any over-dispersion left
in the “could publish” group after accounting for the always-zero group
separately.
- Fit the Zero-Inflated Negative Binomial model
fit.zinb <- zeroinfl(Articles ~ Female + Married + kid5 + Phdprestige + Mentor
| Female + Married + kid5 + Phdprestige + Mentor,
data = dat_art, dist = "negbin")
summary(fit.zinb)##
## Call:
## zeroinfl(formula = Articles ~ Female + Married + kid5 + Phdprestige +
## Mentor | Female + Married + kid5 + Phdprestige + Mentor, data = dat_art,
## dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.2940 -0.7599 -0.2912 0.4452 6.4160
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5144163 0.1289306 3.990 6.61e-05 ***
## FemaleWomen -0.1949255 0.0755492 -2.580 0.00988 **
## MarriedSingle -0.0981557 0.0844325 -1.163 0.24502
## kid5 -0.1518872 0.0541862 -2.803 0.00506 **
## Phdprestige -0.0007336 0.0362614 -0.020 0.98386
## Mentor 0.0247797 0.0034921 7.096 1.29e-12 ***
## Log(theta) 0.9761907 0.1354436 7.207 5.70e-13 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.71006 1.04576 -1.635 0.10200
## FemaleWomen 0.66789 0.85174 0.784 0.43295
## MarriedSingle 1.47875 0.93353 1.584 0.11319
## kid5 0.63321 0.44219 1.432 0.15215
## Phdprestige -0.03395 0.30791 -0.110 0.91221
## Mentor -0.89064 0.32075 -2.777 0.00549 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 2.6543
## Number of iterations in BFGS optimization: 28
## Log-likelihood: -1550 on 13 Df# the estimated dispersion parameter (theta); values well below 1 indicate strong
# over-dispersion remaining in the count part even after modeling excess zeros
fit.zinb$theta## [1] 2.654326- Compare ZIP vs. ZINB
fit.zip (Poisson-based zero-inflation) is nested inside
fit.zinb (Negative-Binomial-based zero-inflation) as the NB
dispersion parameter goes to infinity, so in principle a likelihood
ratio test applies (with the same caveat about boundary testing as in
Lab 6); in practice, most analysts compare AIC or use the Vuong
test.
AIC(fit.zip, fit.zinb)vuong(fit.zip, fit.zinb)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -3.370357 model2 > model1 0.00037535
## AIC-corrected -3.370357 model2 > model1 0.00037535
## BIC-corrected -3.370357 model2 > model1 0.00037535- Compare ZINB vs. a plain Negative Binomial (no zero-inflation)
library(MASS)
fit.nb <- glm.nb(Articles ~ Female + Married + kid5 + Phdprestige + Mentor, data = dat_art)
# Vuong test again, since the two are not nested in the classical sense
vuong(fit.nb, fit.zinb)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -2.243032 model2 > model1 0.012447
## AIC-corrected -1.015746 model2 > model1 0.154875
## BIC-corrected 1.941352 model1 > model2 0.026108AIC(fit.nb, fit.zinb)- Model selection summary
A convenient way to decide among Poisson / Negative Binomial / ZIP / ZINB for a given dataset is to lay out AIC (and, where applicable, Vuong test results) side by side:
data.frame(
Model = c("Poisson", "Negative Binomial", "ZIP", "ZINB"),
AIC = c(AIC(fit.pois), AIC(fit.nb), AIC(fit.zip), AIC(fit.zinb))
)The model with the lowest AIC (supported by a significant Vuong test against its simpler competitors) is generally the preferred model, so long as its two parts (count and zero-inflation) both have a sensible, interpretable set of predictors.
Reading these numbers: a representative run gives AIC ≈ 3314 (Poisson), 3136 (NB), 3234 (ZIP), and 3126 (ZINB) - ZINB has the lowest AIC, with NB a close second, and both comfortably beat Poisson and plain ZIP. However, the ZINB-vs-NB Vuong test tells a more nuanced story: the raw z-statistic favours ZINB, but once you penalise for ZINB’s extra parameters (the AIC-corrected version), the improvement is no longer statistically significant, and the BIC-corrected version (which penalises complexity even more heavily) actually favours the simpler plain NB model. This is a genuinely useful, realistic result to see in teaching: the two most sophisticated models (ZINB and NB) are statistically hard to tell apart here, so a defensible final choice could be either - and if you cannot distinguish them formally, the simpler, easier-to-explain NB model (Lab 6) is usually the better one to report.
- Interpret the final chosen model
exp(coef(fit.zinb, model = "count"))## (Intercept) FemaleWomen MarriedSingle kid5 Phdprestige
## 1.6726618 0.8228959 0.9065077 0.8590852 0.9992667
## Mentor
## 1.0250893exp(coef(fit.zinb, model = "zero"))## (Intercept) FemaleWomen MarriedSingle kid5 Phdprestige
## 0.1808541 1.9501254 4.3874430 1.8836460 0.9666208
## Mentor
## 0.4103926Reading these numbers: the ZINB
count-part rate ratios (Female ≈
0.82, kid5 ≈ 0.86,
Mentor ≈ 1.03) are very close to the plain
NB model’s, because both are describing the same “publication rate among
active researchers.” The zero-inflation part is where
ZINB adds new information, and its odds ratios are considerably larger
than ZIP’s: Female (Women) ≈ 1.95 and
Married (Single) ≈ 4.39 both sharply
increase the odds of being a structural zero, while Mentor
≈ 0.41 sharply decreases it. In a one-sentence summary
for each part: “among biochemists who do publish, having more young
children modestly reduces their output” (count part), and
“single, unmentored biochemists are substantially more likely to
never publish at all” (zero-inflation part) - two distinct
mechanisms that a plain Poisson model would have blurred together into a
single, harder-to-interpret coefficient.
Lab 8 Assignment
Using the articles.csv dataset:
- Fit Poisson, Negative Binomial, ZIP, and ZINB models with the same
set of predictors (
Female,Married,kid5,Phdprestige,Mentor) in every part of every model. Tabulate the AIC of all four.
- Use
odTest()to compare Poisson vs. Negative Binomial, andvuong()to compare Poisson vs. ZIP, and ZIP vs. ZINB. Summarize the sequence of decisions this leads you through.
- Based on your comparisons, which single model would you report as
final? Justify your choice using both statistical criteria (AIC,
Vuong/LR tests) and interpretability.
- For your chosen final model, write one sentence interpreting a
count-part coefficient and one sentence interpreting a
zero-inflation-part coefficient.
- In your own words, explain to a non-statistician colleague why a
Poisson model would have been the wrong choice for this dataset.
Comprehensive comparison: Poisson vs. NB vs. ZIP vs. ZINB
Having built all four count models on the same data, it is worth
stepping back and comparing them side by side, using the same,
final set of predictors
(Female + Married + kid5 + Mentor, as selected by
stepAIC() in Lab 6) in every model. This is exactly the
workflow a data analyst goes through to decide which model to
report.
dt$Female <- factor(dt$Female)
dt$Married <- factor(dt$Married)
fit.poi <- glm(Articles ~ Female+Married+kid5+Mentor, data=dt, family='poisson')
fit.nb <- glm.nb(Articles ~ Female+Married+kid5+Mentor, data=dt)
fit.zip <- zeroinfl(Articles ~ Female+Married+kid5+Mentor, data=dt, dist='poisson')
fit.zinb <- zeroinfl(Articles ~ Female+Married+kid5+Mentor, data=dt, dist='negbin')- Get fitted/predicted values from each model
dt$pred.poi <- predict(fit.poi, newdata=dt, type='response')
dt$pred.nb <- predict(fit.nb, newdata=dt, type='response')
dt$pred.zip <- predict(fit.zip, newdata=dt, type='response')
dt$pred.zinb <- predict(fit.zinb, newdata=dt, type='response')- Compare AIC, BIC, and RMSE across the four models
library(Metrics)
model_comparison <- data.frame(
Model = c("Poisson", "Negative Binomial", "ZIP", "ZINB"),
AIC = c(AIC(fit.poi), AIC(fit.nb), AIC(fit.zip), AIC(fit.zinb)),
BIC = c(BIC(fit.poi), BIC(fit.nb), BIC(fit.zip), BIC(fit.zinb)),
RMSE = c(rmse(dt$Articles, dt$pred.poi),
rmse(dt$Articles, dt$pred.nb),
rmse(dt$Articles, dt$pred.zip),
rmse(dt$Articles, dt$pred.zinb))
)
model_comparisonLower AIC/BIC and lower RMSE both point toward a better-fitting model - but note that RMSE on the mean count rewards models whose average prediction is close to observed values, while AIC/BIC reward the whole fitted distribution (including how well the model captures the excess zeros). A model can have a similar RMSE to a simpler competitor while fitting the zero/near-zero counts far better, which AIC/BIC will reflect but a single RMSE number will not.
- Compare the predicted probability distribution to the observed distribution (a “rootogram” style check)
This is the most direct way to see why Poisson typically loses to the other three: does each model predict the right proportion of scientists with 0, 1, 2, … articles?
# observed relative frequency of each count value
# (force all levels 0..max(Articles) to appear, even ones with zero observations,
# so this lines up exactly with the predprob() columns below)
fr <- table(factor(dt$Articles, levels = 0:max(dt$Articles)))
prob_obs <- prop.table(fr)
# average predicted probability of each count value under each model
pr.poi <- predprob(fit.poi)
pr.nb <- predprob(fit.nb)
pr.zip <- predprob(fit.zip)
pr.zinb <- predprob(fit.zinb)
proby_poi <- apply(pr.poi, 2, mean)
proby_nb <- apply(pr.nb, 2, mean)
proby_zip <- apply(pr.zip, 2, mean)
proby_zinb <- apply(pr.zinb, 2, mean)
xseq <- seq(0, max(dt$Articles), by = 1)
plot(xseq, proby_poi, col = "blue", type = "o", lty = "dashed",
xlab = "Number of articles", ylab = "Probability",
main = "Observed vs. model-implied distribution of Articles")
points(xseq, proby_nb, col = "red", type = "o", lty = "dashed")
points(xseq, proby_zip, col = "purple", type = "o", lty = "dashed")
points(xseq, proby_zinb, col = "green", type = "o", lty = "dashed")
points(xseq, as.vector(prob_obs), col = "black", type = "o", lwd = 2)
legend("topright", legend = c("Poisson","NB","ZIP","ZINB","Observed"),
col = c("blue","red","purple","green","black"), lty = 1, pch = 1,
lwd = c(1,1,1,1,2), cex = 0.8)
The black line (observed) sitting much higher than every colored line
at Articles = 0 is the classic visual signature of
zero-inflation; whichever colored line tracks the black line most
closely across all counts (not just zero) is doing the best
overall job.
- Residual and observed-vs-predicted plots, side by side
dt$res.poi <- resid(fit.poi)
dt$res.nb <- resid(fit.nb)
dt$res.zip <- resid(fit.zip)
dt$res.zinb <- resid(fit.zinb)
par(mfrow=c(2,2))
plot(dt$pred.poi, dt$res.poi, xlab="Fitted", ylab="Residuals", main="Poisson")
plot(dt$pred.nb, dt$res.nb, xlab="Fitted", ylab="Residuals", main="Negative Binomial")
plot(dt$pred.zip, dt$res.zip, xlab="Fitted", ylab="Residuals", main="ZIP")
plot(dt$pred.zinb, dt$res.zinb, xlab="Fitted", ylab="Residuals", main="ZINB")
par(mfrow=c(1,1))- Put everything in long form for a single, faceted
ggplot2comparison
n <- nrow(dt)
dtf <- rbind.data.frame(
data.frame(Model="Observed", Female=dt$Female, Married=dt$Married, kid5=dt$kid5,
Mentor=dt$Mentor, Articles=dt$Articles, yhat=dt$Articles, res=0),
data.frame(Model="Poisson", Female=dt$Female, Married=dt$Married, kid5=dt$kid5,
Mentor=dt$Mentor, Articles=dt$Articles, yhat=dt$pred.poi, res=dt$res.poi),
data.frame(Model="NB", Female=dt$Female, Married=dt$Married, kid5=dt$kid5,
Mentor=dt$Mentor, Articles=dt$Articles, yhat=dt$pred.nb, res=dt$res.nb),
data.frame(Model="ZIP", Female=dt$Female, Married=dt$Married, kid5=dt$kid5,
Mentor=dt$Mentor, Articles=dt$Articles, yhat=dt$pred.zip, res=dt$res.zip),
data.frame(Model="ZINB", Female=dt$Female, Married=dt$Married, kid5=dt$kid5,
Mentor=dt$Mentor, Articles=dt$Articles, yhat=dt$pred.zinb, res=dt$res.zinb)
)
dtf$Model <- factor(dtf$Model, levels=c("Observed","Poisson","NB","ZIP","ZINB"))
# distribution of predicted counts by model, vs the observed distribution
ggplot(dtf, aes(yhat, group=Model, fill=Model)) +
geom_histogram(alpha=0.6, position='identity', bins = 20) +
xlab("Number of articles") +
theme_bw()
# observed vs predicted, by model
ggplot(dtf, aes(x=Articles, y=yhat, group=Model)) +
geom_point(aes(color=Model)) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
theme_bw() +
xlab("Observed values") + ylab("Predicted values")
# predicted values against Mentor, faceted by Married and Female, sized by kid5
ggplot(dtf, aes(x=Mentor, y=yhat, group=Model)) +
geom_point(aes(color=Model, size=kid5)) +
theme_bw() +
xlab("Mentor") + ylab("Predicted values") +
facet_grid(Married ~ Female)
Take-away: for data like this - where over-dispersion and excess zeros are both present - Poisson is essentially always the worst choice, and the decision usually comes down to NB vs. ZINB (ZIP rarely wins once NB is on the table, since NB alone can often absorb much of what looks like “extra zeros” as generic over-dispersion). Use AIC/BIC, the Vuong test, and this kind of predicted-vs-observed distribution plot together, rather than relying on any single number.