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.

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  8

plot(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: 4

confint(fitType)

## Waiting for profiling to be done...

##                  2.5 %     97.5 %
## (Intercept)  0.2275566 2.51346332
## typetube    -2.9292777 0.04379487

anova(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 ' ' 1

Interpretation of the model coefficients:

14/4 # odds on sore throat when wearing mask

## [1] 3.5

exp(fitType$coef[1]) # odds on sore throat when wearing mask

## (Intercept) 
##         3.5

8/9 # odds on sore throat when having tube

## [1] 0.8888889

exp(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.2539683

exp(fitType$coef[2]) # odds ration type

##  typetube 
## 0.2539683

Note 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 (full model, all variable and interaction terms are incorporated).

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.59

Anova(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 ' ' 1

fitTypeDurationInt$coefficients

##       (Intercept)          typetube          duration typetube:duration 
##        0.04978674       -4.47224144        0.02847802        0.07460127

exp(fitTypeDurationInt$coef[1])

## (Intercept) 
##    1.051047

exp(fitTypeDurationInt$coef[2])

##   typetube 
## 0.01142169

exp(fitTypeDurationInt$coef[3])

## duration 
## 1.028887

#############
#odds on sore throat when wearing mask
exp(fitTypeDurationInt$coef[1])

## (Intercept) 
##    1.051047

#odds on sore throat when having tube
exp(fitTypeDurationInt$coef[1]+fitTypeDurationInt$coef[2])

## (Intercept) 
##  0.01200473

#odds ration type
exp(fitTypeDurationInt$coef[2]) 

##   typetube 
## 0.01142169

We 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.16

Anova(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 ' ' 1

confint(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.13216711

We 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.178). 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 main effect of type is not significant (p=0.061). However, since the effect of type is close to be 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 capture a part of the effect of type.

6 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).

---

Home