Case 4 · Positive and negative welfare screens

Check a logistic fit and decide what the observed screen rates support.

All four cases · Diagnostic report

Use this page in class or as a reference later. Complete each formula in your R console or write it down before revealing the solution. Use the displayed output to answer the interpretation questions before checking your answers. To run the analysis, start from the supplied Posit Cloud project root and work through the R blocks in order. If you have yet to fit a model yourself, run the fit in its solution callout before using it in the next block.

Load the diagnostic packages

Run this block in Posit Cloud before starting the worksheet. If RStudio prompts you to install a missing package, install it and run the block again. You can also install both packages from the Console with install.packages(c("DHARMa", "see")).

library("DHARMa")
library("see")

Read the study and prepare the data

Researchers housed 24 animals separately in four rooms, with three animals per diet in each room. Each animal had one welfare screen. In screened_positive, 1 means a positive screen and 0 means a negative screen. The researchers want to compare the probability of a positive screen between diets while allowing the rooms to differ. These data are invented for teaching.

Set Control and R1 as the reference categories.

library("dplyr")

d <- readr::read_csv("data/case-4.csv")
d <- d |>
  mutate(
    diet = factor(diet, levels = c("Control", "High fibre")),
    room = factor(room, levels = c("R1", "R2", "R3", "R4"))
  )

head(d, 3)
# A tibble: 3 × 4
  room  animal diet    screened_positive
  <fct> <chr>  <fct>               <dbl>
1 R1    R1-A01 Control                 1
2 R1    R1-A02 Control                 1
3 R1    R1-A03 Control                 1

Fit a logistic model

A logistic model describes the probability of a positive screen. Its coefficients express changes in log odds.

NoteWrite the starting formula

Which column records the outcome? Which two terms compare diets while allowing different probabilities across rooms? Complete the formula.

Formula template. Fill in the blanks before running it.

original <- glm(___ ~ ___, family = binomial(), data = d)

Diet compares High fibre with Control, and room allows each room to have a different baseline.

original <- glm(
  screened_positive ~ diet + room,
  family = binomial(),
  data = d
)

summary(original)

Call:
glm(formula = screened_positive ~ diet + room, family = binomial(), 
    data = d)

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)     2.557e+01  9.859e+04   0.000        1
dietHigh fibre -5.113e+01  8.818e+04  -0.001        1
roomR2          7.483e-07  1.247e+05   0.000        1
roomR3          7.483e-07  1.247e+05   0.000        1
roomR4          7.483e-07  1.247e+05   0.000        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3.3271e+01  on 23  degrees of freedom
Residual deviance: 3.7848e-10  on 19  degrees of freedom
AIC: 10

Number of Fisher Scoring iterations: 24

Compare the observed screens with the fit

The numerical convergence check asks whether the fitting calculation finished. The range gives the smallest and largest fitted probabilities. The table counts positive and negative screens within each diet.

performance::check_convergence(original)
[1] TRUE
range(fitted(original))
[1] 7.88492e-12 1.00000e+00
screen_counts <- d |>
  summarise(
    positive = sum(screened_positive == 1),
    negative = sum(screened_positive == 0),
    .by = diet
  )
screen_counts
# A tibble: 2 × 3
  diet       positive negative
  <fct>         <int>    <int>
1 Control          12        0
2 High fibre        0       12
NoteDecide whether to trust the diet coefficient

For each diet, how many screens were positive and how many were negative? Look at the fitted probabilities, diet coefficient and its standard error. How could the fit report convergence and still leave the room-adjusted diet comparison unresolved?

All 12 Control animals screened positive and all 12 High fibre animals screened negative. Diet perfectly predicts the observed outcomes, a pattern called complete separation. The fitted probabilities can get closer and closer to the observed zeros and ones as the diet coefficient grows in magnitude. The large coefficient and standard error signal that the ordinary fit has no finite room-adjusted diet estimate, even though the fitting calculation stopped.

Report the observed rates

Pool the rooms and calculate the positive-screen rate for each diet. Each exact binomial interval assumes unrelated outcomes and a common screen probability within that diet across the four rooms.

screen_counts |>
  mutate(rate = positive / (positive + negative))
# A tibble: 2 × 4
  diet       positive negative  rate
  <fct>         <int>    <int> <dbl>
1 Control          12        0     1
2 High fibre        0       12     0
control <- screen_counts |> filter(diet == "Control")
fibre <- screen_counts |> filter(diet == "High fibre")

binom.test(control$positive, control$positive + control$negative)

    Exact binomial test

data:  control$positive and control$positive + control$negative
number of successes = 12, number of trials = 12, p-value = 0.0004883
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.7353515 1.0000000
sample estimates:
probability of success 
                     1 
binom.test(fibre$positive, fibre$positive + fibre$negative)

    Exact binomial test

data:  fibre$positive and fibre$positive + fibre$negative
number of successes = 0, number of trials = 12, p-value = 0.0004883
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.0000000 0.2646485
sample estimates:
probability of success 
                     0
NoteWrite a two-sentence result

Report the observed proportion and 95% interval for each diet. What do those two intervals describe? What part of the original question remains unresolved?

The observed positive-screen rates were 12/12 for Control (100%, 95% interval 73.5% to 100%) and 0/12 for High fibre (0%, 95% interval 0% to 26.5%). Each interval describes a single diet’s screen probability when rooms are pooled. The size of the room-adjusted diet difference remains unresolved under the ordinary logistic model.