Logistic regression

data(diet, package = "Epi")
diet <- diet %>%
  mutate(
    chd = factor(chd, labels = c("No CHD", "CHD"))
  ) %>%
  var_labels(
    chd = "Coronary Heart Disease",
    fibre = "Fibre intake (10 g/day)"
    )

diet %>% estat(~ fibre|chd)

diet %>%
  gf_boxploth(chd ~ fibre, fill = "indianred3", alpha = 0.7) %>%
  axis_labs()
Warning: Removed 4 rows containing non-finite values (stat_boxploth).

model_binom <- glm(chd ~ fibre, data = diet, family = binomial)

model_binom %>%
  glm_coef(labels = model_labels(model_binom))

model_binom %>%
  plot_model("pred", terms = "fibre [all]", title = "")

Matched Case-Control Studies: Conditional Logistic Regression

data(bdendo, package = "Epi") 
bdendo <- bdendo %>%
  mutate(
    cancer = factor(d, labels = c('Control', 'Case')),
    gall = factor(gall, labels = c("No GBD", "GBD")),
    est = factor(est, labels = c("No oestrogen", "Oestrogen"))
  ) %>%
  var_labels(
    cancer = 'Endometrial cancer',
    gall = 'Gall bladder disease',
    est = 'Oestrogen'
  )

bdendo %>%
  mutate(
    cancer = relevel(cancer, ref = "Case"),
    est = relevel(est, ref = "Oestrogen"),
    gall = relevel(gall, ref = "GBD")
  ) %>%
  copy_labels(bdendo) %>%
  cross_tab(cancer ~ est + gall) 

library(survival)
model_clogit <- clogit(cancer == 'Case'  ~ est * gall + strata(set), data = bdendo)

model_clogit %>%
  glm_coef(labels = c("Oestrogen/No oestrogen", "GBD/No GBD",
                      "Oestrogen:GBD Interaction"))

require(ggeffects)
Loading required package: ggeffects
bdendo_pred <- ggemmeans(model_clogit, terms = c('gall', 'est'))

bdendo_pred %>%
  gf_pointrange(predicted + conf.low + conf.high ~ x|group, col = ~ x) %>%
  gf_labs(y = "P(cancer)", x = "", col = get_label(bdendo$gall))

Ordinal Logistic Regression

library(ordinal)

Attaching package: 'ordinal'
The following object is masked from 'package:dplyr':

    slice
data(housing, package = "MASS")
housing <- housing %>%
  var_labels(
    Sat = "Satisfaction",
    Infl = "Perceived influence",
    Type = "Type of rental",
    Cont = "Contact"
    )

model_clm <- clm(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)

model_clm %>%
  glm_coef(labels = model_labels(model_clm, intercept = FALSE))

model_clm %>%
  plot_model(type = "pred", terms = c("Infl", "Cont"), 
           dot.size = 1, title = "") %>%
  gf_theme(axis.text.x = element_text(angle = 45, hjust = 1))

model_clm %>%
  plot_model(type = "pred", terms = c("Infl", "Type"), 
           dot.size = 1, title = "") %>%
  gf_theme(axis.text.x = element_text(angle = 45, hjust = 1))

emmip(model_clm, Infl ~ Type |Cont) %>%
  gf_labs(x = "Type of rental", col = "Perceived influence")

Poisson Regression

data(quine, package = "MASS")
levels(quine$Eth) <- c("Aboriginal", "White")
levels(quine$Sex) <- c("Female", "Male")

quine <- quine %>%
  var_labels(
    Days = "Number of absent days",
    Eth = "Ethnicity",
    Age = "Age group"
    )

EDA

quine %>%
  group_by(Eth, Sex, Age) %>%
  summarise(
    n = n(),
    Mean = mean(Days, na.rm = TRUE),
    SD = sd(Days, na.rm = TRUE),
    Median = median(Days, na.rm = TRUE),
    CV = rel_dis(Days)
  )
`summarise()` regrouping output by 'Eth', 'Sex' (override with `.groups` argument)
# A tibble: 16 x 8
# Groups:   Eth, Sex [4]
   Eth        Sex    Age       n  Mean    SD Median    CV
   <fct>      <fct>  <fct> <int> <dbl> <dbl>  <dbl> <dbl>
 1 Aboriginal Female F0        5 17.6  17.4      11 0.987
 2 Aboriginal Female F1       15 18.9  16.3      13 0.865
 3 Aboriginal Female F2        9 32.6  27.3      20 0.839
 4 Aboriginal Female F3        9 14.6  14.9      10 1.02 
 5 Aboriginal Male   F0        8 11.5   7.23     12 0.629
 6 Aboriginal Male   F1        5  9.6   4.51      7 0.469
 7 Aboriginal Male   F2       11 30.9  17.8      32 0.576
 8 Aboriginal Male   F3        7 27.1  10.4      28 0.382
 9 White      Female F0        5 19.8   9.68     20 0.489
10 White      Female F1       17  7.76  6.48      6 0.834
11 White      Female F2       10  5.7   4.97      4 0.872
12 White      Female F3       10 13.5  11.5      12 0.851
13 White      Male   F0        9 13.6  20.9       7 1.54 
14 White      Male   F1        9  5.56  5.39      5 0.970
15 White      Male   F2       10 15.2  12.9      12 0.848
16 White      Male   F3        7 27.3  22.9      27 0.840

Model

model_pois <- glm(Days ~ Eth + Sex + Age, family = poisson, data = quine)

model_pois %>%
  glm_coef(labels = model_labels(model_pois), se_rob = TRUE)

model_pois %>% glance()

Negative-binomial

deviance(model_pois) / df.residual(model_pois)
[1] 12.44646

Thus, we have over-dispersion. One option is to use a negative binomial distribution.

library(MASS)

Attaching package: 'MASS'
The following objects are masked _by_ '.GlobalEnv':

    birthwt, housing, quine
The following object is masked from 'package:dplyr':

    select
model_negbin <- glm.nb(Days ~ Eth + Sex + Age, data = quine)

model_negbin %>%
  glm_coef(labels = model_labels(model_negbin), se_rob = TRUE)

model_negbin %>% glance()

age group is a factor with more than two levels and is significant:

model_negbin %>% Anova()

model_negbin %>%
  plot_model(type = "pred", terms = c("Age", "Eth"), 
           dot.size = 1.5, title = "")

emmip(model_negbin, Eth ~ Age|Sex) %>%
  gf_labs(y = "Number of absent days", x = "Age group", col = "Ethnicity")

## Adjusting CIs and p-values for multiple comparisons

multiple(model_negbin, ~ Age|Eth)$df

multiple(model_negbin, ~ Age|Eth)$fig_ci %>%
  gf_labs(x = "IRR")

multiple(model_negbin, ~ Age|Eth)$fig_ci %>%
  gf_labs(x = "IRR")

Survival Analysis

effect of thiotepa versus placebo on the development of bladder cancer.

library(survival)
data(bladder)
bladder <- bladder %>%
  mutate(times = stop,
         rx = factor(rx, labels=c("Placebo", "Thiotepa"))
         ) %>%
  var_labels(times = "Survival time",
             rx = "Treatment")

model_surv <- survreg(Surv(times, event) ~ rx, data = bladder)

model_surv %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo", "Scale"), se_rob = TRUE)

model_exp <- survreg(Surv(times, event) ~ rx, data = bladder, dist = "exponential")

model_exp %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo"), se_rob = TRUE)

model_exp %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo"))

model_exp %>%
  plot_model(type = "pred", terms = ~ rx, dot.size = 1.5, title = "") %>%
  gf_labs(y = "Survival time", x = "Treatment", title = "")

model_cox <-  coxph(Surv(times, event) ~ rx, data = bladder)

model_cox %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo"))

model_cox %>%
  plot_model(type = "pred", terms = ~ rx, dot.size = 1.5, 
           title = "") %>%
  gf_labs(x = "Treatment", title = "")

Mixed Linear Regression Models

library(nlme)

Attaching package: 'nlme'
The following objects are masked from 'package:ordinal':

    ranef, VarCorr
The following object is masked from 'package:dplyr':

    collapse
data(Orthodont)

Orthodont <- Orthodont %>%
  var_labels(
    distance = "Pituitary distance (mm)",
    age = "Age (years)"
    )

model_lme <- lme(distance ~ Sex * I(age - mean(age, na.rm = TRUE)), random = ~ 1|Subject, 
                 method = "ML", data = Orthodont)

model_lme %>%
  glm_coef(labels = c(
    "Constant", 
    "Sex: female-male", 
    "Age (years)",
    "Sex:Age interaction"
    ))

model_lme %>%
  plot_model("pred", terms = age ~ Sex, 
           show.data = TRUE, jitter = 0.1, dot.size = 1.5) %>%
  gf_labs(y = get_label(Orthodont$distance), x = "Age (years)", title = "")