Post-hoc test (Tukey) for a robust model (lmrob) with interactions in R

[adriananas]

I am trying to perform post-hoc tests (Tukey) on a robust model with R. The model was found empirically with lmrob since the data violates some normality assumptions and looks like this

my_model <- lmrob(Y ~ A + B + C + A:B + B:C)

Y is a continuous variable, whereas A, B, and C are categorical data. A has 3 levels, B and C have 2 levels each.

In one of my previous projects, I have conducted the post-hoc tests (pairwise comparisons of the estimated means) according to the following procedure (see pages 8ff.): https://cran.r-project.org/web/packages/multcomp/vignettes/multcomp-examples.pdf

There I had the following model, with similar data (Y being continuous; A and B categorical)

old_model <- lmrob(Y ~ A * B)

And I performed the following procedure for the pairwise comparisons:

# First, we compute the means of the response for all combinations of the levels of A and B:
tmp <- expand.grid(A = levels(data$A), B = levels(data$B))
    
# and construct a model matrix
X <- model.matrix(~ A * B, data = tmp)
estimated_means = cbind(tmp, "estimated marginal mean" = predict(old_model, newdata = tmp))

# We now construct a contrast matrix based on Tukey-contrasts for A 
# in a block-diagonal way, i.e., for each level of B

Tukey <- contrMat(table(data$A), "Tukey")

K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = 3*ncol(Tukey)))
rownames(K1) <- paste(levels(data$B)[1], rownames(K1), sep = ":")

K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), 
            Tukey, matrix(0, nrow = nrow(Tukey), ncol = 2*ncol(Tukey)))
rownames(K2) <- paste(levels(data$B)[2], rownames(K2), sep = ":")

K3 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = 2*ncol(Tukey)), 
            Tukey, matrix(0, nrow = nrow(Tukey), ncol = 1*ncol(Tukey)))
rownames(K3) <- paste(levels(data$B)[3], rownames(K3), sep = ":")

K4 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = 3*ncol(Tukey)),
            Tukey)
rownames(K4) <- paste(levels(data$B)[4], rownames(K4), sep = ":")

K <- rbind(K1, K2, K3, K4)
colnames(K) <- c(rep(colnames(Tukey), 4))

# and perform the tests via
summary(glht(old_model, linfct = K %*% X), test = adjusted(type = "holm"))

Then the same procedure was applied vice versa with a Tukey-contrast matrix in B for each level of A.

My question is, how can I conduct the analysis similarly to my previous project? I am not sure how to deal with the interaction terms in the current model and how to construct the model matrix due to the interaction terms.

Any advice would be highly appreciated, thank you.