27
$\begingroup$

In a multi-level model, what are the practical and interpretation-related implications of estimating versus not-estimating random effect correlation parameters? The practical reason for asking this is that in the lmer framework in R, there is no implemented method for estimating p-values via MCMC techniques when estimates are made in the model of the correlations between parameters.

For example, looking at this example (portions quoted below), what are the practical implications of M2 versus M3. Obviously, in one case P5 will not be estimated and in the other it will.

Questions

  1. For practical reasons (the desire to get a p-value through MCMC techniques) one might want to fit a model without correlations between random effects even if P5 is substantially non-zero. If one does this, and then estimates p-values via the MCMC technique, are the results interpretable? (I know @Ben Bolker has previously mentioned that "combining significance testing with MCMC is a little bit incoherent, statistically, although I understand the urge to do so (getting confidence intervals is more supportable)", so if it will make you sleep better at night pretend I said confidence intervals.)
  2. If one fails to estimate P5, is that the same as asserting that it is 0?
  3. If P5 really is non-zero, then in what way are the estimated values of P1-P4 affected?
  4. If P5 really is non-zero, then in what way are the estimates of error for P1-P4 affected?
  5. If P5 really is non-zero, then in what ways are interpretations of a model failing to include P5 flawed?

Borrowing from @Mike Lawrence's answer (those more knowledgeable than I are free to replace this with full model notation, I'm not entirely confident I can do so with reasonable fidelity):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2) (Estimates P1 - P4)

M3: V1 ~ (1+V3|V2) + V3 (Estimates P1-P5)

Parameters that might be estimated:

P1: A global intercept

P2: Random effect intercepts for V2 (i.e. for each level of V2, that level's intercept's deviation from the global intercept)

P3: A single global estimate for the effect (slope) of V3

P4: The effect of V3 within each level of V2 (more specifically, the degree to which the V3 effect within a given level deviates from the global effect of V3), while enforcing a zero correlation between the intercept deviations and V3 effect deviations across levels of V2.

P5: The correlation between intercept deviations and V3 deviations across levels of V2

Answers derived from a sufficiently large and broad simulation along with accompanying code in R using lmer would be acceptable.

$\endgroup$
  • $\begingroup$ Related: stats.stackexchange.com/questions/46610/… $\endgroup$ – Jack Tanner Feb 17 '13 at 22:04
  • $\begingroup$ @JackTanner: It doesn't seem like you got satisfaction there either. It would be great if your concerns were also addressed in the answer to this question. $\endgroup$ – russellpierce Feb 18 '13 at 16:29
  • 4
    $\begingroup$ Giving an exact answer to many of your question - "what happens to _______ when I misspecify the model in _______ way" - is probably impossible without delving into, possibly intractable, theory (although this may be a special case where something is possible - I'm not sure). The strategy I'd probably use is to simulate data when the slope and intercept are highly correlated, fit the model constraining the two to be uncorrelated and compare the results with when the model is correctly specified (i.e. "sensitivity analysis"). $\endgroup$ – Macro Feb 19 '13 at 16:23
  • 4
    $\begingroup$ For your questions, I'm 80 (but not 100) % sure of the following: re. #2, Yes, if you don't estimate the correlation, you force it to be 0; for the rest, if the correlation is actually not exactly 0, then you are mis-specifying the non-independence of your data. The betas can nonetheless be unbiased, but the p-values will be off (& whether they are too high or too low depends & may not be knowable). Thus, interpretations of the betas may be able to proceed as normal, but interpretations of 'significances' will be inaccurate. $\endgroup$ – gung Feb 24 '13 at 21:20
  • 2
    $\begingroup$ @Macro: My hope was that a bounty might knock free a good answer based on theory rather than simulation. With a simulation I'll frequently be concerned I didn't pick up on an appropriate edge case. I am great at running simulations, but always feel a little... uncertain that I am running all of the right simulations (although I suppose I could leave that to journal editors to decide). I may have to ask another question about what scenarios to include. $\endgroup$ – russellpierce Feb 26 '13 at 15:45
15
+50
$\begingroup$

Consider the sleepstudy data, included in lme4. Bates discusses this in his online book about lme4. In chapter 3, he considers two models for the data.

$$M0: \textrm{Reaction} \sim 1 + \textrm{Days} + (1|\textrm{Subject}) +(0+\textrm{Days}|\textrm{Subject}) $$

and

$$MA: \textrm{Reaction} \sim 1 + \textrm{Days} + (\textrm{Days}|\textrm{Subject}) $$

The study involved 18 subjects, studied over a period of 10 sleep deprived days. Reaction times were calculated at baseline and on subsequent days. There is a clear effect between reaction time and the duration of sleep deprivation. There are also significant differences between subjects. Model A allows for the possibility of an interaction between the random intercept and slope effects: imagine, say, that people with poor reaction times suffer more acutely from the effects sleep deprivation. This would imply a positive correlation in the random effects.

In Bates' example, there was no apparent correlation from the Lattice plot and no significant difference between the models. However, to investigate the question posed above, I decided to take the fitted values of the sleepstudy, crank up the correlation and look at the performance of the the two models.

As you can see from the image, long reaction times are associated with greater loss of performance. The correlation used for the simulation was 0.58

enter image description here

I simulated 1000 samples, using the simulate method in lme4, based on the fitted values of my artificial data. I fit M0 and Ma to each and looked at the results. The original data set had 180 observations (10 for each of 18 subjects), and the simulated data has the same structure.

The bottom line is that there is very little difference.

  1. The fixed parameters have exactly the same values under both models.
  2. The random effects are slightly different. There are 18 intercept and 18 slope random effects for each simulated sample. For each sample, these effects are forced to add to 0, which means that the mean difference between the two models is (artificially) 0. But the variances and covariances differ. The median covariance under MA was 104, against 84 under M0 (actual value, 112). The variances of slopes and intercepts were larger under M0 than MA, presumably to get the extra wiggle room needed in the absence of a free covariance parameter.
  3. The ANOVA method for lmer gives an F statistic for comparing the Slope model to a model with only a random intercept (no effect due to sleep deprivation). Clearly, this value was very large under both models, but it was typically (but not always) larger under MA (mean 62 vs mean of 55).
  4. The covariance and variance of the fixed effects are different.
  5. About half the time, it knows that MA is correct. The median p-value for comparing M0 to MA is 0.0442. Despite the presence of a meaningful correlation and 180 balanced observations, the correct model would be chosen only about half the time.
  6. Predicted values differ under the two models, but very slightly. The mean difference between the predictions is 0, with sd of 2.7. The sd of the predicted values themselves is 60.9

So why does this happen? @gung guessed, reasonably, that failure to include the possibility of a correlation forces the random effects to be uncorrelated. Perhaps it should; but in this implementation, the random effects are allowed to be correlated, which means that the data are able to pull the parameters in the right direction, regardless of the model. The wrongness of the wrong model shows up in the likelihood, which is why you can (sometimes) distinguish the two models at that level. The mixed effects model is basically fitting linear regressions to each subject, influenced by what the model thinks they should be. The wrong model forces the fit of less plausible values than you get under the right model. But the parameters, at the end of the day, are governed by the fit to actual data.

enter image description here

Here is my somewhat clunky code. The idea was to fit the sleep study data and then build a simulated data set with the same parameters, but a larger correlation for the random effects. That data set was fed to simulate.lmer() to simulate 1000 samples, each of which was fit both ways. Once I had paired fitted objects, I could pull out different features of the fit and compare them, using t-tests, or whatever.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }
$\endgroup$
  • 1
    $\begingroup$ That is interesting work. Thank you. I do want to see what other comments come up in the next couple days and how things generalize to other cases before I accept the answer. Would you also consider including the relevant R Code in your answer as well as specifying the version of lmer you used? It would be interesting to feed the same simulated cases into PROC MIXED to see how it handles the unspecified random effects correlation. $\endgroup$ – russellpierce Jun 19 '14 at 20:03
  • 1
    $\begingroup$ @rpierce I have added the code sample as requested. I originally wrote it in LaTeX/Sweave, so the lines of code were interwoven with my comments to myself. I used version 1.1-6 of lme4, which is the current version in June 2014. $\endgroup$ – Placidia Jun 24 '14 at 1:14
  • $\begingroup$ @Ben when you say "Model A allows for" in the second paragraph, shouldn't that be MO? $\endgroup$ – nzcoops Apr 20 '16 at 5:14
  • $\begingroup$ I think the text is correct (all I did for this question was prettify the formula a little bit) $\endgroup$ – Ben Bolker Apr 20 '16 at 12:13
  • $\begingroup$ +6. Excellent answer, thanks for your attention to old but worthy questions. $\endgroup$ – amoeba Oct 11 '16 at 20:09
4
$\begingroup$

Placidia has already provided a thorough answer using simulated data based on the sleepstudy dataset. Here is another (less rigorous) answer that also uses the sleepstudy data.

We see that one can affect the estimated correlation between the random intercept and random slope by "shifting" the random predictor variable. Look at the results from models fm1 and fm2 below:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

From the model output, we see that the random variance correlation has changed. However, the slopes (fixed and random) stayed the same, as did the residual variance estimate. The intercepts (fixed and random) estimates changed in response to the shifted variable.

De-correlating random intercept-slope covariance for LMMs is discussed in Dr. Jack Weiss's lecture notes here. Weiss notes that reducing the variance correlation in this fashion can sometimes help with model convergence, among other things.

The above example varies the random correlation (parameter "P5"). Partially addressing the OP's Q3, we see from the above output that:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected
$\endgroup$
  • $\begingroup$ Thanks for adding signal to this long standing question! $\endgroup$ – russellpierce Jul 6 '17 at 3:42
  • $\begingroup$ Note: all of Jack Weiss's excellent lectures and class exercises/notes are linked in this post $\endgroup$ – theforestecologist Jun 5 '18 at 6:07
  • $\begingroup$ How should we then interpret the data in question? What is the "true" correlation? The one from the first or from the second model? Or the ones from BLUPs? $\endgroup$ – User33268 Dec 11 '18 at 7:47

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.