In jepusto/metascore: Model-Based and Robust Score Tests for Meta-Regression

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Vevea & Hedges (1995) developed a selection model for random effects meta-analysis (and meta-regression), which models a potential mechanism for publication bias based on a step function for the probability that significant results are published. The model has two components: a model for the effect size measurement process and a model for the reporting/publication process. Let me start by laying out the model (with some slight tweaks to Vevea and Hedges' notation).

We begin with a set of $n^$ effect size estimates, $T^_1,...,T^_{n^}$ with known sampling variances $\sigma_h^2$, $h = 1,...,n^$. Each effect size estimate also has a set of characteristics captured by a row-vector of $p$ covariates $\mathbf{x}_h$, again for $h = 1,...,n^$. The effect size measurement process is simply a random effects meta-regression, in which: $$ T_h^* \sim N( \mathbf{x}_h \boldsymbol\beta, \ \tau^2 + \sigma_h^2), $$ where $\boldsymbol\beta$ is a $p \times 1$ vector of regression coefficients and $\tau^2$ is between-study variance in the effect size estimands.

The second component is a model for the process by which these effect size estimates are reported or published. We assume that this process is driven by the direction and statistical significance of the estimates, as summarized in the one-sided p-values corresponding to each estimate. The one-sided p-value corresponding to study $h$ is $p_h = 1 - \Phi\left(T_h^* / \sigma_h\right)$, where $\Phi$ is the standard normal cumulative distribution function. The selection model takes the form of a step function in the one-side p-values, with changes in level at a defined set of $q \geq 1$ steps, $a_1,...,a_q$, where $0 < a_s < 1$ for $s = 1,...,q$. Let $\boldsymbol\omega = (\omega_1,...,\omega_q)$ be a $q \times 1$ vector of weights, where $\omega_s \geq 0$ for $s = 1,...,q$. Define

$$ w\left(p_h\right) = \begin{cases} 1 & \text{if} \quad 0 < p_h \leq a_1 \ \omega_j & \text{if} \quad a_j < p_h \leq a_{j+1} \ \omega_q & \text{if} \quad a_q < p_h \leq 1 \end{cases} $$ It is assumed that the effect size $T^*_h$ is observed with probability $w(p_h)$ and censored with probability $1 - w(p_h)$. Note that the step function can be written in terms of the effect size estimates and sampling variances

$$ w\left(T, \ \sigma^2\right) = \begin{cases} 1 & \text{if} \quad -\sigma \Phi^{-1}(a_1) < T \ \omega_s & \text{if} \quad -\sigma \Phi^{-1}(a_{s + 1}) < T \leq -\sigma \Phi^{-1}(a_s) \ \omega_q & \text{if} \quad T \leq -\sigma \Phi^{-1}(a_q) \end{cases} $$ For notational convenience, let $\omega_0 = 1$, $a_0 = 0$, and $a_{q + 1} = 1$.

Let $T_1,...,T_n$ correspond to the observed effect sizes, with corresponding sampling variances $\sigma_1^2,...,\sigma_n^2$. Let $\phi()$ denote the standard normal density. The distribution of the observed effect sizes is a distortion of the normal distribution, given by

$$ f\left(\left. T_i \right| \boldsymbol\beta, \tau^2, \boldsymbol\omega\right) = \frac{w(T_i, \sigma_i^2) \phi\left(\frac{T_i - \mathbf{x}_i \boldsymbol\beta}{\sqrt{\tau^2 + \sigma_i^2}}\right)}{\sqrt{\tau^2 + \sigma_i^2} A(\mathbf{x}_i \boldsymbol\beta, \tau^2, \boldsymbol\omega)}. $$

$A(\mathbf{x}_i \boldsymbol\beta, \tau^2, \boldsymbol\omega)$ is a normalizing constant given by

$$ A(\mu_i, \tau^2, \boldsymbol\omega) = \sum_{s = 0}^q \omega_s B_{is} \left(\mu_i, \tau^2\right), $$

where $\mu_i = \mathbf{x}_i \boldsymbol\beta$ and

$$ B_{is}\left(\mu_i, \tau^2\right) = \Phi\left(\frac{-\sigma_i\Phi^{-1}(a_s) - \mu_i}{\sqrt{\tau^2 + \sigma_i^2}}\right) - \Phi\left(\frac{-\sigma_i\Phi^{-1}(a_{s+1}) - \mu_i}{\sqrt{\tau^2 + \sigma_i^2}}\right). $$

In the simple case where the effect size estimates are mutually independent, we then have joint likelihood

$$ \mathcal{L}\left(\boldsymbol\beta, \tau^2, \boldsymbol\omega\right) = \prod_{i=1}^n \frac{w(T_i, \sigma_i^2) \phi\left(\frac{T_i - \mathbf{x}_i \boldsymbol\beta}{\sqrt{\tau^2 + \sigma_i^2}}\right)}{\sqrt{\tau^2 + \sigma_i^2} A(\mathbf{x}_i \boldsymbol\beta, \tau^2, \boldsymbol\omega)}, $$

and log-likelihood proportional to

$$ l\left(\boldsymbol\beta, \tau^2, \boldsymbol\omega\right) = \sum_{i=1}^n \ln w(T_i, \sigma_i^2) - \frac{1}{2} \sum_{i=1}^n \frac{\left(T_i - \mathbf{x}i \boldsymbol\beta\right)^2}{\tau^2 + \sigma_i^2} - \frac{1}{2}\sum{i=1}^n \ln\left(\tau^2 + \sigma_i^2\right) - \sum_{i=1}^n \ln A(\mathbf{x}_i \boldsymbol\beta, \tau^2, \boldsymbol\omega). $$

Score function

Vevea and Hedges (1995) provide expressions for the score function, which is the derivative of the log likelihood with respect to the parameter vector. Let $n_s = \sum_{i=1}^n I(a_s < p_i < a_{s+1})$ denote the number of observed p-values falling within each level of the step function for $s = 1,...,q$.

To calculate the scores, we will need the derivatives of $A_i = A(\mu_i, \tau^2, \boldsymbol\omega)$ with respect to $\mu_i$ and $\tau^2$. Let $e_i = T_i - \mathbf{x}i \boldsymbol\beta$, $\eta_i^2 = \tau^2 + \sigma_i^2$, and $c{is} = \left(-\sigma_i \Phi^{-1}(a_s) - \mathbf{x}_i \boldsymbol\beta\right) / \eta_i$. We then have

$$ \begin{aligned} \frac{\partial A_i}{\partial \mu_i} &= \frac{1}{\eta_i}\sum_{s=0}^q \omega_s \left[\phi\left(c_{i(s+1)}\right) - \phi\left(c_{is}\right) \right] \ \frac{\partial A_i}{\partial \tau^2} &= \frac{1}{2 \eta_i^2} \sum_{s=0}^q \omega_s \left[ c_{i(s+1)} \phi\left(c_{i(s+1)}\right) - c_{is}\phi\left(c_{is}\right) \right] \end{aligned} $$

The scores are then as follows:

$$ \begin{aligned} \frac{\partial l}{\partial \boldsymbol\beta} &= \sum_{i=1}^n \mathbf{x}i' \left(\frac{e_i}{\eta_i^2} - \frac{\partial A_i / \partial \mu_i}{A_i}\right) \ \frac{\partial l}{\partial \tau^2} &= \frac{1}{2}\sum{i=1}^n \frac{e_i^2}{\eta_i^4} - \frac{1}{2}\sum_{i=1}^n \frac{1}{\eta_i^2} - \sum_{i=1}^n \frac{\partial A_i/ \partial \tau^2}{A_i} \ \frac{\partial l}{\partial \omega_s} &= \frac{n_s}{\omega_s} - \sum_{i=1}^n \frac{B_{is}}{A_i} \end{aligned} $$

If $\boldsymbol\omega = \mathbf{1}$, then $A_i = 1$ for $i = 1,...,n$ and the scores simplify to

$$ \begin{aligned} \frac{\partial l}{\partial \boldsymbol\beta} &= \sum_{i=1}^n \frac{\mathbf{x}i' e_i}{\eta_i^2} \ \frac{\partial l}{\partial \tau^2} &= \frac{1}{2}\sum{i=1}^n \frac{e_i^2}{\eta_i^4} - \frac{1}{2}\sum_{i=1}^n \frac{1}{\eta_i^2} \ \frac{\partial l}{\partial \omega_s} &= n_s - \sum_{i=1}^n B_{is} \end{aligned} $$

Hessian

The score test makes use of the Hessian matrix, which is the matrix of partial second derivatives of the log likelihood with respect to the parameters. To calculate the Hessian, we will need the partial second derivatives of $A_i$, which are given by:

$$ \begin{aligned} \frac{\partial^2 A_i}{\partial \mu_i \partial \mu_i} &= \frac{1}{\eta_i^2} \sum_{s=0}^q \omega_s \left[c_{i(s+1)} \phi\left(c_{i(s+1)}\right) - c_{is} \phi\left(c_{is}\right)\right] \ \frac{\partial^2 A_i}{\partial (\tau^2)^2} &= \frac{1}{4 \eta_i^4} \sum_{s=0}^q \omega_s \left[c_{i(s+1)}\left(c_{i(s+1)}^2 - 3\right) \phi\left(c_{i(s+1)}\right) - c_{is} \left(c_{is}^2 - 3\right) \phi\left(c_{is}\right)\right] \ \frac{\partial^2 A_i}{\partial \mu_i \partial \tau^2} &= \frac{1}{2\eta_i^3} \sum_{s=0}^q \omega_s \left[\left(c_{i(s+1)}^2 - 1\right) \phi\left(c_{i(s+1)}\right) - \left(c_{is}^2 - 1\right) \phi\left(c_{is}\right)\right] \ \end{aligned} $$

The entries of the Hessian matrix are then as follows:

$$ \begin{aligned} \frac{\partial^2 l}{\partial \boldsymbol\beta \partial \boldsymbol\beta'} &= \sum_{i=1}^n \mathbf{x}i \mathbf{x}_i' \left[\left(\frac{\partial A_i / \partial \mu_i}{A_i}\right)^2 - \frac{1}{A_i} \frac{\partial^2 A_i}{\partial \mu_i \partial \mu_i} - \frac{1}{\eta_i^2}\right] \ \frac{\partial^2 l}{\partial (\tau^2)^2} &= \frac{1}{2} \sum{i=1}^n \frac{1}{\eta_i^4} - \sum_{i=1}^n \frac{e_i^2}{\eta_i^6} - \sum_{i=1}^n \frac{1}{A_i} \frac{\partial^2 A_i}{\partial (\tau^2)^2} + \sum_{i=1}^n \left(\frac{\partial A_i / \partial \tau^2}{A_i}\right)^2 \ \frac{\partial^2 l}{\partial \boldsymbol\beta \partial \tau^2} &= \sum_{i=1}^n \mathbf{x}i' \left[\frac{1}{A_i^2} \frac{\partial A_i}{\partial \mu_i} \frac{\partial A_i}{\partial \tau^2} - \frac{1}{A_i} \frac{\partial^2 A_i}{\partial \mu_i \partial \tau^2} - \frac{e_i}{\eta_i^4}\right] \ \frac{\partial^2 l}{\partial \omega_s \partial \omega_t} &= \sum{i=1}^n \frac{B_{is} B_{it}}{A_i^2} - I(s = t) \frac{n_s}{\omega_s^2} \ \frac{\partial^2 l}{\partial \boldsymbol\beta \partial \omega_s} &= \sum_{i=1}^n \frac{\mathbf{x}i'}{A_i} \left[ \frac{B{is}}{A_i} \frac{\partial A_i}{\partial \mu_i} - \frac{\partial B_{is}}{\partial \mu_i}\right] \ \frac{\partial^2 l}{\partial \tau^2 \partial \omega_s} &= \sum_{i=1}^n \frac{1}{A_i} \left[ \frac{B_{is}}{A_i} \frac{\partial A_i}{\partial \tau^2} - \frac{\partial B_{is}}{\partial \tau^2}\right] \end{aligned} $$ If $\boldsymbol\omega = \mathbf{1}$, then the Hessian simplifies to

$$ \begin{aligned} \frac{\partial^2 l}{\partial \boldsymbol\beta \partial \boldsymbol\beta'} &= - \sum_{i=1}^n \frac{\mathbf{x}i \mathbf{x}_i'}{\eta_i^2} \ \frac{\partial^2 l}{\partial (\tau^2)^2} &= \frac{1}{2} \sum{i=1}^n \frac{1}{\eta_i^4} - \sum_{i=1}^n \frac{e_i^2}{\eta_i^6} \ \frac{\partial^2 l}{\partial \boldsymbol\beta \partial \tau^2} &= - \sum_{i=1}^n \frac{\mathbf{x}i'e_i}{\eta_i^4} \ \frac{\partial^2 l}{\partial \omega_s \partial \omega_t} &= \sum{i=1}^n B_{is} B_{it} - I(s = t) n_s \ \frac{\partial^2 l}{\partial \boldsymbol\beta \partial \omega_s} &= - \sum_{i=1}^n \mathbf{x}i' \frac{\partial B{is}}{\partial \mu_i} \ \frac{\partial^2 l}{\partial \tau^2 \partial \omega_s} &= - \sum_{i=1}^n \frac{\partial B_{is}}{\partial \tau^2} \end{aligned} $$

Fisher Information

The Fisher Information matrix is given by the expectation if the negative Hessian matrix: $\mathcal{I}_f = \mathbb{E}\left(-\mathcal{H}\right)$. Under the null hypothesis, the entries of $\mathcal{I}_f$ are as follows:

$$ \begin{aligned} \mathcal{I}{\boldsymbol\beta\boldsymbol\beta} &= \sum{i=1}^n \frac{\mathbf{x}i \mathbf{x}_i'}{\eta_i^2} \ \mathcal{I}{\tau^2 \tau^2} &= \frac{1}{2} \sum_{i=1}^n \frac{1}{\eta_i^4} \ \mathcal{I}{\boldsymbol\beta \tau^2} &= \mathbf{0} \ \mathcal{I}{\omega_s \omega_t} &= I(s = t) \sum_{i=1}^n B_{is} - \sum_{i=1}^n B_{is} B_{it} \ \mathcal{I}{\boldsymbol\beta \omega_s} &= \sum{i=1}^n \mathbf{x}i' \frac{\partial B{is}}{\partial \mu_i} \ \mathcal{I}{\tau^2 \omega_s} &= \sum{i=1}^n \frac{\partial B_{is}}{\partial \tau^2} \end{aligned} $$

Score tests

Let $\boldsymbol{\hat\beta}_0$ and $\hat\tau^2_0$ denote the maximum likelihood estimates of $\boldsymbol\beta$ and $\tau^2$ under the null hypothesis restriction that $\boldsymbol\omega = \mathbf{1}$. Let $\mathbf{\tilde{S}} = \mathbf{S}(\boldsymbol{\hat\beta}_0, \hat\tau^2_0, \mathbf{1})$ denote the score vector and let $\mathbf{\tilde{\mathcal{I}}}_f = \mathbf{\mathcal{I}}_f(\boldsymbol{\hat\beta}_0, \hat\tau^2_0, \mathbf{1})$ denote the Fisher information matrix, both evaluated at the maximum likelihood estimates under the null hypothesis.

Parametric score test

The parametric score test (Rao, 1948; Boos, 1992) uses the statistic

$$ Q^p = \mathbf{\tilde{S}}' \mathbf{\tilde{\mathcal{I}}}f^{-1} \mathbf{\tilde{S}} = \mathbf{\tilde{S}}{\boldsymbol\omega}'\left(\mathbf{\tilde{\mathcal{I}}}{\boldsymbol\omega \boldsymbol\omega} - \mathbf{\tilde{\mathcal{I}}}{(\boldsymbol\beta, \tau^2) \boldsymbol\omega} \mathbf{\tilde{\mathcal{I}}}{(\boldsymbol\beta, \tau^2)(\boldsymbol\beta, \tau^2)}^{-1} \mathbf{\tilde{\mathcal{I}}}{ \boldsymbol\omega (\boldsymbol\beta, \tau^2)}\right)^{-1} \mathbf{\tilde{S}}_{\boldsymbol\omega}. $$

The null hypothesis is rejected at level $\alpha$ if $Q^p$ exceeds the $\alpha$ critical value from a $\chi^2$ distribution with $q$ degrees of freedom.

Robust score test

Boos (1992) described a "generalized" score test, which can be applied in situations where the model is defined in terms of estimating equations rather than a full likelihood. Let

$$ \mathbf{B} = \left[-\mathbf{\tilde{\mathcal{I}}}{ \boldsymbol\omega (\boldsymbol\beta, \tau^2)} \mathbf{\tilde{\mathcal{I}}}{(\boldsymbol\beta, \tau^2)(\boldsymbol\beta, \tau^2)}^{-1}, \ \ \mathbf{I}_q \right]. $$

Then the asymptotic covariance of $\mathbf{S}_{\boldsymbol\omega}$ is

$$ \mathbf{V}{S{\boldsymbol\omega}} = \mathbf{B} \left(\sum_{i=1}^n \mathbf{S}_i \mathbf{S}_i'\right) \mathbf{B}', $$

where $\mathbf{S}_i$ is the vector of scores for observation $i$ for $i = 1,...,n$. The generalized or robust score test statistic is

$$ Q^r = \mathbf{S}{\boldsymbol\omega}' \mathbf{V}{S_{\boldsymbol\omega}}^{-1} \mathbf{S}_{\boldsymbol\omega}, $$

which has an (asymptotic) $\chi^2_q$ distribution under the null hypothesis that $\boldsymbol\omega = \mathbf{1}$.

jepusto/metascore documentation built on June 25, 2019, 8:15 p.m.