Logistic Regression with the Plasma and Student Survey Data
Min Wang
2026-11-09
Binary Logistic Regression: Plasma Data
The plasma dataset contains measurements of fibrinogen
and globulin, along with an erythrocyte sedimentation rate (ESR)
category. The response indicates whether ESR is greater than 20
mm/hour.
Load the data
# Install these packages once, if needed:
# install.packages("HSAUR")
# install.packages("ResourceSelection")
# install.packages("DescTools")
library(HSAUR)
library(ResourceSelection)
library(DescTools)
data("plasma", package = "HSAUR")Prepare the response variable
Code ESR greater than 20 as 1 and ESR less than 20 as 0.
plasma$ESR_high <- as.integer(plasma$ESR == "ESR > 20")
# Check the response coding
table(plasma$ESR, plasma$ESR_high)##
## 0 1
## ESR < 20 26 0
## ESR > 20 0 6Model 1: Fibrinogen as the predictor
plasma_glm_1 <- glm(
ESR_high ~ fibrinogen,
data = plasma,
family = binomial()
)
summary(plasma_glm_1)##
## Call:
## glm(formula = ESR_high ~ fibrinogen, family = binomial(), data = plasma)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.8451 2.7703 -2.471 0.0135 *
## fibrinogen 1.8271 0.9009 2.028 0.0425 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 30.885 on 31 degrees of freedom
## Residual deviance: 24.840 on 30 degrees of freedom
## AIC: 28.84
##
## Number of Fisher Scoring iterations: 5Exponentiating the fibrinogen coefficient gives the estimated odds ratio for a one-unit increase in fibrinogen.
round(exp(coef(plasma_glm_1)), 3)## (Intercept) fibrinogen
## 0.001 6.216The exponentiated intercept is the estimated odds when fibrinogen is zero. The exponentiated fibrinogen coefficient is the odds ratio of interest.
Goodness of fit and pseudo-\(R^2\)
The Hosmer–Lemeshow test can be used as an exploratory check. Interpret it cautiously because this dataset is small.
hoslem.test(
x = plasma$ESR_high,
y = fitted(plasma_glm_1),
g = 5
)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: plasma$ESR_high, fitted(plasma_glm_1)
## X-squared = 4.4293, df = 3, p-value = 0.2187PseudoR2(
plasma_glm_1,
which = c("McFadden", "Nagel", "CoxSnell")
)## McFadden Nagelkerke CoxSnell
## 0.1957136 0.2780366 0.1721257Diagnostic plots
old_par <- par(mfrow = c(2, 2))
plot(plasma_glm_1)
par(old_par)Cook’s distance can help identify observations that may have a large influence on the fitted model.
plot(plasma_glm_1, which = 4)
The residuals-versus-leverage plot can help identify observations with high leverage or influence.
plot(plasma_glm_1, which = 5)
Predictions and classification
The following uses a cutoff of 0.5 to classify an observation as having ESR greater than 20.
# Estimated probability that ESR is greater than 20
prob_1 <- predict(plasma_glm_1, type = "response")
# Predicted class: 1 if probability > 0.5, and 0 otherwise
predicted_class_1 <- as.integer(prob_1 > 0.5)
prediction_results <- data.frame(
fibrinogen = plasma$fibrinogen,
observed_ESR_high = plasma$ESR_high,
predicted_probability = prob_1,
predicted_class = predicted_class_1
)
head(prediction_results)## fibrinogen observed_ESR_high predicted_probability predicted_class
## 1 2.52 0 0.09614329 0
## 2 2.56 0 0.10268435 0
## 3 2.19 0 0.05500436 0
## 4 2.18 0 0.05406235 0
## 5 3.41 0 0.35098115 0
## 6 2.46 0 0.08702953 0# Apparent accuracy on the data used to fit the model
mean(predicted_class_1 == plasma$ESR_high)## [1] 0.875# Apparent classification error rate
mean(predicted_class_1 != plasma$ESR_high)## [1] 0.125These are training-set results because the model is evaluated on the same observations used to fit it.
Model 2: Fibrinogen and globulin as predictors
plasma_glm_2 <- glm(
ESR_high ~ fibrinogen + globulin,
data = plasma,
family = binomial()
)
summary(plasma_glm_2)##
## Call:
## glm(formula = ESR_high ~ fibrinogen + globulin, family = binomial(),
## data = plasma)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.7921 5.7963 -2.207 0.0273 *
## fibrinogen 1.9104 0.9710 1.967 0.0491 *
## globulin 0.1558 0.1195 1.303 0.1925
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 30.885 on 31 degrees of freedom
## Residual deviance: 22.971 on 29 degrees of freedom
## AIC: 28.971
##
## Number of Fisher Scoring iterations: 5Compare the one-predictor and two-predictor models with a likelihood-ratio test.
anova(plasma_glm_1, plasma_glm_2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: ESR_high ~ fibrinogen
## Model 2: ESR_high ~ fibrinogen + globulin
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 30 24.840
## 2 29 22.971 1 1.8692 0.1716Multinomial Logistic Regression: Student Survey Data
The survey dataset in MASS contains
responses from a student survey. The original Smoke
variable has four categories: Never, Occas, Regul, and Heavy. Here,
Regul and Heavy are combined to create three response categories.
Load the data
# Install these packages once, if needed:
# install.packages("MASS")
# install.packages("nnet")
library(MASS)
library(nnet)
data("survey", package = "MASS")Prepare the survey data
Keep observations with complete data for the response and predictors used in the model.
survey3 <- subset(
survey,
!is.na(Smoke) &
!is.na(Sex) &
!is.na(Exer) &
!is.na(Pulse)
)
# Combine regular and heavy smokers
survey3$Smoke3 <- ifelse(
survey3$Smoke %in% c("Regul", "Heavy"),
"Regular/Heavy",
as.character(survey3$Smoke)
)
# Use "Never" as the reference category
survey3$Smoke3 <- factor(
survey3$Smoke3,
levels = c("Never", "Occas", "Regular/Heavy")
)
# Review the variables and response-category counts
head(survey3[, c("Smoke", "Smoke3", "Sex", "Exer", "Pulse")])## Smoke Smoke3 Sex Exer Pulse
## 1 Never Never Female Some 92
## 2 Regul Regular/Heavy Male None 104
## 3 Occas Occas Male None 87
## 5 Never Never Male Some 35
## 6 Never Never Female Some 64
## 7 Never Never Male Freq 83table(survey3$Smoke3)##
## Never Occas Regular/Heavy
## 151 16 23Fit the multinomial logistic regression model
This model predicts smoking category from sex, exercise frequency, and pulse rate.
survey_multinom <- multinom(
Smoke3 ~ Sex + Exer + Pulse,
data = survey3,
trace = FALSE
)
summary(survey_multinom)## Call:
## multinom(formula = Smoke3 ~ Sex + Exer + Pulse, data = survey3,
## trace = FALSE)
##
## Coefficients:
## (Intercept) SexMale ExerNone ExerSome Pulse
## Occas -2.204321 -0.02475553 0.6994517 -0.72116987 0.001730015
## Regular/Heavy -3.204543 0.76465362 -0.7754063 -0.02373879 0.012696221
##
## Std. Errors:
## (Intercept) SexMale ExerNone ExerSome Pulse
## Occas 1.704585 0.5360362 0.7428197 0.6366569 0.02268376
## Regular/Heavy 1.548706 0.4754776 1.1039182 0.4769802 0.01999040
##
## Residual Deviance: 238.988
## AIC: 258.988A convergence value of 0 indicates that the fitting algorithm converged.
survey_multinom$convergence## [1] 0Predict probabilities and classes
The first prediction gives the estimated probability for each smoking category for one student. The second gives the category with the largest estimated probability.
# Estimated class probabilities for the first student
predict(
survey_multinom,
newdata = survey3[1, ],
type = "probs"
)## Never Occas Regular/Heavy
## 0.84011223 0.05283672 0.10705105# Predicted smoking category for the first student
predict(
survey_multinom,
newdata = survey3[1, ],
type = "class"
)## [1] Never
## Levels: Never Occas Regular/HeavyThis survey example demonstrates multinomial logistic regression with three classes. It is not a medical diagnosis dataset.