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

Question

我可以在以下方面获得一些帮助吗？我有一个数据框，它有多个组，我想运行一个线性模型。作为测试，我仅对其中一个组和运行函数 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.value 的减弱？我认为 p.value 应该与我对轨迹进行子集化和运行 lm() 函数时相同。

谢谢

Answer 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

Answer 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)

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

Different P.values when using group_by followed by lm() compared to just lm() only

r

group-by

statistics

lm