---
title: Notes on regression adjustment for selectively reported results in meta-analysis
date: '2026-01-17'
categories:
- effect-size
- distribution-theory
- selective-reporting
code-fold: true
code-tools: true
toc: true
bibliography: "../selection-references.bib"
csl: "../apa.csl"
link-citations: true
---
Several of the most prominent methods for investigating selective reporting bias in meta-analyis take the form of simple linear regressions, estimated by weighted least squares.
Methods of this form include the Precision Effect Test (PET) and Precision Effect Estimator with Standard Error (PEESE) adjustments proposed by @stanley2008metaregression and @stanley2014metaregression as well as the endogenous kink (EK) meta-regression proposed by @bom2019kinked.
These methods are widely used (PET and PEESE especially) and have been studied pretty extensively using Monte Carlo simulations [e.g., @bramley2021examining; @carter2019correcting; @hong2021using; @mcshane2016adjusting; @moreno2009assessment; @stanley2017limitations].
By my read, the simulation evidence is mixed: these regression adjustments seem to be less biased than regular random effects models but they are pretty crude corrections for the bias created by selective reporting, and they don't always work very well.
In this post, I want to look at how the methods work _analytically_. I'll derive expressions for their bias and accuracy and consider how one might use these expressions to say something about the conditions under which regression adjustment might be expected to perform reasonably well.
## Regression adjustments for publication bias
For a meta-analysis of $K$ studies, let $T_i$ be an effect size estimates and $\sigma_i$ be its standard error for $i = 1,...,K$.[^fixed-sigma]
Many bias adjustment methods for this context involve the regression
$$
T_i = \alpha + \beta \ x_i + \epsilon_i,
$$ {#eq-linear-regression}
where $x_i$ is a predictor that is related to the precision of the effect size estimate and $\epsilon_i$ is a residual error. @eq-linear-regression is estimated by some form of weighted least squares with weights $w_i$, which might be the inverses of the sampling variances, $w_i = \sigma_i^{-2}$ (i.e., fixed effects weights) or the inverses of the total residual variance $w_i = \left(\hat\tau^2 + \sigma_i^2\right)^{-1}$, where $\hat\tau^2$ is an estimate of between-study variance in the effect size parameters.
[^fixed-sigma]: I will treat the standard errors as known quantities that are independent of the effect size estimates, even though the standard errors are actually stochastic and, for some effect size metrics, correlated with the effect size estimates.
The intuition behind this sort of estimator is that selective reporting on the basis of statistical significance tends to produce larger biases for smaller primary studies (which have larger standard errors), leading to what looks like "small-study" effects.
PET, PEESE, and EK adjust for the small-study effects using some form of linear regression, treating the estimate of the intercept $\alpha$ in @eq-linear-regression as the average effect size in a population of infinitely large studies (so that their standard errors are $\sigma_i = 0$).[^small-study-effects]
The main difference between the methods is in the predictor used to extrapolate to this hypothetical population:
- PET uses $x_i = \sigma_i$,
- PEESE uses $x_i = \sigma_i^2$, and
- EK uses $x_i = (\sigma_i - c) \times I(\sigma_i > c)$, where $c$ is the so-called "kink" of the meta-regression.[^kinks-with-EK]
[^small-study-effects]: All of the adjustment methods will run into trouble if small study effects are the result of some mechanism other than selective reporting. For instance, if some primary study feature is systematically, positively associated with effect size and tends to be less common in larger studies, then it will act as a confounder, inducing correlation between effect size and standard error even in the absence of selective reporting.
[^kinks-with-EK]: With the EK estimator, the kink is an estimated quantity, defined by a preliminary estimator of $\alpha$ and $\tau^2$. This makes is quite tricky to analyze, but I'm going to gloss over that detail entirely (at least for the remainder of this post).
For now, let's work with the generic predictors $x_1,...,x_K$.
Basic results from linear regression tell us that the estimator of the intercept of a linear regression fitted by weighted least squares has the form
$$
\hat\alpha = \sum_{i=1}^K \left[ \frac{w_i}{W} - \frac{\bar{x} w_i (x_i - \bar{x})}{SS_x} \right] T_i,
$$ {#eq-alpha-hat}
where $W = \sum_{i=1}^K w_i$, $\bar{x} = \frac{1}{W} \sum_{i=1}^K w_i x_i$, and $SS_x = \sum_{i=1}^K w_i (x_i - \bar{x})^2$.
Denoting the full vector of standard errors as $\boldsymbol\sigma = \left(\sigma_1,...,\sigma_K\right)'$ and observing that the predictors are all functions of the standard errors, the conditional expectation and conditional variance of $\hat\alpha$ can be written as
$$
\mathbb{E}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) = \sum_{i=1}^K \left[\frac{w_i}{W} - \frac{\bar{x} w_i (x_i - \bar{x})}{SS_x} \right] \times \mathbb{E}\left(\left. T_i \right| \sigma_i \right)
$$ {#eq-alpha-expectation}
and
$$
\mathbb{V}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) = \sum_{i=1}^K \left[\frac{w_i}{W} - \frac{\bar{x} w_i (x_i - \bar{x})}{SS_x} \right]^2 \times \mathbb{V}\left(\left. T_i \right| \sigma_i \right).
$$ {#eq-alpha-variance}
If the real regression is linear in $x_i$, then $\mathbb{E}\left(\left. T_i \right| \sigma_i \right) = \alpha + \beta \ x_i$.
Substitute this into @eq-alpha-expectation and you'll find $\mathbb{E}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) = \alpha$, so the regression estimator is indeed unbiased.
But what if the real conditional expectation is not linear?
## Data-generating processes
Regression adjustment methods are a bit hard to pin down because they're not full generative models of the process by which the observed set of effect sizes arises.
Instead, they're just different forms of approximations, which don't really commit to any particular model for how selective reporting occurs.
To say anything about how well the regression adjustments work, I will need to look beyond these specific methods to other models for the selective reporting process.
@stanley2014metaregression motivated the PEESE adjustment through a Taylor series approximation to the bias of effect size estimates under a simple form of a [step-function selection model](/posts/step-function-selection-models/), in which the real form of the bias is nonlinear in $\sigma_i$.
There are a number of other models where the strength of selection depends on the statistical significance of the effect size estimate (e.g., the [beta density model](/posts/beta-density-selection-models/)), as well as models that exhibit small-study effects produced by other selection processes (e.g., the [Copas selection model](/posts/Copas-selection-models)).
Here I'll look at a step-function model and the Copas model, mostly because these are fairly tractable to analyze.
Both models posit that effect size estimates follow a conventional random effects model prior to being put through the selection process.
Denoting an effect size estimate prior to selection as $T_i^*$, we assume
$$
T_i^* = \mu + v_i^* + e_i^*
$$ {#eq-random-effects}
where $\mu$ is the overall average effect size, $v_i^*$ is a mean-zero random effect with variance $\tau^2$, which captures heterogeneity in the effect size parameters, and $e_i^*$ is a mean zero sampling error with known standard error $\sigma_i^*$.
In the ideal, a selective reporting correction would recover an estimate of the mean of the effect size distribution $\mu$.
### Step function models
In the step function model, we assume that the generated effect size estimate $T_i^*$ is observed with some probability that depends on its $p$-value for the one-sided hypothesis test of $H_0: \delta \leq 0$.
Letting $p_i^* = 1 - \Phi\left(T_i^* / \sigma_i^*\right)$, the selective reporting process of the step-function model is then
$$
\Pr\left(T_i^* \text{ is observed}| T_i^*, \sigma_i^*\right) = \begin{cases}
1 & \text{if} \quad p_i^* < \alpha_1 \\
\lambda_1 & \text{if} \quad \alpha_h \leq p_i^* < \alpha_{h+1}, \quad h = 1,...,H-1 \\
\lambda_H & \text{if} \quad \alpha_H \leq p_i^* \\
\end{cases}
$$ {#eq-step-selection-process}
for some set of steps $\alpha_1,...,\alpha_H$ and selection weights $\boldsymbol\lambda = \left(\lambda_1,...,\lambda_H\right)$.
In my [previous post on the step-function model](/posts/step-function-selection-models/), I showed that the conditional expectation of the observed effect size estimates is
$$
\mathbb{E}(T_i | \sigma_i) = \mu + \sqrt{\tau^2 + \sigma_i^2} \times \psi(\mu, \tau, \boldsymbol\lambda, \sigma_i),
$$ {#eq-step-expectation}
where $\psi()$ is a certain, rather complicated non-linear function of the model parameters given by
$$
\psi(\mu,\tau, \boldsymbol\lambda, \sigma_i) = \frac{\sum_{h=0}^H \lambda_h \left(\phi(c_{h+1,i}) - \phi(c_{hi})\right)}{\sum_{h=0}^H \lambda_h \left(\Phi(c_{h+1,i}) - \Phi(c_{hi})\right)},
$$ {#eq-psi-bar}
where
$$
c_{hi} = \frac{\sigma_i \Phi^{-1}(1 - \alpha_h) - \mu}{\sqrt{\tau^2+ \sigma_i^2}},
$$
with $c_{0i} = \infty$ and $c_{H+1,i} = -\infty$.[^single-step]
Letting
$$
\kappa_i = \kappa(\mu,\tau, \boldsymbol\lambda, \sigma_i) = \frac{\sum_{h=0}^H \lambda_h \left[c_{hi} \phi\left(c_{hi}\right) - c_{h+1,i} \phi\left(c_{h+1,i}\right)\right]}{\sum_{h=0}^H \lambda_h \left(\Phi(c_{h+1,i}) - \Phi(c_{hi})\right)}
$$
and $\psi_i = \psi(\mu,\tau, \boldsymbol\lambda, \sigma_i)$,
the conditional variance is
$$
\mathbb{V}(T_i | \sigma_i) = \left(\tau^2 + \sigma_i^2\right) \left(1 - \kappa_i - \psi_i^2\right).
$$ {#eq-step-variance}
[^single-step]: For a model with a single step at $\alpha_1$ and a single selection parameter $\lambda_1$,
$$
\psi(\mu,\tau,\lambda_1, \sigma_i) = \frac{(1 - \lambda_1) \phi(c_{1i})}{1 - (1 - \lambda_1)\Phi(c_{1i})}.
$$
Substituting @eq-step-expectation into @eq-alpha-expectation, the bias of a regression-adjustment estimator under the step-function selection model is
$$
\mathbb{E}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) - \mu = \sum_{i=1}^K \left[\frac{w_i}{W} - \frac{\bar{x} w_i (x_i - \bar{x})}{SS_x} \right] \times \sqrt{\tau^2 + \sigma_i^2} \times \psi(\mu, \tau, \boldsymbol\lambda, \sigma_i).
$$ {#eq-alpha-expectation-step}
Thus, the bias of a regression adjustment depends on how well the predictor $x_i$ approximates the non-linear function $\sqrt{\tau^2 + \sigma_i^2} \times \psi(\mu, \tau, \boldsymbol\lambda, \sigma_i)$.
Substituting @eq-step-variance into @eq-alpha-variance, an expression for the sampling variance of the regression adjustment estimator is
$$
\mathbb{V}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) = \sum_{i=1}^K \left[\frac{w_i}{W} - \frac{\bar{x} w_i (x_i - \bar{x})}{SS_x} \right]^2 \times \left(\tau^2 + \sigma_i^2\right) \left(1 - \kappa_i - \psi_i^2\right).
$$ {#eq-alpha-variance-step}
Combining the bias and the variance in the usual way leads to an expression for RMSE:
$$
\text{RMSE}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) = \sqrt{\left(\mathbb{E}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) - \mu\right)^2 + \mathbb{V}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right)}.
$$ {#eq-alpha-RMSE}
In general, the magnitude of the bias and RMSE will depend on:
- the parameters of the true distribution of effects, $\mu$ and $\tau$,
- the strength of selective reporting, measured by $\boldsymbol\lambda$,
- the choice of predictors $x_i$ and weights $w_i$, and
- the distribution of $\sigma_i$ values (i.e., the average and spread of the standard errors).
That's a lot of moving pieces. However, the bias and RMSE are at least _computable_, and analytic expressions are much quicker to calculate than running a couple of thousand replications of a simulation.
### Copas selection models
The Copas selection model posits a different selection process than the step-function model, in which the probability that a generated effect size estimate is reported is expressed in terms of an unobserved, normally distributed variable that is correlated with the sampling error of the effect size estimate (i.e., correlated with $e_i^*$ in @eq-random-effects).
See [my previous post on the Copas model](/posts/Copas-selection-models/) for a more detailed description, interactive illustration, and some spicy takes on why this model is pretty odd.
The selection process here has three parameters, which I will denote as $a$, $b$, and $\rho$ (my previous post used $\alpha$ and $\beta$, but in this post I need these to define the regression equations). The $a$ parameter controls the overall fraction of generated effect size estimates that are observed, the $b$ parameter controls how strongly the probability of observation depends on the precision of the effect size estimate, and the $\rho$ parameter controls the degree to which the probability depends on the magnitude of the observed effect size estimate relative to the overall average effect.
The probability that a generated effect size estimate is observed is given by
$$
\Pr\left(O_i^* = 1 | T_i^*, \sigma_i^*\right) = \Phi\left(\sqrt{\frac{\tau^2 + \sigma_i^{*2}}{\tau^2 + (1 - \rho^2) \sigma_i^{*2}}}\left[\alpha + \frac{\beta}{\sigma_i^*} + \frac{\rho \sigma_i^* \left(T_i^* - \mu\right)}{\tau^2 + \sigma_i^{*2}}\right]\right).
$$ {#eq-Copas-selection-probability}
In [my previous post](/posts/Copas-selection-models/), I gave expressions for the mean and variance of observed effect size estimate conditional on $\sigma_i$.
The conditional expectation is
$$
\mathbb{E}\left(T_i | \sigma_i \right) = \mu + \rho \sigma_i \xi\left(\alpha + \frac{\beta}{\sigma_i}\right)
$$ {#eq-Copas-mean}
where $\xi(u) = \phi(u) / \Phi(u)$.
The conditional variance is
$$
\mathbb{V}\left(T_i | \sigma_i \right) = \tau^2 + \sigma_i^2\left(1 - \rho^2 \xi\left(\alpha + \frac{\beta}{\sigma_i}\right) \times \left[\alpha + \frac{\beta}{\sigma_i} + \xi\left(\alpha + \frac{\beta}{\sigma_i}\right)\right]\right).
$$ {#eq-Copas-variance}
Substituting into @eq-alpha-expectation, the bias of a regression-adjustment estimator under the Copas selection model is
$$
\mathbb{E}\left(\left.\hat\alpha\right| \boldsymbol\sigma \right) - \mu = \rho \sum_{i=1}^K \left[\frac{w_i}{W} - \frac{\bar{x} w_i (x_i - \bar{x})}{SS_x} \right] \times \sigma_i \xi\left(\alpha + \frac{\beta}{\sigma_i}\right).
$$ {#eq-alpha-expectation}
This expression is a bit simpler than the bias under the step function model. For one, the bias depends neither on $\mu$ nor on $\tau$, and is multiplicative in the $\rho$ parameter.
An implication is that, if you were running a simulation where the Copas model was the true data-generating process, the bias of the regression adjustment estimator should be invariant across $\mu$ and $\tau$ and should be strictly proportional to $\rho$.
Just as with the step-function selection model, the bias will depend on how well the predictor $x_i$ approximates the functional form of the bias curve $\sigma_i \xi\left(\alpha + \frac{\beta}{\sigma_i}\right)$, which will in turn depend on the distribution of the standard errors. In other words, the regression adjustments will perform better or worse, depending on the average $\sigma_i$ and shape of the $\sigma_i$ distribution across the included studies.
# Why care?
I've just thrown down a lot of rather ugly looking formulas, and one might well ask why anyone should care about this mess of math.
What are these formulas good for?
To be honest, I worked through this stuff mostly for myself,^[Sometimes, I get a theoretical itch that I just can't help but scratch.]
but I do think there might be some useful implications to be drawn from this.
For one, the expressions for bias and RMSE show how the performance of the regression adjustments depend on the model parameters, which is useful to know if you're trying to determine how to set up a simulation study involving PET or PEESE, choose which parameters to vary, or fit a meta-model to the results.
Another, perhaps more pragmatic application would be to use the expressions to compare the bias or accuracy of different forms of regression adjustment.
For a given application, should we use PET or PEESE? And should we weight by inverse sampling variances (fixed effects weights) or inverse total variances (random effects weights) or something else (like equal weighting)?
We might be able to get some purchase on these questions by computing bias or RMSE for some hypothetical scenarios.
Under either the step function model or the Copas model, the performance of a regression adjustment will depend on the distribution of sampling standard errors, $\sigma_1,...,\sigma_K$.
When conducting a meta-analysis, we will know this sample distribution but we won't know the true values of $\mu$, $\tau$, or the parameters of the selection process.
Still, we could evaluate the expressions above under some plausible hypothetical scenarios and see whether one regression adjustment is consistently better than the others.
This will take some further fleshing out, which I will leave for another day.