Lecture 10: Correlated data

library(tidyverse)
library(tidymodels)
library(knitr)
library(patchwork)

Scenario 1

Scenario 1a

Below is the code for the binomial logistic regression model and confidence intervals calculated using the profile likelihood approach (confidence intervals calculated based on likelihood function with no assumptions on the distribution of the parameter of interest, \(\beta_0\)).

model_1a <- glm(phat_1a ~ 1, family = binomial, weight = rep(10,24), 
                data = scenario_1)

tidy(model_1a) |>
  kable(digits = 3)

term	estimate	std.error	statistic	p.value
(Intercept)	0.067	0.129	0.516	0.606

95% CI of the odds of a pup with deformity

exp(confint(model_1a))

    2.5 %    97.5 % 
0.8298899 1.3778985 

The 95% CI of the probability of a pup with deformity

exp(confint(model_1a)) / (1 + exp(confint(model_1a)))

    2.5 %    97.5 % 
0.4535190 0.5794606 

Next we fit a quasibinomial model

quasi_model_1a <- glm(phat_1a ~ 1, family = quasibinomial,
                      weight = rep(10,24), data = scenario_1)

tidy(quasi_model_1a)

# A tibble: 1 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   0.0667     0.122     0.546   0.590

95% CI of the odds of a pup with deformity

exp(confint(quasi_model_1a))

    2.5 %    97.5 % 
0.8414573 1.3588535 

The 95% CI of the probability of a pup with deformity

exp(confint(quasi_model_1a)) / (1 + exp(confint(quasi_model_1a)))

    2.5 %    97.5 % 
0.4569518 0.5760652 

Scenario 1b

Fit models for Scenario 1b

Fit the binomial and quasibinomial models for the data generated under Scenario 1b. Calculate the 95% CI for \(p\) for each model.

# add code here

# add code here

Questions

Question 1

How do the quasibinomial analysis for Scenario 1b differ from the binomial analysis for Scenario 1b? Consider the coefficient estimates, standard error, predicted probabilities and associated 95% CI.

Question 2

Why are differences between the quasibinomial and binomial models of Scenario 1a less noticeable than the differences in Scenario 1b?

Scenario 2: Random effects model

scenario2_raw <- read_csv("data/scenario-2-raw-data.csv") %>%
  mutate(deformity = factor(deformity),
         dam = factor(dam))

You will need the lme4 R package for random effects models and the broom.mixed package to tidy the output. Once the package is installed, remove eval = F from the code chunk options.

library(lme4)
library(broom.mixed)

random_effects_model <- glmer(deformity ~ dose + (1|dam), 
                          family = binomial, data = scenario2_raw)
tidy(random_effects_model, conf.int = 3) |>
  kable(digits = 3)