# R - Post-hoc test on mixed design ANOVA with variable n's

I am attempting to perform post-hoc tests on a mixed-effects ANOVA, and I am running into issues with the typical R-based solutions, which I think stem from the fact that the groups in my between-subjects manipulation have different n's.

This is a snippet example of what my data looks like:

DT <- as.data.table(read.csv('~/data.csv))
print(DT)

Exp.Condition    Subject    time          y
1:             2            1      0    0.49777778
2:             2            1      1    0.39377778
3:             3            2      0    0.77155556
4:             3            2      1    0.32311111
5:             4            3      0    0.52733333
6:             4            3      1    0.53756565
7:             5            4      0    0.47688889
8:             5            4      1    0.44333333
9:             6            5      0    0.35555556
10:             6            5      1    0.33866667
11:             7            6      0    0.40000000
12:             7            6      1    0.39333333
13:       control            7      0    0.48355556
14:       control            7      1    0.08355556
.........


The size of each group is as follows:

DT[time==0,.N, Exp.Condition]

Exp.Condition     N
1:             2    17
2:             3    17
3:             4    15
4:             5    16
5:             6    13
6:             7    17
7:       control    60


After some searching I found that a common way to run post-hoc tests on the between-subjects factors of a mixed-effects ANOVA's in R is to to use lme(), so I created this model:

model <- lme(y ~ Exp.Condition * time, random = ~ 1 | Subject/Exp.Condition, data=DT)
print(anova(model))

>                   numDF denDF   F-value p-value
> (Intercept)            1   148 1899.4409  <.0001
> Exp.Condition          6   148    5.3273  0.0001
> time                   1   148   31.0840  <.0001
> Exp.Condition:time     6   148    3.2913  0.0045


and attempted this:

summary(glht(model, linfct=mcp(Exp.Condition="Tukey")))


but got this output:

Error in glht.matrix(model = list(modelStruct = list(reStruct = list(Exp.Condition = -9.60134899880588,  :
‘ncol(linfct)’ is not equal to ‘length(coef(model))’


I did some more searching and found another potential solution for making pairwise comparisons in mixed-effects models:

lsmeans(model, pairwise ~ Exp.Condition, adjust='tukey')


But that gave this output:

Error in adjustSigma && object\$method == "ML" :
invalid 'x' type in 'x && y'


I also found a manual solution in this blog post (under '1. Univariate approach using aov()'), but it seems to imply that n should be the same in each group, and I couldn't figure out a workaround. Is it possible to perform pairwise comparisons on the different levels of a between-subjects variable in a mixed-effects ANOVA with unequal n's? Is there a test that I'm missing, or am I headed off into a deep rabbit hole with no end in sight?

• I think the problem in lsmeans is confusion between your adjust='tukey' and the optional adjustSigma argument specific to its support for lme models. If you leave out the adjust argument altogether, it should work. The Tukey method is already the default anyway. – rvl Mar 2 '16 at 15:09
• However, since you have an interaction, doing marginal comparisons of your experimental conditions could be misleading. You ought to plot the predictions at each combination of the two factors (try the lsmip function) and see what's really going in. – rvl Mar 2 '16 at 15:14
• have you read these comments? – rvl Mar 3 '16 at 18:50
• yes, thanks for checking in. both were helpful in finding a solution, thank you for your help! I examined the interaction and determined that marginal comparisons accurately represented the data, so I just removed the adjust='tukey' and got my pairwise comparison p-values. – grrothman Mar 8 '16 at 0:25