knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This vignette details why the LOS_model
dataset was created, how to load it, and gives examples of use for learning/teaching regression modelling using Generalized Linear Models (GLMs), and related techniques.
The data were created specifically for regression tutorials, simulating a small set of data on hospital in-patient spells. The data are 10 sets of 30 simulated patient records, representing 10 different hospitals ("Trusts"). The dataset contains:
library(NHSRdatasets) library(dplyr) library(ggplot2) data("LOS_model") head(LOS_model) summary(LOS_model) # 82.3% survived prop.table(table(LOS_model$Death))
Now lets look at the distributions of Age and LOS. We might expect Age to be normally distributed, but this would be unrealistic, as elderly patient represent higher proportions of the total in-patients. It could be uniform, it could be bi-modal (or have more peaks), but nothing is particularly clear. LOS is highly right-skewed. This means that the high values will affect the mean, dragging it out, and the median would be a better estimation of the centre of the distribution. This shows that most patients don't stay very long, but a few patients stay notably longer.
```{R VisageLOS} ggplot(LOS_model, aes(x=Age)) + geom_histogram(alpha=0.5, col=1, fill=12, bins=20)+ ggtitle("Distribution of Age") ggplot(LOS_model, aes(x=LOS)) + geom_histogram(alpha=0.5, col=1, fill=13, bins=20)+ ggtitle("Distribution of Length-of-Stay")
## Modelling LOS and Death We will attempt to model LOS, using Age and Death, and also to model Death using LOS and Age. This vignette does not describe linear models per se, but assumes you are familiar with linear regression. If not, pause here and do a little reading before you continue. Our data are not continuous, linear, or normally distributed. Death is binary, LOS is a count. We can extend the linear model idea for this as a Generalized Linear Model (GLM)<sup>[1]</sup> , by fitting our model via a 'link function.' This function transforms data and allows our response, and it's variance, to relate to a linear model through the link function: $$\large{g(\mu)= \alpha + \beta x}$$ Where $\mu$ is the expectation of Y, and $g$ is the link function <br> The Ordinary Least Squares (OLS) fitting technique used in linear regression, is not a good fit here as our data are not necessarily distributed like that. Here we can use 'maximum-likelihood estimation' (MLE). MLE is related to probability and can be used to identify the parameters that were most likely to have produced our data (the largest likelihood value, or 'maximum likelihood'). MLE is an iterative process, and is actually performed on the log-likelihood. The model is iteratively refitted to maximise the log-likelihood. Our model is said to 'converge' when we can't improve it any more (usually changes in the log-likelihood < 1e-8). __If you get convergence-related error messages, it means the modelling function can't settle on, or reach the MLE. This often means that the scale of you model predictors is mismatched (e.g. one is binary (0,1), and another is on log scale), or there correlations in your data.__ <br> Although many of the methods for `lm` are common to `glm`, inspecting plots of residuals is trickier. We also can't use R<sup>2</sup> either, as the assumptions reflect OLS. Various methods other methods exist, depend on what you want to compare, including: - <b>AUC:</b> For binary (or other multiple categorical models), the 'area under the receiver operator characteristic curve' ('AUC', 'ROC' or 'C-statistic') can be used and interpreted like R<sup>2</sup>, as the proportion of variation in $y$ explained by our model. - <b>Likelihood Ratio test (LRT):</b> Two `glm`s can be compared using likelihood ratio tests, that compare the log-likelihoods between two models. We test a reduced model against the larger model, with the null hypothesis that the two likelihoods are not significantly different. If our test is positive (usually p<0.05, but beware of p-values...), our models are significantly different, with the one with higher log-likelihood the better model. - <b>Akaike Information Criterion (AIC) <sup>[2]</sup>:</b> is a measure of relative information loss. The value cannot be directly interpreted, but can be compared between models, with lower AICs suggesting less information loss, and being 'better' models. - <b>Prediction error:</b> If your model is primarily about prediction, and less about explanation<sup>[3]</sup>, measures of fit such as Mean-Square Error of Mean Absolute error can be helpful measures of predictive accuracy. ## Generalized Linear Models Let's fit a generalized linear model, firstly for `Death`, and secondly for `LOS`. ### Modelling death: `Death` is a binary variable recorded in our `LOS_model` dataset. This cannot be linear and is usually modelled using the `binomial` family argument, and the `logit` link function. This link function is the default in `R` for `binomial` regressions, and is commonly referred to as __Logistic Regression__. We'll fit this using the `glm` function. Rather than fitting `Death` as such, the link function makes this the log-odds of death. ```r glm_binomial <- glm(Death ~ Age + LOS, data=LOS_model, family="binomial") summary(glm_binomial) ModelMetrics::auc(glm_binomial)
Using the summary
function, we can see our model output, including:
*
, to indicate the significance levels.ModelMetrics
package to calculate the auc
. This is suggests that ~66% of the variation in our data can be explained by our model.If we consider the logic of the model above, it is telling us that LOS is a significant predictor for the log-odds of death, but that age is not. This seems unlikely, unless Age and LOS are not independent. If this is the case, e.g. elderly patients stay longer or that shorter staying elderly patients are more likely to have died, the effects of LOS may be obscuring the effects of Age. This is referred to as confounding. We can examine this using an interaction term.
'Interactions' are where predictor variables affect each other. For example, the Hospital Standardised Mortality Ratio (HSMR) ^{[4]} uses an interaction term between co-morbidity and age. So co-morbidities score may have different effects related to age, as well as the independent effects of age and co-morbidity score.
We can add these to our models with *
to include the interaction and the independent terms, or :
if you just want the interaction without the independent terms:
glm_binomial2<- glm(Death ~ Age + LOS + Age*LOS, data=LOS_model, family="binomial") summary(glm_binomial2) ModelMetrics::auc(glm_binomial2)
Now our model is suggesting that Age, LOS and their interaction are significant. This suggests that, as well as their own effects, their combination is important, and their effects confound each other.
Now we have built a model, we can use it to predict our expected $y$, the probability of death per patient in this case. We can supply a new dataset to fit it to, using newdata
, or omitting this fits our model to the original dataset. Using a glm
, we have built a model via a link function, so we need to decide what scale to make our predictions on: the link
scale (the log-odds of death in this case) or the response
(the probability of death in this case).
LOS_model$preds <- predict(glm_binomial2, type="response") head(LOS_model,5)
We have added our predictions back into the original data.frame
using the column name preds
, and patient ID 1 has a risk of death of 0.136 or 13.6% probability of death. These predictions can be used in many different applications, but common uses of them are identifying the highest risk patients for clinical audit, comparing ratios of the observed deaths and expected deaths (sum of the probability), or risk-adjusted control chart monitoring.
LOS, the time a patient is in hospital, is another variable to examine in the dataset. This is a count, and is therefore bounded at zero (can't have counts < 0), can only take integer values (e.g. 2 or 3, but not 2.5) and is heavily right-skewed, as most patient are only admitted for one or two days, but rarely patients are in for much longer. We'll apply a similar modelling approach to the last example, but the family
argument and link function must be different. We will use the Poisson
distribution and the natural logarithm link function, the standard choices for count data^{[5]}.
glm_poisson <- glm(LOS ~ Age * Death, data=LOS_model, family="poisson" ) summary(glm_poisson)
We can't use the AUC
value here, as that is related to categories. We can only make comparisons between models or in terms of prediction error. Lets fit the same model without the interaction and see if it is different using the AIC
or the likelihood ratio test (that tell you the same thing):
glm_poisson2 <- glm(LOS ~ Age + Death, data=LOS_model, family="poisson" ) AIC(glm_poisson) AIC(glm_poisson2) # anova(glm_poisson2, glm_poisson, test="Chisq") will do the same thing without using the lmtest package lmtest::lrtest(glm_poisson, glm_poisson2)
Our AIC
is lower for the interaction model, and the likelihood ratio test suggest the reduced model is significantly worse.
One of the downfalls of a Poisson models that it's variance is fixed. This is always proportional to the conditional mean. I will spare you the maths at this point, but what it means is that the variance doesn't scale to the data. This would be fine if our data were perfectly Poisson distributed, but that is rare in practice. Data commonly display more variance than we would expect, and this is referred to as overdispersion (OD). OD causes us to underestimate the variance in our model, and our assessments of the significance of parameters are too optimistic because our error is too small. It also affects overall assessment of model fit.
We can test for this, by comparing the ratio of the sum of the squared Pearson residuals over the degrees of freedom. If this value is 1, then the variance is as expected ('equidispersion'), but if it is >1, then we have OD.
sum(residuals(glm_poisson,type="pearson") ^2)/ df.residual(glm_poisson)
Our model is overdispersed, with residual variance around 1.4 times what the Poisson distribution expects.
There are various different ways to deal with this, and this vignette demonstrates three of the following list. Our options include:
unbiased
model, so if we only care about our estimates, or predicted values (and not the error, or 'significance') we might do this.So lets try some of them out. To fit a quasi-Poisson model, we simply change the family. Note the dispersion parameter matches the calculation above:
quasi<-glm(LOS ~ Age * Death, data=LOS_model, family="quasipoisson" ) summary(quasi)
For a negative binomial model, two options are glm.nb
from the MASS
package, or glmmTMB
function from the package of the same name. Note that the function from MASS
does not take a family
argument. The dispersion factor here is quadratic to the mean, so don't compare it directly to the scale factor above.
library(MASS) nb <- glm.nb(LOS ~ Age * Death, data=LOS_model,) summary(nb)
Finally, using the random-intercept model, fitted with glmer
from the lme4
package. The extra formula argument specifies the random-intercept, where 1 means the intercept, and the |
reads as 'intercept varying by organisation.'
library(lme4) glmm <- glmer(LOS ~ scale(Age) * Death + (1|Organisation), data=LOS_model, family="poisson") summary(glmm)
confint(glmm)
These models are significantly more complicated. Confint
can be use to calculate the confidence intervals above by profiling the likelihood (provided you have the MASS
package installed). This allow you to assess the significance of the random-effects. Our CI butts up to zero here, and the AIC is slightly worse than the glm
. This suggests it is not well-modelled with the normally distributed random-intercept. The NB model appears to perform better.
Predictions from GLMMS also requires more thought. You can predict with ("conditional") or without ("marginal") the random effects. Conditional predictions will have the adjustment for their cluster (or whatever random effect you have fitted), but marginal could be considered the model average, subject to the fixed effects, removing the effects of cluster.
The LOS_model
dataset is a fabricated patient dataset with Age, Length of Stay and Death status for 300 patients across 10 hospitals ('Trusts'). It can be used to learn GLM, GLMM and other related modelling techniques, with examples above.
NELDER, J. A. & WEDDERBURN, R. W. M. 1972. Generalized Linear Models. Journal of the Royal Statistical Society. Series A (General), 135, 370-384
AKAIKE, H. 1998. Information Theory and an Extension of the Maximum Likelihood Principle. In: PARZEN, E., TANABE, K. & KITAGAWA, G. (eds.) Selected Papers of Hirotugu Akaike. New York, NY: Springer New York
SHMUELI, G. 2010. To Explain or to Predict? Statist. Sci. , 25, 289-310.
JARMAN, B., GAULT, S., ALVES, B., HIDER, A., DOLAN, S., COOK, A., HURWITZ, B. & IEZZONI, L. I. 1999. Explaining differences in English hospital death rates using routinely collected data. BMJ, 318, 1515-1520
CAMERON, A. C. & TRIVEDI, P. K. 2013. Regression Analysis of Count Data, Cambridge University Press
GOLDSTEIN, H. 2010. Multilevel Statistical Models, John Wiley & Sons Inc.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.