Re: [R] mixed model random interaction term log likelihood ratio test

From: Ben Bolker <bbolker_at_gmail.com>
Date: Thu, 14 Apr 2011 22:14:05 +0000

seatales <ssphadke <at> uh.edu> writes:

>
> Hello,
> I am using the following model
>
> model1=lmer(PairFrequency~MatingPair+(1|DrugPair)+(1|DrugPair:MatingPair),
> data=MateChoice, REML=F)
>
> 1. After reading around through the R help, I have learned that the above
> code is the right way to analyze a mixed model with the MatingPair as the
> fixed effect, DrugPair as the random effect and the interaction between
> these two as the random effect as well. Please confirm if that seems
> correct.

  You should probably send this sort of question to the r-sig-mixed-models mailing list ...

  You probably want (MatingPair|DrugPair) rather than (1|DrugPair:MatingPair).
Whether REML=FALSE or REML=TRUE depends what you want to do next.

>
> 2. Assuming the above code is correct, I have model2 in which I remove the
> interaction term, model3 in which I remove the DrugPair term and model4 in
> which I only keep the fixed effect of MatingPair.
>
> 3. I want to perform the log likelihood ratio test to compare these models
> and that's why I have REML=F. However the code anova(model1, model2, model3,
> model4) gives me a chisq estimate and a p-value, not the LRT values. How do
> I get LRT (L.Ratio) while using lmer?

  The chi-squared values are the differences in deviance (-2 log likelihood) between the respective papers of models, which under the null hypothesis of the LRT will be chi-squared distributed. In other words, these *are* the LRT test statistics.
>
> 4. I am under the impression after reading a few posts that LRT is not
> usually obtained with lmer but it is given if I use lme (the old model).

   I don't know what you mean by this.
   The main difference between lmer and lme in the testing/inference context is that lme is willing to guess at "denominator degrees of freedom" to perform conditional F-tests.

>
> 5. I could not find how to input the random interaction term while using
> lme? Is it the following way? Would someone please guide me to some already
> existing posts or help here?

  See above.

> Thanks a lot. I hope someone can help. Most posts I have found deal with
> nesting but there is absolutely no nesting in my data.
>
> Sujal P.
> p.s: If it matters how data is arranged, then I have one vector called
> MatingPair which has 3 levels and another vector DrugPair which also has 3
> levels. The PairFrequency data is a count data and is normally distributed.
> The data are huge, hence I am not able to post it here.

  It is probably unwise to estimate DrugPair as a random effect if it only has three levels.

>
> Also, here is what I mean by getting chisq value rather than L.Ratio:

  See above.

> Data: MateChoice
> Models:
> model2: PairFrequency ~ MatingPair + (1 | DrugPair)
> model3: PairFrequency ~ MatingPair + (1 | DrugPair:MatingPair)
> model1: PairFrequency ~ MatingPair + (1 | DrugPair) + (1 |
> DrugPair:MatingPair)
>
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> model2 5 274.90 282.82 -132.45
> model3 5 282.44 290.36 -136.22 0.0000 0 1.00000
> model1 6 276.90 286.40 -132.45 7.5443 1 0.00602 **
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1
> ‘ ’ 1
>
> --
> View this message in context:

http://r.789695.n4.nabble.com/mixed-model-random-interaction-term-log-likelihood-ratio-test-tp3448718p3448718.html
> Sent from the R help mailing list archive at Nabble.com.
> [[alternative HTML version deleted]]
>
>



R-help_at_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. Received on Fri 15 Apr 2011 - 01:12:59 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 15 Apr 2011 - 06:20:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive