What happens if we permute after fitting?

permutation_test.R

library(tidyverse)

# Arguments: Sigma of old predictor, Sigma of new predictor, N cases
compare_methods = function(s0, s1, n){
  print(c(s0, s1, n))
  y = rnorm(n, 0, 1)
  x0 = y + rnorm(n, 0, s0) # Old predictor
  x1 = y + rnorm(n, 0, s1) # New predictor to evaluate

  m0 = lm(y ~ x0)
  m1 = lm(y ~ x0 + x1)
  ss1  = sum((predict(m1)-y)^2)
  true_p = anova(m0, m1)$`Pr(>F)`[[2]] # Ground truth p-value
  
  # Permuate data multiple times are refit model ("permute.before")
  permutation_ss1s_a = map_dbl(1:1000, function(i){
    x1.surrogate = sample(x1, n, replace=T)
    m1.surrogate = lm(y ~ x0 + x1.surrogate) 
    ss1  = sum((predict(m1.surrogate)-y)^2)
    ss1
  })
  pval1 = mean(ss1 > permutation_ss1s_a)
  
  # Fit model once, then permuate data multiple times ("permuate.after")
  permutation_ss1s_b = map_dbl(1:1000, function(i){
    x1.surrogate = sample(x1, n, replace=T)
    preds = predict(m1, newdata = list(x1=x1.surrogate))
    ss  = sum((preds-y)^2)
    ss
  })
  pval2 = mean(ss1 > permutation_ss1s_b)
  return(data.frame(s0, s1, n, true_p, permute.before=pval1, permute.after=pval2))
}

# Adjust these values to try different parameters
s0 = c(1)
s1 = c(2)
n = c(10)
n.rep = 1:10

df = expand.grid(n.rep=n.rep, s0=s0, s1=s1, n=n) %>% select(-n.rep)
result = pmap(df, compare_methods)
result = do.call(rbind, result)

# Plot results
result %>%
  gather(Method, Estimate, permute.before, permute.after) %>%
  ggplot(aes(true_p, Estimate, color=Method)) +
  geom_point() +
  geom_abline(intercept=0, slope=1, linetype='dotted') +
  theme_classic() +
  coord_fixed() +
  labs(x='True p value')
# Result: https://imgur.com/K957F9C