library(here)
library(MASS)
library(lme4)
library(lmerTest)
library(patchwork)
library(survival)
library(tidyverse)
source(here("materials/1-workshop1/5-modeling/helpers.R"))
knitr::opts_chunk$set(dev.args = list(bg = "transparent"))
theme_set(
theme_classic(base_size=18) +
theme(
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA),
legend.background = element_rect(fill='transparent'),
legend.box.background = element_rect(fill='transparent')
)
)
set.seed(202606)Modeling and Uncertainty
Learning Outcomes
Explain models as a description of how we think the data were generated, and distinguish between the associated “large” and “small” worlds.
Differentiate between the members of the biostatistical model zoo and explain the types of outcomes and questions to which they are well-suited.
Use R to implement these models, using formula syntax to specify the relationship between covariates and outcomes, and interpret the estimated model parameters in the original problem context.
Diagnose model limitations using visualization and simulation.
Describe how to improve model fit by modifying assumptions (adding covariates, random effects, interactions).
Part 1: What is a Model?
Below are a few views of IL-1a concentration from the workshop dataset. Notice the variation across people and a possible treatment effect at the post-treatment visit.
ids <- read_csv(here("datasets/instructional_dataset/00_sample_ids_UKZN_workshop_2023.csv")) il1 <- read_csv(here("datasets/instructional_dataset/03_elisa_cytokines_UKZN_workshop_2023.csv")) |> filter(cytokine == "IL-1a") |> left_join(ids) head(il1)# A tibble: 6 × 7 sample_id cytokine conc limits pid time_point arm <chr> <chr> <dbl> <chr> <chr> <chr> <chr> 1 SAMP094 IL-1a 174. within limits pid_01 baseline placebo 2 SAMP052 IL-1a 65.7 within limits pid_01 week_1 placebo 3 SAMP017 IL-1a 197. within limits pid_01 week_7 placebo 4 SAMP046 IL-1a 19.7 within limits pid_02 baseline placebo 5 SAMP025 IL-1a 29.9 within limits pid_02 week_1 placebo 6 SAMP050 IL-1a 36.5 within limits pid_02 week_7 placebop1 <- ggplot(il1) + geom_line(aes(time_point, log(conc), group = pid, col = arm)) + scale_y_continuous(expand = c(0, 0)) + scale_x_discrete(expand = c(0, 0)) p2 <- ggplot(il1) + geom_histogram(aes(log(conc), fill = arm), bins = 10, position = "identity", alpha = 0.6) + facet_wrap(~ time_point) (p1 / p2) + plot_layout(guides = "collect")
Imagine running the study again, keeping the recruitment and experimental protocols as similar as possible. Make two lists,
- Aspects of these plots that you think will continue to be true in the new experiment.
- Aspects of these plots that you think will be different.
Take 5 minutes to prepare your description.
A statistical model is a simplified description of how we think the data were generated (the “data story”).
Large vs. small worlds. We should distinguish the “large world” where the data were gathered from the “small world” of the statistical model. The small world makes its assumptions precise, allows hypothetical interventions, and quantifies uncertainty. But there will always be a gap between the two worlds, and that is itself a source of uncertainty not captured by the model.
Modeling Checklist. Organize the work of modeling into these steps.
- Biological Question. A good model should clarify our thinking about a real phenomenon and not be made just for the sake of modeling.
- Plot the Data. What does the outcome look like? Do you already notice any relationships with other variables?
- Working Assumptions. Make a small-world story for how the data were generated. This is the initial statistical model.
- Model Checks. Use simulation and diagnostics to check whether the model captures patterns we see in the real data.
- Revised Models. Revise the working assumptions until the model provides a good match. This model can now be used to answer the original biological question.
Part 2: A Tour through the Model Zoo
Each type of statistical model is designed for a different kind of data or biological question. Knowing all the “species” of model will let you adapt to whatever data or questions you encounter.
Bernoulli and Binomial Models
Bernoulli and binomial models apply when outcomes fall into two groups. The Bernoulli models individual outcomes while the binomial applies to their totals.
- Infected vs. not infected
- Responder vs. non-responder
- Alive vs. dead
A Bernoulli model can simulate new data with frequencies of blue and orange circles similar to (but not exactly the same as) those we originally observed.
A Bernoulli model expects each sample to be associated with a 0/1 label, while the binomial models counts of each group.
The data generating mechanism can let the frequency depend on sample characteristics, e.g., large circles are more likely to be blue.
The model can learn the relationship between the size covariates and probability of the blue outcome and summarizes it as a smooth, increasing curve.
Imagine modeling a disease outbreak in a county. Each person can be infected or not. The probability of infection can depend on the degree of social interaction. This data story is analogous to the large-size and blue circles data story.
infections <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/infections.csv") |> mutate(infected = ifelse(infected == "infected", 1, 0))infections# A tibble: 100 × 5 sample_id social_contact county prob_infection infected <chr> <dbl> <chr> <dbl> <dbl> 1 sample_1 0.748 county_4 0.561 1 2 sample_2 0.379 county_6 0.297 0 3 sample_3 0.483 county_2 0.366 0 4 sample_4 0.813 county_2 0.608 0 5 sample_5 0.670 county_3 0.502 0 6 sample_6 0.400 county_2 0.310 0 7 sample_7 0.0899 county_3 0.151 0 8 sample_8 0.611 county_10 0.458 1 9 sample_9 0.935 county_2 0.691 0 10 sample_10 0.849 county_5 0.634 1 # ℹ 90 more rowsggplot(infections) + geom_jitter(aes(social_contact, infected, col = factor(infected)), height = 0.05) + geom_line(aes(social_contact, prob_infection), linewidth = 1, col = "#232323") + scale_x_continuous(expand = c(0, 0)) + labs( y = "Probability of Infection", x = "Level of Social Contact", col = "Infection Status" )
The binomial model applies when we see totals within categories, but not person-level outcomes.
infections_props <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/infections_props.csv")infections_props# A tibble: 3 × 4 social_level not_infected infected total <chr> <dbl> <dbl> <dbl> 1 (0.0129,0.337] 29 6 35 2 (0.337,0.661] 20 9 29 3 (0.661,0.985] 17 19 36infections_props |> select(-total) |> pivot_longer(-social_level) |> ggplot() + geom_col(aes(social_level, value, fill = name), position = "dodge") + scale_y_continuous(expand = c(0, 0, 0.1, 0)) + labs( x = "Level of Social Contact", y = "Count", fill = "Infection Status" )
To fit these models, we can use the
glmfunction. This function has three arguments,formula: A formula likey ~ xspecifies the relationship between the model outcomeyand contextual variablesx, which we will call “covariates.”family: The type of model to use for the given dataset. This depends on the data type for the outcome (binary, count, etc.) and it differs across the examples in this section.data: Adata.frameobject whose columns give the response and predictor variables to use in the model. Any terms appearing in the formula must be column names in thisdata.frame.
For example, for a Bernoulli model to have probability of infection depend on
social_contact, we can use the command below. Even though we setfamily = binomial(), the fact thatinfectedis a 0/1 variable means that internally the command uses a Bernoulli model1.fit <- glm(infected ~ social_contact, family = binomial(), data = infections) summary(fit)Call: glm(formula = infected ~ social_contact, family = binomial(), data = infections) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.3702 0.5739 -4.130 3.63e-05 *** social_contact 3.0933 0.8908 3.473 0.000515 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 128.21 on 99 degrees of freedom Residual deviance: 113.69 on 98 degrees of freedom AIC: 117.69 Number of Fisher Scoring iterations: 4The
Estimatecolumn describes how infection probability changes across levels of social contact. As given, they are not easily interpretable2. Theplogisfunction converts them to the probability scale. The (Intercept) estimate sets the predicted probability at zero social contact, while thesocial_contactestimate sets how steeply that probability increases with social contact.
plogis(-2.37) # when social contact = 0[1] 0.08548914plogis(-2.37 + 0.5 * 3.09) # when social contact = 0.5[1] 0.3047033plogis(-2.37 + 1 * 3.09) # when social contact = 1[1] 0.672607This Bernoulli model is often called logistic regression, after the inverse logistic function (
plogis) that converts model estimates to the probability scale.The
infectiondata were simulated with a truesocial_contacteffect of 3. The estimate 3.09 is close, but even with 100 samples, slightly off. Nonetheless, it lies within one standard error (0.89) of the truth – the model gives both an estimate of the effect and a measure of its uncertainty.infections$predicted_prob <- plogis(predict(fit)) ggplot(infections) + geom_jitter(aes(social_contact, infected, col = factor(infected)), height = 0.05) + geom_line(aes(social_contact, prob_infection), linewidth = 1, col = "#bbb6b6") + geom_line(aes(social_contact, predicted_prob), linewidth = 1, col = "#414ca1") + scale_x_continuous(expand = c(0, 0)) + labs( y = "Probability of Infection", x = "Level of Social Contact", col = "Infection Status", title = "Original (gray) vs. fitted (blue) infection probability" )
The same function fits the binomial model, but the formula’s left-hand side needs to account for group-level counts, not person-level outcomes. Also, the
weightsargument is added to reflect the different totals in each category.fit <- glm(infected / total ~ social_level, family = binomial(), weights = total, data = infections_props) summary(fit)Call: glm(formula = infected/total ~ social_level, family = binomial(), data = infections_props, weights = total) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.5755 0.4485 -3.513 0.000443 *** social_level(0.337,0.661] 0.7770 0.6019 1.291 0.196703 social_level(0.661,0.985] 1.6868 0.5591 3.017 0.002554 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1.0418e+01 on 2 degrees of freedom Residual deviance: -9.5479e-15 on 0 degrees of freedom AIC: 17.201 Number of Fisher Scoring iterations: 3plogisconverts these estimates to the probability scale. The low category is absorbed into the intercept. For the other two, add their coefficient to the intercept before converting to probabilities for that category.
Poisson and Negative Binomial Models
Poisson and negative binomial models apply to count data. The Poisson arises whenever we count events spread uniformly across space and time, for example, the number of raindrops landing on a window during a minute of rain.
Making a histogram of these counts over time gives us a simple data story. We get a random number of raindrops each minute. Most of the time around five drops land, but a value of zero or ten would not be unusual. Values much larger than that are less likely.
Many biological counts have this flavor.
- Transcripts counts in a cell
- Mutation count along a segment of the genome
- Infection count within a county
Exercise: The last example is again about infections, which we saw in the Bernoulli and binomial examples. What did we know there we don’t know from infection counts alone?
In the data story, we can allow the typical count to depend on sample characteristics. For example, the rain may have been most intense earlier in the night. This figure shows how the probabilities of different number of raindrops progress.
For a more realistic example, note that the infections in counties closer to the epicenter likely have more infections.
outbreak <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/outbreak.csv")outbreak# A tibble: 100 × 5 sample_id epicenter_distance mean county infected <chr> <dbl> <dbl> <chr> <dbl> 1 sample_1 0.952 0.270 county_1 0 2 sample_2 0.213 5.19 county_2 3 3 sample_3 0.498 1.66 county_3 3 4 sample_4 0.307 3.56 county_4 1 5 sample_5 0.888 0.349 county_5 1 6 sample_6 0.605 1.08 county_6 1 7 sample_7 0.930 0.296 county_7 1 8 sample_8 0.737 0.639 county_8 0 9 sample_9 0.399 2.47 county_9 0 10 sample_10 0.104 8.05 county_10 5 # ℹ 90 more rowsggplot(outbreak, aes(epicenter_distance)) + geom_line(aes(y = mean), col = "#414ca1", linewidth = 1) + geom_point(aes(y = infected)) + labs( x = "Epicenter Distance", y = "Infection Count" )
As before,
glmapplies, but we set thefamilyargument topoisson().fit <- glm(infected ~ epicenter_distance, family = poisson(), data = outbreak) summary(fit)Call: glm(formula = infected ~ epicenter_distance, family = poisson(), data = outbreak) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.41503 0.09458 25.54 <2e-16 *** epicenter_distance -3.89306 0.27251 -14.29 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 364.28 on 99 degrees of freedom Residual deviance: 108.18 on 98 degrees of freedom AIC: 304.64 Number of Fisher Scoring iterations: 5The model output is reported on a log scale, so it should be interpreted after exponentiation. The fitted mean at distance 0 is
exp(2.5) = 12.8, matching the plot. Similarly, larger distances are interpreted after plugging into a linear formula and exponentiating. For example, at distance 0.5, the mean isexp(2.5 - .5 * 4.3) = 1.4.outbreak$predicted_mean <- predict(fit, type = "response") ggplot(outbreak, aes(epicenter_distance)) + geom_point(aes(y = infected)) + geom_line(aes(y = mean), col = "#bbb6b6", linewidth = 1) + geom_line(aes(y = predicted_mean), col = "#414ca1", linewidth = 1) + labs( x = "Epicenter Distance", y = "Infection Count", title = "Original (gray) vs. fitted (blue) mean count" )
Exercise: What is the level of uncertainty in these estimates? How would you interpret 0.08 in the
S.E.column? What about 0.29?Often unmeasured factors influence the rates, making the data more heterogeneous than the Poisson can handle. For example, infection might depend on not only distance to the epicenter but also unknown, county-specific population densities. The negative binomial distribution addresses this – it arises as a mixture of Poissons with similar but slightly different rates. In the data simulated below, the counts at a given distance are more spread out even though the blue mean line is unchanged.
outbreak_nb <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/outbreak_nb.csv")head(outbreak_nb)# A tibble: 6 × 5 sample_id epicenter_distance mean county infected <chr> <dbl> <dbl> <chr> <dbl> 1 sample1 0.308 3.55 county_1 6 2 sample2 0.183 5.85 county_2 5 3 sample3 0.860 0.391 county_3 0 4 sample4 0.646 0.920 county_4 0 5 sample5 0.611 1.06 county_5 0 6 sample6 0.473 1.83 county_6 6ggplot(outbreak_nb, aes(epicenter_distance)) + geom_line(aes(y = mean), col = "#414ca1", linewidth = 1) + geom_point(aes(y = infected)) + labs( x = "Epicenter Distance", y = "Infection Count" )
glmdoesn’t provide a negative binomial model, butglm.nbfrom theMASSpackage does, with the same syntax asglm.fit <- glm.nb(infected ~ epicenter_distance, data = outbreak_nb) summary(fit)Call: glm.nb(formula = infected ~ epicenter_distance, data = outbreak_nb, init.theta = 1.888724811, link = log) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.1657 0.2010 10.776 < 2e-16 *** epicenter_distance -3.6485 0.4469 -8.164 3.25e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for Negative Binomial(1.8887) family taken to be 1) Null deviance: 188.04 on 99 degrees of freedom Residual deviance: 105.36 on 98 degrees of freedom AIC: 334.41 Number of Fisher Scoring iterations: 1 Theta: 1.889 Std. Err.: 0.652 2 x log-likelihood: -328.409
Linear Regression Models
When the response is continuous, linear regression applies. For example, it can be used for,
- Antibody titre
- CRP concentration
The model assumes symmetric, random error around a mean value that depends on the inputs. For example, antibody titre might vary as a linear function of age.
antibody <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/antibody.csv")head(antibody)# A tibble: 6 × 5 sample_id age blood_pressure mean antibody <chr> <dbl> <dbl> <dbl> <dbl> 1 sample_1 31.9 99.5 21.6 22.7 2 sample_2 78.2 113. 23.9 27.5 3 sample_3 25.0 118. 21.3 20.2 4 sample_4 31.5 71.0 21.6 17.6 5 sample_5 21.4 120. 21.1 20.3 6 sample_6 45.2 131. 22.3 22.4ggplot(antibody) + geom_line(aes(age, mean), col = "#d53636", linewidth = 1.5) + geom_point(aes(age, antibody)) + labs(x = "Age", y = "Antibody Titre")
Underlying these observations is a linear function with random Gaussian errors.
The
lmfunction fits this model, with the same inputs asglm, but nofamily. The estimates give predicted values at age \(x\) \(19.6 + 0.058 x\). For example, at age 10, the prediction is \(19.6 + 0.58 \approx 20.18\).fit <- lm(antibody ~ age, antibody) summary(fit)Call: lm(formula = antibody ~ age, data = antibody) Residuals: Min 1Q Median 3Q Max -5.345 -1.076 -0.076 1.325 6.698 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 19.57687 0.49545 39.514 < 2e-16 *** age 0.05834 0.01041 5.606 1.91e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.205 on 98 degrees of freedom Multiple R-squared: 0.2428, Adjusted R-squared: 0.2351 F-statistic: 31.42 on 1 and 98 DF, p-value: 1.911e-07What makes Gaussian errors a reasonable assumption? This distribution appears across statistics because of the central limit theorem: the average of many small independent effects tends to a Gaussian distribution. You can find a nice animation of this result here. In this example, we’re implicitly assuming many small, independent factors besides age that perturb the observed antibody measurement.
Survival Models
Survival models describe whether and when an event occurs. For example,
- Time until heart attack
- Time until loan default
- Time until HIV acquisition
If the event doesn’t occur over the study or if they drop out, the participant’s response is called “censored.”
Unlike the other models, whose outcome is a single number, the survival outcome has two parts.
- The
timethe event was observed. Or, if the event was never observed, the time the study ended. - Whether the event occurred (
status = 1) or not (status = 0) before the study end.
For example, imagine modeling HIV acquisition over drug treatments, tracking participants for two years. Some have an event during the study period (
status = 1, triangles) while the rest are censored.surv_data <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/surv_data.csv")head(surv_data)# A tibble: 6 × 4 treatment hazard time status <chr> <dbl> <dbl> <dbl> 1 drug 0.202 2 0 2 control 0.449 0.152 1 3 drug 0.202 0.386 1 4 control 0.449 0.629 1 5 drug 0.202 0.105 1 6 control 0.449 0.996 1surv_data |> mutate(id = row_number()) |> ggplot(aes(time, reorder(id, time), color = treatment, shape = factor(status))) + geom_point(size = 4) + labs( x = "Time", y = "Participant", color = "Treatment", shape = "Event?" ) + theme(panel.grid.major.y = element_line(linewidth = 0.1))
- The
The most popular form of survival model is the Cox proportional hazards model. During any time interval, anyone who hasn’t yet had an event is at risk for it. Some covariates can change the relative risk (“hazard”) of the event. Because a covariate shifts the hazard by the same factor at all times, the model is called a proportional hazards model.
fit <- coxph(Surv(time, status) ~ treatment, data = surv_data) summary(fit)Call: coxph(formula = Surv(time, status) ~ treatment, data = surv_data) n= 50, number of events= 31 coef exp(coef) se(coef) z Pr(>|z|) treatmentdrug -0.8380 0.4326 0.3781 -2.216 0.0267 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 treatmentdrug 0.4326 2.312 0.2062 0.9076 Concordance= 0.589 (se = 0.047 ) Likelihood ratio test= 5.18 on 1 df, p=0.02 Wald test = 4.91 on 1 df, p=0.03 Score (logrank) test = 5.19 on 1 df, p=0.02The model output shows how each covariate rescales the hazard.
Exercise: Can you think of an example where the proportional hazards assumption does not hold? Hint: What if a drug was not effective immediately after treatment start but required several weeks to reach its full effect?
Variations: Multiple Predictors, Mixed Effects, and Interactions
Multiple Predictors. The examples above had a single covariate. We can add more covariates with
+ variableon the right-hand side of the formula. Each variable’s estimate is the effect of that variable with all others held fixed.For example, if we wanted to model antibody titre as a function of both age and blood pressure, we could use the formula
antibody ~ age + blood_pressure. The interpretation of the estimate below is that, fixing blood pressure level, the average antibody level increases by 0.058 units for every year increase in age.fit <- lm(antibody ~ age + blood_pressure, data = antibody) summary(fit)Call: lm(formula = antibody ~ age + blood_pressure, data = antibody) Residuals: Min 1Q Median 3Q Max -5.014 -1.044 -0.175 1.288 6.525 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 21.10179 1.10078 19.170 < 2e-16 *** age 0.05839 0.01033 5.650 1.61e-07 *** blood_pressure -0.01553 0.01003 -1.549 0.125 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.19 on 97 degrees of freedom Multiple R-squared: 0.2611, Adjusted R-squared: 0.2458 F-statistic: 17.14 on 2 and 97 DF, p-value: 4.239e-07
The surface can be visualized in 2D by associating each slice in 3D with a curve.
Mixed Effects. One common situation is that baseline effects (overall probabilities in Bernoulli model, rates in the Poisson, intercepts in the linear regression, …) vary across contexts. For example, infection probability might always increase with social contact, but the baseline rate may vary from county to county, due to unmeasured factors like population density and mobility.
mixed_effect_data <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/mixed_effect_data.csv") |> mutate(infected = ifelse(infected == "infected", 1, 0))head(mixed_effect_data)# A tibble: 6 × 5 social_contact county county_re prob_infection infected <dbl> <chr> <dbl> <dbl> <dbl> 1 0.536 county_4 0.774 0.595 0 2 0.532 county_4 0.774 0.591 1 3 0.854 county_6 -0.358 0.551 1 4 0.127 county_10 0.834 0.313 0 5 0.316 county_2 -0.482 0.177 0 6 0.346 county_8 0.212 0.321 0ggplot(mixed_effect_data) + geom_jitter(aes(social_contact, infected, col = factor(infected)), height = 0.05) + geom_line(aes(social_contact, prob_infection), linewidth = 1, col = "#232323") + facet_wrap(~ reorder(county, prob_infection), nrow = 2) + scale_y_continuous(expand = c(0, 0)) + scale_x_continuous(expand = c(0, 0)) + labs(y = "Probability of Infection")
One solution is to include county as an ordinary covariate. But with many counties this becomes unwieldy – we would have a separate estimate for each. A better approach is to use a random effect (also called a hierarchical model), where the county-level effects are imagined to be drawn from a larger source of randomness. Rather than estimate each county’s effect separately, we estimate the overall variation across counties. Statistically, this is an easier task, and it yields higher power. Because they include both fixed and random effects, these models are called mixed effects models.
To implement this, we replace
glmwithglmerand add(1 | county)to the formula, indicating that county is a random rather than fixed effect.fit <- glmer(infected ~ social_contact + (1 | county), family = binomial(), data = mixed_effect_data) summary(fit)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: infected ~ social_contact + (1 | county) Data: mixed_effect_data AIC BIC logLik -2*log(L) df.resid 230.8 240.7 -112.4 224.8 197 Scaled residuals: Min 1Q Median 3Q Max -2.2664 -0.6116 -0.3655 0.7618 2.5619 Random effects: Groups Name Variance Std.Dev. county (Intercept) 0.5015 0.7082 Number of obs: 200, groups: county, 10 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.4680 0.4783 -5.160 2.47e-07 *** social_contact 3.5811 0.6886 5.201 1.99e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) socil_cntct -0.808To interpret the result, look at the
SDrow. The 0.708 means that estimates vary by roughly this much across counties, at any level ofsocial_contact. Even at social contact 0, the infection probabilities below would be typical.plogis(-2.28 - 0.43) # lower range across counties, social contact = 0[1] 0.06238585plogis(-2.28 + 0.43) # upper range across counties, social contact = 0[1] 0.1358729Interactions. In some data, the effect of a variable on a response might depend on other variables. For example, the age effect might differ between healthy and sick people.
interaction_data <- read_csv("https://huggingface.co/datasets/krisrs1128/durban_workshop_data/resolve/main/interaction_data.csv")head(interaction_data)# A tibble: 6 × 4 age group mean antibody <dbl> <chr> <dbl> <dbl> 1 66.0 disease 4.62 2.20 2 4.19 disease 0.293 3.19 3 28.0 healthy 2.56 5.25 4 79.2 disease 5.54 6.62 5 5.19 healthy 2.10 -0.690 6 33.3 healthy 2.67 3.19ggplot(interaction_data) + geom_line(aes(age, mean, col = group), linewidth = 1.5) + geom_point(aes(age, antibody, col = group)) + labs(x = "Age", y = "Titre")
This type of pattern can be modeled using an interaction term. The R formula syntax has the form
y ~ x1 + x1 : x2, meaning that thex1effect depends on the level ofx2.fit <- lm(antibody ~ age + age : group, data = interaction_data) summary(fit)Call: lm(formula = antibody ~ age + age:group, data = interaction_data) Residuals: Min 1Q Median 3Q Max -6.4765 -1.6597 -0.0768 1.5193 6.1728 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.543777 0.300945 1.807 0.0723 . age 0.061801 0.007332 8.428 7.31e-15 *** age:grouphealthy -0.002380 0.006885 -0.346 0.7300 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.158 on 197 degrees of freedom Multiple R-squared: 0.293, Adjusted R-squared: 0.2858 F-statistic: 40.82 on 2 and 197 DF, p-value: 1.47e-15In the disease group, antibody titre rises by 0.061 units per year of age, while in the healthy group, it changes by 0.061 - 0.002 units (adding the two rows).
Although we introduced these variations with specific Bernoulli and linear regression models, they apply to any model above. They come in handy when revising a model’s working assumptions to improve its fit.
Part 3: Case Studies
Longitudinal Cytokines Analysis
Biological Question. In our study of interest, does treatment with the anti-inflammatory drug change cytokine concentrations post-treatment (week 1), and if so does the effect persist into the long-run (week 7)?
Plot the Data. Roughly, there’s a strong treatment effect one week post-intervention that weakens by week 7. We’ll try capturing this pattern in a small world model. We focus on IL-1a, the cytokine with the fewest measurements below the limit of detection.
ids <- read_csv(here("datasets/instructional_dataset/00_sample_ids_UKZN_workshop_2023.csv")) cytokines <- read_csv(here("datasets/instructional_dataset/03_elisa_cytokines_UKZN_workshop_2023.csv")) |> left_join(ids) ggplot(cytokines) + geom_line(aes(time_point, log(conc), group = pid, col = arm)) + scale_y_continuous(expand = c(0, 0)) + scale_x_discrete(expand = c(0, 0)) + facet_wrap(~ reorder(cytokine, -conc), nrow = 2)
Working Assumptions. Luminex gives a continuous readout of each protein’s concentration, which suggests a linear regression model.
The Gaussian error assumption suits Luminex well: many small, unrelated factors in the measurement or biology could make the fluorescence slightly brighter or dimmer. This also suggests that we model log-concentration instead of raw concentration, since errors are likely multiplicative rather than additive.
Model Checks. The code below fits a regression model of log-concentration on treatment status. The treatment effect is significant, but we know to interpret a model only after checking its fidelity.
cytokines <- cytokines |> filter(cytokine == "IL-1a") fit <- lm(log(conc) ~ arm, data = cytokines) summary(fit)Call: lm(formula = log(conc) ~ arm, data = cytokines) Residuals: Min 1Q Median 3Q Max -2.8298 -0.7570 -0.0689 0.6711 3.7742 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.4001 0.1276 34.486 < 2e-16 *** armtreatment -1.2410 0.1847 -6.719 5.18e-10 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.06 on 130 degrees of freedom Multiple R-squared: 0.2578, Adjusted R-squared: 0.2521 F-statistic: 45.15 on 1 and 130 DF, p-value: 5.176e-10To this end, let’s simulate from the fitted model (full in
helpers.R).sim_data <- simulate_cytokine(fit) plot_cytokine(sim_data)
Reassuringly, the overall concentrations are in the correct range, and the control group has larger concentrations than the treatment group, matching the real data. But there is an issue – in the real data, most blue curves decrease at week 1, while in our simulation about half increase.
Revised Models. We shouldn’t be surprised the model missed the week 1 decrease: we gave it no temporal information. One potential fix is the formula
~ arm + time_point.But considering the large world data story, this model isn’t quite right. We don’t expect IL-1a concentration to change just because the clock reads week 1 or 7. We want the treatment effect itself to depend on the week. This is an interaction,
~ arm + arm : timepoint.The fitted model is printed below. The parameters give the expected log-concentrations,
- Placebo at baseline: 4.66
- Treatment at baseline: 4.66 - 0.29
- Placebo at week 1: 4.66 - 0.48
- Treatment at week 1: 4.66 - 2.34 …
In particular the large negative coefficient for the
armtreatment:time_pointweek_1term captures our earlier observation that the treatment effect seems strongest at week 1.fit <- lm(log(conc) ~ arm + arm : time_point, data = cytokines) summary(fit)Call: lm(formula = log(conc) ~ arm + arm:time_point, data = cytokines) Residuals: Min 1Q Median 3Q Max -2.36776 -0.57246 0.00535 0.47567 2.64101 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.6636 0.1718 27.147 < 2e-16 *** armtreatment -0.2855 0.2487 -1.148 0.2531 armplacebo:time_pointweek_1 -0.4784 0.2429 -1.969 0.0511 . armtreatment:time_pointweek_1 -2.3398 0.2543 -9.203 9.50e-16 *** armplacebo:time_pointweek_7 -0.3121 0.2429 -1.285 0.2013 armtreatment:time_pointweek_7 -1.3173 0.2543 -5.181 8.51e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.8239 on 126 degrees of freedom Multiple R-squared: 0.5653, Adjusted R-squared: 0.5481 F-statistic: 32.77 on 5 and 126 DF, p-value: < 2.2e-16Exercise: Does this model allow for a treatment effect at time zero? Based on what you know about the “large world” experiment, decide whether this is appropriate or if the model could be further improved.
The simulated data look better than before.
sim_data <- simulate_cytokine(fit) plot_cytokine(sim_data)
One more revision. The model doesn’t account for longitudinality, which could inflate our confidence. It treats each sample as independent, when measurements from the same participant may be dependent. The effective sample size lies somewhere between the number of study participants (44) and measurements (132).
We have the machinery to deal with this: a participant-level random effect, allowing a different baseline for each participant.
fit <- lmer(log(conc) ~ arm + arm : time_point + (1 | pid), data = cytokines) summary(fit)Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmerModLmerTest] Formula: log(conc) ~ arm + arm:time_point + (1 | pid) Data: cytokines REML criterion at convergence: 327.2 Scaled residuals: Min 1Q Median 3Q Max -2.8408 -0.6889 -0.0017 0.5540 3.1092 Random effects: Groups Name Variance Std.Dev. pid (Intercept) 0.02198 0.1483 Residual 0.65679 0.8104 Number of obs: 132, groups: pid, 44 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 4.6636 0.1718 125.7363 27.147 < 2e-16 *** armtreatment -0.2855 0.2487 125.7363 -1.148 0.2531 armplacebo:time_pointweek_1 -0.4784 0.2390 83.9999 -2.002 0.0485 * armtreatment:time_pointweek_1 -2.3398 0.2501 83.9999 -9.355 1.14e-14 *** armplacebo:time_pointweek_7 -0.3121 0.2390 83.9999 -1.306 0.1952 armtreatment:time_pointweek_7 -1.3173 0.2501 83.9999 -5.267 1.05e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) armtrt armp:__1 armt:__1 armp:__7 armtreatmnt -0.691 armplcb:__1 -0.696 0.481 armtrtm:__1 0.000 -0.503 0.000 armplcb:__7 -0.696 0.481 0.500 0.000 armtrtm:__7 0.000 -0.503 0.000 0.500 0.000Visually, the model doesn’t seem much better. Indeed, the participant-level variation (SD 0.15) is small compared to the within-sample variation (0.81). Nonetheless, it’s the safer model in practice, because by accounting for within-participant dependence, we are being honest about the effective sample size.
sim_data <- simulate_cytokine(fit) plot_cytokine(sim_data)
Flow Cytometry Cell Compositions
Biological Question. Does the anti-inflammatory drug influence immune cell composition at week 1 post-treatment?
Let’s read in the data. Since our question focuses on
week_1, we keep only these sample IDs. Theright_jointhen drops the other timepoints.ids <- read_csv(here("datasets/instructional_dataset/00_sample_ids_UKZN_workshop_2023.csv")) |> filter(time_point == "week_1") flow <- read_csv(here("datasets/instructional_dataset/04_flow_cytometry_UKZN_workshop_2023.csv")) |> right_join(ids) |> mutate(arm = as.factor(arm)) head(flow)# A tibble: 6 × 13 sample_id live_cd19_negative cd45_negative cd45_positive neutrophils <chr> <dbl> <dbl> <dbl> <dbl> 1 SAMP052 347091 259075 86122 53132 2 SAMP025 115716 27779 83689 56681 3 SAMP040 576299 538725 36244 23491 4 SAMP101 526335 445947 78615 51813 5 SAMP096 239933 187379 49636 41053 6 SAMP087 110029 73529 35053 29168 # ℹ 8 more variables: non_neutrophils <dbl>, cd3_negative <dbl>, # cd3_positive <dbl>, cd4_t_cells <dbl>, cd8_t_cells <dbl>, pid <chr>, # time_point <chr>, arm <fct>Plot the data. We’ve reshaped the flow cytometry data into each sample’s cell type composition, shown as a stacked bar plot below. Samples are sorted so that those with lower cell type diversity appear on the left.
reshaped_flow <- flow |> pivot_longer(matches("cd|neutro"), names_to = "cell_type", values_to = "count") |> group_by(sample_id) |> mutate(proportion = count / sum(count)) ggplot(reshaped_flow) + geom_col(aes(reorder(sample_id, proportion, entropy), proportion, fill = cell_type), col = "black", width = 1, linewidth = 0.3) + facet_wrap(~ arm, scales = "free_x") + scale_y_continuous(expand = c(0, 0)) + labs( x = "Sample ID", fill = "Cell Type", y = "Proportion" ) + theme(axis.text.x = element_text(angle = 90))
Exercise: Replace
proportionwithcountin thegeom_colcommand above. What do you see? What could be some corresponding “large world” explanations?The main question is whether there is a systematic difference in cell compositions between the treatment and control at week 1. If there is a difference, which cell types are responsible?
Before we build our small world model to answer this, it’s worth considering the “large world” data story,
- There are on the order of a trillion immune cells (of various cell types) circulating through a typical body at any moment. Some cell types are much rarer than others, and there are more cell types than we can assay with our flow cytometer. Cells are being created and destroyed at all times and they are not evenly spatially distributed, with concentrations changing in response to signaling molecules and random drift.
- At a doctor’s office, a CVL is collected, catching some of the cells from the previous step that just happened to be present at the right place and time. A technician takes the collected sample and freezes it, killing all the cells present.
- After a long time in the freezer, the sample is prepared for the flow cytometer. Most (but not all) cell-surface markers are bound by labeled antibodies, and the device records the associated profiles. All cells are destroyed in the process; from here on out it’s just data.
- A computational biologist gates the raw cytometer output into distinct cell types. Some don’t fall anywhere we recognize and so are dropped from subsequent analysis. The data are saved to a CSV file and uploaded to the workshop website, living on a server in a data center in Greenland. Three years later, you download a copy of the data to your computer.
This is a big world, and our final file is an imperfect view of the participants’ original immune landscapes. Sampling variation and measurement error enter at each step.
Working Assumptions. Our small world is quite simple. Consider the population of immune cells in one participant. When we draw a sample, we “catch” a mix of CD3- and other cell types. With only two groups (CD3- and “other”) and only totals observed per participant, this is a natural setting for a binomial model.
flow <- flow |> mutate( total = rowSums(across(2:10)), proportion = cd3_negative / total ) fit <- glm(proportion ~ 1, family = binomial(), weights = total, data = flow) summary(fit)Call: glm(formula = proportion ~ 1, family = binomial(), data = flow, weights = total) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.263195 0.001355 -3146 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 387427 on 43 degrees of freedom Residual deviance: 387427 on 43 degrees of freedom AIC: 387907 Number of Fisher Scoring iterations: 5After the inverse logistic transformation, the estimate is just the fraction of observed cells assigned to CD3-.
plogis(-4.26)[1] 0.01392564sum(flow$cd3_negative) / sum(flow$total)[1] 0.01388184Model Checks. We simulate new data using the fitted model. The proportions are the CD3- compositions of hypothetical samples.
y_sim <- simulate(fit)$sim_1 head(y_sim)[1] 0.01386941 0.01351601 0.01383336 0.01397566 0.01386085 0.01430988In the simulated data, the CD3- proportions are close to but not exactly those in the real data, which reflects the randomness of the binomial model.
We’ve written a helper,
fit_cytof, that loops this model over all cell types and simulates counts for each. When we visualize the simulated data, two issues stand out,- The simulated proportions are very close to the estimated proportions across all the samples. The binomial model has much lower variability than observed in the actual data.
- Our model doesn’t allow the fitted cell-type proportions to vary across treatment and control.
We’ll start with the second issue, because it’s easier to address. Of course, we’ll need to answer the first before trusting this model.
m <- fit_cytof(proportion ~ 1) counts <- simulate_counts(m) plot_counts(counts)
Revised Models. To let cell-type proportions depend on treatment, we modify the formula from
proportion ~ 1toproportion ~ arm.m <- fit_cytof(proportion ~ arm)The estimated coefficients with treatment status as a covariate are printed below. The row for
cd3_negativemeans that placebo samples have aboutplogis(-3.98) = 0.018of their composition from CD3- cells. Under treatment this decreases to aboutplogis(-3.98 - 0.578) = 0.010.coef_df <- map_dfr(m, coef, .id = "cell_type") coef_df# A tibble: 9 × 3 cell_type `(Intercept)` armtreatment <chr> <dbl> <dbl> 1 live_cd19_negative -0.283 0.0207 2 cd45_negative -0.773 -0.00320 3 cd45_positive -2.11 0.0531 4 neutrophils -2.48 0.221 5 non_neutrophils -3.45 -0.513 6 cd3_negative -3.98 -0.578 7 cd3_positive -4.54 -0.444 8 cd4_t_cells -5.18 -0.482 9 cd8_t_cells -5.66 -0.338We could in principle interpret changes in composition due to treatment from this model. But we shouldn’t trust this model yet because it still fails to account for sample-to-sample variability.
counts <- simulate_counts(m) plot_counts(counts)
To account for sample-to-sample differences, we add a random effect. Since the sample labels are kept in the
sample_idcolumn, we modify the formula toproportion ~ arm + (1 | sample_id)and switch from
glmtoglmer.The new simulated data are much more realistic. They still don’t quite have the smooth variation in composition across samples seen in the real data – this could be introduced through latent factors or modeling correlation across cell types. Even so, the current model is often sufficient for real data analyses, and underlies the
diffcytBioconductor package Weber et al. (2019).fits <- fit_cytof(proportion ~ arm + (1 | sample_id), glmer) counts <- simulate_counts(fits) plot_counts(counts)
The model output shows how the treatment affects each cell type (
armtreatment) and how much it varies across sample (sample_idStd. Dev). For example, under treatment the CD3- composition tends to decrease relative to the other cell types.summary(fits[["cd3_negative"]])Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: proportion ~ arm + (1 | sample_id) Data: flow2 Weights: total AIC BIC logLik -2*log(L) df.resid 928.7 934.0 -461.3 922.7 41 Scaled residuals: Min 1Q Median 3Q Max -0.129244 -0.003768 0.002199 0.006118 0.011028 Random effects: Groups Name Variance Std.Dev. sample_id (Intercept) 1.117 1.057 Number of obs: 44, groups: sample_id, 44 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.2288 0.2204 -19.188 <2e-16 *** armtreatment -0.6108 0.3190 -1.915 0.0555 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) armtreatmnt -0.691Exercise: According to the fitted model, what is the typical proportion of CD3- cells in the control and treatment groups? Hint: Use the
plogisfunction applied to the right choice of values underEstimate.
References
Weber, Lukas M., Malgorzata Nowicka, Charlotte Soneson, and Mark D. Robinson. 2019. “Diffcyt: Differential Discovery in High-Dimensional Cytometry via High-Resolution Clustering.” Communications Biology 2 (1). https://doi.org/10.1038/s42003-019-0415-5.
Footnotes
The technical reason is that the Bernoulli distribution is a special case of the binomial when the total count parameter is set to 1.↩︎
Because they are given on the log-odds scale. The model estimation algorithm works on this scale because it isn’t constrained to \(\left[0, 1\right]\), and so general-purpose optimization methods can be used.↩︎