0
$\begingroup$

I have various datasets I need to analyse regarding soil properties, all in the same fashion, with one fixed effect (which is a position along a transect, indicating different land uses). Now my main level of replication is across different transects, which will obviously have some form of random variance associated with them, and so I want to account for this in my statistical analyses.

So, in lme4 I specified a mixed model to this specification

model <- lmer(variable.of.interest ~ transect.position + (1|transect))

Now, when I analyse the above model against a model without the transect position term, I get exactly the results I was expecting, and then plugging the above model into anova(), I get the F values, d.f. etc. that I need.

However, I can't figure out how to say, for my overall report, that the random effect of transect does not make any difference to the overall analyses (i.e. I can't get a p-value, F value, d.f. etc.).

Help?

$\endgroup$
  • $\begingroup$ Hi Ewan, welcome. Could you be a bit more precise about your question? What exactly do you want to test/show? $\endgroup$ – COOLSerdash May 16 '13 at 9:41
  • $\begingroup$ Ok, so, Im looking at different soil properties within agricultural fields, with soil samples taken from a seciton of woodland and at various points into the field. On each of these samples, I have performed a variety of analyses, e.g. organic carbon, nitrogen, water holding etc etc. So for example, I want to show that organic carbon declines from the woodland and into the field, which is the "transect.position" term. And all I want to do with the mixed model is acount for the random effect of the potentially different soil properties along each transect. $\endgroup$ – user2037072 May 16 '13 at 9:49
  • $\begingroup$ Okay, I understand that. But you already managed to calculate the mixed effect model to account for the random effects. Is your question "How can I show that the random effects model is an improvement over the normal model with no random effect"? $\endgroup$ – COOLSerdash May 16 '13 at 9:53
  • $\begingroup$ Sorry, no. I want to be able to quote for a statistical output that the random effect of "transect" has no significant effect on the overall model output. $\endgroup$ – user2037072 May 16 '13 at 9:56
  • $\begingroup$ What do you mean by "no significant effect on the overall model output."? $\endgroup$ – COOLSerdash May 16 '13 at 10:00
3
$\begingroup$

You can use a likelihood ratio test (LRT) to test whether a random effect is significant. First, fit the random effects model. Then fit the model without the random effect. Extract and store the log-likelihood for each model using logLik and calculate the twice difference between the log-likelhood of the mixed effects model and the normal model. Use a $\chi^{2}_{1}$ distribution to calculate the $p$-value. Let me give an example:

library(lme4)

fm1 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy, REML=FALSE) # ranodm effects model (random intercept for each Subject)

fm2 <- lm(Reaction ~ Days, data=sleepstudy) # model without a random effect

D <- 2*as.numeric(logLik(fm1) - logLik(fm2))  # likelihood ratio statistic
D
[1] 106.2144                                # the chi2-value is very large

2*pchisq(D, 1, lower.tail=FALSE)            # p-value
[1] 1.323472e-24

In this case, the random effect of "Subject" clearly has an influence and should be included.

Alternatively, use the nlme package and the function lme to run the mixed-effect models. Then you can the anova function directly to calculate the likelihood ratio test:

library(nlme)
library(lme4)

fm1 <- lme(Reaction ~ Days, random=~1|Subject, data=sleepstudy, method="ML")

fm2 <- lm(Reaction ~ Days, data=sleepstudy)

anova(fm1, fm2)
    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm1     1  4 1802.079 1814.851 -897.0393                        
fm2     2  3 1906.293 1915.872 -950.1465 1 vs 2 106.2144  <.0001

Which provides the same result as the one we've manually calculated.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.