[R-sig-ME] Advise on mixed effect model for a longitudinal study (with lme4)
Thierry Onkelinx
thierry.onkelinx at inbo.be
Fri Oct 9 10:05:31 CEST 2015
Dear Sylvain,
You can look at ratio of the number of effective observations and the
number of parameters. The number of effective observations depend on the
distribution: Gaussian = number of observations, Poisson = number of non
zero observations, binomial = min(number absences, number presences). A
common rule of thumb is to have at least 10 effective observations per
parameter.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
2015-10-08 16:12 GMT+02:00 wphantomfr <wphantomfr op gmail.com>:
> Dear Thierry,
>
> First thanks a lot for taking time to answer me.
>
> Thierry Onkelinx wrote:
>
>> Dear anonymous.
>>
> oups.. I forgot to Sign
>
>
>
>> Your random effects seem a bit to complicated for your data. A random
>> slope with 6 observations per subject is pushing the limits. A second
>> degree polynomial on 6 observations is nonsense. I can understand that
>> is does make sense on a conceptual level, but you just don't have enough
>> data to get sensible estimates on such a complex model.
>>
>> My advice would be to restrict the random effect to just a random
>> intercept. The benefit is that you gain 5 parameters and their
>> associated degrees of freedom. Note that a second order polynomial on
>> the fixed effect might be feasible.
>>
>
> I think I understand the problem but did not realize I was pushing the
> limits doing this because (1) it was making sense to use random slopes when
> looking at the data and (2) I was probably too much focused on the results
> of model comparisons that indicated a better fit for this complex model.
>
> But is there a way to know that we are pushing too much the limits ? Would
> only the random intercept + linear slope be acceptable (however the linear
> slope is probably optimal with our design)
>
>
>> You probably want to translate your hypotheses in a set of linear
>> combination of model parameters. Then you can test those hypotheses.
>> Have a look at the examples in the glht() function fom the mulcomp
>> package.
>>
>
> Well that's exactly that... I just didn't realise that glht could work
> with lme models
>
>>
>> PS Please don't post in HTML. It mangles up the output and code.
>>
> sorry. I think I have now disabled the default HTML mode.
>
>
> Thank you very much for your help !
>
> Sylvain CLEMENT
> PSITEC (EA 4072) laboratory
> "Neuropsychologie : Audition, Cognition, Action (NACA)" team
> University of Lille, France.
>
>
>
>> Best regards,
>>
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>> and Forest
>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>> Kliniekstraat 25
>> 1070 Anderlecht
>> Belgium
>>
>> To call in the statistician after the experiment is done may be no more
>> than asking him to perform a post-mortem examination: he may be able to
>> say what the experiment died of. ~ Sir Ronald Aylmer Fisher
>> The plural of anecdote is not data. ~ Roger Brinner
>> The combination of some data and an aching desire for an answer does not
>> ensure that a reasonable answer can be extracted from a given body of
>> data. ~ John Tukey
>>
>> 2015-10-08 10:59 GMT+02:00 wphantomfr <wphantomfr op gmail.com
>> <mailto:wphantomfr op gmail.com>>:
>>
>>
>> Dear list members,
>>
>> I currently trying to apply mixed effects models to a longitudinal
>> studies with a treatment period.
>>
>> The study is :
>> Several evaluation of a participants characteristic (let's call it DV)
>> at different times (in week) : 0 (1 evaluaition), 1 , 3, 6, 8 10.
>>
>> We have 3 different groups of participants (randomized controlled
>> study)
>> : a control group without treatment and groupA and groupB which
>> received
>> a different treatment.
>>
>> The treatment is applied to the groupA & B during weeks 2-5 therefore
>> :
>> - evaluations at w0 and w1 are "baseline measures"
>> (pre-treatment)
>> - evaluation at w3 is during treatment
>> - evaluation at w6 is post treatment
>> - evaluation at w8 & w10 are follow-up evaluation.
>>
>>
>>
>> After reading a lot on mixed effects model I arrived to the following
>> model (using lme4 package):
>>
>> best_model<-lmer(DV ~ (time+I(time^2))*GROUP
>> +(1+(time+I(time^2))|SUBJECT),data=brut)
>>
>> It includes a quadratic effect of time (continuous variable) and both
>> intercept and "slopes" (not a good word for the quadratic par) of the
>> effects of time (both linear and quadratic) are estimated for the
>> random
>> effects of SUBJECT.
>>
>> This model seems to be the best compared to models including only
>> parts
>> of these terms (model comparison with anova ). The summary report of
>> the
>> model is pasted at the end of this message
>>
>>
>> My hypothesis are :
>> - that I should see a superior effect of time in group A & B (if
>> the treatment increases DV, I should have a positive coeficient for
>> time
>> and a negative coeficient for time^2)
>> - that the effect of time may be superior in group B compared to
>> group A
>> - I would also like to see if one of the treatment have a more
>> long-term effect (sustained effect in the foloup evaluations)
>>
>> I'm not sure how to test/answer my questions from my model.
>>
>>
>>
>>
>> My guesses :
>> - the fact that GROUPA et GROUPB have low t-values is a sign that my
>> groups are quite similar at the begining of the study
>> - The coeficient of time:GROUP(A & B) and I(time^2):group(A&B) are in
>> line with my hypothesis
>> - I don't know how, in this context, I can answer to the question of
>> long-term effects...
>>
>>
>> My questions :
>>
>> Am I totally wrong ?
>>
>> How to best assess my hypothesis ?
>>
>> Any more advises ?
>>
>>
>> Thanks in advance
>>
>>
>>
>>
>> Here is the model report :
>>
>>
>> summary(bestM)
>> Linear mixed model fit by REML
>> t-tests use Satterthwaite approximations to degrees of freedom
>> ['lmerMod']
>> Formula: DV ~ (time + I(time^2)) * GROUP + (1 + (time + I(time^2))
>> | SUBJECT)
>> Data: brut
>>
>> REML criterion at convergence: -458.9
>>
>> Scaled residuals:
>> Min 1Q Median 3Q Max
>> -3.1714 -0.3972 0.0660 0.4385 2.7059
>>
>> Random effects:
>> Groups Name Variance Std.Dev. Corr
>> SUBJECT (Intercept) 6.201e-02 0.24902
>> time 2.042e-03 0.04519 -0.45
>> I(time^2) 1.162e-05 0.00341 0.41 -0.95
>> Residual 1.314e-02 0.11462
>> Number of obs: 710, groups: CODE, 117
>>
>> Fixed effects:
>> Estimate Std. Error df t value
>> Pr(>|t|)
>> (Intercept) 1.858e-01 4.557e-02 1.140e+02 4.078
>> 8.44e-05 ***
>> time 1.960e-03 1.159e-02 1.138e+02 0.169
>> 0.86599
>> I(time^2) 2.149e-04 1.032e-03 1.137e+02 0.208
>> 0.83547
>> GROUPA -2.006e-02 6.198e-02 1.140e+02 -0.324
>> 0.74681
>> GROUPB -2.898e-02 6.094e-02 1.138e+02 -0.476
>> 0.63531
>> time:GROUPA 3.837e-02 1.576e-02 1.138e+02 2.435
>> 0.01645 *
>> time:GROUPB 4.212e-02 1.546e-02 1.120e+02 2.724
>> 0.00748 **
>> I(time^2):GROUPA -2.749e-03 1.404e-03 1.137e+02 -1.957
>> 0.05277 .
>> I(time^2):GROUPB -3.222e-03 1.377e-03 1.113e+02 -2.340
>> 0.02108 *
>> ---
>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>> Correlation of Fixed Effects:
>> (Intr) time I(t^2) GROUPA GROUPB t:GROUPA t:GROUPB
>> I(^2):GROUPA
>> time -0.467
>> I(time^2) 0.379 -0.951
>> GROUPA -0.735 0.343 -0.279
>> GROUPB -0.748 0.349 -0.283 0.550
>> t:GROUPA 0.343 -0.735 0.699 -0.467 -0.257
>> t:GROUPB 0.350 -0.749 0.712 -0.257 -0.467 0.551
>> I(^2):GROUPA -0.279 0.699 -0.735 0.379 0.208 -0.951 -0.524
>> I(^2):GROUPB -0.284 0.713 -0.750 0.209 0.379 -0.524 -0.951
>> 0.551
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models op r-project.org
>> <mailto:R-sig-mixed-models op r-project.org> mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list