r - 基于置换法的方差分析中 F 统计量的 Monte Carlo 估计

Question

我正在尝试使用组标识符 g 在 (y₁,...,y_N) 上对 ANOVA 进行排列检验.我应该使用 (1)/(g-1) (muhat_j - muhat)^2 的总和作为测试统计量，而 muhat_j 是第 j 组样本均值，和 muhat=(1/g)summation muhat_j.

## data
y <- c(6.59491, 6.564573, 6.696147, 6.321552, 6.588449, 6.853832, 
6.370895, 6.441823, 6.227591, 6.675492, 6.255462, 6.919716, 6.837458, 
6.41374, 6.543782, 6.562947, 6.570343, 6.993634, 6.666261, 7.082319, 
7.210933, 6.547977, 6.330553, 6.309289, 6.913492, 6.597188, 6.247285, 
6.644366, 6.534671, 6.885325, 6.577568, 6.499041, 6.827574, 6.198853, 
6.965038, 6.58837, 6.498529, 6.449476, 6.544842, 6.496817, 6.499526, 
6.709674, 6.946934, 6.23884, 6.517018, 6.206692, 6.491935, 6.039925, 
6.166948, 6.160605, 6.428338, 6.564948, 6.446658, 6.566979, 7.17546, 
6.45031, 6.612242, 6.559798, 6.568082, 6.44193, 6.295211, 6.446384, 
6.658321, 6.369639, 6.066747, 6.345537, 6.727513, 6.677873, 6.889841, 
6.724438, 6.379956, 6.380779, 6.50096, 6.676555, 6.463236, 6.239091, 
6.797642, 6.608025)

## group
g <- structure(c(2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 2L, 
3L, 2L, 3L, 2L, 3L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 1L, 3L, 1L, 2L, 2L, 1L, 3L, 2L, 2L, 3L, 1L, 2L, 2L, 2L, 1L, 
2L), .Label = c("B1", "B2", "B3"), class = "factor")

这是我现在拥有的，但是当我更改它以测试样本均值而不是 F 统计量时，它不起作用。我很确定我需要将 T.obs 和 T.perm 更改为类似于 by(y, g, mean) 的内容但我认为我还缺少更多。

n <- length(y) #sample size n
T.obs<- anova(lm(y ~ g))$F[1]   #Observed statistic 
n.perm <- 2000   # we will do 2000 permutations
T.perm <- rep(NA, n.perm)   #A vector to save permutated statistic
for(i in 1:n.perm) {
  y.perm <- sample(y, n, replace=F)   #permute data
  T.perm[i] <- anova(lm(y.perm ~ g))$F[1]   #Permuted statistic
  }
mean(T.perm >= T.obs)   #p-value

Answer 1

我真的不知道“它不工作”是什么意思。据我所知，它工作正常，只是有点慢。

set.seed(0)
n <- length(y) #sample size n
T.obs <- anova(lm(y ~ g))$F[1]   #Observed statistic 
n.perm <- 2000   # we will do 2000 permutations
T.perm <- rep(NA, n.perm)   #A vector to save permutated statistic
for(i in 1:n.perm) {
  y.perm <- sample(y, n, replace=F)   #permute data
  T.perm[i] <- anova(lm(y.perm ~ g))$F[1]   #Permuted statistic
  }
mean(T.perm >= T.obs)
# [1] 0.4915

这与理论值相当接近

anova(lm(y ~ g))$Pr[1]
# [1] 0.4823429

所以，是的，你做的都是正确的!

---

从您问题的第一段来看，我们似乎想自己计算 F 统计量，所以下面的函数就是这样做的。有一个开关 "use_lm"。如果设置 TRUE，它使用 anova(lm(y ~ g)) 作为您在原始代码中所做的。 此函数旨在使 F 统计量和 p 值的计算变得透明。此外，手动计算比调用 lm 和 anova 快 15 倍(这是显而易见的...)。

fstat <- function (y, g, use_lm = FALSE) {
  if (!use_lm) {
    ## group mean (like we are fitting a linear model A: `y ~ g`)
    mu_g <- ave(y, g, FUN = mean)
    ## overall mean (like we are fitting a linear model B: `y ~ 1`)
    mu <- mean(y)
    ## RSS (residual sum of squares) for model A
    RSS_A <- drop(crossprod(y - mu_g))
    ## RSS (residual sum of squares) for model B
    RSS_B <- drop(crossprod(y - mu))
    ## increase of RSS from model A to model B
    RSS_inc <- RSS_B - RSS_A
    ## note, according to "partition of squares", we can also compute `RSS_inc` as
    ## RSS_inc <- drop(crossprod(mu_g - mu))
    ## `sigma2` (estimated residual variance) of model A
    sigma2 <- RSS_A / (length(y) - nlevels(g))
    ## F-statistic
    fstatistic <- ( RSS_inc / (nlevels(g) - 1) ) / sigma2
    ## p-value
    pval <- pf(fstatistic, nlevels(g) - 1, length(y) - nlevels(g), lower.tail = FALSE)
    ## retern
    return(c(F = fstatistic, pval = pval))
    }
  else {
    anovalm <- anova(lm(y ~ g))
    return(c(F = anovalm$F[1L], pval = anovalm$Pr[1L]))
    }
  }

让我们先检查一下这个函数的有效性:

F_obs <- fstat(y, g)
#        F      pval 
#0.7362340 0.4823429 

F_obs <- fstat(y, g, TRUE)
#        F      pval 
#0.7362340 0.4823429

不要因为它微不足道而感到惊讶。您的数据并没有真正表明存在显着的群体差异。看看箱线图:

boxplot(y ~ g)    ## or use "factor" method of `plot` function: `plot(g, y)`

现在我们继续排列。我们为此编写了另一个函数 perm。这实际上很容易，因为我们有一个很好定义的 fstat。我们需要做的就是使用replicate 来包装sample + fstat。

lm 实际上很慢:

library(microbenchmark)
microbenchmark(fstat(y, g), fstat(y, g, TRUE), times = 200)

#Unit: microseconds
#              expr     min      lq      mean  median      uq      max neval cld
#       fstat(y, g)  228.44  235.32  272.1204  275.34  290.20   388.84   200  a 
# fstat(y, g, TRUE) 4090.00 4136.72 4424.0470 4181.02 4450.12 16460.72   200   b

所以我们使用 f(..., use_lm = FALSE) 编写此函数:

perm <- function (y, g, n) replicate(n, fstat(sample(y), g)[[1L]])

现在让我们用 n = 2000 来运行它(为再现性设置随机种子):

set.seed(0)
F_perm <- perm(y, g, 2000)

## estimated p-value based on permutation
mean(F_perm > F_obs[[1L]])
# [1] 0.4915

请注意它与理论 p 值的接近程度:

F_obs[[2L]]
# [1] 0.4823429

如您所见，结果与您的原始代码一致。