[R-sig-ME] Regarding mcmcsamp in the most recent alpha release
Dave Warren
dwarren1 at cyrus.psych.uiuc.edu
Sat Mar 29 21:15:50 CET 2008
Hi Simon, Doug, and list,
Thanks for the quick responses. Regarding the dataset, Trial does
indeed vary within Patient; I've got different numbers of trials
(observations, really) for each. The simpler fit you recommended does
the trick (ML here to improve the AIC estimate as suggested):
RT.lmer = lmer( log( RT ) ~ Trial * Status + (1 | Patient), data =
rt_agg, method = "ML" )
> summary( RT.lmer )
Linear mixed model fit by maximum likelihood
Formula: log(RT) ~ Trial * Status + (Trial | Patient)
Data: rt_agg
AIC BIC logLik deviance REMLdev
-71.24 -37.49 43.62 -87.24 -55.81
Random effects:
Groups Name Variance Std.Dev. Corr
Patient (Intercept) 4.0945e-02 0.20234777
Trial 4.4598e-07 0.00066781 0.432
Residual 4.5319e-02 0.21288152
Number of obs: 502, groups: Patient, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 10.9447820 0.0947941 115.46
Trial -0.0038163 0.0006622 -5.76
StatusNormal -0.2812827 0.1337640 -2.10
Trial:StatusNormal 0.0001416 0.0009196 0.15
Correlation of Fixed Effects:
(Intr) Trial SttsNr
Trial -0.046
StatusNorml -0.709 0.033
Trl:SttsNrm 0.033 -0.720 -0.036
Running mcmcsamp on this model yields no exceptions, which is
terrific. I should take a moment to praise the default xyplot output
that mcmcsamp objects now provide. It's very clear and concise.
HPDinterval works just fine with that output as well.
My original model formulation with a random coefficient for Trial
yields an AIC value of -71.24, slightly greater than the suggested
model. I naively modeled a random coefficient for Trial following a
suggestion of Doug's from earlier this month:
"""
Re: Repeated measures using lme
I'm not sure that those are equivalent specifications. If I read the
SAS code correctly (and I don't have a lot of experience with SAS) the
equivalent call to lme would have random = ~ 1 | individual
I would be more inclined to start with a model that does have a random
effect for time but does not have the additional correlation
structure. In lme this would be as in your specification but omitting
the correlation argument. In lmer it would be
lmer(IL6 ~ dust + time * company + (time | individual), data)
"""
It's arguable whether my data are truly repeated measures in the
strictest sense, but they are collected from unique individuals
performing very similar if not identical tasks at several dozen time
points ("Trial"s) each. At any rate, it seemed easy enough to try.
Thanks very much for the feedback!
Dave
Douglas Bates wrote:
> The first thing to notice is that the estimated correlation of the
> random effects is -1.000 which, as Simon indicates, calls into
> question the model that you are fitting. Also, does Trial vary within
> Patient? You can only fit a random Trial:Patient interaction if Trial
> varies within Patient.
>
> You only have 10 different Patients so it is unlikely that you will be
> able to estimate many variance parameters for the random effects. It
> may not be obvious but the model that you have fit is somewhat
> complicated. I would start with simpler models. First fit a model
> with the random effect term (1|Patient). The next model to attempt to
> fit is (1|Patient/Trial), which is what I think that Simon meant to
> write. This provides a random effect for Patient and a random effect
> for the Patient:Trial interaction. It is easier to fit this model
> than to fit the other version of a "random interaction" which is
> (1+Trial|Patient).
>
> Compare the AIC for those models with the AIC for the model that you
> fit. (To be more confident of the results you should fit all three
> models with method = "ML")
>
> That error message can occur even when the model is well-defined. It
> can be caused by the MCMC sampler getting stuck in regions with low,
> but very flat, posterior probability. I have seen it on samples like
>
> mcmcsamp(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy), 10000)
>
> but it is a low probability event for such models and I haven't found
> a seed where I can reproduce the problem repeatably.
>
>
> On Sat, Mar 29, 2008 at 8:21 AM, Simon Blomberg <s.blomberg1 at uq.edu.au> wrote:
>
>> Why do you have a random coefficient for Trial? Perhaps (1|Trial/Patient) might work?
>>
>
>
>> -----Original Message-----
>> From: r-sig-mixed-models-bounces at r-project.org on behalf of Dave Warren
>> Sent: Sat 29/03/2008 2:39 PM
>> To: r-sig-mixed-models at r-project.org
>> Subject: [R-sig-ME] Regarding mcmcsamp in the most recent alpha release
>>
>> Hi all,
>>
>> I've been trying out the most recent alpha release and keep running
>> into the following exception when I try using mcmcsamp.
>>
>> > path.lmer = lmer( log( PathLength ) ~ Trial * Status + (Trial |
>> Patient), data = path_agg )
>> > summary( path.lmer )
>> Linear mixed model fit by REML
>> Formula: log(PathLength) ~ Trial * Status + (Trial | Patient)
>> Data: path_agg
>> AIC BIC logLik deviance REMLdev
>> 640.4 674.1 -312.2 594.9 624.4
>> Random effects:
>> Groups Name Variance Std.Dev. Corr
>> Patient (Intercept) 1.3567e-03 0.0368331
>> Trial 2.1394e-05 0.0046254 -1.000
>> Residual 1.9524e-01 0.4418560
>> Number of obs: 502, groups: Patient, 10
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 2.995689 0.060844 49.24
>> Trial -0.001977 0.002406 -0.82
>> StatusNormal 0.245692 0.084061 2.92
>> Trial:StatusNormal 0.008709 0.003382 2.58
>>
>> Correlation of Fixed Effects:
>> (Intr) Trial SttsNr
>> Trial -0.661
>> StatusNorml -0.724 0.479
>> Trl:SttsNrm 0.470 -0.711 -0.660
>>
>> > a = mcmcsamp( path.lmer, 10000 )
>> Error in .local(object, n, verbose, ...) :
>> crossproduct matrix 1 is not positive definite
>>
>> Have I broken something? Is my model revealed to be invalid?
>> Thanks for any thoughts,
>>
>> Dave
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>> [[alternative HTML version deleted]]
>>
>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>
>
More information about the R-sig-mixed-models
mailing list