---
title: "GRIM in depth"
output:
  rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{GRIM in depth}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
# Print enough significant digits that values like `10.46` are not
# abbreviated to `10.5` in the output tibbles:
options(pillar.sigfig = 4)
```

```{r setup, include = FALSE}
library(scrutiny)
```

This vignette builds on `vignette("grim")`, which covers the basics of GRIM testing with `grim()` and `grim_map()`. Here are the numbers behind a test result, the settings recorded with it, the values around the reported ones with `grim_map_seq()` and `grim_map_total_n()`, and the GRIM statistics behind the `exp_incon` column.

## Behind a test result

### Reconstructed values

With `show_rec = TRUE`, the output includes the numbers on which the test result is based:

```{r}
pigs1 |>
  grim_map(digits_x = 2, items_averaged = 1, data_type = "integer", show_rec = TRUE) |>
  dplyr::select(x, n, consistency, rec_sum:rec_x_lower)
```

- `rec_sum` is `x * n * items_averaged`, the sum implied by the reported mean. The sum of integer data is a whole number, which `rec_sum` generally is not.
- `sum_lower` and `sum_upper` are the smallest and largest whole-number sums that would have been reported as `x`. A value set is consistent if and only if `sum_lower <= sum_upper`. If it is inconsistent, they are the whole numbers just above and below `rec_sum`, so their difference is always 1 and says nothing about how far off the mean is.
- `rec_x_upper` and `rec_x_lower` are the two achievable means closest to `x`: `rec_sum` rounded up and down, divided by `n * items_averaged`. They show how far the reported mean is from an achievable one. With `x_meaning = "percentage"`, they are percentages like `x`.

### Recorded settings

A mapper records the settings of its call on its output, so that a report can state the conditions of each result. `scrutiny_meta()` returns them:

```{r}
pigs1 |>
  grim_map(digits_x = 2, items_averaged = 1, data_type = "integer", rounding = "up") |>
  scrutiny_meta()
```

`default` says whether a setting was left at its default. `explain_map()` and `grim_plot_grid()` read the settings from there, so they apply the same rounding method as the call. For the same reason, `dplyr::bind_rows()` refuses to combine the output of calls with different settings, whose rows would otherwise be read under the first call's settings.

### Achievable means

`grim_values()` returns all means that a sample of `n` integers can have and that would be reported as `x`, each with the whole-number sum behind it:

```{r}
grim_values(x = 5.19, n = 32, digits_x = 2, items_averaged = 1, data_type = "integer")
grim_values(x = 5.19, n = 300, digits_x = 2, items_averaged = 1, data_type = "integer")
grim_values(x = 5.19, n = 28, digits_x = 2, items_averaged = 1, data_type = "integer") # inconsistent: none
```

`grim_closest()` returns the achievable mean nearest to `x`, with the sum behind it, in one row per value set. If `x` is consistent, that is the nearest one that would be reported as `x`:

```{r}
grim_closest(x = 5.19, n = 28, digits_x = 2, items_averaged = 1, data_type = "integer")
```

Both functions are vectorized. A reported mean may have one achievable mean, several, or none, so `grim_values()` returns one row per achievable mean, and `case` says which value set it belongs to. An inconsistent value set has no rows:

```{r}
grim_values(
  x = c(5.19, 5.19, 5.19),
  n = c(32, 300, 28),
  digits_x = 2,
  items_averaged = 1,
  data_type = "integer"
)
```

To line the values up with the data they came from, number its rows as `case` and join. An inner join leaves out the inconsistent rows, and `suffix` keeps the achievable `x` apart from the reported one:

```{r}
pigs1 |>
  dplyr::mutate(case = dplyr::row_number()) |>
  dplyr::inner_join(
    grim_values(pigs1$x, pigs1$n, digits_x = 2, items_averaged = 1, data_type = "integer"),
    by = "case",
    suffix = c("", "_achievable")
  )
```

Once `n * items_averaged` is at least `10^digits_x`, every mean is achievable and GRIM cannot detect anything. The number of achievable means can then become very large, which is why `grim_map()` does not return them.

### A single value set: `grim_plot_line()`

`grim_plot_line()` draws the test for one value set, which helps to explain a result. Of all means that could have been reported with `digits_x` decimal places, it draws a square for each one that a sample of size `n` can have and leaves the others out. The reported value is marked above the row, and the strip behind the squares is the range of unrounded means that would have been reported as `x`. The value set is consistent if a square stands within the strip:

```{r, fig.width=7, fig.height=2}
grim_plot_line(x = 5.19, n = 28, digits_x = 2, items_averaged = 1, data_type = "integer")
```

With `n = 26`, the same mean is consistent:

```{r, fig.width=7, fig.height=2}
grim_plot_line(x = 5.19, n = 26, digits_x = 2, items_averaged = 1, data_type = "integer")
```

By default, the plot shows the reported value and the nearest achievable one. `zoom = "medium"` widens it to the nearest achievable value on either side, and `zoom = "wide"` to a whole unit.

`grimmer_plot_line()` is the same figure over standard deviations; see `vignette("grimmer-in-depth")`.

### The grid behind `grim_plot_grid()`

`grim_grid()` returns the grid behind the raster of `grim_plot_grid()` as a tibble, with GRIM's result for every cell. It is useful for studying GRIM itself, e.g., the share of means it rules out at each sample size:

```{r}
grim_grid(digits_x = 2, data_type = "integer") |>
  dplyr::summarise(share_inconsistent = mean(!consistency), .by = n)
```

## Values around the reported ones: `grim_map_seq()`

An inconsistency may come from a small error in reporting or computing. `grim_map_seq()` varies `x` and `n` a few steps up and down, one at a time, and tests each resulting value set:

```{r}
out_seq1 <- grim_map_seq(pigs1, digits_x = 2, items_averaged = 1, data_type = "integer")
out_seq1
```

By default, only inconsistent value sets are dispersed. `case` is the row in the input data, `var` the variable that was varied, and `diff_var` the number of steps from the reported value.

### Summaries with `audit()`

`audit()` summarizes this output per case:

```{r}
audit(out_seq1)
```

- `x` and `n` are the reported values, and `consistency` their test result.
- `hits_total`, `hits_x`, `hits_n`: the number of consistent value sets found, in total and by varying each variable.
- `diff_x`: the distance from `x` to the nearest consistent value, in dispersion steps. `diff_x_up` and `diff_x_down` are the distances to the nearest consistent value above and below. The `diff_n*` columns are the same for `n`.

The default `dispersion` is `1:5`, for five steps up and down. With more steps, more hits are typically found:

```{r}
out_seq2 <- grim_map_seq(
  pigs1,
  digits_x = 2,
  items_averaged = 1,
  data_type = "integer",
  dispersion = 1:10
)
audit(out_seq2)
```

`explain_print()` explains the rows of this output like those of `grim_map()`, so you can see why a nearby value passes:

```{r}
out_seq1 |>
  dplyr::filter(consistency) |>
  explain_print(1)
```

### Plotting dispersed values

`grim_plot_grid()` accepts this output. Each reported value set appears as a cross: the horizontal arm consists of the dispersed `n` values, the vertical arm of the dispersed `x` values:

```{r, fig.width=6, fig.height=5.5}
grim_plot_grid(out_seq1)
```

Restrict the dispersion to one variable with `var`:

```{r, fig.width=6, fig.height=5.5}
out_seq1_only_n <- grim_map_seq(
  pigs1,
  digits_x = 2,
  items_averaged = 1,
  data_type = "integer",
  var = "n"
)
grim_plot_grid(out_seq1_only_n)
```

## Unknown group sizes: `grim_map_total_n()`

Some studies report group means but only the total sample size. GRIM can then test hypothetical group sizes: start from an even split, move one participant at a time from one group to the other, and record which splits make both means consistent. `grim_map_total_n()` does this, following an analysis by @bauer_expression_2021.

The study they examined reported means of 4.71 and 5.3 for a total of 40 participants. With 20 in each group, only 5.3 is consistent. Bauer and Francis moved to 19/21, 18/22, and 17/23, and found that both means are consistent at 17/23 with 4.71 in the smaller group.

```{r}
df <- tibble::tibble(x1 = 4.71, x2 = 5.3, n = 40)

df_tested <- grim_map_total_n(df, digits_x = c(2, 1), items_averaged = 1, data_type = "integer")
df_tested

audit(df_tested)
```

Each row is one mean with one hypothetical group size. `n_change` is how far that size is from half the total, `both_consistent` says whether both means of the split are consistent, and `dir` says which way the means were paired with the groups (see below). A hit is a split in which both means are consistent with their hypothetical group sizes. The default `dispersion = 0:5` goes five steps away from the even split.

The two means were reported with different numbers of decimal places, so `digits_x` has one value per group. Testing 5.3 as if it were 5.30 would be stricter than the test Bauer and Francis conducted.

`grim_map_total_n()` tests both ways of pairing the means with the group sizes. Bauer and Francis only paired 4.71 with the smaller group; the reverse pairing yields two more hits, at 19/21 and 16/24. In the summary, `hits_forth` counts the hits from the first pairing, `hits_back` those from the reverse one, and `hits_total` their sum. `scenarios_total` is the number of scenarios tested, each a split with one pairing of the means, and `hit_rate` the share of them that are hits. The even split is the same scenario under both pairings, so it counts once: six splits make 11 scenarios here.

## GRIM statistics

### The probability of GRIM inconsistency

The `exp_incon` column of `grim_map()` is, for the simplest rounding methods,

$$
P = \max\left(0, \frac{10^D - NL}{10^D}\right)
$$

where $D$ is the number of decimal places, $N$ the sample size, and $L$ the number of items. This is the probability that a mean with $D$ decimal places but random decimal digits is GRIM-inconsistent with $N$ and $L$. If it is 0, GRIM cannot detect any inconsistency.

The formula is exact for rounding methods in which exactly one range of unrounded values, one rounding step wide, is rounded to each reported value, such as `"up"`. The default, `"up_or_down"`, accepts a mean exactly halfway between two reportable values in both directions, so the formula overstates the probability whenever $NL < 10^D$ and $NL / \gcd(NL, 10^D)$ is even: at $N = 40$ with two decimal places, the probability is $0.4$, not $0.6$. Under `"ceiling_or_floor"`, the range is two steps wide. `grim_exp_incon()` and `grim_total()` therefore take `rounding` and `threshold` arguments and count exactly the values that `grim()` would reject. `grim_map()` passes its own `rounding` on.

```{r}
grim_exp_incon(x = 1.40, n = 72, digits_x = 2, items_averaged = 1, data_type = "integer")

grim_exp_incon(x = 5.93, n = 80, digits_x = 2, items_averaged = 3, data_type = "integer")

grim_exp_incon(
  x = 84.27,
  n = 40,
  digits_x = 2,
  items_averaged = 1,
  data_type = "integer",
  x_meaning = "percentage"
)
```

`grim_total()` returns a count instead of a proportion: the number of the $10^D$ possible fractional parts that are inconsistent:

```{r}
grim_total(x = 1.40, n = 72, digits_x = 2, items_averaged = 1, data_type = "integer")
```

The formula is implicit in @brown_grim_2017, who note (p. 367) that with two decimal places and $N < 100$, a random mean is consistent in about $N$% of cases.

### A study's rate of inconsistency

Correctly reported means are always consistent, so a study's rate of inconsistent value sets should be zero. If the decimal digits of its means were random, that rate would be expected to be the mean of the `exp_incon` column. `grim_binomial()` tests whether the rate of *consistent* value sets is above what random digits would produce. It uses the Poisson binomial distribution, since the probability differs between value sets, and returns the observed consistency rate with a Clopper-Pearson confidence interval:

```{r}
pigs1 |>
  grim_map(digits_x = 2, items_averaged = 1, data_type = "integer") |>
  grim_binomial()
```

`grim_binomial_power()` computes the power of this test for a given number of value sets, so you can judge whether the test is informative before interpreting its p-value:

```{r}
grim_binomial_power(p0 = 0.5, p1 = 0.9, k = 10)
```

Both functions are experimental. `grim_plot_rates()` shows the same comparison for several studies. Each point is a study, placed at its observed rate of inconsistency against the rate that random means would show; the point area is the number of testable value sets. A study on the dotted line at zero is as expected, and one on the dashed diagonal is as inconsistent as random digits:

```{r, eval = rlang::is_installed("ggrepel"), fig.width=6, fig.height=5}
# fmt: skip
codes <- c(
  "B", "B", "C", "C", "A", "C",
  "A", "D", "A", "A", "A", "A"
)

pigs1 |>
  grim_map(digits_x = 2, items_averaged = 1, data_type = "integer") |>
  dplyr::mutate(study = codes) |>
  grim_plot_rates(study)
```

### Granularity

The granularity of a scale is the smallest difference between two means it can produce, $G = 1 / (NL)$. `grim_granularity()` computes it:

```{r}
grim_granularity(n = 80, items_averaged = 5)
```

# References