library(tidyverse)
library(haven)
library(dosresmeta)
library(rms)
library(Epi)
theme_set(theme_light())Orsini N. Weighted mixed-effects dose–response models for tables of correlated contrasts
Summary
This html page reproduces the main examples presented by Orsini in the article Weighted mixed-effects dose-response models for tables of correlated contrasts, Stata Journal, which can be found at the following link.
Reference: Orsini, N. Weighted mixed-effects dose-response models for tables of correlated contrasts. Stata Journal. 2021, Vol.21 (2), p.320-347.
Packages and data
Mean Differences
md_drm <- haven::read_dta(url("http://www.stats4life.se/data/md_drm.dta")) %>%
mutate(
md_low = md - qnorm(.975)*semd,
md_upp = md + qnorm(.975)*semd
)
head(md_drm)# A tibble: 6 × 13
id dose md semd sd n meany mindose maxdose dosesq type_md
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 2.09 0 0 10.2 667 6.14 0.108 3.32 4.35 4
2 1 4.42 -1.83 0.543 9.63 667 4.31 3.32 5.71 19.5 4
3 1 8.50 -0.713 0.571 10.6 666 5.43 5.71 24.4 72.3 4
4 2 2.09 0 0 10.0 334 7.44 0.129 3.21 4.37 4
5 2 4.35 -2.24 0.790 10.4 333 5.20 3.21 5.81 18.9 4
6 2 8.57 4.00 0.838 11.6 333 11.4 5.82 20.1 73.5 4
# ℹ 2 more variables: md_low <dbl>, md_upp <dbl>A plot of the empirical estimates
ggplot(data = filter(md_drm, semd != 0), aes(dose, md, group = id)) +
geom_point(aes(size = 1/semd^2)) +
geom_line(linetype = "dotted") +
geom_errorbar(aes(ymin = md_low, ymax = md_upp), width = .25) +
geom_point(data = filter(md_drm, semd == 0), col = "red") +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_x_continuous("Dose", breaks = seq(0, 20, 5), limits = c(0, 20)) +
scale_y_continuous("Mean Difference", breaks = seq(-15, 25, 5)) +
guides(size = "none")
Mixed-effects model with a linear function
lin_dm <- dosresmeta(meany ~ dose, id = id, sd = sd, covariance = "md",
n = n, data = md_drm, proc = "1stage", method = "ml")Warning in optim(par = par, fn = fn, gr = gr, Xlist = Xlist, Zlist = Zlist, : one-dimensional optimization by Nelder-Mead is unreliable:
use "Brent" or optimize() directlysummary(lin_dm)Call: dosresmeta(formula = meany ~ dose, id = id, n = n, sd = sd, data = md_drm,
covariance = "md", method = "ml", proc = "1stage")
One-stage random-effects meta-analysis
Estimation method: ML
Covariance approximation: Mean Differences
Chi2 model: X2 = 0.7258 (df = 1), p-value = 0.3943
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
dose 0.3138 0.3684 0.8519 0.3943 -0.4082 1.0359
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev
1.1550
10 studies, 19 values, 1 fixed and 1 random-effects parameters
logLik AIC BIC
-118.0202 240.0404 241.9293 Mixed-effects model with restricted cubic splines
k_md <- c(2.086257, 4.350319, 8.50176) # to match Stata
spl_dm <- dosresmeta(meany ~ rcs(dose, k_md), id = id, sd = sd, covariance = "md",
n = n, data = md_drm, proc = "1stage", method = "ml")
summary(spl_dm)Call: dosresmeta(formula = meany ~ rcs(dose, k_md), id = id, n = n,
sd = sd, data = md_drm, covariance = "md", method = "ml",
proc = "1stage")
One-stage random-effects meta-analysis
Estimation method: ML
Covariance approximation: Mean Differences
Chi2 model: X2 = 86.2043 (df = 2), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb
rcs(dose, k_md)dose -1.2623 0.1895 -6.6605 0.0000 -1.6338
rcs(dose, k_md)dose' 2.7680 0.5285 5.2379 0.0000 1.7322
95%ci.ub
rcs(dose, k_md)dose -0.8909 ***
rcs(dose, k_md)dose' 3.8038 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
rcs(dose, k_md)dose 0.3656 rcs(dose, k_md)dose
rcs(dose, k_md)dose' 1.4623 1
10 studies, 19 values, 2 fixed and 3 random-effects parameters
logLik AIC BIC
-48.6696 107.3392 112.0614 A plot of the estimated summary dose-response relationship
doseref <- 5
data.frame(dose = c(doseref, seq(2, 10, length.out = 101))) %>%
bind_cols(predict(spl_dm, newdata = .)) %>%
mutate(true_pred = -2*(dose - doseref) + .2*(dose^2 - doseref^2)) %>%
ggplot(aes(dose, pred)) +
geom_line() +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = .1) +
geom_line(aes(y = true_pred), col = "blue") +
scale_y_continuous(breaks = seq(0, 12, 2)) +
labs(x = "Dose", y = "Mean difference")
Odds Ratios
or_drm <- haven::read_dta(url("http://www.stats4life.se/data/or_drm.dta")) %>%
mutate(
b_low = b - qnorm(.975)*seb,
b_upp = b + qnorm(.975)*seb
)
head(or_drm)# A tibble: 6 × 9
id bmi b seb case n type b_low b_upp
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 22.0 -0.389 0.0947 362 1000 1 -0.575 -0.203
2 1 26.9 0 0 480 1000 1 0 0
3 2 22.1 -0.483 0.194 74 250 1 -0.864 -0.102
4 2 26.9 0 0 107 250 1 0 0
5 3 22.4 -0.928 0.137 148 500 1 -1.20 -0.658
6 3 27.1 0 0 270 500 1 0 0 A plot of the empirical estimates
ggplot(data = filter(or_drm, seb != 0), aes(bmi, exp(b), group = id)) +
geom_point(aes(size = 1/(seb^2))) +
geom_line(linetype = "dotted") +
geom_errorbar(aes(ymin = exp(b_low), ymax = exp(b_upp)), width = .25) +
geom_point(data = filter(or_drm, seb == 0), col = "red") +
geom_hline(yintercept = 1, linetype = "dashed") +
scale_x_continuous(expression(paste("Body Mass Index (kg/", m^2, ")")),
breaks = seq(15, 35, 5), limits = c(15, 35)) +
scale_y_continuous(trans = "log", "Adjusted Odds Ratio", breaks = c(.5, 1, 2, 4)) +
guides(size = "none")
Mixed-effects model with restricted cubic splines
k_or <- c(22.03318, 24.51042, 27.0792) # to match Stata
spl_or <- dosresmeta(b ~ rcs(bmi, k_or), id = id, se = seb, type = "cc",
n = n, cases = case, data = or_drm, proc = "1stage", method = "ml")
summary(spl_or)Call: dosresmeta(formula = b ~ rcs(bmi, k_or), id = id, type = "cc",
cases = case, n = n, data = or_drm, se = seb, method = "ml",
proc = "1stage")
One-stage random-effects meta-analysis
Estimation method: ML
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 236.9607 (df = 2), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
rcs(bmi, k_or)bmi -0.1191 0.0343 -3.4717 0.0005 -0.1863 -0.0518
rcs(bmi, k_or)bmi' 0.3413 0.0442 7.7137 0.0000 0.2546 0.4280
rcs(bmi, k_or)bmi ***
rcs(bmi, k_or)bmi' ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
rcs(bmi, k_or)bmi 0.0242 rcs(bmi, k_or)bmi
rcs(bmi, k_or)bmi' 0.0089 -1
10 studies, 13 values, 2 fixed and 3 random-effects parameters
logLik AIC BIC
6.8770 -3.7541 -0.9293 A plot of the estimated summary dose-response relationship
bmiref <- 24
data.frame(bmi = c(bmiref, seq(21, 28, length.out = 101))) %>%
bind_cols(predict(spl_or, newdata = ., expo = TRUE)) %>%
mutate(true_pred = exp(-2.3*(bmi - bmiref) + 0.05*(bmi^2 - bmiref^2))) %>%
ggplot(aes(bmi, pred)) +
geom_line() +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = .1) +
geom_line(aes(y = true_pred), col = "blue") +
scale_y_continuous(trans = "log", breaks = c(1, 1.2, 1.5, 2, 3, 4)) +
labs(x = expression(paste("Body Mass Index (kg/", m^2, ")")), y = "Adjusted Odds Ratio")
Hazard Ratios
hr_drm <- haven::read_dta(url("http://www.stats4life.se/data/hr_drm.dta")) %>%
mutate(
b_low = b - qnorm(.975)*seb,
b_upp = b + qnorm(.975)*seb
)
head(hr_drm)# A tibble: 6 × 9
id walk b seb case n type b_low b_upp
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.279 1.13 0.110 229 777. 2 0.911 1.34
2 1 2.40 0 0 137 1704. 2 0 0
3 2 0.0882 -0.316 0.0946 245 643. 2 -0.501 -0.130
4 2 0.488 -0.448 0.0941 230 782. 2 -0.632 -0.263
5 2 1.34 -0.493 0.0933 225 920. 2 -0.676 -0.310
6 2 3.77 0 0 236 625. 2 0 0 A plot of the empirical estimates
ggplot(data = filter(hr_drm, seb != 0), aes(walk, exp(b), group = id)) +
geom_point(aes(size = 1/(seb^4))) +
geom_line(linetype = "dotted") +
geom_errorbar(aes(ymin = exp(b_low), ymax = exp(b_upp)), width = .25) +
geom_point(data = filter(hr_drm, seb == 0), col = "red") +
geom_hline(yintercept = 1, linetype = "dashed") +
scale_x_continuous("Brisk walking (hours/week)", breaks = 0:6, limits = c(0, 6)) +
scale_y_continuous(trans = "log", "Adjusted Hazard Ratio", breaks = c(.125, .25, .5, 1, 2, 4, 8)) +
guides(size = "none")
Mixed-effects model with piecewise linear splines
linspl_hr <- dosresmeta(b ~ walk + I((walk > 2)*(walk - 2)), id = id, se = seb, type = "ir",
n = n, cases = case, data = hr_drm, proc = "1stage", method = "ml")
summary(linspl_hr)Call: dosresmeta(formula = b ~ walk + I((walk > 2) * (walk - 2)), id = id,
type = "ir", cases = case, n = n, data = hr_drm, se = seb,
method = "ml", proc = "1stage")
One-stage random-effects meta-analysis
Estimation method: ML
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 110.2668 (df = 2), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb
walk -0.4679 0.0537 -8.7168 0.0000 -0.5731
I((walk > 2) * (walk - 2)) 0.5433 0.0626 8.6741 0.0000 0.4205
95%ci.ub
walk -0.3627 ***
I((walk > 2) * (walk - 2)) 0.6660 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
walk 0.2769 walk
I((walk > 2) * (walk - 2)) 0.2253 -0.2193
30 studies, 61 values, 2 fixed and 3 random-effects parameters
logLik AIC BIC
-13.7733 37.5466 48.1010 Contrasts of interest
# Linear trend per 30 min/wk between 0 and 2 h/wk
Epi::ci.exp(linspl_hr, ctr.mat = cbind(.5, 0)) exp(Est.) 2.5% 97.5%
[1,] 0.7914144 0.7508619 0.834157# Linear trend per 30 min/wk between 2 and 4 h/wk
Epi::ci.exp(linspl_hr, ctr.mat = cbind(.5, .5)) exp(Est.) 2.5% 97.5%
[1,] 1.038426 0.9737464 1.107401A plot of the estimated summary dose-response relationship
walkref <- 2
data.frame(walk = c(walkref, seq(0, 4, length.out = 101))) %>%
bind_cols(select(predict(linspl_hr, newdata = ., expo = TRUE), -walk)) %>%
mutate(true_pred = exp(-.5*(walk - walkref) + .5*(walk > walkref)*(walk - walkref))) %>%
ggplot(aes(walk, pred)) +
geom_line() +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = .1) +
geom_line(aes(y = true_pred), col = "blue") +
scale_y_continuous(trans = "log", breaks = c(1, 1.5, 2, 3)) +
labs(x = "Brisk walking (hours/week)", y = "Adjusted Hazard Ratio")