edit this please > mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)', data = indvar_FODs, + family=binomial) boundary (singular) fit: see help('isSingular') > mod Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod'] Family: binomial ( logit ) Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs) Data: indvar_FODs AIC BIC logLik deviance df.resid 4616.503 4671.476 -2299.252 4598.503 3312 Random effects: Groups Name Std.Dev. subject_FODs (Intercept) 1.403e-07 Number of obs: 3321, groups: subject_FODs, 44 Fixed Effects: (Intercept) Condition_FODsB Condition_FODsC Condition_FODsD 0.02715 -0.17148 -0.14466 -0.12076 Language_used_FODsGerman Language_used_FODsHungarian Language_used_FODsItalian Language_used_FODsTurkish 0.02854 0.16117 0.08189 0.06950 optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

library(multcomp)

mod_glht <- glht(mod, linfct = mcp(Condition_FODs = "Tukey"))
summary(mod_glht)

library(lsmeans)
lsmeans(mod, ~Condition_FODs)