# Is a mixed model appropriate to compare (continous) outcomes between (categorical) groups, with no other parameters?

Design

I have hierarchical data with three categorical levels: group, person, and cell. In group A, B and C I have some unequal number of people. For each person a to zzz I have some unequal number of cells. For each cell I have a continuous measurement x.

I want to compare this continuous measurement between groups, and allow for a random difference between people (and implicitly also between cells of course). As you can see I have no continuous explaining variable like time or concentration or such. I really only have that single value, there is no slope or line that can be fitted to this, it would be just an 'intercept' for each group.

Question

I wanted to use a mixed model because it can handle hierarchy. Is it appropriate to use a mixed model in my case (where the 'line' is a point at the intercept)?

I have found some other people in similar situations, but the basic question of is this even valid does not get treated, and I dont think the answer is trivial. Especially after having tried to perform the fit, which I show below.

Code

If I assume this approach is valid, and I perform it the way I think it should work, the results dont make sense to me. Here I try to explain cell numeric value x (intercept only) with factor group, and allow a random effect (intercept) per person:

fit = lme(cellx ~ group, random = ~1|person, data = longformatmeasurements)

coef(fit)

(Intercept)      groupB     groupC
a  -0.8272526 -0.04870979 0.02885113
b  -0.5515026 -0.04870979 0.02885113
c  -0.7333176 -0.04870979 0.02885113
d  -0.5902652 -0.04870979 0.02885113
e  -0.5966660 -0.04870979 0.02885113
...
zzz-0.4644939 -0.04870979 0.02885113


I already find this strange,why does every person have a coefficient for every group? An where is group A?

Then if I look at the summary:

summary(fit)

Linear mixed-effects model fit by REML
Data: longformatmeasurements
AIC      BIC    logLik
1254.006 1282.229 -622.0032

Random effects:
Formula: ~1 | person
(Intercept)  Residual
StdDev:   0.1592464 0.3174893

Fixed effects: cellx ~ group
Value  Std.Error   DF    t-value p-value
(Intercept) -0.6814135 0.05680689 2054 -11.995260  0.0000
groupB      -0.0487098 0.06878922   35  -0.708102  0.4836
groupC       0.0288511 0.07607308   35   0.379255  0.7068
Correlation:
(Intr) groupB
groupB -0.826
groupC -0.747  0.617

Standardized Within-Group Residuals:
Min         Q1        Med         Q3        Max
-6.5791412 -0.3009672  0.2738938  0.6191442  1.7675351

Number of Observations: 2092
Number of Groups: 38


Where is group A? Why do groups B and C have the number of people minus the number of groups as degrees of freedom, but is group A neglected and does 'intercept' get a degree of freedom for every single cell I measured? Do I just use the wrong syntax or is this just not possible? It looks like a slope is fitted between groups, but no other syntax that I tried is accepted by the function, except lme(cellx ~ group+0, random = ~1|person, data = longformatmeasurements), which does give a random intercept per person, but also a slope for every group and person.

• The coef output makes sense because you have asked the model to determine an intercept for each person. This is what the code in the random= option means. You'll note that the coefficients for Group A and Group B in this output are the same as in the table for summary. – Sal Mangiafico Jul 20 '17 at 15:14
• Ah right I get it, thats where those come from. Beforehand I expected the intercept to be the 'fitted mean', but I guess that only holds if there is no slope and no horizontal axis value under water. The code in your answer seems to calculate the 'fitted mean' per group, thanks. Maybe unde the hood each group has its mean at 'value 1' instead of 0 where the intercept lies. – Leo Jul 20 '17 at 15:43

As far as I can tell, the model you are using is appropriate for the situation you describe.

It sounds like the problem may be that you are looking for output that is not given by the summary. The following code may be helpful. It uses Anova from the car package to produce a traditional anova table, and lsmeans to conduct multiple comparisons among groups. Using least-square means in unbalanced designs can be helpful. The cld output is a compact form of the output from lsmeans. The standard output of the function would be available in the code below with simply printing the output: leastsquare. The code also tests if there is a significant effect of the random terms in the model by fitting a fixed-effects-only model with gls to compare to the mixed model.

To my knowledge, lsmeans and Anova handle objects from lme correctly. More information is available with ?car::Anova and ?lsmeans::models.

As always the caveat: There are no guarantees that the code below or the results are appropriate for your particular data. And, always understand the uses, assumptions, and limitations of functions and techniques before using them.

### Adopted from:
### http://rcompanion.org/handbook/I_07.html

### Install packages

if(!require(multcompView)){install.packages("multcompView")}
if(!require(lsmeans)){install.packages("lsmeans")}
if(!require(nlme)){install.packages("nlme")}
if(!require(car)){install.packages("car")}

### Create data

group = c(rep("A", 5), rep("B", 5), rep("C", 5))
person = rep(c("a","b","c","d","e"), 3)
cellx = c(1,2,3,4,5,5,6,7,8,9,2,4,6,8,10)

longformatmeasurements = data.frame(group, person, cellx)

### Model

fit = lme(cellx ~ group, random = ~1|person, data = longformatmeasurements)

### Anova

library(car)

Anova(fit)

### Effect of random variables

model.fixed = gls(cellx ~ group,
data=longformatmeasurements)

anova(fit,
model.fixed)

### Produce least square mean estimates for groups

library(multcompView)

library(lsmeans)

leastsquare = lsmeans(fit,
pairwise ~ group,

CLD = cld(leastsquare,
alpha   = 0.05,
Letters = letters,