Binary Response (ELMR Chapter 2, IGLM Chapter 7)

Biostat 200C

Dr. Zhe Fei @ UCLA

Apr 6, 2020

## -- Attaching packages -------------------------------------------------------------- tidyverse 1.3.0

## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.1
## v tidyr 1.1.1 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## -- Conflicts ----------------------------------------------------------------- tidyverse_conflicts()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()

Heart disease example

The data frame wcgs in faraway package contains data from the Western Collaborative Group Study.
data(wcgs)
help(wcgs)
wcgs %>% head(10)

## age height weight sdp dbp chol behave cigs dibep chd typechd timechd
## 2001 49 73 150 110 76 225 A2 25 B no none 1664
## 2002 42 70 160 154 84 177 A2 20 B no none 3071
## 2003 42 69 160 110 78 181 B3 0 A no none 3071
## 2004 41 68 152 124 78 132 B4 20 A no none 3064
## 2005 59 70 150 144 86 255 B3 20 A yes infdeath 1885
## 2006 44 72 204 150 90 182 B4 0 A no none 3102
## 2007 44 72 164 130 84 155 B4 0 A no none 3074
## 2008 40 71 150 138 60 140 A2 0 B no none 3071
## 2009 43 72 190 146 76 149 B3 25 A no none 3064
## 2010 42 70 175 132 90 325 A2 0 B no none 1032
## arcus
## 2001 absent
## 2002 present
## 2003 absent
## 2004 absent
## 2005 present
## 2006 absent
## 2007 absent
## 2008 absent
## 2009 absent
## 2010 present
We convert the data frame into a tibble for compatibility with tidyverse.

1
wcgs <- wcgs %>%
as_tibble() %>%
print(width = Inf)

## # A tibble: 3,154 x 13
## age height weight sdp dbp chol behave cigs dibep chd typechd
## <int> <int> <int> <int> <int> <int> <fct> <int> <fct> <fct> <fct>
## 1 49 73 150 110 76 225 A2 25 B no none
## 2 42 70 160 154 84 177 A2 20 B no none
## 3 42 69 160 110 78 181 B3 0 A no none
## 4 41 68 152 124 78 132 B4 20 A no none
## 5 59 70 150 144 86 255 B3 20 A yes infdeath
## 6 44 72 204 150 90 182 B4 0 A no none
## 7 44 72 164 130 84 155 B4 0 A no none
## 8 40 71 150 138 60 140 A2 0 B no none
## 9 43 72 190 146 76 149 B3 25 A no none
## 10 42 70 175 132 90 325 A2 0 B no none
## timechd arcus
## <int> <fct>
## 1 1664 absent
## 2 3071 present
## 3 3071 absent
## 4 3064 absent
## 5 1885 present
## 6 3102 absent
## 7 3074 absent
## 8 3071 absent
## 9 3064 absent
## 10 1032 present
## # ... with 3,144 more rows
For now, we focus just on variables
- chd, whether the person develops coronary heard disease or not,
- height, height of the person in inches,
- cigs, number of cigarettes smoked per day.
wcgs %>%
select(chd, height, cigs) %>%
summary()

## chd height cigs

## no :2897 Min. :60.00 Min. : 0.0
## yes: 257 1st Qu.:68.00 1st Qu.: 0.0
## Median :70.00 Median : 0.0
## Mean :69.78 Mean :11.6
## 3rd Qu.:72.00 3rd Qu.:20.0
## Max. :78.00 Max. :99.0
We use side-by-side boxplots to summarize the qulitative variable chd and quantitative variables height and
cigs
ggplot(data = wcgs) +
geom_boxplot(mapping = aes(x = chd, y = height))

2
 75

70
height

65

60

no yes
chd
and number of cigarretes smoked per day
ggplot(data = wcgs) +
geom_boxplot(mapping = aes(x = chd, y = cigs))

3
 100

75
cigs

50

25

no yes
chd
It seems more cigarettes is associated with heard disease, but not height. How can we formally analyze this?
If we use linear regression (straight line) for the anlaysis, the line will eventually extends beyond the [0, 1]
range, making interpretation hard.
ggplot(data = wcgs) +
geom_point(mapping = aes(x = cigs, y = chd))

4
 yes
chd

no

0 25 50 75 100
cigs

Logistic regression
• Bernoulli model for a binary response
(
1 with probability pi
Yi =
0 with probability 1 − pi

• The parameter pi = E(Yi ) will be related to the predictors X1 , . . . , Xq via an inverse link function

e ηi
pi = ,
1 + e ηi
where ηi is the linear predictor or systematic component

ηi = β0 + β1 xi1 + · · · + βq xiq = xiT β

with   
β0 1
 β1  xi1 
β =  . , xi =  .  .
   
 ..   .. 
βq xiq

• The function  
p
η = g(p) = log
1−p
that links E(Y ) to the systematic component is called the link function. This particular link function
is also called the logit function.

5
 • The function
eη
p = g −1 (η) =
1 + eη
is called the inverse link function. This particular function (inverse logit) is also called the logistic
function. A graph of the logistic function:
ggplot(data = tibble(x = 0), mapping = aes(x = x)) + # null data
stat_function(fun = ilogit) + # ilogit is from faraway
xlim(-6, 6) +
labs(x = expression(eta), y = "p", title = "Logistic function (inverse link)")

Logistic function (inverse link)

1.00

0.75

0.50
p

0.25

0.00

−6 −3 0 3 6
η
# curve(ilogit(x), -6, 6, xlab = expression(eta), ylab = "p")

• Therefore the above model is called the logistic regression.

Fitting logistic regression

• We use maximum likelihood estimation (MLE) to estimate the parameters β0 , . . . , βq .

6
 • Given n data points (yi , xi ), i = 1, . . . , n, the log-likelihood is
X
log pyi i (1 − pi )1−yi
 
`(β) =
i
X
= [yi log pi + (1 − yi ) log(1 − pi )]
i
X e ηi

1
= yi log + (1 − yi ) log
i
1 + e ηi 1 + e ηi
X
= [yi ηi − log(1 + eηi )]
i
Xh T
i
= yi · xiT β − log(1 + exi β ) .
i

Exercise: show that the log-likelihood function of logistic regression is a concave function in β. If you need a
refresher how to take derivatives with respect to a vector or matrix, see Biostat 216 notes.
• Maximization of this log-likelihood function can be carried out by the Newton-Raphson (also known as
Fisher scoring) algorithm.

`(βt )
βt+1 = βt − .
`0 (βt )
lmod <- glm(chd ~ height + cigs, family = binomial, wcgs)

Summary of the result:

(lmod_sm <- summary(lmod))

##
## Call:
## glm(formula = chd ~ height + cigs, family = binomial, data = wcgs)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0041 -0.4425 -0.3630 -0.3499 2.4357
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.50161 1.84186 -2.444 0.0145 *
## height 0.02521 0.02633 0.957 0.3383
## cigs 0.02313 0.00404 5.724 1.04e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1781.2 on 3153 degrees of freedom
## Residual deviance: 1749.0 on 3151 degrees of freedom
## AIC: 1755
##
## Number of Fisher Scoring iterations: 5
## str(lmod_sm)

7
Interpretation
• Exercise: Before we attempt to interpret the results from logistic regression, we first need to understand
how the data is transformed to (yi , xi ).
# dataframe
wcgs %>%
select(chd, height, cigs) %>%
head(10)

## # A tibble: 10 x 3
## chd height cigs
## <fct> <int> <int>
## 1 no 73 25
## 2 no 70 20
## 3 no 69 0
## 4 no 68 20
## 5 yes 70 20
## 6 no 72 0
## 7 no 72 0
## 8 no 71 0
## 9 no 72 25
## 10 no 70 0
# response
lmod$y %>% head(10)

## 1 2 3 4 5 6 7 8 9 10
## 0 0 0 0 1 0 0 0 0 0
# predictors
[Link](lmod) %>% head(10)

## (Intercept) height cigs

## 1 1 73 25
## 2 1 70 20
## 3 1 69 0
## 4 1 68 20
## 5 1 70 20
## 6 1 72 0
## 7 1 72 0
## 8 1 71 0
## 9 1 72 25
## 10 1 70 0
• How to interpret the regression coefficients in logistic regression? Remember
 
p
log = β0 + β1 · height + β2 · cigs.
1−p
The quantity
p
o=
1−p
is called odds (of an event).
Therefore β1 can be interpreted as a unit increase in x1 with x2 held fixed increases the log-odds of
success by β1 , or increase the odds of success by a factor of eβ1 .
• The gtsummary package presents regression results in a much nicer way that facilitates interpretation.
Summarize the log-odds:

8
library(gtsummary)
lmod %>%
tbl_regression() %>%
bold_labels() %>%
bold_p(t = 0.05)

Characteristic log(OR) 95% CI p-value

height 0.03 -0.03, 0.08 0.3
cigs 0.02 0.02, 0.03 <0.001

Summarize the odds:

lmod %>%
tbl_regression(intercept = TRUE, exponentiate = TRUE) %>%
bold_labels() %>%
bold_p(t = 0.05)

Characteristic OR 95% CI p-value

(Intercept) 0.01 0.00, 0.40 0.015
height 1.03 0.97, 1.08 0.3
cigs 1.02 1.02, 1.03 <0.001

• Exercise: Interpret the regression coefficients from wcgs fit.

# same as lmod$coefficients
# coef(lmod) is a named numeric vector
(beta_hat <- unname(coef(lmod)))

## [1] -4.50161397 0.02520779 0.02312740

exp(beta_hat)

## [1] 0.01109108 1.02552819 1.02339691

How to interpret the effect of a pack a day (20 cigarettes) on heart disease?
exp(beta_hat[3] * 20)

## [1] 1.588115
• Suppose the probability of success in the presence of some condition is p1 and p2 in its absence. The
relative risk or risk ratio is p1 /p2 . For example, the predicted probability of a 68in tall person who
smokes a pack (20 cigarettes) a day and who does not smoke are, respectively
(p1 <- ilogit(sum(beta_hat * c(1, 68, 20))))

## [1] 0.08907868
and
(p2 <- ilogit(sum(beta_hat * c(1, 68, 0))))

## [1] 0.05800425
Then the relative risk is

9
p1 / p2

## [1] 1.535727
• When the probability of event is very small (rare disease assumption), i.e., p1 , p2 ≈ 0, then the odds
ratio is approximately equal to the risk ratio

o1 p1 /(1 − p1 ) p1
= ≈ .
o2 p2 /(1 − p2 ) p2

Inference (analysis of deviance)

• The deviance of a logistic regression fit is
X
D = 2 [yi log yi + (1 − yi ) log(1 − yi )]
i
X
−2 [yi log pbi + (1 − yi ) log(1 − pbi )]
i
X
= −2 [yi log pbi + (1 − yi ) log(1 − pbi )] .
i

It comes from the likelihood ratio test (LRT) statistic

LΩ
2 log ,
Lω
where Ω is the full/saturated model (same number of parameters as observations) and ω is the smaller
model.
• In the model output, the residual deviance, denoted DL , is the deviance of the current model and
the null deviance, denoted DS , is the deviance of the model with just an intercept term. Assuming the
null model, the test statistic DS − DL is asymptotically distributed χ2`−s . In our case, the test statistic
is
lmod$[Link] - lmod$deviance

## [1] 32.19451
giving p-value
pchisq(lmod$[Link] - lmod$deviance, 2, [Link] = FALSE)

## [1] 1.02106e-07
Therefore our model gives a significantly better fit than the null (intercept-only) model.
• We can also test the significance of individual predictor using analysis of deviance (anova function).
For example, is height necessary in the model?
# fit a model without height
lmodc <- glm(chd ~ cigs, family = binomial, wcgs)
anova(lmodc, lmod, test = "Chi")

## Analysis of Deviance Table

##
## Model 1: chd ~ cigs
## Model 2: chd ~ height + cigs
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 3152 1750
## 2 3151 1749 1 0.92025 0.3374

10
 • Similar to linear regression, the convenience function drop1 tests each individual predictor in one shot.
drop1(lmod, test = "Chi")

## Single term deletions

##
## Model:
## chd ~ height + cigs
## Df Deviance AIC LRT Pr(>Chi)
## <none> 1749.0 1755.0
## height 1 1750.0 1754.0 0.9202 0.3374
## cigs 1 1780.1 1784.1 31.0695 2.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
• The coefficient test from summary is based on the z-value β̂j /se(β̂j ).
lmod_sm$coefficients

## Estimate Std. Error z value Pr(>|z|)

## (Intercept) -4.50161397 1.841862743 -2.444055 1.452321e-02
## height 0.02520779 0.026327385 0.957474 3.383281e-01
## cigs 0.02312740 0.004040183 5.724344 1.038339e-08
In general, deviance-based test is preferred over the z-test.

Confidence intervals
• Confidence interval can be constructed either from normal approximation

β̂j ± z α/2 se(β̂j )

tibble(
`coef` = beta_hat,
`2.5%` = beta_hat - 1.96 * lmod_sm$coefficients[, 2],
`97.5%` = beta_hat + 1.96 * lmod_sm$coefficients[, 2])

## # A tibble: 3 x 3
## coef `2.5%` `97.5%`
## <dbl> <dbl> <dbl>
## 1 -4.50 -8.11 -0.892
## 2 0.0252 -0.0264 0.0768
## 3 0.0231 0.0152 0.0310
or from profile-likelihood, reference
confint(lmod)

## Waiting for profiling to be done...

## 2.5 % 97.5 %
## (Intercept) -8.13475465 -0.91297018
## height -0.02619902 0.07702835
## cigs 0.01514949 0.03100534

Diagnostics
• There are two kinds of fitted values (or predicted values). The first is on the scale of the linear predictor,
η,

11
linpred <- predict(lmod)
linpred %>% head(10)

## 1 2 3 4 5 6 7 8
## -2.083261 -2.274521 -2.762277 -2.324936 -2.274521 -2.686653 -2.686653 -2.711861
## 9 10
## -2.108468 -2.737069
The second on the scale of response, p = logit−1 (η),
predprob <- predict(lmod, type = "response")
predprob %>% head(10)

## 1 2 3 4 5 6 7
## 0.11073449 0.09325523 0.05939705 0.08907868 0.09325523 0.06376553 0.06376553
## 8 9 10
## 0.06227708 0.10827647 0.06082112
• We compute the raw residuals
y − pb
# same as residuals(lmod, type = "response")
rawres <- lmod$y - predprob

The plot of raw residuals against the fitted values is not very informative.
wcgs %>%
mutate(rawres = rawres, linpred = linpred) %>%
ggplot() +
geom_point(mapping = aes(x = linpred, y = rawres)) +
labs(x = "Linear predictor", y = "Raw residuals")

12
 1.0

0.5
Raw residuals

0.0

−3 −2 −1
Linear predictor
We do not expect the raw residuals to have equal variance because the binary variance is p(1 − p).
• The deviance residuals are standardized residuals defined by
p
di = sign(yi − pbi ) −2[yi log pbi + (1 − yi ) log(1 − pbi )].

Note X
d2i = deviance
i

2i
P
is analogy to ib = RSS in linear regression. The term sign(yi − pbi ) ensures that di has the same
sign as the raw residual yi − pbi .
devres <- residuals(lmod)
devres %>% head(10)

## 1 2 3 4 5 6 7
## -0.4844779 -0.4424800 -0.3499548 -0.4319693 2.1782631 -0.3630133 -0.3630133
## 8 9 10
## -0.3586106 -0.4787466 -0.3542579
Sanity check:
sqrt(-2 * (lmod$y * log(predprob) + (1 - lmod$y) * log(1 - predprob))) %>%
head(10)

## 1 2 3 4 5 6 7 8
## 0.4844779 0.4424800 0.3499548 0.4319693 2.1782631 0.3630133 0.3630133 0.3586106
## 9 10
## 0.4787466 0.3542579

13
The plot of deviance residuals against the fitted values.
wcgs %>%
mutate(devres = devres, linpred = linpred) %>%
ggplot() +
geom_point(mapping = aes(x = linpred, y = devres)) +
labs(x = "Linear predictor", y = "Deviance residuals")

2
Deviance residuals

−1

−3 −2 −1
Linear predictor
Again we see the residuals are clustered into two lines: the upper one corresponding to yi = 1 and the lower
one to yi = 0. We can improve this plot by binning: divide the range of linear predictor into 100 bins of
roughly equal points and plot average residual against average linear predictors per bin.
wcgs %>%
mutate(devres = devres, linpred = linpred) %>%
group_by(cut(linpred, breaks = unique(quantile(linpred, (1:100)/101)))) %>%
summarize(devres = mean(devres),
linpred = mean(linpred)) %>%
ggplot() +
geom_point(mapping = aes(x = linpred, y = devres)) +
labs(x = "Linear predictor", y = "Binned deviance residual")

## `summarise()` ungrouping output (override with `.groups` argument)

14
 0.25
Binned deviance residual

0.00

−0.25

−0.50
−2.5 −2.0
Linear predictor
Question: is there a concern that the deviance residuals are not centered around 0? No.
• Exercise: Do similar binned plots for deviance residuals vs. height and deviance residuals vs. cigs to
check the linearity assumption.
• QQ plot is not helpful since there is no reason these deviance residuals are approximately standard
normal.
qqnorm(devres)

15
 Normal Q−Q Plot

2.0
Sample Quantiles

1.0
0.0
−1.0

−2 0 2

Theoretical Quantiles
• Half-normal plot (hat values against half-normal quantiles) can help detect unusual cases in predictor
space (high leverage cases). For logistic regression, we use the generalized hat matrix

H=W c −1 XW
c 1/2 XT (XT WX) c 1/2 ,

where  
pb1 (1 − pb1 )
W=
 .. .

c .
pbn (1 − pbn )
halfnorm(hatvalues(lmod))

16
 2527

0.020
2695
Sorted Data

0.010
0.000

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Half−normal quantiles
We see two high leverage cases, who smoke an unusual number of cigarettes per day!
wcgs %>%
slice(c(2527, 2695)) %>%
print(width = Inf)

## # A tibble: 2 x 13
## age height weight sdp dbp chol behave cigs dibep chd typechd timechd
## <int> <int> <int> <int> <int> <int> <fct> <int> <fct> <fct> <fct> <int>
## 1 52 71 168 120 80 251 A1 99 B no none 2956
## 2 47 64 158 116 76 206 A1 80 B no none 2114
## arcus
## <fct>
## 1 present
## 2 absent
• Plot of Cook distance against half-normal quantiles may reveal high influential points (high residual
combined with high leverage) cases.
halfnorm([Link](lmod))

17
 0.012
9532082

0.008
Sorted Data

0.004
0.000

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Half−normal quantiles
wcgs %>%
slice(c(953, 2082)) %>%
print(width = Inf)

## # A tibble: 2 x 13
## age height weight sdp dbp chol behave cigs dibep chd typechd timechd
## <int> <int> <int> <int> <int> <int> <fct> <int> <fct> <fct> <fct> <int>
## 1 57 63 155 128 88 196 A2 0 B yes silent 2349
## 2 49 77 210 138 86 235 B3 0 A yes angina 3048
## arcus
## <fct>
## 1 absent
## 2 absent

Goodness of fit
Hosmer-Lemeshow statistic
• The usual goodness of fit test using χ2n−q−1 asymptotic null distribution can not be applied here since
we only have a single observation for each predictor pattern. This is different from the binomial model
in next chapter. The Hosmer-Lemeshow test partitions the predicted probabilities into J bins and then
carries out a Pearson X 2 type test to assess the goodness of fit (ELMR 2.6).
• Intuitively if we divide observations into J bins according to linear predictors η, then yj /nj (observed
proportion of “successes”) for j-th bin should be close to the average predicted probabilities in that bin.
wcgs_binned <- wcgs %>%
mutate(predprob = predict(lmod, type = "response"),
linpred = predict(lmod, type = "link"),
bin = cut(linpred, breaks = unique(quantile(linpred, (1:100) / 101)))) %>%
group_by(bin) %>%
summarize(y = sum(ifelse(chd == "yes", 1, 0)),
avgpred = mean(predprob),
count = n()) %>%

18
 mutate(se_fit = sqrt(avgpred * (1 - avgpred) / count))

## `summarise()` ungrouping output (override with `.groups` argument)

wcgs_binned %>%
ggplot(mapping = aes(x = avgpred, y = y / count)) +
geom_point() +
geom_linerange(mapping = aes(ymin = y / count - 2 * se_fit,
ymax = y / count + 2 * se_fit), alpha = 0.5) +
geom_abline(intercept = 0, slope = 1) +
labs(x = "Predicted probability", y = "Observed proportion")

0.4
Observed proportion

0.2

0.0

0.050 0.075 0.100 0.125 0.150

Predicted probability
• The Hosmer-Lemeshow test formalizes this idea by testing the statistic
J
2
X (yi − mj pbj )2
XHL =
j=1
mj pbj (1 − pbj )

against the χ2 distribution with J − 1 degrees of freedom.

# Hosmer-Lemeshow test statistic
(hlstat <- with(wcgs_binned, sum((y - count * avgpred)ˆ2 / (count * avgpred * (1 - avgpred)))))

## [1] 64.87001
# J
nrow(wcgs_binned)

## [1] 52

19
# p-value
pchisq(hlstat, nrow(wcgs_binned) - 1, [Link] = FALSE)

## [1] 0.09172918
We see a moderate p-value, which indicates no lack of fit.

ROC curve
• Logistic regression is often used as a tool for classification.
• If we choose a threshold, say 0.2, then the predicted probabilities give a classification rule
(
"success" if pbi ≥ 0.2
case i is a .
"failure" if pbi < 0.2

wcgs %>%
mutate(predprob = predict(lmod, type = "response")) %>%
mutate(predout = ifelse(predprob >= 0.2, "yes", "no")) %>%
xtabs(~ chd + predout, data = .)

## predout
## chd no yes
## no 2886 11
## yes 254 3
• With this classification rule, we see the error rate is about
(11 + 254) / (2886 + 254 + 11 + 3)

## [1] 0.08402029

The sensitivity is
TP 3
= = 1.17%
TP + FN 257
and the specificity is
TN 2886
= = 99.62%
FP + TN 11 + 2886
• If we lower the threshold, then we increase the sensitivity but decrease the specificity. If we plot sensitivity
against 1-specificity by varying the threshold, then we get the receiver operating characteristic
(ROC) curve.
library(pROC)

## Type 'citation("pROC")' for a citation.

##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
lmod_roc <- roc(chd ~ predprob, wcgs)

## Setting levels: control = no, case = yes

20
## Setting direction: controls < cases
ggroc(lmod_roc)

1.00

0.75
sensitivity

0.50

0.25

0.00

1.00 0.75 0.50 0.25 0.00

specificity
The area under the curve (AUC) is a measure of the overal classification of the classifier. Larger AUC means
better classification performance.
auc(lmod_roc)

## Area under the curve: 0.6007

Alternative, Detection error trade-off (DET) curve

Model selection by AIC

• We can perform sequential search using the Akaike information criterion (AIC)

AIC = deviance + 2q.

We start from a rich enough model

wcgs <- wcgs %>%
# 1 in = 0.0254 m, 1 lb = 0.4536 kg
mutate(bmi = 703 * weight / height)
biglm <- glm(chd ~ age + height + weight + bmi + sdp + dbp + chol + dibep + cigs + arcus,
family = binomial, data = wcgs)

and then do sequential backward search using the step function

21
 Figure 1: ROC vs. DET

step(biglm, trace = TRUE, direction = "both") %>%

tbl_regression() %>%
bold_labels()

## Start: AIC=1591.05
## chd ~ age + height + weight + bmi + sdp + dbp + chol + dibep +
## cigs + arcus
##
## Df Deviance AIC
## - dbp 1 1569.1 1589.1
## - weight 1 1569.3 1589.3
## - bmi 1 1569.5 1589.5
## - height 1 1569.5 1589.5
## <none> 1569.0 1591.0
## - arcus 1 1571.1 1591.1
## - sdp 1 1576.8 1596.8
## - dibep 1 1590.4 1610.4
## - cigs 1 1592.0 1612.0
## - age 1 1593.7 1613.7
## - chol 1 1619.8 1639.8
##
## Step: AIC=1589.06
## chd ~ age + height + weight + bmi + sdp + chol + dibep + cigs +
## arcus
##
## Df Deviance AIC
## - weight 1 1569.4 1587.4
## - bmi 1 1569.5 1587.5
## - height 1 1569.5 1587.5
## <none> 1569.1 1589.1
## - arcus 1 1571.2 1589.2

22
## + dbp 1 1569.0 1591.0
## - sdp 1 1586.2 1604.2
## - dibep 1 1590.4 1608.4
## - cigs 1 1592.5 1610.5
## - age 1 1593.7 1611.7
## - chol 1 1619.9 1637.9
##
## Step: AIC=1587.36
## chd ~ age + height + bmi + sdp + chol + dibep + cigs + arcus
##
## Df Deviance AIC
## - height 1 1570.3 1586.3
## <none> 1569.4 1587.4
## - arcus 1 1571.5 1587.5
## + weight 1 1569.1 1589.1
## + dbp 1 1569.3 1589.3
## - bmi 1 1574.4 1590.4
## - sdp 1 1586.7 1602.7
## - dibep 1 1590.7 1606.7
## - cigs 1 1592.8 1608.8
## - age 1 1594.0 1610.0
## - chol 1 1620.0 1636.0
##
## Step: AIC=1586.32
## chd ~ age + bmi + sdp + chol + dibep + cigs + arcus
##
## Df Deviance AIC
## <none> 1570.3 1586.3
## - arcus 1 1572.6 1586.6
## + height 1 1569.4 1587.4
## + weight 1 1569.5 1587.5
## + dbp 1 1570.3 1588.3
## - bmi 1 1577.3 1591.3
## - sdp 1 1587.2 1601.2
## - dibep 1 1591.9 1605.9
## - age 1 1594.2 1608.2
## - cigs 1 1594.6 1608.6
## - chol 1 1620.0 1634.0
## Table printed with `knitr::kable()`, not {gt}. Learn why at
## [Link]
## To suppress this message, include `message = FALSE` in code chunk header.

Characteristic log(OR) 95% CI p-value

age 0.06 0.04, 0.08 <0.001
bmi 0.00 0.00, 0.00 0.008
sdp 0.02 0.01, 0.03 <0.001
chol 0.01 0.01, 0.01 <0.001
dibep
A
B 0.66 0.38, 1.0 <0.001
cigs 0.02 0.01, 0.03 <0.001
arcus
absent

23
 Characteristic log(OR) 95% CI p-value
present 0.22 -0.06, 0.50 0.13

Model selection by lasso

• A well-known approach for model selection is the lasso, which minimizes the function
q
X
−n−1 `(β) + λ |βj |,
j=1

where `(β) is the log-likelihood and λ > 0 is a tuning parameter. We notice that
– when λ = ∞, all non-intercept regression coefficients will be pushed to 0, and

– when λ = 0, the regression coefficients are same as those from regular logistic regression.
If we vary λ from 0 to larger values, we will obtain intermediate models with lesser and lesser non-zero
predictors. This way we are achieving continuous model selection, i.e. solution path.
• For details of glmnet package, see the vignette at [Link]
[Link]
• How do we choose λ, which determines the model size? One natural idea is to split the data into a
training set and a validation set. The training set is used to fit the logistic regression at different λ
values. Then the validation set is used to evaluate and compare the performance of different models.
We will choose the model that gives the best performance on the validation set.
• First let’s remove cases with missing values
(wcgs <- wcgs %>%
select(-c(behave, typechd, timechd)) %>%
drop_na())

## # A tibble: 3,140 x 11
## age height weight sdp dbp chol cigs dibep chd arcus bmi
## <int> <int> <int> <int> <int> <int> <int> <fct> <fct> <fct> <dbl>
## 1 49 73 150 110 76 225 25 B no absent 1445.
## 2 42 70 160 154 84 177 20 B no present 1607.
## 3 42 69 160 110 78 181 0 A no absent 1630.
## 4 41 68 152 124 78 132 20 A no absent 1571.
## 5 59 70 150 144 86 255 20 A yes present 1506.
## 6 44 72 204 150 90 182 0 A no absent 1992.
## 7 44 72 164 130 84 155 0 A no absent 1601.
## 8 40 71 150 138 60 140 0 B no absent 1485.
## 9 43 72 190 146 76 149 25 A no absent 1855.
## 10 42 70 175 132 90 325 0 B no present 1758.
## # ... with 3,130 more rows
split data into 80% training cases and 20% validation cases.
library(glmnet)

## Loading required package: Matrix

##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':

24
##
## expand, pack, unpack
## Loaded glmnet 4.0-2
library(caret)

## Loading required package: lattice

##
## Attaching package: 'lattice'
## The following object is masked from 'package:faraway':
##
## melanoma
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
# set seed for reproducibility
[Link](200)

# list = FALSE, request result to be in a matrix (of row position) not a list
training_samples <- wcgs$chd %>%
createDataPartition(p = 0.8, list = FALSE)
# (train_data <- wcgs %>%
# slice(training_samples[, 1]))
# (val_data <- wcgs %>%
# slice(-training_samples[, 1]))

The glmnet package takes a matrix (of predictors) and a vector (of responses) as input. We use [Link]
function to create them. glmnet will add intercept by default, so we drop intercept term when forming x
matrix.
# X and y from original data
x_all <- [Link](
chd ~ - 1 + age + height + weight + bmi + sdp + dbp + chol + dibep + cigs + arcus,
data = wcgs)
y_all <- ifelse(wcgs$chd == "yes", 1, 0)
# training X and y
x_train <- x_all[training_samples[, 1], ]
y_train <- y_all[training_samples[, 1]]
# validation X and y
x_val <- x_all[-training_samples[, 1], ]
y_val <- y_all[-training_samples[, 1]]

Fit lasso regression and plot solution path:

lasso_fit <- glmnet(x_train, y_train, family = "binomial", alpha = 1)
summary(lasso_fit)

## Length Class Mode

## a0 49 -none- numeric
## beta 539 dgCMatrix S4
## df 49 -none- numeric

25
## dim 2 -none- numeric
## lambda 49 -none- numeric
## [Link] 49 -none- numeric
## nulldev 1 -none- numeric
## npasses 1 -none- numeric
## jerr 1 -none- numeric
## offset 1 -none- logical
## classnames 2 -none- character
## call 5 -none- call
## nobs 1 -none- numeric
plot(lasso_fit, xvar = "lambda", label = TRUE)

10 9 8 6 0
0.2

11

1
0.0

10
2
5
7
3
6
9
Coefficients

−0.2
−0.4
−0.6

−7 −6 −5 −4 −3

Log Lambda
• Here
– α = 1 corresponds to the lasso regression, in the family of elastic net penalties

(1 − α)kβk22 /2 + αkβk1 .

– Other choices for xvar are "lambda" for log lambda value, "norm" for the `1 -norm of the coefficients
(default), and "dev" for the percentage of deviance explained.
• Now we can evaluate the performance of the models (corresponding to different λ values) on the
validation set.
# predict validation case probabilities at different \lambda values and calculate test deviance
pred_val <- predict(lasso_fit, newx = x_val, type = "response", s = lasso_fit$lambda)
dev_val <- -2 * colSums(y_val * log(pred_val) + (1 - y_val) * log(1 - pred_val))
tibble(lambda = lasso_fit$lambda, dev_val = dev_val) %>%
ggplot() +
geom_point(mapping = aes(x = lambda, y = dev_val)) +
scale_x_log10() +

26
 labs(y = "Binomial deviance on validation set", x = "Lambda")

350
Binomial deviance on validation set

340

330

0.001 0.003 0.010 0.030

Lambda
From the graph, it seems that λ = 0.005 yields a model that performs best on the validation set.
• Now the question is whether we should make decision on this single training-validation split? A common
strategy is to use cross validation.
[Link](200)
cv_lasso <- [Link](x_all, y_all, alpha = 1, family = "binomial",
[Link] = "auc")
plot(cv_lasso)

27
 10 9 9 9 8 8 8 8 8 8 8 7 6 6 6 2 0

0.75
0.70
AUC

0.65
0.60
0.55

−7 −6 −5 −4 −3

Log(λ)
The plot displays the cross-validation error (devience by default, we chose AUC here) according to the log of
lambda. The left dashed vertical line indicates that the log of the optimal value of log lambda is approximately
-5.2, which is the one that minimizes the average AUC. This lambda value will give the most accurate model.
The exact value of lambda and corresponding model can be viewed as follow:
cv_lasso$[Link]

## [1] 0.005696827
coef(cv_lasso, cv_lasso$[Link])

## 12 x 1 sparse Matrix of class "dgCMatrix"

## 1
## (Intercept) -9.893281e+00
## age 4.954011e-02
## height .
## weight 5.334625e-03
## bmi .
## sdp 1.537576e-02
## dbp .
## chol 9.228292e-03
## dibepA -5.162308e-01
## dibepB 3.010844e-14
## cigs 1.656021e-02
## arcuspresent 9.149817e-02
cv_lasso$lambda.1se

## [1] 0.0190935
coef(cv_lasso, cv_lasso$lambda.1se)

## 12 x 1 sparse Matrix of class "dgCMatrix"

## 1

28
## (Intercept) -6.202330e+00
## age 2.267635e-02
## height .
## weight .
## bmi .
## sdp 9.964743e-03
## dbp .
## chol 6.236502e-03
## dibepA -2.217104e-01
## dibepB 2.545112e-14
## cigs 6.193633e-03
## arcuspresent .
We see that this model differs from the best model chosen by AIC by replacing the predictor bmi by weight.

29