[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