1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Sal Mangiafico Jul 20 '17 at 15:14
  • $\begingroup$ 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. $\endgroup$ – Leo Jul 20 '17 at 15:43
1
$\begingroup$

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,
                      adjust="tukey")

CLD = cld(leastsquare,
          alpha   = 0.05,
          Letters = letters,
          adjust  = "tukey")

CLD
$\endgroup$
  • $\begingroup$ Hi Sal, thanks for your answer. It might be that the output I`m not looking for in summary indeed. I dont quite understand what your code does though. You make fit with lme, after that Im not sure what lsmeans and cld do. Lsmeans calculates the mean using the fitted model from fit? And the standard error reported by cld also still takes into account the lme fit? $\endgroup$ – Leo Jul 20 '17 at 15:26
  • $\begingroup$ Ah yes I see that those functions can handle nlme objects, also in the reference you gave: rcompanion.org/handbook/I_07.html . Now the important bit of your answer is in the comments, might as well copy paste it into the answer for other people to find. Maybe also add an answer to the main question, then I'll accept it. Thanks for your help. $\endgroup$ – Leo Jul 20 '17 at 16:16
1
$\begingroup$

The basic answer to your question is "Yes. They are appropriate".

The question about group A is the same in regular regression: You have to have one group to compare to; you can also compare to the grand mean, and R has options for that.

I am not sure about the R output, but the coefficients are going to be used only as appropriate, unless there are mistakes in your code (but coding questions are off topic here).

$\endgroup$
  • 1
    $\begingroup$ Hi Peter, thanks for your answer. So if I perform the fit and interpret the results correctly it should be okay. Thats a step in the good direction :) $\endgroup$ – Leo Jul 20 '17 at 15:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.