6  测量误差和元分析

Measurement Error and Meta-Analysis

本节译者:马桢

初次校审:李君竹

二次校审:邱怡轩(DeepSeek 辅助)

Most quantities used in statistical models arise from measurements. Most of these measurements are taken with some error. When the measurement error is small relative to the quantity being measured, its effect on a model is usually small. When measurement error is large relative to the quantity being measured, or when precise relations can be estimated being measured quantities, it is useful to introduce an explicit model of measurement error. One kind of measurement error is rounding.

统计模型中使用的多数量都是通过测量获得的,而大多数的测量都有一定的误差。当测量误差相对于被测量量较小时,其对模型的影响通常较小。当测量误差相对于被测量量较大,或者可以估计被测量量之间的精确关系时,引入显式的测量误差模型是有用的。一种测量误差就是舍入(rounding)。

Meta-analysis plays out statistically much like measurement error models, where the inferences drawn from multiple data sets are combined to do inference over all of them. Inferences for each data set are treated as providing a kind of measurement error with respect to true parameter values.

从统计学角度来看,元分析与测量误差模型非常相似,其中从多个数据集中得出的推断被结合,来整体的对它们进行推断。对每个数据集的推断都被当作是提供了一种相对于真实参数值的测量误差。

6.1 Bayesian measurement error model

贝叶斯测量误差模型

A Bayesian approach to measurement error can be formulated directly by treating the true quantities being measured as missing data (Clayton 1992; Richardson and Gilks 1993). This requires a model of how the measurements are derived from the true values.

处理测量误差的贝叶斯方法可以直接将真实测量值视为缺失数据来构建模型 (Clayton 1992; Richardson and Gilks 1993)。这需要一个关于如何从真实值中获取测量值的模型。

Regression with measurement error

含测量误差的回归模型

Before considering regression with measurement error, first consider a linear regression model where the observed data for \(N\) cases includes a predictor \(x_n\) and outcome \(y_n\). In Stan, a linear regression for \(y\) based on \(x\) with a slope and intercept is modeled as follows.

在考虑含测量误差的回归之前,首先考虑一个具有预测变量 \(x_n\) 和结果 \(y_n\) 的线性回归模型,其中观测到的案例数量为 \(N\)。在 Stan 中,一个基于 \(x\)\(y\) 进行带有斜率和截距项的线性回归建模如下。

data {
  int<lower=0> N;       // number of cases
  vector[N] x;          // predictor (covariate)
  vector[N] y;          // outcome (variate)
}
parameters {
  real alpha;           // intercept
  real beta;            // slope
  real<lower=0> sigma;  // outcome noise
}
model {
  y ~ normal(alpha + beta * x, sigma);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ cauchy(0, 5);
}

Now suppose that the true values of the predictors \(x_n\) are not known, but for each \(n\), a measurement \(x^{\textrm{meas}}_n\) of \(x_n\) is available. If the error in measurement can be modeled, the measured value \(x^{\textrm{meas}}_n\) can be modeled in terms of the true value \(x_n\) plus measurement noise. The true value \(x_n\) is treated as missing data and estimated along with other quantities in the model. A simple approach is to assume the measurement error is normal with known deviation \(\tau\). This leads to the following regression model with constant measurement error.

假设预测变量 \(x_n\) 的真实值未知,但是对于每个 \(n\),有一个 \(x_n\) 的测量值 \(x^{\textrm{meas}}_n\) 是可用的。如果测量的误差可以被建模,那么测量值 \(x^{\textrm{meas}}_n\) 可以用真实值 \(x_n\) 加上测量噪声来建模。真实值 \(x_n\) 被视为缺失数据,与模型中的其他量一起被估计。一种简单的方法是假设测量误差服从正态分布,且标准差 \(\tau\) 已知。由此得到下面含常测量误差的回归模型。

data {
  // ...
  array[N] real x_meas;   // measurement of x
  real<lower=0> tau;     // measurement noise
}
parameters {
  array[N] real x;    // unknown true value
  real mu_x;          // prior location
  real sigma_x;       // prior scale
  // ...
}
model {
  x ~ normal(mu_x, sigma_x);  // prior
  x_meas ~ normal(x, tau);    // measurement model
  y ~ normal(alpha + beta * x, sigma);
  // ...
}

The regression coefficients alpha and beta and regression noise scale sigma are the same as before, but now x is declared as a parameter rather than as data. The data are now x_meas, which is a measurement of the true x value with noise scale tau. The model then specifies that the measurement error for x_meas[n] given true value x[n] is normal with deviation tau. Furthermore, the true values x are given a hierarchical prior here.

回归系数 alphabeta 以及回归噪声尺度 sigma 和之前一样,但现在 x 被当作参数而不是数据。现在的数据是 x_meas,它是对真实 x 值的测量,噪声尺度为 tau。然后模型规定,给定真值 x[n]x_meas[n] 的测量误差是正态的,标准差为 tau。此外,真实值 x 在这里具有分层先验。

In cases where the measurement errors are not normal, richer measurement error models may be specified. The prior on the true values may also be enriched. For instance, Clayton (1992) introduces an exposure model for the unknown (but noisily measured) risk factors \(x\) in terms of known (without measurement error) risk factors \(c\). A simple model would regress \(x_n\) on the covariates \(c_n\) with noise term \(\upsilon\),

当测量误差不服从正态分布时,有更丰富的测量误差模型可以使用。对真实值的先验信息也可以得到丰富。例如, Clayton (1992) 引入了一个曝露(exposure)模型,用已知(没有测量误差)的风险因素 \(c\) 来表示未知(但受噪声干扰的)风险因素 \(x\)。一个简单的模型是利用协变量 \(c_n\) 回归 \(x_n\),并加上噪声项 \(\upsilon\)

\[ x_n \sim \textsf{normal}(\gamma^{\top}c, \upsilon). \] This can be coded in Stan just like any other regression. And, of course, other exposure models can be provided.

这可以像任何其他回归一样在 Stan 中编码。当然,Stan 也可以提供其他曝露模型。

Rounding

舍入

A common form of measurement error arises from rounding measurements. Rounding may be done in many ways, such as rounding weights to the nearest milligram, or to the nearest pound; rounding may even be done by rounding down to the nearest integer.

舍入误差是一种常见的测量误差形式。舍入可以用许多方式进行,例如将重量舍入到最接近的毫克或最接近的磅,甚至可以通过向下舍入到最接近的整数来进行。

Exercise 3.5(b) by Gelman et al. (2013) provides an example.

Gelman et al. (2013) 的练习3.5(b)提供了一个例子。

3.5. Suppose we weigh an object five times and measure weights, rounded to the nearest pound, of 10, 10, 12, 11, 9. Assume the unrounded measurements are normally distributed with a noninformative prior distribution on \(\mu\) and \(\sigma^2\).

  1. Give the correct posterior distribution for \((\mu, \sigma^2)\), treating the measurements as rounded.

3.5. 假设我们对一个物体进行了五次称重,并测得如下舍入到最近整数磅的 质量:10、10、12、11、9。假设未舍入的测量值服从正态分布,并对 \(\mu\)\(\sigma^2\) 采用无信息的先验分布。

  1. 把测量值视为取整后,给出 \((\mu, \sigma^2)\) 正确的后验分布。

Letting \(z_n\) be the unrounded measurement for \(y_n\), the problem as stated assumes

\(z_n\)\(y_n\) 未取整的测量值,按照问题陈述,假设

\[ z_n \sim \textsf{normal}(\mu, \sigma). \]

The rounding process entails that \(z_n \in (y_n - 0.5, y_n + 0.5)\)1. The probability mass function for the discrete observation \(y\) is then given by marginalizing out the unrounded measurement, producing the likelihood

舍入过程意味着 \(z_n \in (y_n - 0.5, y_n + 0.5)\)2。离散观测值 \(y\) 的概率质量函数可以通过边缘化消去未舍入的测量值得到,从而生成似然函数。

\[\begin{align*} p(y_n \mid \mu, \sigma) &= \int_{y_n - 0.5}^{y_n + 0.5} \textsf{normal}(z_n \mid \mu, \sigma) \,\textsf{d}z_n \\ &= \Phi\!\left(\frac{y_n + 0.5 - \mu}{\sigma}\right) -\Phi\!\left(\frac{y_n - 0.5 - \mu}{\sigma}\right). \end{align*}\] Gelman’s answer for this problem took the noninformative prior to be uniform in the variance \(\sigma^2\) on the log scale, but we replace it with more recently recommended half-normal prior on \(\sigma\)

Gelman 解决这个问题的方法是令无信息的先验在对数尺度上关于方差 \(\sigma^2\) 是均匀的,但是我们将其替换为最近更被推荐使用的半正态先验:

\[ \sigma \sim \textsf{normal}^+(0, 1). \] The posterior after observing \(y = (10, 10, 12, 11, 9)\) can be calculated by Bayes’s rule as

观测到 \(y = (10, 10, 12, 11, 9)\) 后的后验分布可以通过贝叶斯定理计算得出

\[\begin{align*} p(\mu, \sigma \mid y) &\propto p(\mu, \sigma) \ p(y \mid \mu, \sigma) \\ &\propto \textsf{normal}^+(\sigma \mid 0, 1)\prod_{n=1}^5 \left( \Phi\!\left(\frac{y_n + 0.5 - \mu}{\sigma}\right) -\Phi\!\left(\frac{y_n - 0.5 - \mu}{\sigma}\right) \right). \end{align*}\]

The Stan code simply follows the mathematical definition, providing an example of the direct definition of a probability function up to a proportion.

Stan 代码简单地遵循数学定义,再次提供了一个直接定义正比概率函数的例子。

data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  sigma ~ normal(0, 1);
  for (n in 1:N) {
    target += log_diff_exp(normal_lcdf(y[n] + 0.5 | mu, sigma),
                           normal_lcdf(y[n] - 0.5 | mu, sigma));
  }
}

where normal_lcdf(y[n]+0.5 | mu, sigma) is equal to log(Phi((y[n] + 0.5 - mu) / sigma)), and log_diff_exp(a, b) computes log(exp(a) - exp(b)) in numerically more stable way.

其中 normal_lcdf(y[n]+0.5 | mu, sigma) 等价于 log(Phi((y[n] + 0.5 - mu) / sigma)),而 log_diff_exp(a, b) 使用了一种数值更稳定的方法来计算 log(exp(a) - exp(b))

Alternatively, the model may be defined with latent parameters for the unrounded measurements \(z_n\). The Stan code in this case uses a distribution statement for \(z_n\) directly while respecting the constraint \(z_n \in (y_n - 0.5, y_n + 0.5)\).

或者,模型也可以使用潜参数(latent parameters)来定义未经舍入的测量值 \(z_n\)。在这种情况下,Stan 代码直接使用 \(z_n\) 的似然函数来保证 \(z_n \in (y_n - 0.5, y_n + 0.5)\) 这一约束条件。

data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
  vector<lower=y-0.5, upper=y+0.5>[N] z;
}
model {
  sigma ~ normal(0, 1);
  z ~ normal(mu, sigma);
}

This explicit model for the unrounded measurements \(z\) produces the same posterior for \(\mu\) and \(\sigma\) as the previous model that marginalizes \(z\) out. Both approaches mix well, but the latent parameter version is about twice as efficient in terms of effective samples per iteration, as well as providing a posterior for the unrounded parameters.

这个针对未舍入测量值 z 的显式模型会产生与之前将 \(z\) 边缘化掉的模型相同的 \(\mu\)\(\sigma\) 后验分布。两种方法十分相似,但就每次迭代的有效样本而言,潜参数版本的效率约为前者的两倍,并且为未舍入参数提供了后验分布。

6.2 Meta-analysis

元分析

Meta-analysis aims to pool the data from several studies, such as the application of a tutoring program in several schools or treatment using a drug in several clinical trials.

元分析的目标是整合来自多项研究的数据,例如在多所学校实施辅导计划或在多项临床试验中使用药物治疗。

The Bayesian framework is particularly convenient for meta-analysis, because each previous study can be treated as providing a noisy measurement of some underlying quantity of interest. The model then follows directly from two components, a prior on the underlying quantities of interest and a measurement-error style model for each of the studies being analyzed.

贝叶斯方法特别适合用于元分析,因为每个先前的研究都可以被看作是对某些感兴趣的潜在量的有误差的测量。然后,模型直接来源于两个组成部分,即对感兴趣的潜在量的先验和对每个被分析的研究的测量误差式模型。

Treatment effects in controlled studies

对照研究中的治疗效果

Suppose the data in question arise from a total of \(M\) studies providing paired binomial data for a treatment and control group. For instance, the data might be post-surgical pain reduction under a treatment of ibuprofen (Warn, Thompson, and Spiegelhalter 2002) or mortality after myocardial infarction under a treatment of beta blockers (Gelman et al. 2013, sec. 5.6).

假设所涉及的数据来自总共 \(M\) 项研究,这些研究提供了治疗组和对照组的配对二项数据。例如,这些数据可能是在应用布洛芬进行手术后疼痛减轻方面的数据 (Warn, Thompson, and Spiegelhalter 2002) 或在应用 beta 受体阻滞剂治疗心肌梗死后的死亡率数据 (Gelman et al. 2013, sec. 5.6)

Data

数据

The clinical data consists of \(J\) trials, each with \(n^t\) treatment cases, \(n^c\) control cases, \(r^t\) successful outcomes among those treated and \(r^c\) successful outcomes among those in the control group. This data can be declared in Stan as follows.3

临床数据由 \(J\) 个试验组成,每个试验包含 \(n^t\) 个治疗组病例,\(n^c\) 个对照组病例,治疗组中有 \(r^t\) 个成功结果,对照组中有 \(r^c\) 个成功结果。可以在 Stan 中如下声明这些数据。4

data {
  int<lower=0> J;
  array[J] int<lower=0> n_t;  // num cases, treatment
  array[J] int<lower=0> r_t;  // num successes, treatment
  array[J] int<lower=0> n_c;  // num cases, control
  array[J] int<lower=0> r_c;  // num successes, control
}

Converting to log odds and standard error

转换为对数几率和标准差

Although the clinical trial data are binomial in its raw format, it may be transformed to an unbounded scale by considering the log odds ratio

尽管临床试验数据在原始形式下是二项分布的,但可以将其转换到无穷尺度上,通过考虑对数几率(log odds ratio)

\[\begin{align*} y_j &= \log \left( \frac{r^t_j / (n^t_j - r^t_j)} {r^c_j / (n^c_j - r^c_j)} \right) \\ &= \log \left( \frac{r^t_j}{n^t_j - r^t_j} \right) -\log \left( \frac{r^c_j}{n^c_j - r^c_j} \right) \end{align*}\] and corresponding standard errors

以及对应的标准误

\[ \sigma_j = \sqrt{ \frac{1}{r^T_i} + \frac{1}{n^T_i - r^T_i} + \frac{1}{r^C_i} + \frac{1}{n^C_i - r^C_i} }. \]

The log odds and standard errors can be defined in a transformed data block, though care must be taken not to use integer division.5

对数几率和标准误可以在变换数据模块中定义,但需要注意不要使用整数除法。6

transformed data {
  array[J] real y;
  array[J] real<lower=0> sigma;
  for (j in 1:J) {
    y[j] = log(r_t[j]) - log(n_t[j] - r_t[j])
            - (log(r_c[j]) - log(n_c[j] - r_c[j]));
  }
  for (j in 1:J) {
    sigma[j] = sqrt(1 / r_t[j] + 1 / (n_t[j] - r_t[j])
                     + 1 / r_c[j] + 1 / (n_c[j] - r_c[j]));
  }
}

This definition will be problematic if any of the success counts is zero or equal to the number of trials. If that arises, a direct binomial model will be required or other transforms must be used than the unregularized sample log odds.

如果任何一次试验的成功次数为零或等于试验次数,那么这个定义就会有问题。如果出现这种情况,就需要使用直接的二项模型或应该使用其他变换来代替未经正则化的样本对数几率。

Non-hierarchical model

非分层模型

With the transformed data in hand, two standard forms of meta-analysis can be applied. The first is a so-called “fixed effects” model, which assumes a single parameter for the global odds ratio. This model is coded in Stan as follows.

有转换后的数据之后,有两种标准的元分析形式可以用。第一种是所谓的“固定效应”模型,它假设全局几率(global odds ratio)只有一个参数。该模型在 Stan 中的代码如下。

parameters {
  real theta;  // global treatment effect, log odds
}
model {
  y ~ normal(theta, sigma);
}

The distribution statement for y is vectorized; it has the same effect as the following.

y 的抽样语句是向量化的;它的效果与下面代码相同。

  for (j in 1:J) {
    y[j] ~ normal(theta, sigma[j]);
  }

It is common to include a prior for theta in this model, but it is not strictly necessary for the model to be proper because y is fixed and \(\textsf{normal}(y \mid \mu,\sigma) = \textsf{normal}(\mu \mid y,\sigma)\).

在这个模型中通常会对 theta 设定一个先验分布,但是这并不是必需的,因为 y 是固定的,因此 \(\textsf{normal}(y \mid \mu,\sigma) = \textsf{normal}(\mu \mid y,\sigma)\)

Hierarchical model

分层模型

To model so-called “random effects,” where the treatment effect may vary by clinical trial, a hierarchical model can be used. The parameters include per-trial treatment effects and the hierarchical prior parameters, which will be estimated along with other unknown quantities.

为了模拟所谓的“随机效应”,即治疗效果可能因临床试验而异,可以使用分层模型。参数包括每个试验的治疗效果和分层先验参数,这些参数将与其他未知量一起被估计。

parameters {
  array[J] real theta;  // per-trial treatment effect
  real mu;              // mean treatment effect
  real<lower=0> tau;    // deviation of treatment effects
}
model {
  y ~ normal(theta, sigma);
  theta ~ normal(mu, tau);
  mu ~ normal(0, 10);
  tau ~ cauchy(0, 5);
}

Although the vectorized distribution statement for y appears unchanged, the parameter theta is now a vector. The distribution statement for theta is also vectorized, with the hyperparameters mu and tau themselves being given wide priors compared to the scale of the data.

虽然 y 的向量化抽样语句似乎没有变化,但是参数 theta 现在是一个向量。theta 的抽样语句也是向量化的,其中超参数 mutau 相对于数据的尺度给出了“宽”的先验分布。

Rubin (1981) provided a hierarchical Bayesian meta-analysis of the treatment effect of Scholastic Aptitude Test (SAT) coaching in eight schools based on the sample treatment effect and standard error in each school.

Rubin (1981) 提供了一种基于每所学校样本处理效果和标准差的层次贝叶斯元分析,用于研究八所学校中 Scholastic Aptitude Test (SAT) 辅导的处理效果。

Extensions and alternatives

扩展和替代方案

Smith, Spiegelhalter, and Thomas (1995) and Gelman et al. (2013, sec. 19.4) provide meta-analyses based directly on binomial data. Warn, Thompson, and Spiegelhalter (2002) consider the modeling implications of using alternatives to the log-odds ratio in transforming the binomial data.

Smith, Spiegelhalter, and Thomas (1995)Gelman et al. (2013, sec. 19.4) 直接基于二项数据进行元分析。Warn, Thompson, and Spiegelhalter (2002) 考虑了在转换二项数据时使用对数几率(log-odds ratio)之外的替代方法所产生的建模影响。

If trial-specific predictors are available, these can be included directly in a regression model for the per-trial treatment effects \(\theta_j\).

Clayton, D. G. 1992. “Models for the Analysis of Cohort and Case-Control Studies with Inaccurately Measured Exposures.” In Statistical Models for Longitudinal Studies of Exposure and Health, edited by James H. Dwyer, Manning Feinleib, Peter Lippert, and Hans Hoffmeister, 301–31. New York: Oxford University Press.
Gelman, Andrew, J. B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. Third Edition. London: Chapman & Hall / CRC Press.
Richardson, Sylvia, and Walter R. Gilks. 1993. “A Bayesian Approach to Measurement Error Problems in Epidemiology Using Conditional Independence Models.” American Journal of Epidemiology 138 (6): 430–42.
Rubin, Donald B. 1981. “Estimation in Parallel Randomized Experiments.” Journal of Educational Statistics 6: 377–401.
Smith, Teresa C., David J. Spiegelhalter, and Andrew Thomas. 1995. Bayesian Approaches to Random-Effects Meta-Analysis: A Comparative Study.” Statistics in Medicine 14 (24): 2685–99.
Warn, David E., S. G. Thompson, and David J. Spiegelhalter. 2002. Bayesian Random Effects Meta-Analysis of Trials with Binary Outcomes: Methods for the Absolute Risk Difference and Relative Risk Scales.” Statistics in Medicine 21: 1601–23.

  1. There are several different rounding rules (see, e.g., Wikipedia: Rounding), which affect which interval ends are open and which are closed, but these do not matter here as for continuous \(z_n\) \(p(z_n=y_n-0.5)=p(z_n=y_n+0.5)=0\).↩︎

  2. 存在多种不同的舍入规则(参见 Wikipedia: Rounding),这些规则会影响区间端点的开闭性,但由于 \(z_n\) 是连续的,\(p(z_n=y_n-0.5)=p(z_n=y_n+0.5)=0\),因此这些规则在此并不重要。↩︎

  3. Stan’s integer constraints are not powerful enough to express the constraint that \(\texttt{r}\mathtt{\_}\texttt{t[j]} \leq \texttt{n}\mathtt{\_}\texttt{t[j]}\), but this constraint could be checked in the transformed data block.↩︎

  4. Stan的整数约束不足以表达约束条件 \(\texttt{r}\mathtt{\_}\texttt{t[j]} \leq \texttt{n}\mathtt{\_}\texttt{t[j]}\),但这个约束条件可以在变换数据模块中被写入。↩︎

  5. When dividing two integers, the result type is an integer and rounding will ensue if the result is not exact. See the discussion of primitive arithmetic types in the reference manual for more information.↩︎

  6. 当对两个整数进行除法运算时,结果类型为整数,并且如果结果不精确,则会进行舍入。更多相关信息,请参阅参考手册中原始算术类型的讨论。↩︎