---
title: "Methodology"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Methodology}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(respondeR)
```

This vignette sets out the statistics behind respondeR: the cut-point approach,
how responder probabilities and their variances are computed, each pooling
method, the relative measures, the threshold-free common-language effect size,
the standardized-mean-difference bridge, random effects, threshold uncertainty,
multi-arm trials, and the assumptions and their limits. It closes with a guide to
choosing a method and a note on how version 0.2.0 differs from 0.1.0.

## The cut-point approach

For one study arm with mean change $\mu$, standard deviation $\sigma$ and a
minimal important difference (MID) threshold $m$, assume the patient-level change
$X$ is Normally distributed. A *responder* is a patient whose change crosses the
threshold. The responder probability is

$$
p = \Pr(X > m) = \Phi\!\left(\frac{\mu - m}{\sigma}\right)
\quad\text{(higher change is better),}
$$

or $p = \Phi\!\left(\frac{m - \mu}{\sigma}\right)$ when a *lower* change is
better. This is the cut-point ("dichotomization") method reviewed by Thorlund and
colleagues (2011) and detailed by Anzures-Cabrera, Sarpatwari & Higgins (2011).
The between-arm contrast is then a familiar binary effect measure.

respondeR evaluates the requested tail directly on the log scale
(`pnorm(..., log.p = TRUE)`), so $\log p$ and $\log(1 - p)$ keep full precision
even when $p$ is far below machine precision. Risk ratios and odds ratios are
computed from these log probabilities, never from a clamped $p$.

## Sampling variance of a responder probability

The arm mean $\bar x$ and SD $s$ are estimates, so $p$ is too. By the delta
method, with $a = (\mu - m)/\sigma$ for the Normal model,

$$
\mathrm{Var}(\hat p) \approx \phi(a)^2\left[\frac{1}{n} + \frac{a^2}{2(n-1)}\right],
\qquad
\mathrm{Var}(\log \hat p) \approx \lambda(a)^2\left[\frac{1}{n} + \frac{a^2}{2(n-1)}\right],
$$

where $\lambda(a) = \phi(a)/\Phi(a)$ is the inverse Mills ratio, computed as
$\exp\{\log\phi(a) - \log\Phi(a)\}$ so it stays finite in the tails. The logit
scale uses $\lambda(a) + \lambda(-a)$. For the Student-$t$ and lognormal models
the gradients are numerical (central differences of $\log p$) and the joint
sampling distribution of $(\bar x, s)$ uses the model's own third and fourth
central moments $\mu_3, \mu_4$:

$$
\mathrm{Var}(s) \approx \frac{\mu_4 - \sigma^4 (n-3)/(n-1)}{4 n \sigma^2},
\qquad
\mathrm{Cov}(\bar x, s) \approx \frac{\mu_3}{2 n \sigma}.
$$

For the Normal model these reduce to $\sigma^2 / (2(n-1))$ and 0. These are
first-order (large-sample) formulas. In finite samples the sampling variance of
an SD converges slowly under heavy tails or skew, so for the $t$ ($\nu > 4$) and
lognormal models the intervals tend to be somewhat too wide in the tails of the
distribution. A $t$ model with $\nu \le 4$ degrees of freedom has no finite
fourth moment; the Normal-theory value is then used, which makes the intervals
much too narrow, and respondeR warns. Treat the non-Normal models as
sensitivity analyses.

Each per-arm variance is floored at half a patient's resolution,
$p_0(1 - p_0)/n$ with $p_0 = 0.5/n$, so that a probability deep in a tail (whose
model-based variance underflows) can never receive an infinite weight.

## The pooling methods

Studies report per-arm summaries; the methods differ in how those are combined.
Study $i$ contributes $(\bar d_{e,i}, s_{e,i}, n_{e,i})$ and
$(\bar d_{c,i}, s_{c,i}, n_{c,i})$.

### Individual (recommended)

Dichotomize each study, form its contrast, then pool. Every comparison stays
within a randomized study. With $p_{e,i}$ and $p_{c,i}$ as above,
$\mathrm{RD}_i = p_{e,i} - p_{c,i}$, and the log risk ratio and log odds ratio
are differences of exact log probabilities.

The per-study variance follows `se_method`:

* `"delta"` (default): the delta-method variances above, which carry the
  uncertainty in the reported mean and SD.
* `"binomial"`: the pseudo-binomial $p(1-p)/n$, which treats the model-based
  probability as if it were an observed proportion of dichotomized patients
  (the version 0.1.0 default).

**Fixed effects.** The risk difference is pooled with sample-size weights
$a_i = n_{e,i} n_{c,i} / (n_{e,i} + n_{c,i})$ (the Mantel-Haenszel risk
difference weights), with variance $\sum_i a_i^2 V_i / (\sum_i a_i)^2$
(`rd_weights = "sample_size"`, the default). Inverse-variance weights
(`rd_weights = "inverse_variance"`) let any study whose responder probabilities
sit near 0 or 1, and so whose model-based variance is tiny, dominate the pooled
risk difference, which is exactly what happens in an ordinary MID sweep. The log
risk ratio and log odds ratio are pooled by inverse variance.

```{r}
responder_analysis(sample_responder_data, mid = 1, method = "individual")
```

respondeR warns whenever one study holds more than half of the weight for a
pooled measure (with three or more studies).

### Weighted mean

Pool the summaries, then dichotomize. respondeR pools the **control mean** and
the **within-study mean difference** separately by inverse variance (fixed or
random effects, following `pooling`), and sets the experimental mean to their
sum:

$$
\hat\mu_c = \frac{\sum_i w_{c,i} \bar d_{c,i}}{\sum_i w_{c,i}},\quad
\widehat{\mathrm{MD}} = \frac{\sum_i w_{d,i} (\bar d_{e,i} - \bar d_{c,i})}{\sum_i w_{d,i}},\quad
\hat\mu_e = \hat\mu_c + \widehat{\mathrm{MD}},
$$

with $w_{c,i} = n_{c,i}/s_{c,i}^2$ and
$w_{d,i} = 1/(s_{e,i}^2/n_{e,i} + s_{c,i}^2/n_{c,i})$ (plus $\tau^2$ under random
effects). Each arm's SD is the within-study pooled SD
$S = \sqrt{\sum_i (n_i - 1)s_i^2 / \sum_i (n_i - 1)}$.

The estimator depends on every study's arm means and SDs, including through the
SD-dependent weights and the pooled SDs. Its sampling covariance is the full
Jacobian form $\sum J V_{\text{arm}} J^\top$ over all arms, where $J$ holds the
derivatives of the arm responder probabilities (on the probability, log and
logit scales) with respect to that arm's mean and SD, and $V_{\text{arm}}$ is
the joint sampling covariance of the arm mean and SD under the chosen
distribution (including their covariance for the skewed lognormal model).
Because the control arm enters both pooled means, the two are negatively
correlated, and the Jacobian carries that automatically. Under random effects
the between-study variances of the two pooled means,
$\hat\tau^2 \sum_i \tilde w_i^2$ with normalized weights $\tilde w$, are added;
with $\hat\tau^2 = 0$ and a Wald interval the result equals the fixed-effect
one. `ci_method = "hksj"` or `"hksj_adhoc"` multiplies the variance path of each
pooled mean by its HKSJ factor $q$ (scaling its gradient by $\sqrt q$, which
keeps the joint covariance positive semidefinite) and uses a $t_{k-1}$ critical
value. The between-study variances are held at their estimates.

Anchoring on the within-study difference matters. Pooling the two arms
separately (as version 0.1.0 did) compares patients across studies and can show
a Simpson's-paradox bias when allocation ratios differ between studies (Altman
and Deeks, 2002):

```{r}
paradox <- data.frame(
  study = c("A", "B"), change_e = c(2.5, 0.5), sd_e = 1, n_e = c(400, 20),
  change_c = c(2.0, 0.0), sd_c = 1, n_c = c(20, 400)
)
responder_rd_individual(paradox, mid = 1)[, c("study", "rd")]
as.data.frame(responder_analysis(paradox, mid = 1,
  method = c("individual", "weighted")))[, c("method", "rd")]
```

Separate arm pooling gave a risk difference of about 0.74 here, although neither
trial shows more than 0.15.

In a simulation (1000 replicates; 10 studies of 100 patients per arm; Normal
change scores with SD 1; control means drawn from $N(0.3, 0.5^2)$; a mean
difference of 0.5, fixed or drawn from $N(0.5, 0.2^2)$; MID 1), the 95 percent
weighted interval covered the realized-study target at 0.94 to 0.96 with fixed
effects and the population target at 0.93 to 0.96 with random effects. With 5
studies the random-effects Wald interval covered at 0.91, and with
`ci_method = "hksj_adhoc"` at 0.96.

### Unweighted mean and median

Dichotomize the arithmetic mean or the median of the study means and SDs. These
are descriptive robustness summaries with **no variance model**, so respondeR
reports point estimates only.

### Baseline risk: matched or median control

By default (`control = "matched"`) the control proportion comes from the same
pooled summaries. The simulation study that motivated this package (Sofi-Mahmudi,
2024) held the baseline risk at the **median control arm**; that is available as
`control = "median"`. For the weighted method the experimental mean is then the
median control mean plus the pooled within-study difference. Point estimates
only.

```{r}
responder_analysis(sample_responder_data, mid = 1, control = "median")
```

## Relative effect measures and the NNT

$$
\mathrm{RR} = \frac{p_e}{p_c}, \quad
\mathrm{OR} = \frac{p_e/(1 - p_e)}{p_c/(1 - p_c)}, \quad
\mathrm{NNT} = \frac{1}{\mathrm{RD}}.
$$

Confidence intervals for RR and OR are formed on the log scale. Following Altman
(1998), a positive risk difference is reported as a number needed to treat for
one more responder (NNTB) and a negative one as a number needed to treat for one
fewer (NNTH); when the risk-difference interval includes zero the NNT interval
passes through infinity and is printed in that form.

The responder risk ratio deserves care. When responders are common in the
control arms, $\log p_c$ is close to 0 with a small variance, so the pooled risk
ratio is compressed toward 1. The odds ratio is symmetric in responders and
non-responders and, under the cut-point model, closer to constant across
baseline risks; it is the default route for summary-of-findings tables.

## Summary-of-findings rows

`responder_sof()` applies the pooled relative effect to an assumed control risk
and reports, per 1000 patients, the control risk, the intervention risk and the
difference with their intervals (Guyatt and colleagues, 2013). The default
control risk is the median control-arm responder probability:

```{r}
ind <- responder_analysis(sample_responder_data, mid = 1, method = "individual")
responder_sof(ind, control_risk = c(0.1, 0.3))
```

## Common-language effect size (threshold-free)

The common-language effect size (CLES, the probabilistic index) is the
probability that a randomly chosen treated patient has a better change than a
randomly chosen control. Under a Normal model it is exact:

$$
\mathrm{CLES} = \Phi(\delta), \qquad
\delta = \frac{\mu_e - \mu_c}{\sqrt{\sigma_e^2 + \sigma_c^2}}.
$$

Per-study $\delta_i$ are pooled by inverse variance and back-transformed. It is
a population-level probability, not the chance that treatment helps a particular
patient.

```{r}
cles <- responder_cles(sample_responder_data)
c(cles = cles$cles, lb = cles$cles_lb, ub = cles$cles_ub)
```

## The SMD bridge (`method = "smd"`)

The second approach of Anzures-Cabrera et al. (2011) pools the standardized mean
difference and maps it to an odds ratio. respondeR pools Hedges' $g$ and uses the
logistic-distribution conversion $\ln\mathrm{OR} = \frac{\pi}{\sqrt 3}\, g$
(Hasselblad and Hedges, 1995; Chinn, 2000). (The Cox conversion uses a factor of
about 1.65 instead of $\pi/\sqrt 3 \approx 1.81$.) The odds ratio interval is the
pooled $g$ interval transformed, so `ci_method` applies. Risks follow from the
pooled control risk $p_c$:
$p_e = \mathrm{logit}^{-1}(\mathrm{logit}\, p_c + \ln\mathrm{OR})$.
Because $p_c$ is itself estimated, the risk-difference and risk-ratio intervals
use the delta method with
$\partial\mathrm{RD}/\partial \ln\mathrm{OR} = p_e(1 - p_e)$ and
$\partial\mathrm{RD}/\partial p_c = p_e(1-p_e)/\{p_c(1-p_c)\} - 1$, including
the covariance between $g$ and the pooled control mean. That covariance comes
from the full Jacobian of the within-study sampling errors; under random effects
the same correlation is assumed between the pooled $g$ and the pooled control
mean (an assumption with a visible effect on the risk-difference interval), and
the control mean is pooled with a Wald variance while `ci_method` applies to
$g$. The prediction interval is transformed to the risk-difference scale at
$p_c$.

```{r}
responder_analysis(sample_responder_data, mid = 1, method = "smd")
```

## Converting a pooled effect fitted elsewhere

`responder_from_effect()` applies the same bridge to a single pooled effect, for
example from a published meta-analysis, given a baseline control responder risk;
`responder_analysis()` calls it for you when you pass a `metafor` `rma` object,
and reads the arms directly when you pass a `meta::metacont` object. The effect
is taken as experimental minus control and oriented by `direction`, so a pain
reduction with `direction = "lower"` raises the responder risk. For a mean
difference the experimental proportion is the control probability shifted along
the Normal cut-point,

$$
p_e = \Phi\!\left(\Phi^{-1}(p_c) + \frac{\mathrm{MD}}{\sigma}\right),
$$

which needs a reference within-arm SD $\sigma$ (`sd_ref`). With a fixed baseline
probability the intervals are conditional on it, as in a GRADE
summary-of-findings table; with an estimated baseline (a mean, SD and $n$) its
uncertainty is propagated jointly with that of the effect. See
`vignette("meta-analysis-objects")` for worked examples.

## Random effects and heterogeneity

The individual, weighted and SMD methods can use random effects
(`pooling = "random"`) with DerSimonian-Laird or REML (`tau_method = "REML"`,
via *metafor*) between-study variance. Cochran's $Q$, $I^2$ and $\tau^2$ are
inverse-variance based and refer to the scale that is pooled: the risk
difference for `individual`, Hedges' $g$ for `smd`.

Three interval methods are available for the pooled estimate:

* `"wald"`: Normal critical value; can under-cover with few studies.
* `"hksj"`: Hartung-Knapp-Sidik-Jonkman, a $t_{k-1}$ interval scaled by the
  observed dispersion. When the studies agree more closely than their
  within-study variances predict, it can be **much narrower** than Wald.
* `"hksj_adhoc"`: HKSJ with the ad hoc floor $\max(1, q)$ of Rover, Knapp and
  Friede (2015), never narrower than Wald. Recommended when there are few
  studies.

```{r}
near_identical <- data.frame(
  study = 1:4, change_e = c(1, 1.01, 0.99, 1), sd_e = 1, n_e = 100,
  change_c = 0, sd_c = 1, n_c = 100
)
sapply(c("wald", "hksj", "hksj_adhoc"), function(ci) {
  r <- responder_analysis(near_identical, mid = 0.5, method = "individual",
                          pooling = "random", ci_method = ci)
  r$rd_ub - r$rd_lb
})
```

The prediction interval is $\hat\mu \pm t_{k-2}\sqrt{\hat\tau^2 + V}$, where $V$
is the variance that built the confidence interval, for every `tau_method`
(Higgins, Thompson and Spiegelhalter, 2009). It is not reported for fewer than
three studies and is unstable for small $k$.

## Uncertainty in the MID

If the threshold is itself estimated, `mid_sd` treats it as uncertain. The MID is
**one** threshold shared by every study and both arms, so its uncertainty is
fully correlated across studies. respondeR pools without it and then adds it
once, through the derivative of the whole pooled estimator with respect to the
MID:

$$
\mathrm{Var}_{\text{total}}(\hat\theta) = \mathrm{Var}_{\text{sampling}}(\hat\theta)
+ \left(\frac{\partial \hat\theta}{\partial m}\right)^2 \mathrm{SD}(m)^2 .
$$

It therefore does not change the weights, $\tau^2$ or the point estimate, and it
does not shrink as studies are added. A complementary approach is to report the
results across a range of thresholds with `responder_curve()`:

```{r, fig.width = 6, fig.height = 3.6, fig.alt = "Pooled risk difference across MID thresholds"}
plot(responder_curve(sample_responder_data, mids = seq(0, 2, by = 0.25)), at = 1)
```

The bands are pointwise.

## Multi-arm trials

A trial with two active arms and one control contributes two comparisons that
reuse the same control patients. Entering them as independent studies counts
those patients twice and overstates the precision. Mark such comparisons with a
`trial` column: when comparisons of one trial share an identical control arm,
respondeR divides the control sample size equally among them (Cochrane Handbook,
section 23.3.4). The bundled `vas_pain` data use this for two three-arm trials.

## Final values and absolute thresholds

A responder can also be defined by the outcome reaching an absolute state, for
example a pain score of 3 or less at follow-up, rather than by a change. The same
cut-point model applies to the final-value distribution: put the final means and
SDs in the `change_e`/`change_c` columns and the absolute threshold in `mid`.
The interpretation differs (the proportion in the target state, not the
proportion who improved by the MID), so label the analysis accordingly.

```{r}
final <- data.frame(
  study = c("A", "B"), change_e = c(2.6, 3.1), sd_e = c(1.9, 2.1), n_e = c(60, 45),
  change_c = c(3.8, 4.2), sd_c = c(2.0, 2.2), n_c = c(58, 47)
)
# Proportion with a final pain score of 3 or less (lower is better).
responder_analysis(final, mid = 3, direction = "lower", method = "individual")
```

## Interval options

* `ci_type = "logit"` builds the weighted method's risk-difference interval by
  the MOVER method from logit-scale arm intervals, allowing for their
  correlation, so it stays within $[-1, 1]$. It also applies to
  `responder_rd_individual()`. Pooled `individual` and `smd` intervals are
  formed on the scale that is pooled. A nonzero `mid_sd` correlates the arms in a
  way MOVER cannot carry, so Wald intervals are used, with a warning.
* `dist = "lognormal"` or `"t"` changes the change-score model as a sensitivity
  analysis for skewed or heavy-tailed data.

## Assumptions and limitations

* **Change-score distribution.** The probabilities assume a parametric model
  for the patient-level change within each arm. Means and SDs alone cannot
  establish that the model fits; skewed outcomes can bias the responder
  proportions. Use `dist` as a sensitivity analysis.
* **Summary-statistic input.** Only means, SDs and sample sizes are used.
* **Choice of MID.** Results depend on the threshold. Report the MID, show the
  curve across plausible thresholds, and consider the threshold-free CLES.
* **Normal-approximation intervals.** Wald intervals can misbehave for extreme
  proportions or tiny samples, and the delta method is a first-order
  approximation.

## Choosing a method

| If you want… | Use |
|--------------|-----|
| A defensible default that keeps comparisons within trials | `individual` |
| A pool-then-dichotomize estimator with an interval | `weighted` |
| A descriptive robustness summary | `median` / `unweighted` (point estimates) |
| A cross-check through a different bridge to risks | `smd` |
| Few studies, random effects | `pooling = "random", ci_method = "hksj_adhoc"` |
| To show threshold dependence | `responder_curve()` |
| A GRADE summary-of-findings row | `responder_sof()` |
| To avoid choosing a threshold | `responder_cles()` |

## Changes from version 0.1.0

Version 0.2.0 changes several results. The per-study variance defaults to the
delta method; fixed-effect risk differences use sample-size weights; ratio
measures use exact log probabilities; MID uncertainty is added after pooling;
the weighted method anchors on the within-study difference; the SMD intervals
keep the chosen interval method and the control-risk uncertainty; prediction
intervals follow one rule; and shared control arms of multi-arm trials are
split. `se_method = "binomial", rd_weights = "inverse_variance"` reproduces the
version 0.1.0 individual point estimate.

## Relation to the thesis

The thesis behind this package (Sofi-Mahmudi, 2024) evaluated the median,
unweighted, weighted and individual strategies by simulation, pooling means and
SDs by inverse variance and converting a pooled risk ratio to a risk difference
with a representative baseline risk. The estimators here follow its framework
but differ in detail (the pooled SD, the anchored weighted method, the variance
models and the interval options), so the thesis simulations do not validate
every option in this package.

## References

Sofi-Mahmudi, A. (2024). Identifying an optimal strategy for converting pain as
a continuous outcome to a responder analysis [Master's thesis, McMaster
University]. MacSphere. https://hdl.handle.net/11375/30210

Thorlund, K., Walter, S. D., Johnston, B. C., Furukawa, T. A., & Guyatt, G. H.
(2011). Pooling health-related quality of life outcomes in meta-analysis: a
tutorial and review of methods for enhancing interpretability. *Research
Synthesis Methods*, 2(3), 188 to 203. doi:10.1002/jrsm.46

Anzures-Cabrera, J., Sarpatwari, A., & Higgins, J. P. T. (2011). Expressing
findings from meta-analyses of continuous outcomes in terms of risks.
*Statistics in Medicine*, 30(25), 2967 to 2985. doi:10.1002/sim.4298

Altman, D. G. (1998). Confidence intervals for the number needed to treat.
*BMJ*, 317(7168), 1309 to 1312.

Altman, D. G., & Deeks, J. J. (2002). Meta-analysis, Simpson's paradox, and the
number needed to treat. *BMC Medical Research Methodology*, 2, 3.

Chinn, S. (2000). A simple method for converting an odds ratio to effect size for
use in meta-analysis. *Statistics in Medicine*, 19(22), 3127 to 3131.

Guyatt, G. H., Oxman, A. D., Santesso, N., et al. (2013). GRADE guidelines: 12.
Preparing summary of findings tables: binary outcomes. *Journal of Clinical
Epidemiology*, 66(2), 158 to 172.

Hasselblad, V., & Hedges, L. V. (1995). Meta-analysis of screening and
diagnostic tests. *Psychological Bulletin*, 117(1), 167 to 178.

Higgins, J. P. T., Thompson, S. G., & Spiegelhalter, D. J. (2009). A
re-evaluation of random-effects meta-analysis. *Journal of the Royal Statistical
Society: Series A*, 172(1), 137 to 159.

McGraw, K. O., & Wong, S. P. (1992). A common language effect size statistic.
*Psychological Bulletin*, 111(2), 361 to 365.

Rover, C., Knapp, G., & Friede, T. (2015). Hartung-Knapp-Sidik-Jonkman approach
and its modification for random-effects meta-analysis with few studies. *BMC
Medical Research Methodology*, 15, 99. doi:10.1186/s12874-015-0091-1