我在方差分析lsmeans中没有得到显著性差异,但有一些真正有意义的数据。我的脚本出了什么问题?
df <- structure(list(value = c(1.693086732, 0.25167691, 1.100272527, 1.60428654, 0.908237338, 1.449864567, 1.06604818, 0.596785144, 0.652925021, 0.453697295, 0.544252785, 1.464221767, 1.043720641, 0.735035158, 0.938875327, 0.712832947, 1.701854524, 1.021094251, 0.564349482, 2.326316679, 1.10170484, 1.075217638, 1.397442796, 0.501086703, 0.675502908, 0.846651623, 1.578086856, 1.857360967, 1.194232629, 1.875837087, 1.106882408, 1.112407609, 1.30479321, 0.637491754, 1.281566883, 1.103115742, 1.895286629, 1.623933836, 0.941989812, 1.30636425, 0.69977606, 1.937055334, 0.666069131, 0.829396619, 0.892844633, 0.573255443, 1.27370148, 0.531593222, 2.782899244, 0.972928201, 0.729463812, 1.121965821, 2.55117084, 0.999302442, 1.138902544, 1.656807007, 0.545349299, 0.550315908, 2.346074577, 0.637551271), gene = c("WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "WT", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox2", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox5", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox7", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9", "aox9"), time = c("0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t", "0t", "0t", "0t", "1t", "1t", "1t", "3t", "3t", "3t", "5t", "5t", "5t")), row.names = c(NA, -60L), class = c("data.table", "data.frame"))
我用过的脚本:
library(emmeans)
fit <- aov(value ~ gene*time, df)
summary(fit)
em <- emmeans(fit, ~ gene | time)
pairs(em)
pairs(em, adjust=NULL)
当我用Std Err绘制条形图中的数据时,我看到一些样本确实有意义,特别是在3天的时候。我已经在配对t检验中对此进行了测试,这些样本在p0.05处有显着差异。
这是配对测试,即使excel在时间点3给出了sig结果。
你能试着修正我的方差分析和最小二乘均值吗?
发布于 2018-08-16 11:10:47
不确定为什么你期望Tukey HSD成对比较是显着的;Anova中的主要影响都不是,在情节中也没有任何明显的东西。
为了完整性,即使主要效果不是很明显,你也可以使用Tukey HSD,但这种情况很少见,而且提前决定哪个测试与你的问题更相关是很重要的,而不是选择给出“更好”结果的测试。
summary(fit)
## Df Sum Sq Mean Sq F value Pr(>F)
## gene 4 0.873 0.2184 0.625 0.648
## time 3 1.670 0.5568 1.593 0.206
## gene:time 12 1.577 0.1314 0.376 0.965
## Residuals 40 13.984 0.3496
library(ggplot2)
sem <- summary(em)
sem$gene <- relevel(factor(sem$gene), 'WT')
ggplot(sem) + aes(gene, emmean, ymin=lower.CL, ymax=upper.CL) + geom_pointrange() +
facet_wrap(~time, nrow=1) + ggtitle("emmeans with 95% CIs")
该图以红色显示原始数据(这一点很重要)。
顺便说一句,没有一个成对t检验在.05水平上是显着的,即使没有对多次比较进行校正。就像在你的代码中一样,我每次都会成对地测试基因。
library(broom)
library(tidyr)
library(dplyr)
options(digits=2)
lapply(split(df, df$time), function(x) {
data.frame(time=x$time[1], tidy(pairwise.t.test(x$value, x$gene, p.adjust.method="none")))
}) %>% bind_rows() %>% spread(time, p.value)
## group1 group2 0t 1t 3t 5t
## 1 aox5 aox2 0.82 0.32 0.71 0.97
## 2 aox7 aox2 0.32 0.73 0.21 0.70
## 3 aox7 aox5 0.43 0.50 0.37 0.67
## 4 aox9 aox2 0.31 0.40 0.60 0.71
## 5 aox9 aox5 0.42 0.86 0.88 0.74
## 6 aox9 aox7 0.99 0.62 0.45 0.45
## 7 WT aox2 0.85 0.72 0.20 0.74
## 8 WT aox5 0.97 0.51 0.34 0.71
## 9 WT aox7 0.41 0.99 0.95 0.96
## 10 WT aox9 0.40 0.63 0.42 0.49
但是,如果在不合并方差的情况下运行此函数(或者再次运行时不校正多个比较),则确实会得到预期的结果,其中aox5和aox7的p值为3t = 0.016。
然而,由于这两个原因,这样做是错误的;首先,您确实应该纠正多次比较,其次,您没有足够的数据来很好地估计方差,因此池化很重要。您可以在原始数据中看到这一点;这两个数据恰好比其他数据更紧密,但这几乎可以肯定是随机机会造成的。
## group1 group2 0t 1t 3t 5t
## 1 aox5 aox2 0.70 0.25 0.794 0.96
## 2 aox7 aox2 0.17 0.73 0.413 0.61
## 3 aox7 aox5 0.32 0.49 0.016 0.53
## 4 aox9 aox2 0.46 0.52 0.744 0.79
## 5 aox9 aox5 0.56 0.89 0.868 0.80
## 6 aox9 aox7 0.99 0.71 0.428 0.59
## 7 WT aox2 0.82 0.65 0.398 0.70
## 8 WT aox5 0.97 0.36 0.096 0.65
## 9 WT aox7 0.41 0.99 0.892 0.95
## 10 WT aox9 0.57 0.69 0.409 0.63
所以最简短的答案是,带成对测试的ANOVA与你的成对t-test不匹配,因为ANOVA汇集了SD,但你的t-test没有,但更大的答案是,你应该同时汇集和纠正多次比较,所以你在开始时的代码,即:
pairs(em)
才是正确的选择。
https://stackoverflow.com/questions/51868951
复制相似问题