In this tutorial, we analyze the sore throat dataset.

1 The sore throat dataset

Dataset sore contains results of a study on throat age 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 Goal

  • Explore the dataset
  • Fit a logit model using “type” and interpret the parameter estimate
  • Fit a logit model considering all predictors
  • Fit an optimal logit model, interpret the parameter estimates, and conduct inference about the effects.

Load the libraries:

library(tidyverse)
## ── Attaching packages ──────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ─────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
#library(MASS)
#library(boot)
#library(car)

sore <- read_csv("https://raw.githubusercontent.com/GTPB/PSLS20/master/data/sore.csv")
## Parsed with column specification:
## cols(
##   duration = col_double(),
##   type = col_character(),
##   sore = col_double()
## )
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
sore$type <- as.factor(sore$type)
sore$sore <- as.factor(sore$sore)

3 Data Exploration

#####################
#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))

4 Analysis

4.1 Fit for “type” only

Note: We aill assess 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).

#fit <- glm(...)

#inference/interpretation

Interpretation:

# odds on sore throat when wearing mask
14/4 # in our dataset, out of 18 patients that wore a mask, 14 had a sour throat
exp(fit$coef[1]) # when exponentiating the intercept, we obtain the odds of a sore throat for the reference type (mask) 

# odds on sore throat when wearing a tube
8/9 # in our dataset, out of 17 patients that wore a mask, 8 had a sour throat
exp(fit$coef[1]+fit$coef[2]) # when exponentiating the intercept and add the effect of changing the ttype to tube, we obtain the odds of a sore throat for the tube type 

# odds ratio for type 
exp(fit$coef[2]) 

Here, we constructed a model with type as the sole explanatory variable to explain soarness of the throat. Can you see a problem with that?

4.2 Fit for type, duration and type-duration interaction

Do we need the interaction term?

4.3 Fit an “optimal” model

4.4 Interpret the model parameters of the final model

5 Conclusion

Formulate a conclusion that answers the research question (see “goal”).

s

LS0tCnRpdGxlOiAiVHV0b3JpYWwgMTAuMjogTG9naXN0aWMgcmVncmVzc2lvbiBvbiB0aGUgc29yZSB0aHJvYXQgZGF0YXNldCIgICAKb3V0cHV0OgogICAgaHRtbF9kb2N1bWVudDoKICAgICAgY29kZV9kb3dubG9hZDogdHJ1ZSAgICAKICAgICAgdGhlbWU6IGNvc21vCiAgICAgIHRvYzogdHJ1ZQogICAgICB0b2NfZmxvYXQ6IHRydWUKICAgICAgaGlnaGxpZ2h0OiB0YW5nbwogICAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKLS0tCgpJbiB0aGlzIHR1dG9yaWFsLCB3ZSBhbmFseXplIHRoZSBzb3JlIHRocm9hdCBkYXRhc2V0LgoKIyBUaGUgc29yZSB0aHJvYXQgZGF0YXNldAoKRGF0YXNldCBzb3JlIGNvbnRhaW5zIHJlc3VsdHMgb2YgYSBzdHVkeSBvbiB0aHJvYXQgYWdlIHVwb24gc3VyZ2VyeSB3aXRoIGdlbmVyYWwgYW5hZXN0aGVzaWEuICAKVGhlIHZhcmlhYmxlIHNvcmUgZW5jb2RlcyB3aGV0aGVyIGEgcGF0aWVudCBleHBlcmllbmNlZCBhIHNvcmUgdGhyb2F0IG9uCndha2luZyAoMDogbm8sIDE6IHllcykgIGFzIGEgZnVuY3Rpb24gb2YgdGhlIGR1cmF0aW9uIG9mIHRoZSBzdXJnZXJ5CihkdXJhdGlvbikgaW4gbWludXRlcyBhbmQgdGhlIHR5cGUgb2YgZGV2aWNlICh2YXJpYWJsZSB0eXBlKSB1c2VkIHRvIHNlY3VyZQp0aGUgYWlyd2F5OiBsYXJ5bmdlYWwgbWFzayAoIm1hc2siKSBvciB0cmFjaGVhbCB0dWJlICgidHViZSIpIChEYXRhIHRha2VuCmZyb20gQWdyZXN0aSAyMDAyKS4KCiMgR29hbAoKLSBFeHBsb3JlIHRoZSBkYXRhc2V0Ci0gRml0IGEgbG9naXQgbW9kZWwgdXNpbmcgInR5cGUiIGFuZCBpbnRlcnByZXQgdGhlIHBhcmFtZXRlciBlc3RpbWF0ZQotIEZpdCBhIGxvZ2l0IG1vZGVsIGNvbnNpZGVyaW5nIGFsbCBwcmVkaWN0b3JzCi0gRml0IGFuIG9wdGltYWwgbG9naXQgbW9kZWwsIGludGVycHJldCB0aGUgcGFyYW1ldGVyIGVzdGltYXRlcywgCmFuZCBjb25kdWN0IGluZmVyZW5jZSBhYm91dCB0aGUgZWZmZWN0cy4KCkxvYWQgdGhlIGxpYnJhcmllczoKCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKI2xpYnJhcnkoTUFTUykKI2xpYnJhcnkoYm9vdCkKI2xpYnJhcnkoY2FyKQoKc29yZSA8LSByZWFkX2NzdigiaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL0dUUEIvUFNMUzIwL21hc3Rlci9kYXRhL3NvcmUuY3N2IikKaGVhZChzb3JlKQpgYGAKCmBgYHtyfQpzb3JlJHR5cGUgPC0gYXMuZmFjdG9yKHNvcmUkdHlwZSkKc29yZSRzb3JlIDwtIGFzLmZhY3Rvcihzb3JlJHNvcmUpCmBgYAoKCiMgRGF0YSBFeHBsb3JhdGlvbgoKYGBge3J9CiMjIyMjIyMjIyMjIyMjIyMjIyMjIwojZGF0YSBleHBsb3JhdGlvbgojIyMjIyMjIyMjIyMjIyMjIyMjIyMKdGFibGUoc29yZVssYygidHlwZSIsInNvcmUiKV0pCmBgYAoKYGBge3J9CnBsb3QoYXMuZG91YmxlKHNvcmUpfmR1cmF0aW9uLHBjaD1hcy5kb3VibGUodHlwZSksY29sPWFzLmRvdWJsZSh0eXBlKSwgZGF0YT1zb3JlKQpsZWdlbmQoImJvdHRvbXJpZ2h0IixwY2g9MToyLGNvbD0xOjIsbGVnZW5kPWxldmVscyhzb3JlJHR5cGUpKQpgYGAKCiMgQW5hbHlzaXMKCiMjIEZpdCBmb3IgInR5cGUiIG9ubHkKCk5vdGU6IFdlIGFpbGwgYXNzZXNzIHRoZSBzaWduaWZpY2FuY2Ugb2YgdGhlIGRpZmZlcmVudCBjb3ZhcmlhdGVzIGluIHRoZQptb2RlbCBieSB1c2luZyBhIGxpa2VsaWhvb2QgcmF0aW8gdGVzdHMgKEFub3ZhIGZ1bmN0aW9uIGluIFIpIGJlY2F1c2UgdGhlCnAtdmFsdWVzIGFyZSBtb3JlIGFjY3VyYXRlIHRoYW4gdGhlIHAtdmFsdWVzIG9mIHRoZSBXYWxkLXRlc3RzIChzdW1tYXJ5CmZ1bmN0aW9uIGluIFIpLiAKCmBgYHtyfQojZml0IDwtIGdsbSguLi4pCgojaW5mZXJlbmNlL2ludGVycHJldGF0aW9uCgpgYGAKCkludGVycHJldGF0aW9uOgoKYGBgCiMgb2RkcyBvbiBzb3JlIHRocm9hdCB3aGVuIHdlYXJpbmcgbWFzawoxNC80ICMgaW4gb3VyIGRhdGFzZXQsIG91dCBvZiAxOCBwYXRpZW50cyB0aGF0IHdvcmUgYSBtYXNrLCAxNCBoYWQgYSBzb3VyIHRocm9hdApleHAoZml0JGNvZWZbMV0pICMgd2hlbiBleHBvbmVudGlhdGluZyB0aGUgaW50ZXJjZXB0LCB3ZSBvYnRhaW4gdGhlIG9kZHMgb2YgYSBzb3JlIHRocm9hdCBmb3IgdGhlIHJlZmVyZW5jZSB0eXBlIChtYXNrKSAKCiMgb2RkcyBvbiBzb3JlIHRocm9hdCB3aGVuIHdlYXJpbmcgYSB0dWJlCjgvOSAjIGluIG91ciBkYXRhc2V0LCBvdXQgb2YgMTcgcGF0aWVudHMgdGhhdCB3b3JlIGEgbWFzaywgOCBoYWQgYSBzb3VyIHRocm9hdApleHAoZml0JGNvZWZbMV0rZml0JGNvZWZbMl0pICMgd2hlbiBleHBvbmVudGlhdGluZyB0aGUgaW50ZXJjZXB0IGFuZCBhZGQgdGhlIGVmZmVjdCBvZiBjaGFuZ2luZyB0aGUgdHR5cGUgdG8gdHViZSwgd2Ugb2J0YWluIHRoZSBvZGRzIG9mIGEgc29yZSB0aHJvYXQgZm9yIHRoZSB0dWJlIHR5cGUgCgojIG9kZHMgcmF0aW8gZm9yIHR5cGUgCmV4cChmaXQkY29lZlsyXSkgCmBgYAoKSGVyZSwgd2UgY29uc3RydWN0ZWQgYSBtb2RlbCB3aXRoIGB0eXBlYCBhcyB0aGUgc29sZSBleHBsYW5hdG9yeQp2YXJpYWJsZSB0byBleHBsYWluIHNvYXJuZXNzIG9mIHRoZSB0aHJvYXQuIENhbiB5b3Ugc2VlIGEgcHJvYmxlbSB3aXRoIHRoYXQ/CgojIyBGaXQgZm9yIHR5cGUsIGR1cmF0aW9uIGFuZCB0eXBlLWR1cmF0aW9uIGludGVyYWN0aW9uCgpgYGB7cn0KCmBgYAoKRG8gd2UgbmVlZCB0aGUgaW50ZXJhY3Rpb24gdGVybT8KCiMjIEZpdCBhbiAib3B0aW1hbCIgbW9kZWwKCmBgYHtyfQoKYGBgCgojIyBJbnRlcnByZXQgdGhlIG1vZGVsIHBhcmFtZXRlcnMgb2YgdGhlIGZpbmFsIG1vZGVsCgoKIyBDb25jbHVzaW9uCgpGb3JtdWxhdGUgYSBjb25jbHVzaW9uIHRoYXQgYW5zd2VycyB0aGUgcmVzZWFyY2ggcXVlc3Rpb24gKHNlZSAiZ29hbCIpLgoKCgoKCnM=