Как производить Tukey CLD, используя lmer с несколькими факторами

Question

Я пытаюсь проанализировать сельскохозяйственное испытание и хотел бы выполнить многократное сравнение, аналогичное HSD.test. Я предоставил пример кода, демонстрирующий, чего я хотел бы достичь с помощью lmer.

Я пробовал использовать emmeans и cld/glht, но могу получить только половину того, что мне нужно.

library(agricolae)
library(lme4)
library(multcomp)

df <- data.frame(iris,Petal.Colour=c("red","blue"), Country=c("UK","France","Spain"))

#this is what I would usually do to produce means with Tukey letter display

lm.model <- lm(data=df, Petal.Length~ Species + Country)
anova(lm.model)
aov.model <- aov(lm.model)
HSD.test(aov.model, trt = c("Species","Country"), console=T)

#my agricultural trial analysis uses lme4 to control for rows and columns (represented here by 1|Petal.Colour for ease)

lmer.model <- lmer(data=df, Petal.Length ~ Species + Country + (1|Petal.Colour))
anova(lmer.model)

#using emmeans I can produce multiple comparisons with associated p values

emmeans(lmer.model, list(pairwise ~ Species*Country), adjust = "tukey")

#using cld and glht I can get letters but only for one factor

cld(glht(lmer.model,linfct=mcp(Species ="Tukey")))
#or
cld(glht(lmer.model,linfct=mcp(Country ="Tukey")))

#how can I combine these??

Любая помощь с благодарностью получена!