rvlenth / emmeans

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

emtrends in a model with two interacting categorical predictors and one (or more) numeric #423

Closed samyogita-hardikar closed 1 year ago

samyogita-hardikar commented 1 year ago

Hello! I was wondering how to get appropriate emtrends in a model that has 2 interacting categorical predictors and a numeric predictor.

Reproduced below using the obk.long dataset

library(lmerTest)
library(emmeans)

data("OBrienKaiser", package = "carData")
set.seed(1)
OBrienKaiser2 <- within(OBrienKaiser, {
  id <- factor(1:nrow(OBrienKaiser))
  age <- scale(sample(18:35, nrow(OBrienKaiser), replace = TRUE), scale = FALSE)})
attributes(OBrienKaiser2$age) <- NULL # needed or resahpe2::melt throws an error.
OBrienKaiser2$age <- as.numeric(OBrienKaiser2$age)
obk.long <- reshape2::melt(OBrienKaiser2, id.vars = c("id", "treatment", "gender", "age"))
obk.long[,c("phase", "hour")] <- lapply(as.data.frame(do.call(rbind,
                                                              strsplit(as.character(obk.long$variable), "\\."),)), factor)
obk.long <- obk.long[,c("id", "treatment", "gender", "age", "phase", "hour", "value")]
obk.long <- obk.long[order(obk.long$id),]
rownames(obk.long) <- NULL 

My data and model are comparable to the following where "phase" is a within subject variable, and "age" and "gender" are between.

model<-lmer(value~ (1|id) + phase +
                 phase:age + gender  , data = obk.long)

I am interested in the effect of age at every level of phase but also for both genders. The first part is achieved by doing

emtrend_age_phase<-emtrends(model, ~phase,  var="age", adjust = 'bonferroni', infer = TRUE)

emtrend_age_phase

phase age.trend SE df lower.CL upper.CL t.ratio p.value fup 0.1592 0.0845 15.9 -0.0666 0.385 1.885 0.2333 post 0.1810 0.0845 15.9 -0.0449 0.407 2.143 0.1437 pre 0.0254 0.0845 15.9 -0.2005 0.251 0.300 1.0000

Results are averaged over the levels of: gender Degrees-of-freedom method: satterthwaite Confidence level used: 0.95 Conf-level adjustment: bonferroni method for 3 estimates P value adjustment: bonferroni method for 3 tests

but after adding gender as an interaction to emtrends ...

emtrend_age_phase_gender<-emtrends(model, ~phase*gender,  var="age", adjust = 'bonferroni', infer = TRUE)

emtrend_age_phase_gender

phase gender age.trend SE df lower.CL upper.CL t.ratio p.value fup F 0.1592 0.0845 15.9 -0.0950 0.413 1.885 0.4667 post F 0.1810 0.0845 15.9 -0.0732 0.435 2.143 0.2873 pre F 0.0254 0.0845 15.9 -0.2288 0.280 0.300 1.0000 fup M 0.1592 0.0845 15.9 -0.0950 0.413 1.885 0.4667 post M 0.1810 0.0845 15.9 -0.0732 0.435 2.143 0.2873 pre M 0.0254 0.0845 15.9 -0.2288 0.280 0.300 1.0000

Degrees-of-freedom method: satterthwaite Confidence level used: 0.95 Conf-level adjustment: bonferroni method for 6 estimates P value adjustment: bonferroni method for 6 tests

even though 6 "age.trend" rows are generated, only the pvalues are adjusted to reflect 6 tests instead of 3. But the actual age.trend is still only given for each level of phase, while averaging over two levels of gender.

What would be the solution here?

rvlenth commented 1 year ago

What you see is a consequence of the model you fitted. In that model, gender does not interact with the other factors; in so doing, you have specified that the effects of phase, age, and phase:age are the same for both genders.

samyogita-hardikar commented 1 year ago

Aah, I see. Thanks for the quick help!