Search code examples
rp-valueemmeanstukey

Tukey post-hoc test not giving exact p-value

I am running a model in R to look at differences between subareas and age classes over many years.

enter image description here (snippet of what the data looks like)

This is my model: m1<-glmmTMB(Quantity ~ Subarea+ AgeClass, family=truncated_nbinom2(link = "log"), data=df1_m1_52)

After the results came back significant,

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 ***
---

I ran a Tukey's post-hoc test and got the following results:

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.0058

Results 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

Now my issue is that I need the exact p-value for these results. How do I tell R to give them to me rather than "<0.0001"?

Solution

  • This is likely just a function of how it prints outputs. You should be able to do something like this

    library(glmmTMB)

data(Owls)


m1 = glmmTMB(SiblingNegotiation ~ FoodTreatment + Nest + BroodSize, family = truncated_nbinom2(link = 'log'), zi = ~FoodTreatment, data = Owls)
#> dropping columns from rank-deficient conditional model: BroodSize

get_pvals = pairs(emmeans::emmeans(m1, ~FoodTreatment, component = 'cond'), type = 'response', bias.adjust = FALSE, adj = 'Tukey', infer = TRUE) 
  
get_pvals |> 
  broom::tidy()
#> # A tibble: 1 × 9
#>   term         contrast null.value ratio std.error    df  null statistic p.value
#>   <chr>        <chr>         <dbl> <dbl>     <dbl> <dbl> <dbl>     <dbl>   <dbl>
#> 1 FoodTreatme… Deprive…          0  1.27     0.102   Inf     1      2.97 0.00302

summary(get_pvals)$p.value
#> [1] 0.003019641

    Created on 2024-08-30 with reprex v2.1.1