r - 从 data.frame 迭代计算 F 统计数据，特定列具有要比较的类别

Question

我有一个过于简化的 data.frame，看起来像这样(真正的 data.frame 在“类”列中有 > 10 个类和 > 1000 行):

Bin    Class       Var      n
0.1    benign      0.04     15
0.1    damaging    0.3      14
0.1    all         0.0006   16
0.2    benign      0.1      13
0.2    damaging    0.04     16
0.2    all         0.03     10
0.3    benign      0.07     8
0.3    damaging    0.06     12
0.3    all         0.1      10

对于“Bin”中的每个值，我想使用相应的方差(“Var”列)和样本大小计算“所有”、“良性”和“破坏性”(“Class”列)之间的 F 统计量(“n”列)。作为输出，对于“所有与良性”和“所有与破坏性”比较，我会得到一个观察到的 F 统计量 (Obs_F)、一个预期的 F 统计量 (Exp_F) 和一个 p 值。

所有与良性的示例，Bin“0.1”，以及相应的公式:

Obs_F = 0.04/0.0006 # higher Var/lower Var
Exp_F = qf(.95, df1= 15 , df2 = 16) # df1 and df2 = "n" of higher and lower Var, respectively
p-value = pf(Obs_F, df1= 15 , df2 = 16 ) # df1 and df2 = "n" of higher and lower Var, respectively

# I suspect using ifelse() function is a good way of sorting the highest vs lower Var and degrees of freedom (df1 and df2). 

我希望获得如下所示的输出 data.frame:

Bin  Comparison        Obs_F    Exp_F   p-value
0.1  all_vs_benign      …        …        …
0.1  all_vs_damaging    …        …        …
0.2  all_vs_benign      …        …        …
0.2  all_vs_damaging    …        …        …
0.3  all_vs_benign      …        …        …
0.3  all_vs_damaging    …        …        …

我尝试过使用 dplyr、聚合和类似函数，但到目前为止，我一次只能计算 F stats 1。

Answer 1

这是一个基于 R 的想法，使用 Map 和 mapply。首先，我们需要创建一个函数来返回您想要的结果。我们在 Class 上拆分原始数据框并创建第二个函数 (fun2) 以 Map 第一个 fun1 到那个名单。然后下一步是创建一个包含所有兴趣组合的矩阵。最后，使用 mapply 将函数应用于矩阵。

fun1 <- function(d1, d2){
  Obs_F <- pmax(d1$Var, d2$Var)/pmin(d1$Var, d2$Var)
  dd <- rbind(d1, d2)
  n_min <- dd$n[dd$Var == pmin(d1$Var, d2$Var)]
  n_max <- dd$n[dd$Var == pmax(d1$Var, d2$Var)]
  Exp_F <- qf(.95, df1= n_min, df2 = n_max)
  p_value <- pf(Obs_F, df1= n_min, df2 = n_max)
  return(data.frame(Obs_F, Exp_F, p_value, stringsAsFactors = FALSE))
}

l1 <- split(df, df$Class)
fun2 <- function(x, y){ Map(fun1, l1[x], l1[y])}

m1 <- combn(1:length(l1), 2)[,1:length(l1)-1] 

final_list <- mapply(fun2, m1[1,], m1[2,])

#tidy up to required data frame
final_df <- do.call(rbind, c(final_list, make.row.names = FALSE))
final_df$Bin <- rep(unique(df$Bin), 2) 
final_df <- final_df[order(final_df$Bin),]
final_df$Comparison <- rep(c('all_vs_benign', 'all_vs_damaging'), length(unique(df$Bin)))

final_df
#       Obs_F    Exp_F   p_value Bin      Comparison
#1  66.666667 2.352223 1.0000000 0.1   all_vs_benign
#4 500.000000 2.373318 1.0000000 0.1 all_vs_damaging
#2   3.333333 2.671024 0.9772730 0.2   all_vs_benign
#5   1.333333 2.493513 0.7067062 0.2 all_vs_damaging
#3   1.428571 3.071658 0.7068978 0.3   all_vs_benign
#6   1.666667 2.753387 0.8009820 0.3 all_vs_damaging