In a multi-level model, what are the practical implications of estimating versus not-estimating random effect correlation parameters? 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


*

*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.)

*If one fails to estimate P5, is that the same as asserting that it is 0?

*If P5 really is non-zero, then in what way are the estimated values of P1-P4 affected?

*If P5 really is non-zero, then in what way are the estimates of error for P1-P4 affected?

*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.
 A: 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

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.


*

*The fixed parameters have exactly the same values under both models.

*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.

*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).

*The covariance and variance of the fixed effects are different.

*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.

*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.

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)
    }

