使用glmmTMB后,如何在R中获取Tukey事后检验的精确p值?
我在R中使用glmmTMB拟合截断负二项模型,分析不同子区域(Subarea)和年龄组(AgeClass)对数量(Quantity)的影响,模型代码如下:
m1 <- glmmTMB(Quantity ~ Subarea + AgeClass, family = truncated_nbinom2(link = "log"), data = df1_m1_52)
通过Deviance分析(Type II Wald卡方检验),Subarea和AgeClass的效应均显著:
glmmTMB:::Anova.glmmTMB(m1, type = c("II", "III", 2, 3), test.statistic = c("Chisq", "F"), component = "cond", vcov. = vcov(m1))
输出结果:
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Quantity
Chisq Df Pr(>Chisq)
Subarea 65.298 10 3.555e-10 ***
AgeClass 19.858 2 4.873e-05 ***
随后执行Tukey事后检验:
pairs(emmeans(m1, ~AgeClass, component="cond"), type="response", bias.adjust=F, adj="Tukey", infer=(TRUE))
输出结果:
contrast ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
A / SA 4.214 1.370 Inf 1.967 9.028 1 4.425 <.0001
A / Y 1.393 0.325 Inf 0.807 2.405 1 1.423 0.3291
SA / Y 0.331 0.119 Inf 0.142 0.767 1 -3.083 0.0058Results are averaged over the levels of: Subarea
Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 3 estimates
Intervals are back-transformed from the log scale
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log scale
现在需要获取这些检验结果的精确p值,而非<0.0001这类格式化后的结果,请问在R中如何设置实现?
可以通过调整数值打印选项或直接提取结果对象中的原始值来获取精确p值,分两种场景处理:
1. 针对Anova.glmmTMB的精确p值
- 方法一:调整全局打印精度
在运行Anova命令前,设置更高的数值显示精度:options(scipen = 999, digits = 10) glmmTMB:::Anova.glmmTMB(m1, type = "II", test.statistic = "Chisq", component = "cond", vcov. = vcov(m1))scipen=999关闭科学计数法,digits=10设置显示10位有效数字,可根据需求调整位数。
- 方法二:直接提取p值
将Anova结果保存为对象,直接访问p值列:anova_res <- glmmTMB:::Anova.glmmTMB(m1, type = "II", test.statistic = "Chisq", component = "cond", vcov. = vcov(m1)) # 查看精确p值 anova_res$`Pr(>Chisq)`
2. 针对pairs(emmeans)的精确p值
- 方法一:调整
emmeans打印选项
使用emm_options()设置打印精度:emm_options(digits = 10) pairs(emmeans(m1, ~AgeClass, component="cond"), type="response", bias.adjust=F, adj="Tukey", infer=(TRUE)) - 方法二:提取结果对象中的原始p值
把事后检验结果保存为对象,直接提取p值:
也可以转换为数据框查看完整数值:tukey_res <- pairs(emmeans(m1, ~AgeClass, component="cond"), type="response", bias.adjust=F, adj="Tukey", infer=(TRUE)) # 查看精确p值 tukey_res$p.valueas.data.frame(tukey_res)$p.value
内容的提问来源于stack exchange,提问作者YvonneS

