Robust Adjusted Profile Score Estimation

Why Use an Adjusted Profile Score?

Inverse-variance weighted (IVW) estimation can be efficient when instruments are strong and pleiotropic effects do not distort the estimate. With weak instruments, however, the SNP–exposure associations are themselves measured with error. Treating their estimates as exact can bias estimation. Pleiotropic effects can also make the SNP–outcome associations more variable than the causal model alone predicts. Robust Adjusted Profile Score (MR-RAPS) addresses these issues in a summary-data model: it accounts for uncertainty in both association estimates and allows mean-zero pleiotropic effects to have nonzero variance.

The method does not make the instrumental-variable assumptions unnecessary. In particular, it does not make a nonzero average pleiotropic effect harmless, nor can it guarantee that extreme variants are invalid or that their influence should always be removed. Its value is a structured way to combine weak-instrument adjustment with a robust score when the model is a reasonable approximation.

Model and Assumptions

Definition: For variant $j$, let $\hat\gamma_j$ and $s_{Xj}$ be its estimated SNP–exposure association and standard error, and let $\hat\Gamma_j$ and $s_{Yj}$ be its estimated SNP–outcome association and standard error. MR-RAPS uses the summary-data model $$ \hat\gamma_j \mid \gamma_j \sim N(\gamma_j,s_{Xj}^2), \qquad \hat\Gamma_j \mid \gamma_j,\alpha_j \sim N(\beta\gamma_j+\alpha_j,s_{Yj}^2). $$ Here $\gamma_j$ is the true SNP–exposure association, $\beta$ is the causal effect, and $\alpha_j$ is the direct pleiotropic effect. A common random-effects assumption is $\alpha_j\sim N(0,\tau^2)$, independent of $\gamma_j$. Thus $\tau^2$ describes the variance of mean-zero pleiotropic effects, not their mean.

The model assumes the association estimates are approximately normal and that the standard errors are adequately estimated. It usually treats variants as independent, as is reasonable after suitable linkage-disequilibrium pruning; if variants are correlated, their covariance must instead be handled. The standard two-sample formulation also relies on suitable separation of estimation errors. Sample overlap or correlated estimation errors can require a different model.

The random-effects assumption permits heterogeneous direct effects, but it does not establish that their mean is zero. If pleiotropic effects have a systematic nonzero mean, the causal estimate can still be biased. Nor does the normal random-effects model describe every pattern of pleiotropy. A robust score reduces sensitivity to unusually large residuals; it does not repair every violation of the causal assumptions.

Weak-Instrument Adjustment

For a proposed value of $\beta$, define the residual and its model variance by

R_j(\beta)=\hat\Gamma_j-\beta\hat\gamma_j,\qquad V_j(\beta,\tau^2)=s_{Yj}^2+\beta^2s_{Xj}^2+\tau^2.

The $\beta^2s_{Xj}^2$ term appears because uncertainty in $\hat\gamma_j$ contributes to uncertainty in the residual. The $\tau^2$ term allows additional variation from pleiotropic effects. A central feature of RAPS is an adjustment to the score that accounts for measurement error in $\hat\gamma_j$. Without the adjustment, a score based on $\hat\gamma_j R_j$ has a nonzero expectation even at the true causal effect.

Key result: Under the normal measurement-error model, a basic adjusted profile score contribution is $$ U_j(\beta,\tau^2)= \frac{\hat\gamma_j\hat\Gamma_j-\beta(\hat\gamma_j^2-s_{Xj}^2)} {s_{Yj}^2+\beta^2s_{Xj}^2+\tau^2}. $$ The adjustment adds $\beta s_{Xj}^2$ to $\hat\gamma_j R_j(\beta)$. Summing these contributions gives an adjusted score whose root estimates $\beta$.

A short justification makes the correction clear. Write $\hat\gamma_j=\gamma_j+e_{Xj}$ and, at the true $\beta$, $R_j=\alpha_j+e_{Yj}-\beta e_{Xj}$. Under independent, mean-zero measurement errors, with $\alpha_j$ also mean zero and independent of $\gamma_j$,

E\{\hat\gamma_jR_j\} =E\{(\gamma_j+e_{Xj})(\alpha_j+e_{Yj}-\beta e_{Xj})\} =-\beta s_{Xj}^2.

The nonzero term comes from the shared exposure-estimation error: $e_{Xj}$ appears in both $\hat\gamma_j$ and $R_j$. Adding $\beta s_{Xj}^2$ makes the adjusted numerator have expectation zero under the model. This correction is important for weak instruments, for which $s_{Xj}^2$ need not be small compared with $\hat\gamma_j^2$.

Worked Example: Seeing the Weak-Instrument Correction

Suppose a hypothetical variant has $\hat\gamma=0.10$ exposure units per allele, $s_X=0.05$ exposure units per allele, $\hat\Gamma=0.08$ outcome units per allele, and $s_Y=0.04$ outcome units per allele. At the proposed causal effect $\beta=0.50$ outcome units per exposure unit, take $\tau=0$.

The residual is $R=0.08-(0.50)(0.10)=0.03$ outcome units per allele. Its variance is $V=0.04^2+0.50^2(0.05^2)=0.0016+0.000625=0.002225$, in squared outcome units per allele. The unadjusted numerator is $\hat\gamma R=(0.10)(0.03)=0.003$. The adjustment is $\beta s_X^2=(0.50)(0.0025)=0.00125$, so the adjusted numerator is $0.00425$.

Equivalently, it is $\hat\gamma\hat\Gamma-\beta(\hat\gamma^2-s_X^2)=(0.10)(0.08)-(0.50)(0.0100-0.0025)=0.008-0.00375=0.00425$. The score contribution is $0.00425/0.002225=1.910$. This single positive contribution says the score at $\beta=0.50$ is positive for this variant; it does not by itself determine the estimate from a set of variants.

Robust Loss and Overdispersion

The basic adjusted score can still be sensitive to a variant with an extreme standardized residual. A robust version replaces the linear residual contribution with a bounded-influence score. One common choice is the Huber loss. Let $r_j=R_j/\sqrt{V_j}$, and choose a tuning constant $c>0$. Its loss and score are

\rho_c(r)= \begin{cases} r^2/2,& |r|\le c,\\ c|r|-c^2/2,& |r|>c, \end{cases} \qquad \psi_c(r)= \begin{cases} r,& |r|\le c,\\ c\,\operatorname{sign}(r),& |r|>c. \end{cases}

For a normally calibrated adjusted robust score, a corresponding contribution is

U_{R,j}(\beta,\tau^2)= \frac{\hat\gamma_j}{\sqrt{V_j}}\psi_c(r_j) +\frac{\beta s_{Xj}^2}{V_j}a_c, \qquad a_c=E\{\psi_c'(Z)\}=P(|Z|\le c),\quad Z\sim N(0,1).

The second term is the weak-instrument correction for this robust score. It includes $a_c$, the expected derivative of the Huber score under a standard normal residual. In the quadratic region, Huber loss is the usual squared-error loss; beyond $c$, the score is capped rather than growing without bound. As $c$ becomes very large, $a_c$ approaches 1 and the adjusted score reduces to the non-robust form.

The overdispersion parameter $\tau^2$ is estimated along with $\beta$, generally using a profile or estimating procedure. It increases the residual variance used in standardization when pleiotropic effects create extra dispersion. A larger $\tau$ therefore changes both the residual scale and the relative weight of variants. It is not a measure of causal-effect uncertainty, and a large estimate is not proof that a particular variant is invalid.

Worked Example: How Huber Loss Limits an Extreme Residual

Suppose a hypothetical variant has standardized residual $r=3.0$ under a fitted model, and the Huber tuning constant is $c=1.5$. The quadratic loss would be $r^2/2=4.50$, with score $\psi(r)=3.0$. Under Huber loss, the value is instead $\rho_{1.5}(3)=1.5(3)-1.5^2/2=3.375$, and its score is $\psi_{1.5}(3)=1.5$.

Thus this residual still contributes evidence against the fitted model, but its score contribution is half the uncapped value. Huber loss does not delete the variant or declare it an outlier; it limits how rapidly its contribution increases. Whether a robust fit is preferable depends on the scientific model and the pattern of residuals, not on the fact that one residual is inconvenient.

Worked Summary-Data Analysis

Worked Example: An Illustrative MR-RAPS Estimate

Consider three hypothetical, independent variants. Suppose their exposure estimates are $0.20,\ 0.30,\ 0.10$, and their outcome estimates are $0.10,\ 0.15,\ 0.05$, respectively. For a hand calculation, take the common standard errors to be $s_X=0.05$ exposure units per allele and $s_Y=0.04$ outcome units per allele. Suppose a joint robust fit to the full analysis has supplied $\hat\tau=0.03$ outcome units per allele. To keep this calculation transparent, use $c=1.5$ and the normal calibration $a_c=P(|Z|\le1.5)\approx0.8664$.

At the estimate, the standardized residuals are within the Huber quadratic region, so $\psi_c(r_j)=r_j$ for all three variants. Because the standard errors and $V_j$ are common, setting the summed robust score to zero gives

\sum_j\{\hat\gamma_j\hat\Gamma_j-\beta\hat\gamma_j^2 +\beta s_X^2a_c\}=0, \qquad \hat\beta= \frac{\sum_j\hat\gamma_j\hat\Gamma_j} {\sum_j\hat\gamma_j^2-3s_X^2a_c}.

Substitute the values: $\sum_j\hat\gamma_j\hat\Gamma_j=(0.20)(0.10)+(0.30)(0.15)+(0.10)(0.05)=0.0700$. Also, $\sum_j\hat\gamma_j^2=0.20^2+0.30^2+0.10^2=0.1400$, and $3s_X^2a_c=3(0.0025)(0.8664)=0.006498$. Therefore,

\hat\beta=\frac{0.0700}{0.1400-0.006498} =\frac{0.0700}{0.133502} \approx 0.5243.

The fitted residual variance for each variant at this estimate is $V=0.04^2+(0.5243)^2(0.05^2)+0.03^2\approx0.003187$, in squared outcome units per allele. Its residual standard deviation is $\sqrt{V}\approx0.05646$ outcome units per allele. The outcome residuals are approximately $-0.00487,\ -0.00730,$ and $-0.00243$; dividing by $0.05646$ gives standardized residuals about $-0.0863,\ -0.1293,$ and $-0.0431$. These all lie within $c=1.5$, consistent with using the linear part of the Huber score in the calculation.

The estimate is about $0.5243$ outcome units per one-unit increase in the exposure, under the model and its causal assumptions. The supplied $\hat\tau=0.03$ represents estimated dispersion of direct effects on the outcome-association scale; it is not an additional causal effect. This arithmetic illustrates the point estimate, not its uncertainty. In an analysis, use the fitted procedure’s appropriate standard error or profile interval and report it with the estimate. For example, a hypothetical 95% interval of $0.21$ to $0.84$ outcome units per exposure unit would describe uncertainty in the estimated causal effect under the method’s assumptions; it would not mean there is a 95% probability that the fixed causal effect lies in that particular interval.

Diagnostics, Comparisons, and Limitations

A robust fit should be examined rather than treated as an automatic safeguard. Check the instrument-strength information and the association standard errors: an adjustment cannot manufacture information when exposure associations are extremely weak. Examine standardized residuals, the fitted $\hat\tau$, and how much the estimate changes between robust and non-robust fits. A large $\hat\tau$, many large residuals, or substantial estimate changes can signal model mismatch, but none identifies its cause by itself. Check that variants are aligned consistently and that their dependence has been handled; an independence-based score is not justified merely because the data are in summary form.

MR-RAPS addresses a different combination of concerns from several other estimators. Fixed-effect IVW does not allow extra pleiotropic dispersion, while multiplicative random-effects IVW inflates uncertainty for heterogeneity but does not provide the same weak-instrument correction. MR-Egger models an intercept for average directional pleiotropy under its own assumptions, including InSIDE; MR-RAPS instead uses a mean-zero random-effects model for direct effects. The weighted median and weighted mode rely on different assumptions about the distribution of valid instruments. These methods are not interchangeable robustness switches: compare estimates in light of their assumptions.

The robust score limits the influence of extreme standardized residuals, but it can also reduce efficiency when the normal model is correct. Results can depend on the tuning constant and the scale fit, especially with few variants. A near-zero fitted overdispersion does not prove absence of pleiotropy, and a positive value does not establish that pleiotropy is balanced in the relevant sense. Finally, robustification does not solve correlated pleiotropy, selection effects, weak identification so severe that the data contain little causal information, or systematic directional pleiotropy that violates the mean-zero assumption.

Common Pitfall

Wrong reasoning: “MR-RAPS is robust, so large pleiotropic effects and a nonzero average direct effect cannot bias its causal estimate.”
Correction: Huber scoring caps the influence of large standardized residuals, and $\tau^2$ models extra dispersion around a mean-zero direct-effect assumption. Neither feature guarantees protection from directional pleiotropy, correlated pleiotropy, or other violations of the causal model. State the assumptions and examine sensitivity rather than interpreting “robust” as “assumption-free.”

Check Your Understanding

Answer each question, then compare with the solutions.

What does $\tau^2$ represent in the MR-RAPS model, and what does it not represent?
Why does the adjusted score include a term involving $s_{Xj}^2$?
For Huber tuning constant $c=1.5$, what are $\psi_c(2)$ and $\rho_c(2)$?
Why is a robust MR-RAPS fit not guaranteed to remove bias from a nonzero mean pleiotropic effect?

Solutions. (1) It is the variance of mean-zero random direct pleiotropic effects, on the SNP–outcome association scale; it is not the variance of $\hat\beta$ or a causal effect. (2) Exposure-estimation error enters both $\hat\gamma_j$ and the residual $\hat\Gamma_j-\beta\hat\gamma_j$, giving $E\{\hat\gamma_jR_j\}=-\beta s_{Xj}^2$ at the true effect; the adjustment corrects this expectation. (3) Since $2>1.5$, $\psi_{1.5}(2)=1.5$, and $\rho_{1.5}(2)=1.5(2)-1.5^2/2=1.875$. (4) Robust scoring limits the contribution of extreme residuals, but MR-RAPS still relies on a mean-zero direct-effect model; a nonzero mean violates that assumption and can shift the causal estimate.