我有一个数据集,其中测量了多年和多种作物的表型。我寻找关系的第一种方法是双向方差分析,我问“在考虑了抽样结构后,多样性对表型是否有影响?”即
表型〜采样事件+品种
我真正感兴趣的问题是,是否有任何品种明显优于或差于总体人口平均值。
鉴于我发现表型和品种之间存在显着关系,我一直在努力寻找合适的事后测试。 TukeyHSD 对我不起作用,因为我对不同品种之间的比较不感兴趣,而是与平均值比较。
作为一种解决方法,我一直在使用 T 检验,并针对多次测试调整 p 值。但是,我不确定这是否a)有效,b)是否可以在R中以更简化的方式完成。目前我必须手动创建一个新列,其中我感兴趣的品种是“a”并且其他一切都是“b”,因此 T 检验将样本“a”与平均值进行比较。然而,考虑到我有 30 个不同的品种和 8 个不同的表型,这将需要很长时间手动完成!
非常欢迎任何帮助!我附上了一些示例数据。为了简单起见,我只包含了一种表型
variety <- c(A,A,A,A,B,B,B,B,C,C,C,C,D,D,D,D,E,E,E,E,F,F,F,F,G,G,G,G,H,H,H,H)
sampling_event <- c(A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2,A1,A1,A2,A2)
phenotype1 <- c(13.86,14.48,15.5,16.22,14.72,14.72,16.6,16.98,16.98,12.34,15.6,17.82,17.6,9.26,13.46,12.24,13.1,16.22,15.94,10.86,12.44,10.58,17.3,13.38,15.2,13.66,18.2,14.9,15.68,18.8,15.94,13.38)
dummydata <-data.frame(variety, sampling_event, phenotype1)
我按照 https://www.datanovia.com/en/blog/how-to-perform-multiple-t-test-in-r-for- Different-variables/ 上的说明进行操作。由此,我设法为每个品种的所有表型生成 T 检验,这很好,但我宁愿对每个表型的所有品种进行 T 检验。
我不确定这是否是针对您的特定问题的最佳统计数据,但您是否考虑过使用估计边际平均值来寻找组之间的差异?
library(emmeans)
dummydata <- data.frame(
variety = c('A','A','A','A',
'B','B','B','B',
'C','C','C','C',
'D','D','D','D',
'E','E','E','E',
'F','F','F','F',
'G','G','G','G',
'H','H','H','H'),
sampling_event = c('A1','A1','A2','A2','A1','A1','A2','A2','A1','A1','A2','A2',
'A1','A1','A2','A2','A1','A1','A2','A2','A1','A1','A2','A2','A1',
'A1','A2','A2','A1','A1','A2','A2'),
phenotype1 <- c(13.86,14.48,15.5,16.22,14.72,14.72,16.6,16.98,16.98,12.34,15.6,17.82,17.6,9.26,13.46,12.24,13.1,16.22,15.94,10.86,12.44,10.58,17.3,13.38,15.2,13.66,18.2,14.9,15.68,18.8,15.94,13.38)
)
# Create linear model (setup based on what was asked in post)
lm_mod <- lm(phenotype1 ~ sampling_event + variety, data = dummydata)
# Estimated marginal means
emm_lm_mod <- emmeans(lm_mod, pairwise ~ sampling_event | variety)
emm_lm_mod
#> $emmeans
#> variety = A:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 14.6 1.23 23 12.0 17.1
#> A2 15.5 1.23 23 12.9 18.0
#>
#> variety = B:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 15.3 1.23 23 12.8 17.8
#> A2 16.2 1.23 23 13.7 18.8
#>
#> variety = C:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 15.2 1.23 23 12.7 17.8
#> A2 16.1 1.23 23 13.6 18.7
#>
#> variety = D:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 12.7 1.23 23 10.1 15.2
#> A2 13.6 1.23 23 11.1 16.1
#>
#> variety = E:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 13.6 1.23 23 11.0 16.1
#> A2 14.5 1.23 23 11.9 17.0
#>
#> variety = F:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 13.0 1.23 23 10.4 15.5
#> A2 13.9 1.23 23 11.3 16.4
#>
#> variety = G:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 15.0 1.23 23 12.5 17.6
#> A2 15.9 1.23 23 13.4 18.5
#>
#> variety = H:
#> sampling_event emmean SE df lower.CL upper.CL
#> A1 15.5 1.23 23 12.9 18.0
#> A2 16.4 1.23 23 13.9 19.0
#>
#> Confidence level used: 0.95
#>
#> $contrasts
#> variety = A:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = B:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = C:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = D:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = E:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = F:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = G:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
#>
#> variety = H:
#> contrast estimate SE df t.ratio p.value
#> A1 - A2 -0.917 0.82 23 -1.120 0.2745
创建于 2024-05-27,使用 reprex v2.1.0