11
$\begingroup$

I am finishing up some analysis on a large set of data. I would like to take the linear model used in the first part of the work and re-fit it using an linear mixed model (LME). The LME would be very similar with the exception that one of the variables used in the model would be used as a random effect. This data comes from many observations (>1000) in a small group of subjects (~10) and I know that modeling the effect of subject is better done as a random effect (this is a variable that I want to shift). The R code would look like:

my_modelB <- lm(formula = A ~ B + C + D)    
lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')

Everything runs fine and the results are vastly similar. It would be nice if I could use something like RLRsim or an AIC/BIC to compare these two models and decide which is the most appropriate. My colleagues don't want to report the LME because there isn't an easily accessible way of choosing which is "better", even though I think the LME is the more appropriate model. Any suggestions?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
8
$\begingroup$

This is to add to @ocram's answer because it is too long to post as a comment. I would treat A ~ B + C as your null model so you can assess the statistical significance of a D-level random intercept in a nested model setup. As ocram pointed out, regularity conditions are violated when $H_0: \sigma^2 = 0$, and the likelihood ratio test statistic (LRT) will not necessarily be asymptotically distributed $\chi^2$. The solution was I taught was to bootstrap the LRT (whose bootstrap distribution will likely not be $\chi^2$) parametrically and compute a bootstrap p-value like this:

library(lme4)
my_modelB <- lm(formula = A ~ B + C)
lme_model <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
lrt.observed <- as.numeric(2*(logLik(lme_model) - logLik(my_modelB)))
nsim <- 999
lrt.sim <- numeric(nsim)
for (i in 1:nsim) {
    y <- unlist(simulate(mymodlB))
    nullmod <- lm(y ~ B + C)
    altmod <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
    lrt.sim[i] <- as.numeric(2*(logLik(altmod) - logLik(nullmod)))
}
mean(lrt.sim > lrt.observed) #pvalue

The proportion of bootstrapped LRTs more extreme that the observed LRT is the p-value.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you for completing my answer. Also, sometimes people use a mixture of chi-squares instead of a chi-square distribution for the test statistic. $\endgroup$
    – ocram
    Commented Mar 17, 2011 at 20:46
  • $\begingroup$ @ocram +1 for your comment on deciding whether to treat the variable as random or fixed separately from the analysis. @MudPhud If your PI doesn't understand the issue and insists on a p-value, then maybe just show him the result of the test of the random effect (which you would include anyway in the write-up). $\endgroup$
    – lockedoff
    Commented Mar 17, 2011 at 20:55
  • $\begingroup$ Thanks for the code. When I ran it the result is none of the bootstrapped LRTs are greater than the observed, so this means that I can stick to the lm without the random effects or even the original variable thrown in. $\endgroup$
    – MudPhud
    Commented Mar 17, 2011 at 21:05
  • $\begingroup$ @MudPhud: Did you get any errors? Try typing lrt.sim to make sure they're not all zeros, in which case the most likely culprit would be that you don't have the package lme4 installed. $\endgroup$
    – lockedoff
    Commented Mar 17, 2011 at 21:08
  • $\begingroup$ They're not 0, just very small (~1e-6) compared to the observed (63.95). $\endgroup$
    – MudPhud
    Commented Mar 17, 2011 at 21:17
2
$\begingroup$

I am not totally sure to figure out what model is fitted when you use the lme function. (I guess the random effect is supposed to follow a normal distribution with zero mean?). However, the linear model is a special case of the mixed model when the variance of the random effect is zero. Although some technical difficulties exist (because $0$ is in the boundary of the parameter space for the variance) it should be possible to test $H_0:variance = 0$ vs $H_1: variance > 0$...

EDIT

In order to avoid confusion: The test mentioned above is sometimes used to decide whether or not the random effect is significant... but not to decide whether or not it should be transformed into a fixed effect.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The question is: is there are test to decide if the variable should be modeled as a mixed effect or random effect? Otherwise you could do the test you described and then test it with a chi-square dist (I'm not sure what the appropriate test would be). $\endgroup$
    – MudPhud
    Commented Mar 17, 2011 at 20:16
  • 2
    $\begingroup$ @MudPhud: Modelling a variable as a fixed or as a random effect should actually be decided before the analysis, when the study is planned. It depends, in particular, on the scope of your conclusions. Random effects allow more generalisability. It could also avoid some technical difficulties. For example, the asymptotics might break down when the number of parameters grow up, as it is the case when a categorical variable with a lot of levels is considered as a fixed variable. $\endgroup$
    – ocram
    Commented Mar 17, 2011 at 20:24
  • $\begingroup$ I agree, but when I tried explaining this to my PI he just turned around and asked for a p-value of some kind. I want to include this analysis in a manuscript, but he won't put it in if there isn't a more concrete justification. $\endgroup$
    – MudPhud
    Commented Mar 17, 2011 at 20:28
  • 1
    $\begingroup$ @MudPhud: To the best of my knowledge, there is no p-value for such a decision. If interest centers on the effect of the specific levels chosen then it should be considered as fixed. If the available factor levels are seen as a random sample from a larger population and that inferences are wanted for the larger population, the effect should be random. $\endgroup$
    – ocram
    Commented Mar 17, 2011 at 20:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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