t2_lkl: Compute marginal log-likelihood of tau^2

In Replicate: Statistical Metrics for Multisite Replication Studies

Description

Given point estimates and their variances, returns the marginal restricted log-likelihood of a specified tau^2 per Veroniki AA, et al. (2016), Section 3.11. Useful as a diagnostic for p_orig per Mathur VanderWeele (under review).

Usage

t2_lkl(yi, vi, t2)

Arguments

yi	Vector of point estimates
vi	Vector of variances of point estimates Can be a vector for multiple replication studies.
t2	Heterogeneity value (tau^2) at which to compute the marginal log-likelihood

References

1. Veroniki AA et al. (2016). Methods to estimate the between-study variance and its uncertainty in meta-analysis. Research Synthesis Methods.

2. Mathur MB & VanderWeele TJ (under review). New statistical metrics for multisite replication projects.

Examples

# replication estimates (Fisher's z scale) and SEs
# from moral credential example in Mathur & VanderWeele
# (under review)
yir = c(0.303, 0.078, 0.113, -0.055, 0.056, 0.073,
        0.263, 0.056, 0.002, -0.106, 0.09, 0.024, 0.069, 0.074,
        0.107, 0.01, -0.089, -0.187, 0.265, 0.076, 0.082)

seir = c(0.111, 0.092, 0.156, 0.106, 0.105, 0.057,
         0.091, 0.089, 0.081, 0.1, 0.093, 0.086, 0.076,
         0.094, 0.065, 0.087, 0.108, 0.114, 0.073, 0.105, 0.04)

vir = seir^2

# fit meta-analysis
.m = metafor::rma.uni( yi = yir,
              vi = vir,
              knha = TRUE ) 

# vector and list of tau^2 at which to compute the log-likelihood
t2.vec = seq(0, .m$tau2*10, .001)
t2l = as.list(t2.vec)

# compute the likelihood ratio vs. the MLE for each tau^2 in t2l
temp = lapply( t2l,
               FUN = function(t2) {
                 # log-lkl itself
                 t2_lkl( yi = yir,
                         vi = vir,
                         t2 = t2 )
                 
                 # lkl ratio vs. the MLE
                 exp( t2_lkl( yi = yir,
                              vi = vir,
                              t2 = t2 ) ) / exp( t2_lkl( yi = yir,
                                                         vi = vir,
                                                         t2 = .m$tau2 ) )
               })

# plotting dataframe
dp = data.frame( tau = sqrt(t2.vec),
                 V = t2.vec,
                 lkl = unlist(temp) )

# fn: ratio of the plotted tau^2 vs. the actual MLE (for secondary x-axis)
g = function(x) x / .m$tau2

# breaks for main and secondary x-axes
breaks.x1 = seq( 0, max(dp$V), .005 )
breaks.x2 = seq( 0, max( g(dp$V) ), 1 )

p = ggplot2::ggplot( data = dp,
        ggplot2::aes(x = V, 
            y = lkl ) ) +
  
  ggplot2::geom_vline(xintercept = .m$tau2,
             lty = 2,
             color = "red") +  # the actual MLE
  
  ggplot2::geom_line(lwd = 1.2) +
  ggplot2::theme_classic() +
  ggplot2::xlab( bquote( hat(tau)["*"]^2 ) ) +
  ggplot2::ylab( "Marginal likelihood ratio of " ~ hat(tau)["*"]^2 ~ " vs. " ~ hat(tau)^2 ) +
  ggplot2::scale_x_continuous( limits = c(0, max(breaks.x1)),
                      breaks = breaks.x1,
                      sec.axis = ggplot2::sec_axis( ~ g(.),  
                                           name = bquote( hat(tau)["*"]^2 / hat(tau)^2 ),
                                           breaks=breaks.x2 ) ) +
  ggplot2::scale_y_continuous( limits = c(0,1),
                      breaks = seq(0,1,.1) )
graphics::plot(p)

Replicate documentation built on Dec. 1, 2019, 1:14 a.m.