In binary logistic regression, there is a single binary dependent variable, coded by an indicator variable. For example, if we represent a response as 1 and non-response as 0, then the corresponding probability of response, can be between 0 (certainly not a response) and 1 (certainly a response) - hence the labeling !
The logistic model models the log-odds of an event as a linear combination of one or more independent variables (explanatory variables). If we observed \((y_i, x_i),\) where \(y_i\) is a Bernoulli variable and \(x_i\) a vector of explanatory variables, the model for \(\pi_i = P(y_i=1)\) is
The model is especially useful in case-control studies and leads to the effect of risk factors by odds ratios.
Example: Lung Cancer Data
Data source: Loprinzi CL. Laurie JA. Wieand HS. Krook JE. Novotny PJ. Kugler JW. Bartel J. Law M. Bateman M. Klatt NE. et al. Prospective evaluation of prognostic variables from patient-completed questionnaires. North Central Cancer Treatment Group. Journal of Clinical Oncology. 12(3):601-7, 1994.
wgt_catn consists of: 1= patients a weight loss of zero or less, 0= patients with a weight loss of more than zero
trt01pn consists of 1= active treatment, 0 = placebo
Model Fit
We analyze the event of weight gain (or staying the same weight) in lung cancer patients in dependency of treatment (active or placebo), age, sex, ECOG performance score and calories consumed at meals. One of the most important things to remember is to ensure you tell R what your event is and what treatment comparison you are doing Active / Placebo or Placebo/Active! The easiest way to do this is to have event (or non-reference treatment) as 1, and non-event (reference treatment) as 0. Hence we will use the wt_catn and trt01pn variables.
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss trt01p
1 3 306 2 74 1 1 90 100 1175 NA Active
2 3 455 2 68 1 0 90 90 1225 15 Active
3 3 1010 1 56 1 0 90 90 NA 15 Active
4 5 210 2 57 1 1 90 60 1150 11 Active
5 1 883 2 60 1 0 100 90 NA 0 Active
6 12 1022 1 74 1 1 50 80 513 0 Active
trt01pn dose_mg dose_id wt_cat wt_catn
1 1 10 1 <NA> NA
2 1 10 1 loss 0
3 1 10 1 loss 0
4 1 10 1 loss 0
5 1 10 1 gain 1
6 1 10 1 gain 1
m1 <-glm(wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, data = lung, family =binomial(link="logit"))summary(m1)
Call:
glm(formula = wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal,
family = binomial(link = "logit"), data = lung)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.8623998 1.7675776 -2.185 0.0289 *
trt01pn 0.3887667 0.3781566 1.028 0.3039
age 0.0122549 0.0211553 0.579 0.5624
sex 0.8321005 0.3743793 2.223 0.0262 *
ph.ecog -0.3763592 0.2638322 -1.427 0.1537
meal.cal 0.0008500 0.0004486 1.895 0.0581 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 202.36 on 169 degrees of freedom
Residual deviance: 190.46 on 164 degrees of freedom
(58 observations deleted due to missingness)
AIC: 202.46
Number of Fisher Scoring iterations: 4
The model summary contains the parameter estimates \(\beta_j\) for each explanatory variable \(x_j\), corresponding to the log-odds for the response variable to take the value \(1\), conditional on all other explanatory variables remaining constant. For better interpretation, we can exponentiate these estimates, to obtain estimates for the odds instead and provide 95% confidence intervals.
# model coefficients summarysummary(m1)$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.8623998098 1.7675776454 -2.1851373 0.02887878
trt01pn 0.3887666677 0.3781565596 1.0280574 0.30392281
age 0.0122549015 0.0211552875 0.5792831 0.56239813
sex 0.8321005169 0.3743792762 2.2226137 0.02624186
ph.ecog -0.3763592487 0.2638321918 -1.4265100 0.15372119
meal.cal 0.0008499918 0.0004486401 1.8945961 0.05814593
est 2.5 % 97.5 %
(Intercept) 0.0210175 0.000580124 0.6197446
trt01pn 1.4751603 0.696210092 3.0850890
age 1.0123303 0.971670916 1.0561942
sex 2.2981410 1.107651762 4.8367698
ph.ecog 0.6863557 0.405659156 1.1474521
meal.cal 1.0008504 0.999978126 1.0017613
NOTE: that there are two types of confidence intervals that you can calculate. Function confint.default gives the Wald confidence limits, which is the default option in SAS PROC LOGISTIC procedure; whereas confint gives the profile-likelihood limits.
Model Comparison
To compare two logistic models, the residual deviances (-2 * log likelihoods) are compared against a \(\chi^2\)-distribution with degrees of freedom calculated using the difference in the two models’ parameters. Below, the only difference is the inclusion/exclusion of age in the model, hence we test using \(\chi^2\) with 1 df. Here testing at the 5% level.
m2 <-glm(wt_catn ~ trt01pn + sex + ph.ecog + meal.cal, data = lung, family =binomial(link="logit"))summary(m2)
Call:
glm(formula = wt_catn ~ trt01pn + sex + ph.ecog + meal.cal, family = binomial(link = "logit"),
data = lung)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.9951557 0.9245655 -3.240 0.0012 **
trt01pn 0.3681563 0.3761976 0.979 0.3278
sex 0.7919227 0.3674936 2.155 0.0312 *
ph.ecog -0.3312031 0.2527586 -1.310 0.1901
meal.cal 0.0007896 0.0004378 1.803 0.0713 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 202.36 on 169 degrees of freedom
Residual deviance: 190.80 on 165 degrees of freedom
(58 observations deleted due to missingness)
AIC: 200.8
Number of Fisher Scoring iterations: 4
anova(m1, m2, test ="LRT")
Analysis of Deviance Table
Model 1: wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal
Model 2: wt_catn ~ trt01pn + sex + ph.ecog + meal.cal
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 164 190.46
2 165 190.80 -1 -0.33867 0.5606
Stackexchange here has a good article describing this method and the difference between comparing 2 models using the likelihood ratio tests versus using wald tests and Pr>chisq (from the maximum likelihood estimate). Note: anova(m1, m2, test = "Chisq") and using test="LRT" as above are synonymous in this context.
Prediction
Predictions from the model for the log-odds of a patient with new data to experience a weight loss are derived using predict():
# new female, symptomatic but completely ambulatory patient consuming 2500 caloriesnew_pt <-data.frame(trt01pn=1, sex=2, ph.ecog=1, meal.cal=2500)predict(m2, new_pt, type ="response")
1
0.6455793
Contrast statements for 2 or more treatments
To create contrasts, you can use the fit.contrast() function from the gmodels package.
This can be used with lm and glm objects:
Here if we use the 3 level treatment variable (dose_id), we have 1=10mg, 2=20mg doses for active treatment and then 3= placebo which is 0mg.
You would fit the model as above, followed by fit.contrast(). This is effective testing the null hypothesis that 0.5dose10mg + 0.5 dose20mg - placebo = 0.
lung2<-lung |>mutate(dose_id2 =as.factor(lung$dose_id))m3 <-glm (wt_catn ~ dose_id2 + age + sex + ph.ecog + meal.cal, data = lung2, family =binomial(link="logit"))fit.contrast (m3,'dose_id2',c(0.5,0.5,-1),conf.int=0.95)
Estimate Std. Error z value Pr(>|z|) lower CI
dose_id2 c=( 0.5 0.5 -1 ) -0.4096323 0.3801683 -1.077502 0.2812558 -1.160322
upper CI
dose_id2 c=( 0.5 0.5 -1 ) 0.3410574
Reference
#\| echo: false#List all the packages neededsi <- sessioninfo::session_info(c('tidyverse','gmodels')) si