---
title: "Falsification controls and power curves: is the statistic real?"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Falsification controls and power curves: is the statistic real?}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
# the vignette runs in full in GitHub CI and for the package tarball (NOT_CRAN), on r-universe
# (MY_UNIVERSE), on the package website (IN_PKGDOWN) and in an interactive session; when CRAN
# re-builds it in its check, none of these is set and only this chunk runs
if (!interactive() && !any(nzchar(Sys.getenv(c("NOT_CRAN", "MY_UNIVERSE", "IN_PKGDOWN"))))) {
  knitr::knit_exit()
}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(rmoriebricklayer)
```

A capsule can prove that a statistic was computed from exactly these bytes. It cannot, by
itself, say whether the statistic means anything. Two complementary checks address that: a
**falsification control** asks whether the statistic survives deliberate attempts to break
it, and a **power curve** asks whether the procedure could have detected an effect of a given
size at all. Both are borrowed from the refutation tests of causal-inference practice
(placebo treatments, random subsets, permutation) and made to run on any statistic you can
write as a function of the data.

## Falsification controls

`cpsfals()` takes the data, a statistic (a function of the data frame) and the name
of the treatment column, and runs the controls that a real association must pass:

* **permutation**: shuffle the treatment; the statistic should collapse toward its null;
* **placebo**: replace the treatment with noise; the statistic should collapse;
* **subsets**: recompute on random subsets; the statistic should stay stable.

```{r falsification-and-power-falsify}
set.seed(1)
d <- data.frame(x = rnorm(200))
d$y <- 0.8 * d$x + rnorm(200)

real <- cpsfals(d, function(z) cor(z$x, z$y), trt = "x", n = 199, seed = 42)
real$controls[, c("control", "passed")]
```

A statistic that ignores the data fails the controls that can see it. A constant "0.5"
passes the subset check (it is perfectly stable) and fails the permutation check (it does
not collapse when the treatment is destroyed):

```{r falsification-and-power-falsify-fake}
fake <- cpsfals(d, function(z) 0.5, trt = "x", n = 199, seed = 42)
fake$controls[fake$controls$control == "permutation", "passed"]
```

This is the point of running controls rather than looking at a p-value: the p-value of a
constant is undefined, but its behaviour under permutation is not.

## Power curves

A null result is only informative if the procedure could have found something. `cpspwr()`
injects an effect of known size into the data, re-runs the statistic with its permutation
test, and reports the detection rate at each size:

```{r falsification-and-power-power}
set.seed(1)
d <- data.frame(x = rnorm(120), y = rnorm(120))
inj <- function(z, size) {
  z$y <- z$y + size * z$x
  z
}

pw <- cpspwr(d, function(z) cor(z$x, z$y), inj,
                    sizes = c(0, 0.3, 0.6), trt = "x",
                    n = 99, reps = 5, seed = 42)
pw$curve
```

The rate at size 0 estimates the false-positive rate. With few repetitions it is usually 0,
because alpha is small (five draws at 0.05 come up empty about three times in four), so read
it as a sanity check that it is not *large*, not as an estimate of alpha:

```{r falsification-and-power-power-zero}
pw$curve[pw$curve$size == 0, "rate"]
```

Use more `reps` for a real study; the example keeps them small so the vignette builds
quickly. The curve answers the question a reviewer will ask of a null finding: "what effect
size would you have detected with this data?"

## Families of statistics: the multiple-testing correction

When several statistics are examined, the smallest p-value among them is not the finding.
`falsfam()` applies a multiple-testing correction across a family of falsification
results, so one real association is still reported as real and the noise is not:

```{r falsification-and-power-family}
set.seed(1)
d <- data.frame(x = rnorm(150))
d$y <- 0.6 * d$x + rnorm(150)
d$a <- rnorm(150)
d$b <- rnorm(150)
fam <- list(
  real    = cpsfals(d, function(z) cor(z$x, z$y), trt = "x", n = 199, seed = 1),
  noisea = cpsfals(d, function(z) cor(z$a, z$y), trt = "a", n = 199, seed = 2),
  noiseb = cpsfals(d, function(z) cor(z$b, z$y), trt = "b", n = 199, seed = 3)
)
falsfam(fam)
```

It also accepts bare p-values, which is what stops the smallest of several from being read
as the finding:

```{r falsification-and-power-family-p}
falsfam(c(a = 0.01, b = 0.04, c = 0.2, d = 0.5))$results
```

## Sensitivity to an unmeasured confounder

For an observed risk ratio, the E-value (VanderWeele and Ding 2017) is the minimum strength
of association an unmeasured confounder would need with both treatment and outcome to
explain the result away:

```{r falsification-and-power-evalue}
evlrr(2)                       # needs a confounder associated by 3.41 with both
evlrr(0.5)                     # a protective effect is inverted; same answer
evlrr(2, lo = 0.9, hi = 4.4)   # an interval that includes the null needs nothing
evlrr(3, lo = 2.1, hi = 4.3)   # a strong result with a limit well above the null
```

## Where this fits in a capsule

`cpsrpt()` runs the falsification controls alongside the schema, rule and drift
checks when a statistic is supplied, so the capsule records not only that the number was
computed from these bytes but that it survived its controls. A pre-registered hypothesis
(`prrgdecl()`, `prrgchk()`) pins what was going to be tested before the data were
seen, which is the other half of keeping a finding honest.

## Several gaps, one family of tests

A report on Indigenous health usually tests many indicators at once. `falsfam()`
adjusts the family; here the differences from non-Indigenous adults in 2019/2022 for
First Nations adults, with z statistics from the published 95 percent intervals:

```{r falsification-and-power-indigenous-falsify}
hi <- indiex("health")
s <- hi[hi$period == "2019/2022" & hi$age == "Total, 18 years and over" &
          hi$quality == "ok", ]
fn <- s[s$group == "First Nations", ]
ni <- s[s$group == "Non-Indigenous", ]
both <- merge(fn, ni, by = "indicator", suffixes = c("fn", "ni"))
se <- function(l, u) (u - l) / (2 * 1.96)
z <- (both$estimatefn - both$estimateni) /
  sqrt(se(both$lowerfn, both$upperfn)^2 + se(both$lowerni, both$upperni)^2)
p <- stats::setNames(2 * stats::pnorm(-abs(z)), both$indicator)
res <- falsfam(p, meth = "holm")
head(res)
```

## The tests in a picture

`plot()` draws each result. For the real statistic, the observed correlation stands far outside
the correlations the permuted tables give. The fake statistic returns 0.5 whatever the data, so
every permuted table also gives 0.5: the null collapses onto the observed value, and the
permutation control fails, as it should. The last panel is the power of the test across injected
effect sizes.

```{r falsification-and-power-falsification-plots, fig.width = 7, fig.height = 4.5}
plot(real)
plot(fake)
plot(pw)
```


## References

VanderWeele, Ding (2017). Sensitivity analysis in observational research: introducing the
E-value. *Annals of Internal Medicine* 167(4).
Benjamini, Hochberg (1995). Controlling the false discovery rate. *JRSS B* 57(1).