r - How to obtain Tukey compact letter display from a GLM with interactions -

i have set of data i've analyzed generalized linear model has 3 categorical factors in 3-way interaction (factora, factorb, factorc) , fourth continuous factor (factord) added in model. trying obtain set of tukey letter groups (ie, compact letter display) model haven't found way include interaction successfully. i'm not interested in including factord, 3 in interaction.

i have gotten tukey-adjusted pairwise comparisons this:

lsmeans(my.glm, factora*factorb*factorc) 

but not able figure out how produce compact letters display that. can done multcomp package find ways main effects package, not interactions.

so tried agricolae package, post (https://stats.stackexchange.com/questions/31547/how-to-obtain-the-results-of-a-tukey-hsd-post-hoc-test-in-a-table-showing-groupe) discusses that should work. however, following instructions in answer led non-functional response hsd.test. specifically, main effects tests work fine, e.g. hsd.test(my.glm,"factora") not interactions work. tried this:

intxns<-with(my.data, interaction(factora,factorb,factorc)) hsd.test(my.glm,"intxns",group=true) 

but error indicates hsd.test function didn't recognize "intxns" valid object, looks (also, checked intxns object , looks , number of rows matched number of residuals of glm):

name: inxtns  factora factorb factorc factord 

i same error if put nonsense factor field in hsd.test function call. checked inxtns object , looks , number of rows matched number of residu agricolae notes don't cover use of interactions in hsd.test, assume can work.

does know how hsd.test work interactions? or there other function you've gotten work produce compact letter displays glm interactions?

i've been working on number of days , haven't been able find solution, i'm not missing obvious.

thanks!

i don't know how you've specified glm model, hsd.test, it's looking match particular treatment name same name specified in glm formula data frame. why main effect, factora work, not 3-way interaction. multiple comparison tests on interactions, find easiest generate interactions separately , add them data frame additional columns. glm model can specified using new variables code interaction.

for example,

set.seed(42) glm.dat <- data.frame(y = rnorm(1000), factora = sample(letters[1:2],     size = 1000, replace = true),  factorb = sample(letters[1:2], size = 1000, replace = true),  factorc = sample(letters[1:2], size = 1000, replace = true))  # generate interactions explicitly , add them data.frame glm.dat$factorab <- with(glm.dat, interaction(factora, factorb)) glm.dat$factorac <- with(glm.dat, interaction(factora, factorc)) glm.dat$factorbc <- with(glm.dat, interaction(factorb, factorc)) glm.dat$factorabc <- with(glm.dat, interaction(factora, factorb, factorc))  # general linear model   glm.mod <- glm(y ~ factora + factorb + factorc +  factorab + factorac +     factorbc + factorabc, family = 'gaussian', data = glm.dat)   # multiple comparison test  library(agricolae) comp <- hsd.test(glm.mod, trt = "factorabc", group = true) 

giving

comp$groups giving      trt        means m  1 a.a.a  0.070052189  2 a.b.b  0.035684571  3 b.a.a  0.020517535  4 b.b.b -0.008153257  5 a.b.a -0.036136140  6 a.a.b -0.078891136  7 b.a.b -0.080845419  8 b.b.a -0.115808772