Permutation test for paired data

[mine-cetinkaya-rundel]

Permutation test for paired data

[andrewpbray]

I just bonked my head into this gap in infer this week too! This fix seems the most straightfoward.

[hardin47]

1. I'm totally okay with the single sample framework of hypothesize(null = "point", mu = 0). But what happens if someone types in a different value for mu?

2. Fine if the wording is plain language. I think that in the t.test() function there is a paired = TRUE option, so adding that option to generate() would be consistent with the base R language. But I also like type = "permute paired".

[mine-cetinkaya-rundel]

I haven't quite wrapped my head around the right/better approach for (1) yet, but I wanted to chime in for (2): I think I prefer type = "permute" with paired = TRUE to emphasize that this is a variant of a "permute" operation.

[simonpcouch]

Thanks for bringing these up, @hardin47. I hear you that we don't yet have enough information to do anything about (1) at the hypothesize() step—maybe we wait to check mu inside of the generate() internals for permute paired = TRUE, and warn with the "results may not mean what you think they mean" sort of message used elsewhere? Delaying the warning for the mistake made in hypothesize() until generate() is definitely a bit awkward, though there's precedent for that elsewhere in the package.

[andrewpbray]

Hmm, I'm not sure what to do about 1.

hypothesize() has generally held in the information necessary for generate() to know how to generate data under the null. That's most clear to me in hypothesize(null = "independence") and hypothesize(null = "point", p = .2) (or multiple p's).

If we did something like hypothesize(null = "point", mu = 0), that's more specific than what's actually being done in generate(). For that reason I'd be inclined to stick to hypothesize(null = "independence"). Here's a first take at what that implementation could look like:

1. If you use hypothesize(null = "independence") without an explanatory variable, the output data frame has a new column called diff_direction or sign that has a shuffled vector of +1 and -1. It also outputs a warning or message that it added this column because there was no explanatory variable.

2. At this point, generate() would have everything it needed to know to proceed like normal with a permutation test on response and diff_direction. So options would be to a) accept generate(reps = 500, type = "permute") or b) require generate(reps = 500, type = "permute", paired = TRUE). There is precedent for the pedantic / superfluous 2b).

I'm not sure how pipelines get evaluated, but @simonpcouch probably knows: would it be possible/desirable to suppress the warning / message from 1 if followed by a generate() with paired = TRUE? That is, if the user indicates that they know what they're doing, we don't need to warn them, but if they haven't (say, if they just end the pipeline at hypothesize()), then we do.

Two downsides that I see to this approach:

1. There are no other examples where hypothesize() adds a column to the data frame.
2. The output of generate() would then presumably be three cols: replicate, response, and diff_direction. The operation of calculate() then isn't the same stat = "mean" as is normally used, since it'd have to multiply response * diff_direction before taking the mean.

This is a pretty verbose implementation of independence = "null". I can also imagine a much slimmer version that implements @mine-cetinkaya-rundel 's sample code and just adjusts the error from hypothesize() and uses paired = TRUE.

[simonpcouch]

I'm with it! A couple thoughts:

1. This meaning of null = "independence" is a departure from how we've used it in the past:

infer/R/hypothesize.R

Lines 94 to 98 in 5b7ea26

	if ((null == "independence") && !has_explanatory(x)) {
	abort(paste0(
	'Please `specify()` an explanatory and a response variable when ',
	'testing a null hypothesis of `"independence"`.'
	), call = call)

Does this feel okay? We can definitely adjust the conditions for that error, but it's been around in some form since the first release of infer, which maybe speaks to the significance of that change. I'd wonder if hypothesize(null = "point", mu = 0) (or otherwise) would feel more expressive.

2. Since that new type is in quotes, it doesn't need to be a valid column name (i.e. contain an underscore). We've used this rationale in the past to justify more plain-language argument values—do we want this to be "permute paired"? Or possibly just "permute" with a paired = logical(1L) argument?

[hardin47]

I agree with this (which is a better way of expressing my problem of setting mu to a specific number):

If we did something like hypothesize(null = "point", mu = 0), that's more specific than what's actually being done in generate(). For that reason I'd be inclined to stick to hypothesize(null = "independence").

@andrewpbray 's ideas seem good to me.

[simonpcouch]

would it be possible/desirable to suppress the warning / message from 1 if followed by a generate() with paired = TRUE?

It's definitely possible! We can raise that warning lazily depending on generate() input. I'd say that behavior is not R-idiomatic and may cause confusion for users developing pipelines interactively, but maybe is worth it if null = "independence" is more expressive. We can make sure that we hit the nail on the head with the warning message to prevent too much confusion.

I do think it's worth making sure there are no better options before overloading null = "independence" in this way. I'm concerned that:

• The input type for null = "independence" is no longer consistent, in that it used to require an explicit explanatory variable and now does not. Again, this requirement has been around since the first release of infer.
• The function injects code to raise a warning that could be raised from any number of functions. i.e. if a user just writesdf |> specify(response = diff) |> hypothesize(null = "independence"), the warning would be raised in hypothesize(). Appending generate(type = "permute", paired = TRUE) silences the warning. (Note that, if a user is writing their pipeline progressively, they're introducing a warning even when they're on the right track.) Instead appending generate(type = "permute") raises a warning in generate(), but generate() could also be any other function.
• hypothesize() is no longer type-stable (in that it sometimes adds a column). The conditions where it does so is also not visible in the call to hypothesize() itself—i.e. calling the function with the same arguments hypothesize(type = "independence") can result in two different things depending on what's piped in.

Would something like hypothesize(null = "paired") be more principled? I understand this is another intricacy to teach around, but a null = "independence" interface also feels pedagogically difficult.

[mine-cetinkaya-rundel]

Thanks @simonpcouch for outlining the potential drawbacks of null = "independence". I'm convinced that when specify(response = diff) this is not the way to go. Here are two options.

1. Option 1: New null type -- "paired independence". This is basically what you suggested above but emphasizing that we're doing some independence thing, but a special one. Then, proposed HT pipeline would be

df |>
  specify(response = diff) |>
  hypothesize(null = "paired independence") |>
  generate(1000, type = "permute", paired = TRUE) |>
  calculate(stat = "mean")

Pro: Addresses concerns re: making how null = "independence" works inconsistent.
Con: Need to specify pairing in two places (hypothesize() and generate()) and that seems excessive.

2. Option 2: (Clarifying why I suggested null = "independence" over null = "point", mu = 0 in the first place, but opening a whole other can of worms 🪱.)

We would use specify(response = diff) if the data is in wide format, which I think is customary for paired data, e.g., for repeated measurements or same subject being exposed to two treatments. In such cases either diff column is already computed in the raw data or it's trivial to compute it in a mutate() statement.

For non-paired two mean comparisons, infer assumes the data comes in long format, since that's almost always the case. Then, we do something like

df |>
  specify(response = outcome, explanatory = group) |>
  hypothesize(null = "independence") |>
  generate(1000, type = "permute") |>
  calculate(stat = "diff in means", order = c("treatment", "control")

If for paired comparisons we could assume the data come in long format as well, we could do something like

df |>
  specify(response = outcome, explanatory = group) |>
  hypothesize(null = "independence") |>
  generate(1000, type = "permute", paired = TRUE) |>
  calculate(stat = "diff in means", order = c("treatment", "control")

but this wouldn't be sufficient enough, we would also need to specify some id or similar to pair on.

Hence the can of worms...

Pro: Similarity in syntax/setup between paired and non-paired is elegant. Keeping data in long format makes the input data frame consistent across various inference procedures.
Con: Rarely does one get paired data in long format. The id (or similar) variable to pair on might be non-trivial to explain.

[simonpcouch]

No specific thoughts yet on Option 2, but a possibility for Option 1... we could leave out the paired argument in generate() since the mention of pairing appears already in hypothesize(). infer is "smart" about the generate() type in some ways already, so I'd feel okay with it knowing how to treat that null specially:

df |>
  specify(response = diff) |>
  hypothesize(null = "paired independence") |>
  generate(1000, type = "permute") |>
  calculate(stat = "mean")

[hardin47]

I don't like having to specify grouping = id which is what we'd have to do if the data are long. I'd prefer to have the data be set up as differences (wide). Although I can see a situation down the road where the long data flexibility is built in somehow.

I like the following option (from above), and it feels like we are converging on very good syntax!

df |>
  specify(response = diff) |>
  hypothesize(null = "paired independence") |>
  generate(1000, type = "permute") |>
  calculate(stat = "mean")