Hi:
Let's start with the data:
> str(test.data)
'data.frame': 100 obs. of 4 variables:
$ StudentID: num 17370 17370 17370 17370 17379 ...
$ GroupID : num 1 1 1 1 1 1 1 1 1 1 ...
$ Time : num 1 2 3 4 1 2 3 4 1 2 ...
$ Score : num 76.8 81.8 89.8 92.8 75.9 ...
Both StudentID and GroupID are numeric; in the model, they would be treated
as continuous covariates rather than factors, so we need to convert:
test.data$StudentID <- factor(test.data$StudentID)
test.data$GroupID <- factor(test.data$GroupID)
Secondly, I believe there are some flaws in your model. After converting
your variables to factors, I ran
library(lme4)
mlmoded1.lmer <- lmer(Score ~ Time + (Time | GroupID/StudentID), data =
test.data)
You have two groups, so they should be treated as a fixed effect - more
specifically, as a fixed blocking factor. The StudentIDs are certainly
nested within GroupID, and Time is measured on each StudentID, so it is a
repeated measures factor. The output of this model is
> mlmoded1.lmer
Linear mixed model fit by REML
Formula: Score ~ Time + (Time | GroupID/StudentID)
Data: test.data
AIC BIC logLik deviance REMLdev
393.1 416.5 -187.5 376.9 375.1
Random effects:
Groups Name Variance Std.Dev. Corr
StudentID:GroupID (Intercept) 0.504131 0.71002
Time 0.083406 0.28880 1.000
GroupID (Intercept) 12.809567 3.57905
Time 3.897041 1.97409 -1.000
Residual 1.444532 1.20189
Number of obs: 100, groups: StudentID:GroupID, 25; GroupID, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 72.803 2.552 28.530
Time 4.474 1.401 3.193
Correlation of Fixed Effects:
(Intr)
Time -0.994
The high correlations among the random effects and then among the fixed
effects suggests that the model specification may be a bit off.
The above model fits random slopes to GroupIDs and StudentIDs, along with
random intercepts, but GroupID is a between-subject effect and should be at
the top level. Time is a within-subject effect and StudentIDs are the
observational units. I modified the model to provide fixed effects for
GroupIDs, scalar random effects for StudentIDs and random slopes for
StudentIDs.
> mod3 <- lmer(Score ~ 1 + GroupID + Time + (1 | StudentID) +
+ (0 + Time | StudentID), data = test.data)
> mod3
Linear mixed model fit by REML
Formula: Score ~ 1 + GroupID + Time + (1 | StudentID) + (0 + Time |
StudentID)
Data: test.data
AIC BIC logLik deviance REMLdev
430.9 446.5 -209.4 418.4 418.9
Random effects:
Groups Name Variance Std.Dev.
StudentID (Intercept) 4.2186e-13 6.4951e-07
StudentID Time 1.8380e+00 1.3557e+00
Residual 1.6301e+00 1.2768e+00
Number of obs: 100, groups: StudentID, 25
Fixed effects:
Estimate Std. Error t value
(Intercept) 70.7705 0.4204 168.33
GroupID2 4.0248 0.5854 6.88
Time 4.5292 0.2942 15.39
Correlation of Fixed Effects:
(Intr) GrpID2
GroupID2 -0.668
Time -0.264 0.000
I didn't check the quality of the fit, but on the surface it seems to be
more stable, FWIW. Perhaps one could also add a term (GroupID | StudentID),
but I don't know offhand if that would make any sense. Another issue to
consider is whether to fit by REML or ML, but that is secondary to getting
the form of the model equation right. I don't claim this as a final model,
but rather a 're-starting point'. It may well be in need of improvement, so
comments are welcome.
The confusion between subjects nested in time or vice versa has occurred
several times this week with respect to repeated measures/longitudinal
models using lmer(), so perhaps it merits a comment: subjects/experimental
units are NOT nested in time. Measurements taken on an individual at several
time points *entails* that time be nested within subject. Just saying...
This discussion may be better continued on the R-sig-mixed list, so I've
cc-ed to that group as well.
HTH,
Dennis
On Wed, Aug 25, 2010 at 1:27 AM, Bertolt Meyer <bme...@sozpsy.uzh.ch> wrote:
> Ben Bolker <bbolker <at> gmail.com> writes:
>
> Bertolt Meyer <bmeyer <at> sozpsy.uzh.ch> writes:
>>
>>
>>> Hello lmer() - users,
>>>
>>> A call to the lmer() function causes my installation of R (2.11.1 on
>>> Mac OS X 10.5.8) to crash and I am trying to figure out the problem.
>>>
>>
>> [snip snip]
>>
>> detach("package:nlme")
>>> library(lme4)
>>>
>>> mod1 <- lmer(performance ~ time + (time | GroupID/StudentNumber), data
>>> = dataset.long, na.action = na.omit)
>>>
>>> However, this call results in a segfault:
>>>
>>> *** caught segfault ***
>>> address 0x154c3000, cause 'memory not mapped'
>>>
>>> and a lengthy traceback. I can reproduce this error. It also occurs
>>> when I don't load nlme before lme4. Can someone tell me what I am
>>> doing wrong? Any help is greatly appreciated.
>>>
>>
>> This may well be a bug in lmer. There have been a number of
>> fussy computational issues with the lme4 package on the Mac platform.
>>
>
> Ben, thanks for your reply. I tried to replicate this issue with a small
> clean data set on a windows machine. You can find the code for the data
> frame (100 observations from my data) at the end of this mail. Very simple:
> four test scores per student over time, and students are nested in groups.
> On my Windows installation, lmer() throws an error that does not seem to get
> caught on the Mac, resulting in the segfault:
>
> library(lme4)
> mlmoded1.lmer <- lmer(Score ~ Time + (Time | GroupID/StudentID), data =
> test.data)
>
> Error: length(f1) == length(f2) is not TRUE
> Addditional Warnings:
> 1: In StudentID:GroupID :
> numeric expression has 100 elements: only first one is used
> 2: In StudentID:GroupID :
> numeric expression has 100 elements: only first one is used
>
> It seems to me that I am committing a trivial error here and that I am too
> blind to see it. Any ideas?
>
> Regards,
> Bertolt
>
> If it is at all possible, please (1) post the results of sessionInfo()
>> [which will in particular specify which version of lme4 you are using];
>> (2) possibly try this with the latest development version of lme4, from
>> R-forge, if that's feasible (it might be necessary to build the package
>> from source), and most importantly:
>>
>> (3) create a reproducible (for others) example -- most easily by
>> posting your data on the web somewhere, but if that isn't possible
>> by simulating data similar to yours (if it doesn't happen with another
>> data set of similar structure, that's a clue -- it says it's some more
>> particular characteristic of your data that triggers the problem) and
>>
>> (4) post to to *either* the R-sig-mac or the R-sig-mixed-models list,
>> where the post is more likely to come to the attention of those who
>> can help diagnose/fix ...
>>
>> good luck
>> Ben Bolker
>>
>>
>
> test.data <- data.frame(c(17370, 17370, 17370, 17370, 17379, 17379, 17379,
> 17379, 17387, 17387, 17387, 17387, 17391, 17391, 17391, 17391, 17392, 17392,
> 17392, 17392, 17394, 17394, 17394, 17394, 17408, 17408, 17408, 17408, 17419,
> 17419, 17419, 17419, 17429, 17429, 17429, 17429, 17432, 17432, 17432, 17432,
> 17436, 17436, 17436, 17436, 17439, 17439, 17439, 17439, 17470, 17470, 17470,
> 17470, 17220, 17220, 17220, 17220, 17348, 17348, 17348, 17348, 17349, 17349,
> 17349, 17349, 17380, 17380, 17380, 17380, 17398, 17398, 17398, 17398, 17400,
> 17400, 17400, 17400, 17402, 17402, 17402, 17402, 17403, 17403, 17403, 17403,
> 17413, 17413, 17413, 17413, 17416, 17416, 17416, 17416, 17420, 17420, 17420,
> 17420, 17421, 17421, 17421, 17421), c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), c(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4,
> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1,
> 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2,
> 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3,
> 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4), c(76.76, 81.83, 89.78, 92.82, 75.86,
> 81.84, 88.96, 92.28, 75.28, 80.68, 88.62, 92.29, 76.60, 84.59, 92.03, 94.05,
> 75.57, 79.94, 86.11, 90.25, 74.54, 81.42, 87.50, 90.71, 76.02, 83.68, 91.11,
> 94.14, 76.31, 83.76, 90.44, 94.58, 72.29, 80.51, 86.09, 90.41, 74.99, 82.28,
> 88.77, 92.26, 75.28, 81.92, 89.25, 92.64, 76.31, 83.93, 91.00, 94.60, 76.31,
> 82.44, 90.57, 95.17, 76.94, 82.21, 83.81, 85.00, 79.96, 81.92, 86.32, 90.05,
> 82.01, 84.81, 88.79, 93.10, 77.87, 82.94, 86.86, 90.31, 77.87, 79.64, 85.66,
> 86.97, 79.35, 80.44, 84.26, 83.62, 79.06, 81.56, 85.00, 87.43, 79.34, 81.47,
> 83.23, 86.86, 79.44, 80.37, 84.36, 89.11, 78.77, 81.02, 81.60, 87.21, 75.75,
> 79.35, 80.38, 86.87, 76.04, 80.57, 83.36, 86.31))
>
> names(test.data) <- c("StudentID", "GroupID", "Time", "Score")
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
