rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
340 stars 30 forks source link

unadjusted p-value for pairwise emtrends contrasts #450

Closed samyogita-hardikar closed 8 months ago

samyogita-hardikar commented 8 months ago

Is there a way to retreive the raw/ unadjusted p-values for pairwise emtrends contrasts?
(That is it, really. But let me know if I should repost with reproducible example)

My model looks something like: model<-lmerTest::lmer(DV~ (1|sub) + (1|Family) + within_IV + between_IV1 + between_IV2 + within_IV:between_IV1 + within_IV:between_IV2 + age + gender)

I obtain the estimates for interaction effects by doing: interaction<-emtrends(model, pairwise ~ within_IV, var="Between_IV1", adjust="bonferroni", infer = TRUE)

Here my only options seem to be"bonferroni" or "tukey" (default), and the output never gives just the raw p-values

However, the adjustment is only done for the number of pairwise contrasts for one between_IV. If I want to be very conservative and correct for ALL tests I am doing, A) is it safe to assume that the "bonferroni" adjusted p-values are just multiplied by the number of pairwise comparisons, so I would just multiply that number accordingly to account for more than one between_IV that I'm testing? or B) Is it possible to retrieve the raw uncorrected p-values somehow?

samyogita-hardikar commented 8 months ago
library(emmeans)

###  create example data
### subject ID
sub<-rep(1:300, 4)

###within subjects factor 4 levels
within_IV<-rep(c("level1", "level2", "level3", "level4"), each=300) 

### dependent variable
DV<-as.numeric(rnorm(1200))  

### between-subjects continuous independent variables 1 and 2
between_IV1<-as.numeric(rep(rnorm(300),4))  
between_IV2<-as.numeric(rep(rnorm(300),4)) 

### age and gender as covariates of no interest
age<-as.numeric(rep(rnorm(300),4)) 
gender<-rep(sample(c("m", "f"), size=300, replace = T), 4) 

### data.frame
df<-as.data.frame(cbind(sub, within_IV, DV, as.numeric(between_IV1), as.numeric(between_IV2), 
                        as.numeric(age), gender))

###  model 
mod<-lmerTest::lmer(as.numeric(DV)~ (1|sub)+ within_IV + as.numeric(age) + gender + 
                       between_IV1 + between_IV2+
                      within_IV:between_IV1 + within_IV:between_IV2, data = df)

###  followup tests for interaction between within_IV and between_IV1
emm_options( lmer.df = "satterthwaite")
interaction_bonf<-emtrends(mod, pairwise ~ within_IV, var="between_IV1", adjust="bonferroni", infer = TRUE)
interaction_no_correction<-emtrends(mod, pairwise ~ within_IV, var="between_IV1",  infer = TRUE)

In interaction_bonf, "bonferroni" adjusts for 6 comparisons. But if I do followups for both between_IVs, that's 12 comparisons. In interaction_no_correction, the Tukey method is applied by default.

My question is, if I want to correct for pairwise tests done for both between_IV1 and between_IV2 (2*6=12 pairwise comparisons), does it make sense to simply take the "bonferroni" corrected p-values from interaction_bonf, and multiply them by 2? Or, alternatively, is there a way to retrieve the raw "interaction" p-values without any adjustment?

rvlenth commented 8 months ago

The answer to your initial question is well documented: specify adjust = "none".

For the rest,. I'll go step-by-step. For one set of comparisons:

WT <- emtrends(mod,  ~ within_IV, var="between_IV1")
WCMP <- pairs(WT)

I separated out the pairwise stuff because that call produces a list of two results, and it's confusing to understand what arguments go with which.

Simlarly, for the other:

BT <- emtrends(mod,  ~between_IV, var="between_IV1")
BCMP <- pairs(BT)

Now combine them:

summary(rbind(WT, BT), infer = TRUE, adjust = "bonf")
## or just:
WT + BT

There are other adjust possibilities such as "none", "mvt", "sidak", "FDR", etc. But "tukey" will not be allowed because it can only be used with one set of pairwise comparisons.