Two-level model
First, for simplicity, we run a two-level model ignoring the fact that certain effect sizes come from the same study.
dat2l <- datsel
dat2l[["author"]] <- NULLTwo-level model with the lasso prior
We start with running a penalized meta-analysis using the lasso prior. Compared to the horseshoe prior, the lasso is easier to use because it has only two hyperparameters to set. However, the lighter tails of the lasso can result in large coefficients being shrunken too much towards zero thereby leading to potentially more bias compared to the regularized horseshoe prior.
For the lasso prior, we need to specify the degrees of freedom df and the scale. Both default to 1. The degrees of freedom determines the chi-square prior that is specified for the inverse-tuning parameter. Increasing the degrees of freedom will allow larger values for the inverse-tuning parameter, leading to less shrinkage. Increasing the scale parameter will also result in less shrinkage. The influence of these hyperparameters can be visualized through the implemented shiny app, which can be called via shiny_prior().
fit_lasso <- brma(yi ~ .,
data = dat2l,
vi = "vi",
method = "lasso",
prior = c(df = 1, scale = 1),
mute_stan = FALSE)Assessing convergence and interpreting the results
We can request the results using the summary function. Before we interpret the results, we need to ensure that the MCMC chains have converged to the posterior distribution. Two helpful diagnostics provided in the summary are the number of effective posterior samples n_eff and the potential scale reduction factor Rhat. n_eff is an estimate of the number of independent samples from the posterior. Ideally, the ratio n_eff to total samples is as close to 1 as possible. Rhat compares the between- and within-chain estimates and is ideally close to 1 (indicating the chains have mixed well). Should any values for n_eff or Rhat be far from these ideal values, you can try increasing the number of iterations through the iter argument. By default, the brma function runs four MCMC chains with 2000 iterations each, half of which is discarded as burn-in. As a result, a total of 4000 iterations is available on which posterior summaries are based. If this does not help, non-convergence might indicate a problem with the model specification.
sum <- summary(fit_lasso)
sum$coefficients[, c("mean", "sd", "2.5%", "97.5%", "n_eff", "Rhat")]#> mean sd 2.5% 97.5% n_eff Rhat
#> Intercept -21.0631 14.3441 -50.8793 3.1855 2202 1
#> mTimeLength -0.0033 0.0049 -0.0145 0.0050 2525 1
#> year 0.0106 0.0071 -0.0014 0.0254 2205 1
#> modelLG 0.0884 0.1442 -0.1647 0.4143 2414 1
#> modelLNB 0.1436 0.1043 -0.0259 0.3670 1828 1
#> modelM 0.0230 0.0509 -0.0693 0.1368 2548 1
#> modelMD 0.0225 0.0767 -0.1272 0.1873 2779 1
#> ageWeek -0.0073 0.0052 -0.0182 0.0014 1571 1
#> strainGroupedC57Bl6 -0.0198 0.0681 -0.1706 0.1161 2507 1
#> strainGroupedCD1 -0.1256 0.1860 -0.5367 0.1996 2470 1
#> strainGroupedDBA -0.0305 0.1634 -0.3682 0.3089 3241 1
#> strainGroupedlisterHooded -0.0556 0.3242 -0.7567 0.5836 2808 1
#> strainGroupedlongEvans 0.0752 0.1231 -0.1449 0.3406 2496 1
#> strainGroupedlongEvansHooded 0.1744 0.2155 -0.1850 0.6499 2308 1
#> strainGroupedNMRI 0.1186 0.2402 -0.3263 0.6381 2703 1
#> strainGroupedNS 0.0348 0.1482 -0.2655 0.3497 2741 1
#> strainGroupedother -0.0148 0.1220 -0.2692 0.2396 2371 1
#> strainGroupedspragueDawley -0.0049 0.0620 -0.1351 0.1233 2942 1
#> strainGroupedswissWebster 0.5754 0.4996 -0.2536 1.6622 2418 1
#> strainGroupedwistar 0.0245 0.0558 -0.0780 0.1472 2791 1
#> strainGroupedwistarKyoto 0.0347 0.2762 -0.5335 0.6131 2912 1
#> bias -0.0085 0.0274 -0.0681 0.0434 2789 1
#> speciesrat 0.1183 0.0850 -0.0166 0.2939 1973 1
#> domainsLearning 0.0032 0.0490 -0.0925 0.1069 3120 1
#> domainnsLearning 0.1158 0.0807 -0.0195 0.2859 1723 1
#> domainsocial 0.1733 0.1204 -0.0263 0.4318 1631 1
#> domainnoMeta -0.3626 0.1846 -0.7341 -0.0248 1714 1
#> sexM 0.1154 0.0661 -0.0027 0.2554 1864 1
#> tau2 0.4280 0.0388 0.3560 0.5041 1368 1If we are satisfied with the convergence, we can continue looking at the posterior summary statistics. The summary function provides the posterior mean estimate for the effect of each moderator. Since Bayesian penalization does not automatically shrink estimates exactly to zero, some additional criterion is needed to determine which moderators should be selected in the model. Currently, this is done using the 95% credible intervals, with a moderator being selected if zero is excluded in this interval. In the summary this is denoted by an asterisk for that moderator. In this model, the only significant moderator is the dummy variable for sex. All other coefficients are not significant after being shrunken towards zero by the lasso prior.
Also note the summary statistics for tau2, the (unexplained) residual between-studies heterogeneity. The 95% credible interval for this coefficient excludes zero, indicating that there is a non-zero amount of unexplained heterogeneity. It is customary to express this heterogeneity in terms of \(I^2\), the percentage of variation across studies that is due to heterogeneity rather than chance (Higgins and Thompson, 2002; Higgins et al., 2003). The helper function I2() computes the posterior distribution of \(I^2\) based on the MCMC draws of tau2:
I2(fit_lasso)#> mean sd 2.5% 25% 50% 75% 97.5%
#> I2 68 2 64 67 68 69 71Two-level model with the horseshoe prior
Next, we look into the regularized horseshoe prior. The horseshoe prior has five hyperparameters that can be set. Three parameters are degrees of freedom parameters which influence the tails of the distributions in the prior. Generally, it is not needed to specify different values for these hyperparameters. Here, we focus instead on the global scale parameter and the scale of the slab.
# use the default settings
fit_hs1 <- brma(yi ~ .,
data = dat2l,
vi = "vi",
method = "hs",
prior = c(df = 1, df_global = 1, df_slab = 4, scale_global = 1, scale_slab = 1, relevant_pars = NULL),
mute_stan = FALSE)
# reduce the global scale
fit_hs2 <- brma(yi ~ .,
data = dat2l,
vi = "vi",
method = "hs",
prior = c(df = 1, df_global = 1, df_slab = 4, scale_global = 0.1, scale_slab = 1, relevant_pars = NULL),
mute_stan = FALSE)
# increase the scale of the slab
fit_hs3 <- brma(yi ~ .,
data = dat2l,
vi = "vi",
method = "hs",
prior = c(df = 1, df_global = 1, df_slab = 4, scale_global = 1, scale_slab = 5, relevant_pars = NULL),
mute_stan = FALSE)Note that the horseshoe prior results in some divergent transitions. This can be an indication of non-convergence. However, these divergences arise often when using the horseshoe prior and as long as there are not too many of them, the results can still be used.
Next, we plot the posterior mean estimates for a selection of moderators for the different priors.
make_plotdat <- function(fit, prior){
plotdat <- data.frame(fit$coefficients)
plotdat$par <- rownames(plotdat)
plotdat$Prior <- prior
return(plotdat)
}
df0 <- make_plotdat(fit_lasso, prior = "lasso")
df1 <- make_plotdat(fit_hs1, prior = "hs default")
df2 <- make_plotdat(fit_hs2, prior = "hs reduced global scale")
df3 <- make_plotdat(fit_hs3, prior = "hs increased slab scale")
df <- rbind.data.frame(df0, df1, df2, df3)
df <- df[!df$par %in% c("Intercept", "tau2"), ]
pd <- 0.5
ggplot(df, aes(x=mean, y=par, group = Prior)) +
geom_errorbar(aes(xmin=X2.5., xmax=X97.5., colour = Prior), width=.1, position = position_dodge(width = pd)) +
geom_point(aes(colour = Prior), position = position_dodge(width = pd)) +
geom_vline(xintercept = 0) +
theme_bw() + xlab("Posterior mean") + ylab("")
We can see that, in general, the different priors give quite similar results in this application. A notable exception is the estimate for the dummy variable testAuthorGrouped_stepDownAvoidance which is much smaller for the lasso compared to the horseshoe specification. This indicates that the lasso can shrink large coefficients more towards zero whereas the horseshoe is better at keeping them large.