Logistic Regression - II
Hello dear readers and welcome back to my blog with a new post! Although I am on leave, I would like to keep updating my posts to follow-up from the topic introduced in my last post about logistic regression (also because I feel a little guilty for the relatively small introduction I provided last time). Moving from the introductory concepts dealt with in the last post, today I expand the topic to deal with more advanced elements such as hypothesis testing and more complex model specifications.
As per usual, feel free to let me know if you have any questions on the material and/or you have any suggestions to better introduce the concepts or examples. I will use any type of constructive feedback to improve my educational material for my future students!
Hypothesis Testing
There are different types of tests that can be used to assess the possible association between predictor variables and the binary outcome variable of logistic regression. Common examples include:
the Wald Test: for a binary predictor \(X_1\), the tested hypotheses are \(H_0: \beta_1=0\) vs \(H_1: \beta_1 \neq 0\), which are equivalent to testing whether \(\text{exp}(\beta_1)=1\) under \((H_0)\) vs \(\text{exp}(\beta_1)\neq 1\) under \((H_1)\).
the Likelihood Ratio Test (LRT): for a group of indicator predictors \(X_2,X_3,\ldots\) (generated from a categorical predictor), the tested hypotheses are \(H_0:\beta_2=\beta_3=\ldots=0\) vs \(H_1:\) at least one coefficient among \(\{\beta_2,\beta_3,\ldots\}\) is \(\neq 0\).
Wald Test
The test statistic is the ratio of the coefficient to its standard error, which is assumed to be asymptotically Normally distributed under \(H_0\):
\[ \frac{\hat{\beta}_1}{\text{SE}(\hat{\beta}_1)} \sim N(0,1). \] For instance, considering the example of surgical mortality introduced in the last post, where death status (dead/not dead) was regressed on past hart attack (yes/no - HA) and ejection fraction status (good/fair/poor - EFF/EFP), the logistic regression model is:
\[ \log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_0 + \beta_1 \times \text{HA} + \beta_2 \times \text{EFF} + \beta_3 \times \text{EFP}. \] Using the estimates \(\hat{\beta}_1=0.145\) and \(\text{SE}(\hat{\beta}_1)=0.323\), it is possible to derive the test statistic \(\frac{\hat{\beta}_1}{\text{SE}(\hat{\beta}_1)}=0.45\) and associated two-sided p-value of \(0.65\).
Likelihood Ratio Test
The LRT can be used to compare (log) likelihoods between two nested models. For instance, in our example, we wish to test the effect of ejection fraction and fit two competing models to obtain the associated likelihoods \(L_1\) and \(L_2\):
model 1: ejection fraction and past heart attack are included into the model as predictors (for a total of \(p=4\) coefficients)
model 2: only past heart attack is included into the model as predictor (for a total of \(q=2\) coefficients)
The test statistic is calculated as the (log) ratio of the likelihoods of the two models, re-scaled on the so-called deviance scale (i.e. multiplied by \(-2\)) which, under \(H_0\), is assumed to follow a \(\chi^2\) distribution with \(\nu=p-q\) degrees of freedom:
\[ -2 \log\Bigg(\frac{L_1}{L_2}\Bigg) = -2 (\log L_2 - \log L_1) \sim \chi^2_{p-q}. \] For instance, in our example, the log-likelihood associated with the two models are \(\log L_1 = -206.41808\) and \(\log L_2=-217.25855\), with an associated LRT statistic of \(-2(-217.25855 + 206.41808)=21.6809\) following a \(\chi^2\) with \(\nu=2\) degrees of freedom under \(H_0\) and an associated p-value of \(\approx 0\).
The two tests often give the same results for large data sets, but the LRT is more reliable for small data sets with irregular shaped likelihoods.
Continuous Predictors and Fitted Values
It is relatively simple to include continuous predictors into logistic regression models. For instance, in our example, we can expand model 1 to further include age (in years) as a new predictor:
\[ \log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_0 + \beta_1 \times \text{HA} + \beta_2 \times \text{EFF} + \beta_3 \times \text{EFP} + \beta_4 \times \text{age}. \tag{1}\]
The interpretation of the coefficient \(\beta_4\) is a little different with respect to binary predictors: \(\beta_4\) represents the increase in log-odds for a one unit (year) increase in age, while \(\text{exp}(\beta_4)\) is the odds ratio for a one year increase in age, adjusted for past hart attack and ejection fraction.
After fitting Equation 1, the estimated value for the age coefficient is \(\hat{\beta_4}=0.05827\), so that \(\text{exp}(\hat{\beta}_4)=1.060\) with associated \(95\%\) CI of \((1.024;1.098)\). We can interpret this results as follows: there is \(6\%\) increase in the odds of surgical mortality for a one year increase in age. The associated Wald test statistic is \(\frac{\hat{\beta}_4}{\text{SE}(\hat{\beta}_4)}=\frac{1.060}{0.01888}=3.29\) with a two-sided p-value of \(\approx 0\).
It is also possible to use the model to generate fitted values (predictions) for specific individual profiles in terms of the predictor variables. For example, for a \(60\) year old patient with a past heart attack and poor ejection fraction, the predicted log odds is:
\[ \log\Bigg(\frac{\hat{\pi}}{1-\hat{\pi}}\Bigg)=-7.91 + 0.20 \times 1 + 0.80 \times 0 + 1.64 \times 1 + 0.0586 \times 60 = -2.55, \] and a predicted risk of \(\hat{\pi}=\frac{\text{exp}(-2.55)}{1+\text{exp}(-2.55)}=0.073\).
A possible drawback for this type of regression model rests on its assumption that the log-odds of mortality increases linearly with age. If this assumption is difficult to justify, the association between the variables may be incorrectly modelled, the adjustment for other predictors may be only partial, and predictions may be inaccurate.
Goodness of Fit
It is possible to check the fit of the logistic regression model by comparing the observed outcomes to the predicted values using the Hosmer-Lemeshow goodness of fit test. First, we assign each patient to one of \(10\) equal-sized “risk groups” based on their predicted values with the aim to group patients with similar predictions: \(10\%\) of patients (\(n_1\)) with the lowest predictions (\(G_1\)); \(10\%\) of patients (\(n_2\)) with the next lowest predictions (\(G_2\)); \(\ldots\); \(10\%\) of patients (\(n_{10}\)) with the highest predictions (\(G_{10}\)). The test statistic is
\[ \text{HL}=\sum_{j=1}^{10}\frac{(o_j-e_j)^2}{e_j\times (1-e_j/n_j)}, \] where \(o_j=\sum_{i \in G_j}y_i\) is the observed number of surgical deaths, \(e_j=\sum_{i \in G_j} \hat{\pi}_i\) is the expected number of surgical deaths, \(n_j\) is the size of group \(j\). The test statistic is assumed to follow a \(\chi^2_{\nu}\) distribution with \(\nu=J-2=8\) degrees of freedom if the model is correctly specified.
?@tbl-ex1 shows an example of the data needed to perform the test.
| group | obs.y | exp.y | n |
|---|---|---|---|
| 1 | 2 | 0.9 | 152 |
| 2 | 2 | 1.6 | 155 |
| - | - | - | - |
| 9 | 7 | 8.8 | 149 |
| 10 | 19 | 16.9 | 149 |
For instance, group \(2\) contributes \(\frac{(2-1.6)^2}{1.6\times (1-1.6/155)}=0.11\) to the test statistic. The HL test statistic in this example is \(\text{HL}=4.49\) with an associated p-value of \(0.8104\) suggesting no evidence of a poor fit. Alternatively, a scatter plot of the observed vs expected values can also be used to assess a possible misfit of the model.
Interactions
The relationship between a predictor and the outcome may be modified by another predictor, which requires the inclusion into the model of an interaction term, often used in clinical trials to investigate whether the intervention works for different subgroups. Let us consider the hip dysplasia example where infants were categorised by the level of “clinical suspicion” (strong/moderate) regarding the risk of hip dysplasia:
- \(61\) infants had “moderate suspicion” and received ultrasound treatment
- \(73\) infants had “moderate suspicion” and did not receive ultrasound treatment
- \(65\) infants had “strong suspicion” and received ultrasound treatment
- \(86\) infants had “strong suspicion” and did not receive ultrasound treatment
The odds ratios for each “suspicion” group are:
\(\text{OR}_{\text{mod}}=\frac{61/158}{73/153}=0.81\)
\(\text{OR}_{\text{str}}=\frac{65/30}{86/3}=0.08\)
The question becomes on whether ultrasound treatment is more effective in the “strong” vs “moderate” suspicion subgroup. A possible approach to address this problem is to fit a logistic regression model to each subgroup separately, although this makes it difficult to perform a simple test (e.g. Wald test) to test the null hypothesis \(H_0: \beta_{\text{mod}}=\beta_{\text{str}}\). Alternatively, we could use interaction terms to fit a single model to the entire sample. This is achieved by creating an indicator variable for “clinical suspicion” (\(1\) for strong suspicion and \(0\) for moderate suspicion). We then create an interaction between ultrasound and suspicion by multiplying the two associated indicators: \(\text{ultrasound} \times \text{suspicion} = 1\) if an infant received ultrasound treatment and had a strong suspicion or \(0\) otherwise. ?@tbl-ex2 shows an example of the re-coded data.
| id | ultrasound | suspicion | ultra.susp | treatment |
|---|---|---|---|---|
| 1 | 0 | 0 | 0 | 0 |
| 2 | 1 | 1 | 1 | 1 |
| 3 | 1 | 0 | 0 | 0 |
| 4 | 1 | 1 | 1 | 1 |
| 5 | 1 | 0 | 0 | 1 |
| - | - | - | - | - |
We then fit the model
\[ \log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_0 + \beta_1 \times \text{ultrasound} + \beta_2 \times \text{suspicion} + \beta_3 \times \text{ultrasound} \times \text{suspicion}, \] where:
for the “moderate” subgroup (\(\text{suspicion}=0\)), the model is \(\log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_0 + \beta_1 \times \text{ultrasound}\), and \(\beta_1\) represents the ultrasound coefficient.
for the “strong” subgroup (\(\text{suspicion}=1\)), the model is \(\log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_0 + \beta_1 \times \text{ultrasound} + \beta_2 + \beta_3 \times \text{ultrasound}\), and \(\beta_1+\beta_3\) represents the ultrasound coefficient.
After fitting the model, the estimated odds ratios for ultrasound treatment (yes vs no) can be easily derived as:
for the “moderate” subgroup: \(\text{exp}(\hat{\beta}_1)=\text{exp}(-0.21)=0.81\)
for the “strong” subgroup: \(\text{exp}(\hat{\beta}_1 + \hat{\beta}_3)=\text{exp}(-0.21 - 2.37)=0.08\)
Using this approach, we can formally test for a subgroup by testing the interaction term: \(H_0: \beta_3=0\) vs \(\beta_3 \neq 0\). The Wald test statistic is:
\[ \frac{\hat{\beta}_3}{\text{SE}(\hat{\beta}_3)}=\frac{-2.3708}{0.66078}=-3.5879, \] with an associated two-sided p-value of \(\approx 4\times 10^{-4}\), thus suggesting enough counter evidence against \(H_0\).
Case-Control Studies
Logistic regression is suitable to analyse case-control studies. For instance, consider a study aimed at identifying risk factors for low birth weight (\(<2.5\) kg) which collected data on \(189\) women, of whom \(59\) had low birth weight babies (cases) and \(130\) had normal weight babies (controls). Four potential risk factors were considered: age in years, maternal weight before pregnancy in kg (mwt), race (white=1,black=2,3=other), and number of doctor visits in first trimester (visits).
The binary outcome is represented by the case (1) vs control (0) status (i.e. low birth weight), with a possible logistic regression model to consider being:
\[ \log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_0 + \beta_1 \times \text{age} + \beta_2 \times \text{mwt} + \beta_3 \times \text{rBlack} + \beta_4 \times \text{rOther} + \beta_5 \times \text{visits}. \] Estimates of regression coefficients can be interpreted in the usual way. For instance, an estimate of the mwt coefficient of \(\hat{\beta}_2=-0.031286\) corresponds to an OR of \(\text{exp}(\hat{\beta}_2)=0.97\), thus suggesting that the odds of a low birth weight baby decreases by \(\approx 3\%\) for every 1 kg increase in maternal weight. However, we cannot interpret the intercept (\(\hat{\beta}_0\)) nor make predictions \(\hat{\pi}\) based on this model.
Let us consider a new example, where a matched case-control study was conducted to investigate the relationship between presenile dementia and family history: \(109\) patients with dementia (cases) were identified and matched by age and sex with community patients without dementia (controls), while information on family history was recorded for all subjects. ?@tbl-ex3 shows a snapshot of this data set.
| Matched Group | famHisotry | Dementia |
|---|---|---|
| 1 | 0 | 0 |
| 1 | 0 | 1 |
| 2 | 0 | 0 |
| 2 | 1 | 1 |
| - | - | - |
A suitable regression model to analyse these data is:
\[ \log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_{0,1} + \ldots + \beta_{0,109} + \beta_1 \times \text{famHistory}, \] where a separate intercept is required for each matched pair. However, logistic regression should not be used to estimate this model as estimation problems are likely to arise due to the large number of parameters. Instead, we use conditional logistic regression to estimate this model, which relies on conditional likelihood not including the \(\beta_{0,j}\) terms and whose results can be interpreted in the same way as logistic regression (but intercepts are not estimated).
Conclusions
I hope that, following this introduction to logistic regression, those of you who were not familiar with the method have learned a bit more about:
using logistic or conditional logistic regression models to answer different research questions for different study designs.
model binary, categorical and continuous predictors within logistic regression models
model non-linearity and interaction terms within logistic regression models
assess the fit of logistic regression models
calculate the likelihood for a logistic regression model
Appendix
Linear Regression with Binary Outcomes
There are three main problems when using linear regression, e.g. \(\pi = \beta_0 + \beta_1 \times x\), to estimate a risk parameter \(\pi\) associated with a binary outcome variable:
The fitted values (predictions) from the model are not bounded by \([0,1]\) and impossible estimates of \(\pi\) may be obtained
The variance term is not constant as it depends on the mean parameter, i.e. for a binary variable \(\text{E}[Y]=\pi\) and \(\text{Var}(Y)=\pi \times (1-\pi)\).
Normal distribution theory, i.e. \(Y \sim N(\mu,\sigma^2)\), does not apply and CI and p-values obtained may be incorrect
Investigating Linearity: Quadratic Term
There are many approaches to investigate non-linearity, such as adding a quadratic term to the model along with the linear term:
\[ \log\Bigg(\frac{\pi}{1-\pi}\Bigg)=\beta_{0} + \ldots + \beta_{p} \times \text{age} + \beta_{p+1} \times \text{age}^2, \] where a test on \(\beta_{p+1}\) may be used to check if there is evidence of non-linearity (comparing the model with vs without the non-linear term). For instance, the LRT may be used to assess significance since linear and quadratic terms are often strongly correlated. If evidence of non-linearity is present, we should include an appropriate term into the model.
Modelling Continuous Predictors
There are many approaches to modelling continuous predictors and checking for non-linearity. Typical examples include:
Standard predictor transformation (e.g. logarithm, square root, reciprocal), which requires caution in presence of zero values, or fractional polynomials for more complicated transformations (Royston and Altman (1994)).
Categorised predictors, where the categorical terms can then be tested statistically.