7. Comparison of Fixed Weights
Matt Secrest and Isaac Gravestock
weighting.RmdIntroduction
In a psborrow2 analysis it is possible to specify fixed
weights for an observation’s log-likelihood contribution. This is
similar to a weighted regression or a fixed power prior parameter.
This vignette will show how weights can be specified and compare
regression model results with other packages. We will compare a
glm model with weights, a weighted likelihood in Stan with
psborrow2, and BayesPPD::glm.fixed.a0 for
generalized linear models with fixed a0 (power prior
parameter).
Logistic regression
We fit logistic regression models with the external control arm
having weights (or power parameters) equal to 0, 0.25, 0.5, 0.75, 1. The
internal treated and control patients have weight = 1. The model has a
treatment indicator and two covariates,
resp ~ trt + cov1 + cov2.
Results
logistic_res_df <- do.call(
rbind,
c(logistic_glm_reslist, logistic_ppd_reslist, logistic_psb_reslist)
)
logistic_res_df$est_ci <- sprintf(
"%.3f (%.3f, %.3f)",
logistic_res_df$estimate, logistic_res_df$lower, logistic_res_df$upper
)
wide <- reshape(
logistic_res_df[, c("fitter", "borrowing", "variable", "est_ci")],
direction = "wide",
timevar = "fitter",
idvar = c("borrowing", "variable"),
)
new_order <- order(wide$variable, wide$borrowing)
knitr::kable(wide[new_order, ], digits = 3, row.names = FALSE)| borrowing | variable | est_ci.glm | est_ci.BayesPPD | est_ci.psborrow2 |
|---|---|---|---|---|
| 0.00 | (Intercept) | 0.646 (-0.038, 1.357) | 0.661 (-0.028, 1.406) | 0.658 (-0.054, 1.379) |
| 0.25 | (Intercept) | 0.394 (-0.131, 0.931) | 0.401 (-0.123, 0.931) | 0.407 (-0.116, 0.944) |
| 0.50 | (Intercept) | 0.293 (-0.158, 0.751) | 0.302 (-0.160, 0.763) | 0.296 (-0.164, 0.745) |
| 0.75 | (Intercept) | 0.235 (-0.168, 0.642) | 0.239 (-0.179, 0.656) | 0.242 (-0.167, 0.658) |
| 1.00 | (Intercept) | 0.196 (-0.172, 0.567) | 0.209 (-0.155, 0.572) | 0.197 (-0.170, 0.556) |
| 0.00 | cov1 | -0.771 (-1.465, -0.095) | -0.790 (-1.482, -0.089) | -0.790 (-1.491, -0.075) |
| 0.25 | cov1 | -0.781 (-1.340, -0.231) | -0.793 (-1.373, -0.225) | -0.797 (-1.357, -0.246) |
| 0.50 | cov1 | -0.769 (-1.252, -0.291) | -0.786 (-1.281, -0.292) | -0.778 (-1.262, -0.303) |
| 0.75 | cov1 | -0.758 (-1.191, -0.329) | -0.769 (-1.218, -0.328) | -0.768 (-1.211, -0.340) |
| 1.00 | cov1 | -0.749 (-1.145, -0.357) | -0.769 (-1.156, -0.379) | -0.753 (-1.139, -0.357) |
| 0.00 | cov2 | -0.730 (-1.472, -0.008) | -0.755 (-1.505, 0.003) | -0.751 (-1.511, -0.025) |
| 0.25 | cov2 | -0.559 (-1.114, -0.014) | -0.571 (-1.131, -0.032) | -0.572 (-1.128, -0.029) |
| 0.50 | cov2 | -0.459 (-0.926, 0.003) | -0.464 (-0.935, -0.004) | -0.464 (-0.924, -0.003) |
| 0.75 | cov2 | -0.398 (-0.811, 0.011) | -0.403 (-0.819, 0.010) | -0.405 (-0.814, 0.001) |
| 1.00 | cov2 | -0.358 (-0.731, 0.013) | -0.359 (-0.716, 0.006) | -0.362 (-0.735, 0.016) |
| 0.00 | trt | 0.154 (-0.558, 0.871) | 0.162 (-0.611, 0.877) | 0.164 (-0.567, 0.903) |
| 0.25 | trt | 0.349 (-0.183, 0.885) | 0.357 (-0.184, 0.899) | 0.352 (-0.194, 0.888) |
| 0.50 | trt | 0.405 (-0.082, 0.894) | 0.413 (-0.081, 0.897) | 0.411 (-0.080, 0.899) |
| 0.75 | trt | 0.434 (-0.031, 0.900) | 0.442 (-0.018, 0.898) | 0.437 (-0.035, 0.905) |
| 1.00 | trt | 0.452 (-0.000, 0.905) | 0.451 (-0.001, 0.899) | 0.458 (0.000, 0.919) |
logistic_res_df$borrowing_x <- logistic_res_df$borrowing +
(as.numeric(factor(logistic_res_df$fitter)) - 3) / 100
ggplot(logistic_res_df, aes(x = borrowing_x, y = estimate, group = fitter, colour = fitter)) +
geom_errorbar(aes(ymin = lower, ymax = upper)) +
geom_point() +
facet_wrap(~variable, scales = "free")
Exponential models
Now we fit models with an exponentially distributed outcome. There is
no censoring in this data set. For glm we use
family = Gamma(link = "log") and specify fixed
dispersion = 1 to fit a exponential model. As before, the
external control arm having weights (or power parameters) equal to 0,
0.25, 0.5, 0.75, 1. The internal treated and control patients have
weight = 1. The model has a treatment indicator and two covariates,
eventtime ~ trt + cov1 + cov2.
head(sim_data_exp)
# id eventtime status trt cov1 cov2 ext censor
# 1 1 0.182694079 1 0 0 1 0 0
# 2 2 0.009720194 1 0 1 0 0 0
# 3 3 0.039758408 1 0 1 0 0 0
# 4 4 0.064351107 1 0 1 0 0 0
# 5 5 0.106819972 1 0 0 0 0 0
# 6 6 0.003613664 1 0 1 0 0 0Results
Note: Wald confidence intervals are displayed here for
glm for the exponential models.
res_df <- do.call(rbind, c(glm_reslist, ppd_reslist, psb_reslist))
res_df$est_ci <- sprintf(
"%.3f (%.3f, %.3f)",
res_df$estimate, res_df$lower, res_df$upper
)
wide <- reshape(
res_df[, c("fitter", "borrowing", "variable", "est_ci")],
direction = "wide",
timevar = "fitter",
idvar = c("borrowing", "variable"),
)
new_order <- order(wide$variable, wide$borrowing)
knitr::kable(wide[new_order, ], digits = 3, row.names = FALSE)| borrowing | variable | est_ci.glm | est_ci.BayesPPD | est_ci.psborrow2 |
|---|---|---|---|---|
| 0.00 | (Intercept) | 1.925 (1.592, 2.258) | 1.908 (1.574, 2.230) | 1.908 (1.579, 2.229) |
| 0.25 | (Intercept) | 1.455 (1.197, 1.712) | 1.443 (1.172, 1.702) | 1.446 (1.180, 1.687) |
| 0.50 | (Intercept) | 1.343 (1.121, 1.565) | 1.332 (1.105, 1.551) | 1.337 (1.112, 1.554) |
| 0.75 | (Intercept) | 1.287 (1.088, 1.485) | 1.281 (1.082, 1.477) | 1.283 (1.080, 1.481) |
| 1.00 | (Intercept) | 1.252 (1.070, 1.434) | 1.248 (1.064, 1.423) | 1.247 (1.061, 1.421) |
| 0.00 | cov1 | 0.949 (0.620, 1.279) | 0.952 (0.625, 1.291) | 0.955 (0.622, 1.292) |
| 0.25 | cov1 | 0.910 (0.641, 1.179) | 0.914 (0.640, 1.186) | 0.913 (0.646, 1.188) |
| 0.50 | cov1 | 0.931 (0.697, 1.165) | 0.938 (0.697, 1.177) | 0.931 (0.696, 1.162) |
| 0.75 | cov1 | 0.947 (0.737, 1.157) | 0.948 (0.737, 1.159) | 0.949 (0.740, 1.165) |
| 1.00 | cov1 | 0.959 (0.766, 1.151) | 0.962 (0.770, 1.162) | 0.962 (0.771, 1.154) |
| 0.00 | cov2 | 0.282 (-0.070, 0.634) | 0.274 (-0.087, 0.629) | 0.273 (-0.105, 0.635) |
| 0.25 | cov2 | 0.203 (-0.060, 0.467) | 0.202 (-0.069, 0.468) | 0.202 (-0.071, 0.463) |
| 0.50 | cov2 | 0.234 (0.011, 0.457) | 0.234 (0.010, 0.451) | 0.234 (0.013, 0.455) |
| 0.75 | cov2 | 0.255 (0.059, 0.452) | 0.257 (0.063, 0.456) | 0.253 (0.049, 0.451) |
| 1.00 | cov2 | 0.270 (0.092, 0.448) | 0.267 (0.089, 0.446) | 0.269 (0.090, 0.447) |
| 0.00 | trt | 0.951 (0.606, 1.297) | 0.958 (0.621, 1.320) | 0.958 (0.614, 1.322) |
| 0.25 | trt | 1.476 (1.218, 1.734) | 1.478 (1.212, 1.739) | 1.475 (1.223, 1.730) |
| 0.50 | trt | 1.564 (1.327, 1.800) | 1.562 (1.329, 1.803) | 1.563 (1.323, 1.791) |
| 0.75 | trt | 1.602 (1.375, 1.828) | 1.597 (1.368, 1.818) | 1.595 (1.363, 1.817) |
| 1.00 | trt | 1.623 (1.403, 1.843) | 1.619 (1.401, 1.831) | 1.618 (1.393, 1.833) |
res_df$borrowing_x <- res_df$borrowing +
(as.numeric(factor(res_df$fitter)) - 3) / 100
ggplot(res_df, aes(x = borrowing_x, y = estimate, group = fitter, colour = fitter)) +
geom_errorbar(aes(ymin = lower, ymax = upper)) +
geom_point() +
facet_wrap(~variable, scales = "free")