与仅使用 lm() 相比,使用 group_by 后跟 lm() 时的不同 P.values
Different P.values when using group_by followed by lm() compared to just lm() only
我可以在以下方面获得一些帮助吗?我有一个数据框,它有多个组,我想 运行 一个线性模型。作为测试,我仅对其中一个组和 运行 函数 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() 函数时相同。
谢谢
我用 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
正如我在评论中提到的,问题是您没有使用过滤后的数据,而是使用了整个数据集。因此不匹配。
下面是带有示例数据的代码,在其上使用 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() 进行子集化,并得到以下输出:
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() 函数时相同。
谢谢
我用 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
正如我在评论中提到的,问题是您没有使用过滤后的数据,而是使用了整个数据集。因此不匹配。
下面是带有示例数据的代码,在其上使用 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)