Tutorial 10.2: Logistic regression on the sore throat dataset
1 The sore throat dataset
The sore throat dataset contains results of a study on throat soarness upon surgery with general anaesthesia. The variable sore encodes whether a patient experienced a sore throat on waking (0: no, 1: yes) as a function of the duration of the surgery (duration) in minutes and the type of device (variable type) used to secure the airway: laryngeal mask (“mask”) or tracheal tube (“tube”) (Data taken from Agresti 2002).
2 Goals
- Explore the dataset
- Fit a logit model using “type” and interpret the parameter estimate
- Fit a logit model considering all predictors, interpret parameter estimates, and conduct inference about the effects.
This file contains many details for explaining the models to students. We do not expect that reports are this extensive. I.e. a report on the exercise typically would contain as much information as in the soreBase rmarkdown file.
library(tidyverse)
library(car)3 Import data
sore <- read_csv("https://raw.githubusercontent.com/GTPB/PSLS20/master/data/sore.csv")
head(sore)## # A tibble: 6 x 3
## duration type sore
## <dbl> <chr> <dbl>
## 1 45 mask 0
## 2 15 mask 0
## 3 40 mask 1
## 4 83 tube 1
## 5 90 tube 1
## 6 25 tube 1## set the type and sore variables as factors
sore <- sore %>%
mutate(type = as.factor(type), sore = as.factor(sore))4 Data Exploration
table(sore[,c("type","sore")])## sore
## type 0 1
## mask 4 14
## tube 9 8plot(as.double(sore)~duration,pch=as.double(type),col=as.double(type), data=sore)
legend("bottomright",pch=1:2,col=1:2,legend=levels(sore$type))
5 Analysis: logistic regression
5.1 Type variable as only predictor
fitType=glm(sore~type,data=sore,family=binomial)
#inference/interpretation
summary(fitType)##
## Call:
## glm(formula = sore ~ type, family = binomial, data = sore)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.734 -1.128 0.709 0.709 1.228
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2528 0.5669 2.210 0.0271 *
## typetube -1.3705 0.7467 -1.836 0.0664 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 46.180 on 34 degrees of freedom
## Residual deviance: 42.578 on 33 degrees of freedom
## AIC: 46.578
##
## Number of Fisher Scoring iterations: 4confint(fitType)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 0.2275566 2.51346332
## typetube -2.9292777 0.04379487anova(fitType,test="LRT")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: sore
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 34 46.180
## type 1 3.6022 33 42.578 0.0577 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretation of the model coefficients:
14/4 # odds on sore throat when wearing mask## [1] 3.5exp(fitType$coef[1]) # odds on sore throat when wearing mask## (Intercept)
## 3.58/9 # odds on sore throat when having tube## [1] 0.8888889exp(fitType$coef[1]+fitType$coef[2]) # odds on sore throat when having tube## (Intercept)
## 0.8888889(8/9)/(14/4) # odds ration type## [1] 0.2539683exp(fitType$coef[2]) # odds ration type## typetube
## 0.2539683Note that the model above did not account for differences in the duration of the treatment.
5.2 Model Selection
We first will fit a model with a main effects for type, duration and a type x duration interaction.
fitTypeDurationInt=glm(sore~type*duration,data=sore,family=binomial)
summary(fitTypeDurationInt,cor=TRUE)##
## Call:
## glm(formula = sore ~ type * duration, family = binomial, data = sore)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9707 -0.3779 0.3448 0.7292 1.9961
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.04979 1.46940 0.034 0.9730
## typetube -4.47224 2.46707 -1.813 0.0699 .
## duration 0.02848 0.03429 0.831 0.4062
## typetube:duration 0.07460 0.05777 1.291 0.1966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 46.180 on 34 degrees of freedom
## Residual deviance: 28.321 on 31 degrees of freedom
## AIC: 36.321
##
## Number of Fisher Scoring iterations: 6
##
## Correlation of Coefficients:
## (Intercept) typetube duration
## typetube -0.60
## duration -0.92 0.55
## typetube:duration 0.55 -0.91 -0.59Anova(fitTypeDurationInt,test="LR",type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: sore
## LR Chisq Df Pr(>Chisq)
## type 3.9458 1 0.04699 *
## duration 0.8069 1 0.36904
## type:duration 1.8169 1 0.17768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We can also fit the model without the interaction effect:
fitTypeDurationNoInt=glm(sore~type+duration,data=sore,family=binomial)
summary(fitTypeDurationNoInt,cor=TRUE)##
## Call:
## glm(formula = sore ~ type + duration, family = binomial, data = sore)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3802 -0.5358 0.3047 0.7308 1.7821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.41734 1.09457 -1.295 0.19536
## typetube -1.65895 0.92285 -1.798 0.07224 .
## duration 0.06868 0.02641 2.600 0.00931 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 46.180 on 34 degrees of freedom
## Residual deviance: 30.138 on 32 degrees of freedom
## AIC: 36.138
##
## Number of Fisher Scoring iterations: 5
##
## Correlation of Coefficients:
## (Intercept) typetube
## typetube -0.25
## duration -0.83 -0.16Anova(fitTypeDurationNoInt,test="LR",type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: sore
## LR Chisq Df Pr(>Chisq)
## type 3.5134 1 0.0608744 .
## duration 12.4396 1 0.0004203 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(fitTypeDurationNoInt)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -3.80786158 0.61635531
## typetube -3.64627873 0.07434058
## duration 0.02547651 0.13216711We assessed the significance of the different covariates in the model by using a likelihood ratio tests (Anova function in R) because the p-values are more accurate than the p-values of the Wald-tests (summary function in R). The interaction between type and duration is not significant (p=0.177). Hence the interaction term can be removed from the model. In the model with main effects for type and duration, the duration effect is extremely significant (p$<$0.001). The effect of type is not significant (p=0.061). However, since the effect of type is close to being statistically significant we will leave it in the model so as to estimate a “pure” effect for duration. If we remove type from the model and if the effect truly exists but there was not enough power for detecting it, then removing the type from the model might result in a biased parameter estimate for duration as it might captures a part of the effect of type.
6 Interpretation of the different models
6.1 Model Type
# Model 1 only type
plot(sore~duration,pch=as.double(type),col=type,data=sore,main="Type",xlim=c(0,120))## Warning in spineplot.default(x, y, ...): x axis is on a cumulative probability
## scale, 'xlim' must be in [0,1]legend("bottomright",pch=c(1:2,NA,NA),lty=c(0,0,1,1),col=c(1,2,1,2),legend=rep(levels(sore$type),2))
#construct new data to make a nice plot (predict data using the model)
x=seq(0,120,1)
#use type response in order to predict probabilities
newdataMask=data.frame(duration=x,type=factor("mask",levels=levels(sore$type)))
yPredMask=predict(fitType,newdataMask,type="response")
lines(x,yPredMask)
#use type response in order to predict probabilities
newdataTube=data.frame(duration=x,type=factor("tube",levels=levels(sore$type)))
yPredTube=predict(fitType,newdataTube,type="response")
lines(x,yPredTube,col=2)
The model for the log odds on a sore throat after surgery in function of anaesthesia type is given:
\[ \log \text{Odds}_i=\beta_0+\beta_\text{tube}x_\text{tube,i} \]
with \(x_\text{tube,i}=0\) when the patient was anaesthesized with a mask and \(x_\text{tube,i}=1\) when a tube was used. Based on this model we can estimate the log odds on a sore throat after surgery when wearing a mask with
\[ \log \widehat{\text{Odds}}^\text{mask}=\hat \beta_0. \]
Hence, the intercept has the interpretation of the log odds on a sore throat after a surgery with mask anaesthesia. The logg odds on a sore throat after a surgery with a tube is estimated as
\[ \log \widehat{\text{Odds}}^\text{tube}=\hat \beta_0+\hat\beta_{tube}. \]
The log odds ratio on a sore throat between tube and mask can then be estimated as:
\[ \log \widehat{\text{Odds}}^\text{tube}-\log \widehat{\text{Odds}}^\text{mask}=\hat \beta_0+\hat\beta_{tube}-\hat\beta_{0}=\hat\beta_\text{tube}. \]
Hence, the slope for the dummy variable tube has the interpretation of the log odds ratio on a sore throat between tube and mask.
Interpretation based on the parameter estimates: The probability on a sore throat when wearing a mask is on average 3.5 times higher than the probability to wake up from surgery without a sore throat when wearing a mask (95% CI [1.3,12.3]). The odds on a sore throat after surgery with anaesthesia using a tube is on average 4 times smaller than the odds with a mask, but this odds ratio is not statistically significant ( 95% CI [0.05,1.04]).
6.2 Model Duration
fitDuration=glm(sore~duration,data=sore,family=binomial)
summary(fitDuration,cor=TRUE)##
## Call:
## glm(formula = sore ~ duration, family = binomial, data = sore)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0964 -0.7392 0.3020 0.8711 1.3753
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.21358 0.99874 -2.216 0.02667 *
## duration 0.07038 0.02667 2.639 0.00831 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 46.180 on 34 degrees of freedom
## Residual deviance: 33.651 on 33 degrees of freedom
## AIC: 37.651
##
## Number of Fisher Scoring iterations: 5
##
## Correlation of Coefficients:
## (Intercept)
## duration -0.91Anova(fitDuration,test="LR",type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: sore
## LR Chisq Df Pr(>Chisq)
## duration 12.528 1 0.0004008 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(fitDuration)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -4.45076704 -0.4359574
## duration 0.02639986 0.1334656plot(sore~duration,pch=as.double(type),col=type,data=sore,main="Type",xlim=c(0,120))## Warning in spineplot.default(x, y, ...): x axis is on a cumulative probability
## scale, 'xlim' must be in [0,1]legend("bottomright",pch=c(1:2,NA,NA),lty=c(0,0,1,1),col=c(1,2,1,2),legend=rep(levels(sore$type),2))
#construct new data to make a nice plot (predict data using the model)
x=seq(0,120,1)
#use type response in order to predict probabilities
newdataMask=data.frame(duration=x,type=factor("mask",levels=levels(sore$type)))
yPredMask=predict(fitDuration,newdataMask,type="response")
lines(x,yPredMask)
#use type response in order to predict probabilities
newdataTube=data.frame(duration=x,type=factor("tube",levels=levels(sore$type)))
yPredTube=predict(fitDuration,newdataTube,type="response")
lines(x,yPredTube,col=2)
The model for the log odds on a sore throat after surgery in function of the surgery time can be denoted as: \[ \log \text{Odds}_i=\beta_0+\beta_\text{time}t_\text{i} \]
Then the intercept has the interpretation of the odds on a sore throat when the surgery took \(t=0\) minutes. Note, that this parameter does not have a useful interpretation in practice. Moreover, it also involves an extrapolation. The intercept in a model were the time is centered, however, has a different interpretation: it is the odds on a sore throat when the surgery took 48.0 minutes, i.e. the average duration of surgery in the dataset.
The odds on a surgery that took 30 minutes can be estimated as:
\[ \log \text{Odds}^{30min}=\beta_0+\beta_\text{time}30 \]
and the odds on a surgery that took 31 minutes can be estimated as:
\[ \log \widehat{\text{Odds}}^{31min}=\beta_0+\hat \beta_\text{time}31, \]
hence the slope \(\beta_\text{time}\) has the interpretation of the log odds ratio between surgeries with a one minute difference in duration.
6.3 Model Type + Duration
plot(sore~duration,pch=as.double(type),col=type,data=sore,main="Type",xlim=c(0,120))## Warning in spineplot.default(x, y, ...): x axis is on a cumulative probability
## scale, 'xlim' must be in [0,1]legend("bottomright",pch=c(1:2,NA,NA),lty=c(0,0,1,1),col=c(1,2,1,2),legend=rep(levels(sore$type),2))
#construct new data to make a nice plot (predict data using the model)
x=seq(0,120,1)
#use type response in order to predict probabilities
newdataMask=data.frame(duration=x,type=factor("mask",levels=levels(sore$type)))
yPredMask=predict(fitTypeDurationNoInt,newdataMask,type="response")
lines(x,yPredMask)
#use type response in order to predict probabilities
newdataTube=data.frame(duration=x,type=factor("tube",levels=levels(sore$type)))
yPredTube=predict(fitTypeDurationNoInt,newdataTube,type="response")
lines(x,yPredTube,col=2)
summary(fitTypeDurationNoInt)##
## Call:
## glm(formula = sore ~ type + duration, family = binomial, data = sore)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3802 -0.5358 0.3047 0.7308 1.7821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.41734 1.09457 -1.295 0.19536
## typetube -1.65895 0.92285 -1.798 0.07224 .
## duration 0.06868 0.02641 2.600 0.00931 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 46.180 on 34 degrees of freedom
## Residual deviance: 30.138 on 32 degrees of freedom
## AIC: 36.138
##
## Number of Fisher Scoring iterations: 5The model for the log odds on a sore throat after surgery with a main effect for surgery time and type can be written as:
\[ \log \text{Odds}_i=\beta_0+\beta_\text{tube}x_\text{tube}+\beta_\text{time}t_i \]
The intercept has the interpretation of the log odds on a sore throat for patients that were anaesthetised with a mask (\(x_\text{tube}=0\)) and that had surgery that had a duration of 0 minutes (\(t=0\)). Again this log odds does not have a useful interpretation.
The slope for tube again has the interpretation of the log odds ratio on a sore throat between groups of patients that were anaesthetised with tubes vs masks and got a surgery with the same duration. This can be rephrased as the log odds ratio on a sore throat between groups of patients that were anaesthetised with tubes vs groups of patients masks after correcting for the duration of the surgery. (This will result in a shift of the curve between patients with mask and tubes.)
The slope for time has the interpretation of the log odds ratio on a sore throat between groups of patients that had surgery with a one minute duration difference and were anaesthetised with the same type, i.e. the log odds ratio on a sore throat between groups of patients that had surgery with a one minute duration difference after correction for the anaesthesia type.
6.4 Model Type x Duration
plot(sore~duration,pch=as.double(type),col=type,data=sore,main="Type",xlim=c(0,120))## Warning in spineplot.default(x, y, ...): x axis is on a cumulative probability
## scale, 'xlim' must be in [0,1]legend("bottomright",pch=c(1:2,NA,NA),lty=c(0,0,1,1),col=c(1,2,1,2),legend=rep(levels(sore$type),2))
#construct new data to make a nice plot (predict data using the model)
x=seq(0,120,1)
#use type response in order to predict probabilities
newdataMask=data.frame(duration=x,type=factor("mask",levels=levels(sore$type)))
yPredMask=predict(fitTypeDurationInt,newdataMask,type="response")
lines(x,yPredMask)
#use type response in order to predict probabilities
newdataTube=data.frame(duration=x,type=factor("tube",levels=levels(sore$type)))
yPredTube=predict(fitTypeDurationInt,newdataTube,type="response")
lines(x,yPredTube,col=2)
summary(fitTypeDurationInt)##
## Call:
## glm(formula = sore ~ type * duration, family = binomial, data = sore)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9707 -0.3779 0.3448 0.7292 1.9961
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.04979 1.46940 0.034 0.9730
## typetube -4.47224 2.46707 -1.813 0.0699 .
## duration 0.02848 0.03429 0.831 0.4062
## typetube:duration 0.07460 0.05777 1.291 0.1966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 46.180 on 34 degrees of freedom
## Residual deviance: 28.321 on 31 degrees of freedom
## AIC: 36.321
##
## Number of Fisher Scoring iterations: 6The model for the log odds on a sore throat after surgery with a main effect for surgery time and type and a time-type interaction is given by:
\[ \log \text{Odds}_i=\beta_0+\beta_\text{tube}x_\text{tube}+\beta_\text{time}t_i+\beta_\text{time-tube}x_\text{tube,i}t_i \]
Hence, for patients with a mask the odds can be predicted in function of time by
\[ \log \text{Odds}^\text{mask}=\beta_0+\beta_\text{time}t \]
and the odds on a sore throat for patients that were anaesthetised with a tube in function of time by
\[ \log \text{Odds}^\text{tube}=\beta_0+\beta_\text{tube}+(\beta_\text{time}+\beta_\text{time-tube})t \]
The intercept \(\beta_0\) is again the log odds on a sore throat for patients that were anaesthetised with a mask (\(x_\text{tube}=0\)) and that had surgery that had a duration of 0 minutes (\(t=0\)).
The slope \(\beta_\text{tube}\) has the interpretation of the log odds ratio on a sore throat for patients that were anaesthetised with a tube (\(x_\text{tube}=0\)) vs patients that were anaesthetised with a mask and that had surgery that had a duration of 0 minutes (\(t=0\)). Again these two parameters do not have a useful interpretation.
The parameter \(\beta_\text{time}\) has the interpretation of the log odds ratio on a sore throat between groups of patients that had surgery with a one minute duration difference and were anaesthetised with mask.
The parameter \(\beta_\text{time-tube}\) allows for a different log odds ratio on a sore throat for patients that had surgery with a one minute duration difference in groups of patients that were anaesthetised with a tube than that within patient groups with a mask, i.e. it allows for a different odds profile in function of time for patients with a tube than for patients with a mask.
7 Conclusion
The effect of the duration of surgery on the odds on a sore throat is extremely significant (p<0.001). If the time of surgery increases, patients have a higher probability to wake up with a sore throat. The odds ratio for a one minute increase in surgery duration is 1.07 (95% CI [1.03,1.14]). There is no significant interaction between the type of anaesthesia and duration of the surgery (p=0.117). While there is some evidence in the data suggesting that the use of a tube is beneficial, i.e. the odds on a sore throat is more than a factor of 5 lower when using an anaesthesia with a tube than with a mask (OR=0.19, 95% CI [0.02,1.08]), the latter effect is not statistically significant (p=0.061).