thrombectomy <- medisr::thrombectomy |>
mutate(
death_1year_followup = if_else(time_to_death < 365 & death ==1,1,0),
age_at_diagnosis = as.numeric(age_at_diagnosis)
)
th_model <- glm(
death_1year_followup ~ thrombectomy +
prior_diabetes +
age_at_diagnosis +
prior_hypertention +
sex +
thrombectomy:prior_diabetes +
thrombectomy:age_at_diagnosis +
thrombectomy:prior_hypertention +
thrombectomy:sex,
family = binomial(link = "logit"),
data = thrombectomy
)Marginalizing the problem
Case 15: “The odds ratio is 0.65… and so what?”
At the stroke research meeting, the biostatistician puts a slide on the screen.
Thrombectomy vs. no thrombectomy: OR = 0.65
The clinical team pauses.
“That sounds promising,” someone says. “But what does 0.65 actually mean for a patient?”
The analysis comes from an observational study of patients with acute stroke. The outcome is death within one year, and the researchers have fitted a logistic regression model to account for differences between patients who received thrombectomy and those who did not.
The statistician explains that the regression coefficient corresponds to an odds ratio of 0.65.
Another question immediately follows:
“But what we really want to know is whether patients have a lower probability of dying if they receive thrombectomy. Why are we talking about odds?”
Now the room gets filled with discussion.
Suppose 20% of patients in the control group die within one year. What would you actually want to show a clinician—or a patient?
- An odds ratio of 0.65?
- The probability of death under each treatment?
- The absolute difference in probability of death?
- Something else?
And whose probability are we talking about?
- A particular patient with a particular age, diabetes status, and clinical profile?
- An older adult with diabetes?
- The average patient in this study?
- Or the population-level probability we would see if everyone received thrombectomy, compared with what we would see if no one received it?
These are related questions, but they are not the same question. The clinical team decides to work through the problem one scientific question at a time.
- Is thrombectomy associated with a different probability of death?
- Does the effect differ for older adults with diabetes?
- What does the model imply for the average stroke patient?
- What if everyone could participate twice — both in the treatment and the no-treatment group?
The team notices something — the regression model has produced one coefficient, but the researchers have several different scientific questions.
So they ask:
What exactly is the model estimating—and how do we turn that into the quantity we actually want to know?
Your task is to make sure that the quantity we calculate and communicate answers the scientific question we started with.
The unfortunate truth…
Reflection points
Reading Task: 15.4
Communicating results of regression models
Discussion: “What do we actually want to know?
An observational study retrospecively investigates whether a new treatment reduces the risk of a complication after surgery. The researchers fit a logistic regression model and obtain 𝛽1 (odds ratio) as 0.65 for treatment vs. control.
Discuss with our neighbor:
- What is the scientific question? What would you actually want to know about the treatment?
- What is the model giving us? What exactly is the odds ratio of 0.65? Is it itself a clinical outcome or something else?
- What would be easier to interpret? Imagine the complication occurs in 20% of patients in the control group. What kind of result would you rather show a clinician or patient: the odds ratio, a risk/probability, an absolute risk difference, or something else? Why?
- What is the estimand? Try to state, in one sentence, the quantity you would ideally like to learn about.
- If the model coefficient is not the scientific answer, what role does the model play?
Workflow
The dataset
Read for yourself…
Scientific questions
Does the population who got thrombectory have lower probability of death at 1-year follow-up after stroke compared to the population who did not get thrombectomy?
What is the effect of thrombectory in older adults with diabetes, compared to older adults without diabetes?
What is effect of the average person profile in the sample?
What would be the average probability of death, had all participants received thrombectomy, compared to not receiving it?
The estimand
We are interested in the probability of death at 1-year follow-up:
\[ P(death = 1 | X) \]
for an observational study looking at the effect of thrombectomy after stroke:
The estimand
- We need to adjust for age, hypertension and diabetes.
- We will still incl. sex, as the dataset is unbalanced and we want to adjust for that bias.
The estimator
One is better
The estimator
A matter of choice
Our estimator
- follow individuals for one year
- no censoring
- probability of death is relative high
- n can be defined as the sample size
\[ y_i = Binomial(n,p_i) \\ logit(p_i) = α + βx_i \]
where n is the number of observations and p is the probability of death.
Our estimator
The logit function is defined as the log-odds:
\[ logit(p_i) = log(\frac{p_i}{1 - p_i}) \]
The “odds” of an event is just the probability it happens divided by the probability it does not happen. In other words:
\[ log(\frac{p_i}{1 - p_i}) = α + βx_i \] To figure out the definition of p_i, we do some algebra to solve for \(p_i\):
\[ p_i = \frac{exp( α + βx_i)}{1 + exp(α + βx_i)} \] The above function is the logistic
- Note that the log-odds space is linear (what we model)
- The logistic space is the sigmoidal shape (the transformation)
The statistical model
\[ logit( \text{death at 1 year followup}) = α + β1_{thrombectomy} + β2_{diabites} + β3_{age} \\ + β4_{hypertention} + β5_{sex} + β6_{thrombectomy×diabites} \\ + β7_{thrombectomy×age} + β8_{thrombectomy×hyptertension} + β9_{thrombectomy×sex} \\ \]
In R
thrombectomy <- medisr::thrombectomy |>
mutate(death_1year_followup = if_else(time_to_death < 365 & death ==1,1,0),
age_at_diagnosis = as.numeric(age_at_diagnosis)
)
th_model <- glm(
death_1year_followup ~ thrombectomy +
prior_diabetes +
age_at_diagnosis +
prior_hypertention +
sex +
thrombectomy:prior_diabetes +
thrombectomy:age_at_diagnosis +
thrombectomy:prior_hypertention +
thrombectomy:sex,
family = binomial(link = "logit"),
data = thrombectomy
)The coefficients
# Fixed Effects
Parameter | Odds ratio | SE
-------------------------------------------------------
(Intercept) | 1.82e-03 | 1.55e-03
thrombectomy | 9.06 | 10.29
prior_diabetes | 2.70 | 0.74
age_at_diagnosis | 1.07 | 0.01
prior_hypertention | 1.48 | 0.30
sex | 1.04 | 0.21
thrombectomy:prior_diabetes | 1.64 | 0.76
thrombectomy:age_at_diagnosis | 0.97 | 0.01
thrombectomy:prior_hypertention | 0.43 | 0.14
thrombectomy:sex | 1.34 | 0.43
Read Step 3: Stock with regression coefficients…
Mistake #1
Don’t read coefficient, thinking it generalizes to the sample!
Mistake #2
Confusing OR for a probability
Mistake #3
1-unit change OR gives different probability for two individuals
Read Step 4: Clear communication — transforming coefficients
- Run all code examples and
- Make predictions for
- A 70 year-old female without diabetes, hyptersions and thrombectomy
- A 70 year-old male without diabetes but having hyptersions and received thrombectomy
Step 5-7: Model-based transformations
A model is conditional on it’s sample
Regression coefficients let’s us predict any arbitrary sample that is conditional on being explained by the model!
post-hoc weights on certain participant profiles, so that our estimate is targeting a specific aspect of the model inference.
Questions that form such discussions can be:
- A particular patient?
- A clinically meaningful scenario?
- The average population?
- Different levels of an interacting variable?
RQ 1: Our sample as reality
“Does the population who got thrombectory have lower probability of death at 1-year follow-up after stroke compared to the population who did not get thrombectomy?”
Can we use the model coefficient for thrombectomy? 👀
# Fixed Effects
Parameter | Odds ratio | SE
-------------------------------------------------------
(Intercept) | 1.82e-03 | 1.55e-03
thrombectomy | 9.06 | 10.29
prior_diabetes | 2.70 | 0.74
age_at_diagnosis | 1.07 | 0.01
prior_hypertention | 1.48 | 0.30
sex | 1.04 | 0.21
thrombectomy:prior_diabetes | 1.64 | 0.76
thrombectomy:age_at_diagnosis | 0.97 | 0.01
thrombectomy:prior_hypertention | 0.43 | 0.14
thrombectomy:sex | 1.34 | 0.43Here, the odds ratio is 2.20 …
when all interacting variables are at their reference/zero values. In particular, the interpretation depends on how age is coded/centered and on the reference categories for diabetes, hypertension, and sex.
But how intuitive is that? 👀
If we don’t supply a newdata grid, the predictions() function will take all participants in our sample and compute predictions for them, one by one:
Estimate Pr(>|z|) S 2.5 % 97.5 %
0.2677 <0.001 28.9 0.2083 0.337
0.0706 <0.001 40.2 0.0361 0.133
0.3220 <0.001 19.9 0.2606 0.390
0.1979 <0.001 36.8 0.1417 0.269
0.1275 <0.001 40.1 0.0794 0.198
--- 1013 rows omitted. See ?print.marginaleffects ---
0.4588 0.683 0.6 0.2776 0.652
0.6696 0.102 3.3 0.4651 0.825
0.3413 0.137 2.9 0.1788 0.552
0.4494 0.617 0.7 0.2694 0.644
0.6611 0.118 3.1 0.4575 0.819
Type: invlink(link)You can see that the row length is the same as the original dataset:
- A table with 1023 rows is not super intuitive and nice to make inference from.
Why can we do about that?
- Make a plot to illustrate it all visually, or
- take the average of all these predictions
Design is not balanced!
0 1
583 440 about 0.4301075 received thrombectomy.
Marginalizing the problem
We can compute the average for them by using avg_predictions()
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
0.267 0.013 20.5 <0.001 307.1 0.241 0.292
Type: responseSo the average probability of death at the 1-year follow-up is 26.7%.
This is a reasonable estimate for answering the question
“how many Danes die at 1-year follow-up after stroke?”
However, as we now know that the groups are unbalanced, this is not a good estimate for the treatment effect! And it doesn’t have the granularity to address more nuanced questions of sub-population level.
Finally, to compute the average for the patients receiving thrombectomy and those who did not, we can group the average for these conditions
Can we do better?
the risk difference or the risk ratio.
\[ RD = 0.214 - 0.307 = -9.3 \] or a 9.3% points lower predicted risk of death on the absolute scale.
\[ RR = \frac{0.214}{0.307} = 0.697 \] or a 30% lower predicted risk of death on the relative scale.
So we conclude that:
The model-estimated 1-year mortality was 30.7% in the absence of thrombectomy and 21.4% under thrombectomy, corresponding to an un-adjusted absolute risk difference of −9.3 percentage points
RQ 2: Treatment effect of older adults conditional on diabetes status
“What is the effect of thrombectory in older adults with diabetes, compared to older adults without diabetes?”
Here, we’re interested specifically in a sub-population.
we need to consider multiple variables (diabetes and thrombectory) and the proportions who have this profile in our sample.
Making new grids
prior_diabetes thrombectomy age_at_diagnosis Estimate Pr(>|z|) S 2.5 %
1 0 80 0.587 0.2257 2.1 0.446
1 0 81 0.603 0.1539 2.7 0.461
1 0 82 0.619 0.1020 3.3 0.476
1 0 83 0.634 0.0659 3.9 0.491
1 0 84 0.649 0.0417 4.6 0.506
1 0 85 0.664 0.0259 5.3 0.520
97.5 %
0.716
0.730
0.743
0.757
0.770
0.782
Type: invlink(link)Estimated Marginal Mean
thrombectomy prior_diabetes Estimate Std. Error z Pr(>|z|) S 2.5 %
0 0 0.423 0.0455 9.30 <0.001 65.9 0.334
0 1 0.662 0.0674 9.82 <0.001 73.2 0.530
1 0 0.267 0.0490 5.46 <0.001 24.3 0.171
1 1 0.616 0.0889 6.93 <0.001 37.8 0.442
97.5 %
0.513
0.794
0.363
0.790
Type: responsefrom 42.3% probability of death without treatment to 26.7% for non-diabetic.
from 66.2% probability of death without treatment to 61.6% for non-diabetic.
RQ:3 “The average person”
What is effect of the average person profile in the sample?
Our sample is unbalanced
First of all, males are nearly 48% propensity of receiving thrombectory or not, while it’s only 38% of females.
prior_diabetes
thrombectomy 0 1
0 512 71
1 400 40- more males in the treatment group
- few with diabetes — 10-15% in our sample
- 6.4% in DK population
- A balanced design should have 50%
The balanced grid
thrombectomy Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
0 0.334 0.0296 11.29 <0.001 95.8 0.276 0.392
1 0.365 0.0432 8.45 <0.001 54.9 0.280 0.450
Type: response::: {.r-fit-text .} Is “the average person”! 👀 :::
RQ 4: “What if” 🤔
What would be the average probability of death, had all participants received thrombectomy, compared to not receiving it?
thrombectomy Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
0 0.271 0.0168 16.1 <0.001 191.8 0.238 0.304
1 0.239 0.0212 11.3 <0.001 95.6 0.198 0.281
Type: responseFor each individual in the Syntheric DanStroke medisr sample, what is the predicted probability of death in the counterfactual world where they receive thrombectory, and where they do not?
Can we do better?
Risk difference
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
-0.0317 0.0271 -1.17 0.242 2.0 -0.0847 0.0214
Term: thrombectomy
Type: response
Comparison: 1 - 0A risk difference of -0.0317, favoring thrombectomy
The comparisons() function above computed the predicted probability of mortality (death_1year_followup == 1) for each observed row of the data in two counterfactual cases: when thrombectomy == 1, and when thrombectomy == 0. Then, it computed the differences between these two sets of predictions. Finally, it computed the population-average of risk differences.
population-averaged (marginal) adjusted risk ratios:
Estimate Pr(>|z|) S 2.5 % 97.5 %
0.883 0.251 2.0 0.714 1.09
Term: thrombectomy
Type: response
Comparison: ln(mean(1) / mean(0))The model-estimated 1-year mortality was 27.1% in the absence of thrombectomy and 23.9% under thrombectomy, corresponding to an adjusted risk ratio of 0.883.
Time for exercises 🎉
