与仅使用 lm() 相比，使用 group_by 后跟 lm() 时的不同 P.values答案

【问题标题】：Different P.values when using group_by followed by lm() compared to just lm() only与仅使用 lm() 相比，使用 group_by 后跟 lm() 时的不同 P.values
【发布时间】：2017-11-18 20:07:43
【问题描述】：

我可以就以下问题获得一些帮助吗？我有一个数据框，它有多个组，我想在上面运行线性模型。作为测试，我只对其中一个组进行子集化并运行函数 lm() 并得到以下输出：

test <- filter(dat, locus == "ChrX_1")
test.result <- lm(methylation ~ Pheno, dat)

              term estimate  std.error statistic    p.value
1 (Intercept)   56.955      0.9729203 58.540254  9.080525e-250
2      Pheno1    9.015      1.1915791  7.565591  1.464884e-13

然后我使用 dplyr 包中的 group_by 对不同的组执行 lm() 函数。但是轨迹 "ChrX_1" 的 p.value 的输出现在不同并且更弱了。

test.result4 <- group_by(dat, locus) %>%
  do(model.test2 = lm(methylation ~ Pheno, data = .))
tidy(test.result4, model.test2)  

    locus        term estimate std.error statistic      p.value
    <chr>       <chr>    <dbl>     <dbl>     <dbl>        <dbl>
1   ChrX_1 (Intercept)    59.40  4.476666 13.268804 1.342225e-13
2   ChrX_1      Pheno1     9.05  5.482773  1.650624 1.099895e-01
3  ChrX_10 (Intercept)    59.00  4.069398 14.498459 1.522725e-14
4  ChrX_10      Pheno1    11.40  4.983974  2.287331 2.993721e-02
5  ChrX_11 (Intercept)    58.90  4.665565 12.624408 4.460131e-13
6  ChrX_11      Pheno1     9.10  5.714127  1.592544 1.224905e-01
7  ChrX_12 (Intercept)    52.80  3.717022 14.204921 2.526739e-14
8  ChrX_12      Pheno1    10.65  4.552403  2.339424 2.667444e-02
9  ChrX_13 (Intercept)    53.10  3.556734 14.929427 7.343091e-15
10 ChrX_13      Pheno1     7.10  4.356092  1.629901 1.143224e-01
# ... with 30 more rows

因此，我想知道是什么导致 p.values 减弱？我认为 p.value 应该与我对轨迹进行子集化并在其上运行 lm() 函数时相同。

谢谢

【问题讨论】：

当您尝试将lm 用于locus == "ChrX_1" 时，您没有使用test 数据集，而是使用dat。是这个问题吗？
我怀疑您将整个数据集用于第一个 lm。你需要使用test 而不是dat
感谢您指出。不知何故，我完全错过了。

标签： r group-by statistics lm

【解决方案1】：

我用 iris 尝试过，两种方法的结果都是一样的。您的 group_by() 行有问题。试试我的方法。

看：

test <- filter(iris, Species=="setosa")
test.lm <- lm(Sepal.Length ~Sepal.Width, data=test)

  Species        term  estimate  std.error statistic      p.value
   <fctr>       <chr>     <dbl>      <dbl>     <dbl>        <dbl>
1  setosa (Intercept) 2.6390012 0.31001431  8.512514 3.742438e-11
2  setosa Sepal.Width 0.6904897 0.08989888  7.680738 6.709843e-10

然后用 group_by()

iris %>% group_by(Species) %>% do(tidy(lm(Sepal.Length~Sepal.Width, data=.)))

 Species        term  estimate  std.error statistic      p.value
      <fctr>       <chr>     <dbl>      <dbl>     <dbl>        <dbl>
1     setosa (Intercept) 2.6390012 0.31001431  8.512514 3.742438e-11
2     setosa Sepal.Width 0.6904897 0.08989888  7.680738 6.709843e-10
3 versicolor (Intercept) 3.5397347 0.56287357  6.288685 9.069049e-08
4 versicolor Sepal.Width 0.8650777 0.20193757  4.283887 8.771860e-05
5  virginica (Intercept) 3.9068365 0.75706053  5.160534 4.656345e-06
6  virginica Sepal.Width 0.9015345 0.25310551  3.561892 8.434625e-04

【讨论】：

【解决方案2】：

正如我在评论中提到的，问题在于您没有使用过滤后的数据，而是使用了整个数据集。因此错配。

下面是带有示例数据的代码，显示在其上使用 group_by 和 lm 时没有不匹配。

library(dplyr)
library(tidyr)
library(broom)

set.seed(123)
dat <- data.frame(methylation=runif(1000, min=10, max=200), 
  Pheno=runif(1000, min=10, max=200), 
  locus=sample(paste0("ChrX_", 1:10), 1000, replace=TRUE)
  )
dat$locus <- as.character(dat$locus)

test <- filter(dat, locus == "ChrX_1")
test.result <- lm(methylation ~ Pheno, test)
summary(test.result)

test.result4 <- group_by(dat, locus) %>%
  do(model.test2 = lm(methylation ~ Pheno, data = .))
tidy(test.result4, model.test2)

【讨论】：