# nlme estimates near zero variance for the random effects

I am doing various analysis on a small sample. Basically, we have an experiment where 14 subjects (UID 1 ~ 14) used one of the 6 instruments (MID 1 ~ 6) on 3 occasions (Sequence 1 ~ 3). Each time an outcome score was registered (between 1 ~ 100).

The test was double blind. The subjects were told they are measuring 3 different conditions while in reality they were either measuring conditions A, B, A or B, A, B (randomly assigned to the machines and users). The objective was to see if A and B are different or not.

To see if there is any significant difference between the ratings for the conditions A and B, I tried to fit a simple, random intercept model using the nlme package in R. I tried:

f.1 <- lme(Score ~ Condition, random = ~1|UID, data)


However, for some reason lme fails to fit the model: it gives no error or warning but the variance of the fitted random effect is essentially zero:

> summary(f.1)
Linear mixed-effects model fit by REML
Data: data
AIC      BIC   logLik
349.3259 356.0815 -170.663

Random effects:
Formula: ~1 | UID
(Intercept) Residual
StdDev: 0.0009303203 15.98295

Fixed effects: Score ~ Condition
Value Std.Error DF   t-value p-value
(Intercept) 77.47619  3.487766 27 22.213700  0.0000
ConditionA  -0.85714  4.932446 27 -0.173776  0.8633
Correlation:
(Intr)
ConditionA -0.707

Standardized Within-Group Residuals:
Min         Q1        Med         Q3        Max
-2.9704269 -0.4677603  0.2472873  0.7835730  1.4628682

Number of Observations: 42
Number of Groups: 14


I tried doing the same thing using lme4 and got the same results. The estimates for the intercept and the Condition factor is almost identical to a linear model if I use lm.

I am trying hard to understand what lme or lmer fail to estimate the random effect. I generated some data by simulation and both routines had no problem fitting the model so I doubt there is something wrong with the syntax of what I have used.

The data is here:

   UID MID Seq Score Condition
1    1   1   1    90  B
2    1   1   2    85  A
3    1   1   3    75  B
4    2   4   1    75  A
5    2   4   2    95  B
6    2   4   3    85  A
7    3   6   1    60  A
8    3   6   2    82  B
9    3   6   3    85  A
10   4   3   1    60  A
11   4   3   2    70  B
12   4   3   3    75  A
13   5   2   1    85  B
14   5   2   2    85  A
15   5   2   3    85  B
16   6   5   1    90  B
17   6   5   2    95  A
18   6   5   3   100  B
19   7   2   1    90  B
20   7   2   2    70  A
21   7   2   3    50  B
22   8   1   1    70  B
23   8   1   2    75  A
24   8   1   3    80  B
25   9   3   1    90  A
26   9   3   2    30  B
27   9   3   3    90  A
28  10   6   1    50  A
29  10   6   2    85  B
30  10   6   3    92  A
31  11   4   1    50  A
32  11   4   2    85  B
33  11   4   3    92  A
34  12   5   1    65  B
35  12   5   2    50  A
36  12   5   3    90  B
37  13   4   1    65  A
38  13   4   2    70  B
39  13   4   3    80  A
40  14   2   1    60  B
41  14   2   2   100  A
42  14   2   3    80  B

• Not an answer, but some comments on your data. First off I am not sure whether fitting a linear model to the scores is appropiate here. The mean in group A is already 77.5 with a stdev of 16. That means the possible max value is just 1.4 standard deviations from the mean. This means your distribution will be skewed. Also you seem to assume that the score does not depend on where the subject is in his sequence (1,2,3) or if he is in sequence ABA or BAB. – Erik Nov 21 '14 at 10:42
• You might be right. However, as a first step I like to fit a linear model and then take additional steps to make a better/more complex models. The problem is that the fitting procedure fails without warning. – Shapul Nov 21 '14 at 10:46
• Also, there might be a learning effect, true. If I add Seq as a covariate I can have an idea about it. However, with or without that, the procedure fails to fit the random effects. – Shapul Nov 21 '14 at 10:46
• The estimate is just very small. You just get a meaningful random effect if the variation between the individuals is larger than the error that would be expected when considering just inter-intravidual variance. Your stdev is 16, given 3 samples per subject you can estimate the individual means with an error of about 16/sqrt(3). If the standard deviation of the individual means is of a magnitude consistent with that you won't get a meaningful random effect. That's why my comments matter, since if your model does not fit the data well your residual standard deviation might be too high. – Erik Nov 21 '14 at 11:19
• I see. I like to expand the study to measure more subjects. But I don't think that the within subject variance would reduce a lot: the score is about their perception of delay, which is rather hard for individuals to assess. Based on what you said, it seems that increasing the sample size would not help much and we should focus more on possibly increasing the repetitions. – Shapul Nov 21 '14 at 13:48