Dear Ben and Jake,

Thanks a lot for your answer, I appreciate it a lot.

Ben:
Thanks for your quick answer! Ok, this makes sense. To be honest, I did not
really get the idea of the firs suggestion ("fit the model with the random
factor as a fixed effect, get the level coefficients in the sum to zero
form, and then compute the standard deviation of the coefficients"). I am
not sure what is meant by level coefficients or how to calculate them.
Therefore, I tried to estimate the same model using the R library MCMCglmm.
However, I did not succeed in retrieving the correlation coefficient (yet).
It is not in included in the standard summary(model) output and the
function posterior.cor(model$VCV) that calculates correlations from the
covariance matrix only results in an error message.

Luckily, I found JAGS code (
http://people.ucsc.edu/~abrsvn/general_correlated_ranefs_bayes_jags.r, see
third search result for the phrase '# Bayesian analysis using JAGS' to find
the code I used) that could be easily changed to my purposes. This still
results in an estimate of very low variance for the random slopes (0.0003),
a variance of 0.62 for the random intercepts (which is comparable to the
variance estimated by more simple models) and a correlation between random
slopes and intercepts of 0.6. The covariance between the slopes is 0.014.
This seems more plausible. There are however a few things I still do not
understand conceptually (I a m a newbie concerning multilevel modelling):
The correlation of 0.6 seems moderate, but the variance of the random
slopes is very low. How can this be? It seems not logical to me that there
can be a moderate correlation with another variable if there is no almost
no variance at all. Secondly, in the model the correlation of the slopes
and intercepts is estimated as 0.6. However, when I just use the posterior
means for the random intercepts and the posterior means for the random
slopes and correlate them manually (cor(out$mean$rand_intercepts,
out$mean$rand_slopes) I am getting an estimate of 0.98, which 1) seems to
high, 2) is not equal to the estimated correlation in the model. How can
that be? I thought that it would make more sense if this correlation would
be equal to (or less) than as estimated in the model. If my questions are
unclear or lack information/detail, please let me know.

Jake:
Thanks for your remark a lot. When you say 'grand mean', do you mean the
population mean (the average test score of the population)?
I am not sure if I understand your remark correctly. I standardized both
the SES scores as well as the test scores. Therefore, the population mean
is zero, so the SES predictor is already centered around the grand mean.
Therefore, I think the interpretation of my results is already as you are
suggesting. Please let me know if I misunderstood your suggestion.

Kind regards,
Inga


2014-07-17 0:39 GMT+02:00 Jake Westfall <jake987722@hotmail.com>:

> Hi Inga,
>
> Interpreting covariances between the random intercepts and slopes is
> tricky because the covariance depends on the scaling of the predictor.
> Deviating the predictor around different constants can drastically alter
> the intercept/slope covariance. So I think you want to make sure that the
> scale of the predictor that you are using really does make good sense for
> the interpretation you are after. In your case, I wonder if centering the
> SES predictor around the grand mean wouldn't actually make more sense; the
> correlation in that case would be interpreted as: do schools with higher
> predicted outcomes at the overall mean SES also have higher SES slopes?
> Perhaps more relevantly, this centering may also help with the implausible
> correlation estimate, which likely has something to do with the 0-point on
> SES being far outside the range of the observed data. Just an idea.
>
> Jake
>
> > Date: Wed, 16 Jul 2014 09:57:08 -0400
> > From: bbolker@gmail.com
> > To: r-sig-mixed-models@r-project.org
> > Subject: Re: [R-sig-ME] Question regarding the estimation of the
> correlation between random intercepts and random slopes in the lme4 package
> >
> > On 14-07-16 08:37 AM, Inga Schwabe wrote:
> > > Dear list members,
> > >
> > > I posted a question regarding the lme4 package in R on
> stackoverflow.com
> > > (see
> > >
> http://stackoverflow.com/questions/24778642/estimating-correlation-between-random-slopes-and-random-intercepts-using-the-lme
> )
> > > and a stackoverflow-user kindly  referred me to this mailing list. I
> hope
> > > that my question is appropriate for this mailing list. If it is not,
> please
> > > let me know!
> > >
> > > This is my question:
> > >
> > > For answering my research question I am interested in the correlation
> > > between the random slopes and random intercepts in a multilevel model,
> > > estimated using the R library lme4.
> > >
> > > The data I have is: Y (test-scores from students), SES (socio-economic
> > > status for each student) and schoolid (ID for each school).
> > >
> > > I am using the following syntax to estimate random intercepts and
> slopes
> > > for the schools:
> > >
> > > library(lme4)
> > > model3 <- lmer(Y ~ SES + (1 + SES | schoolid))
> > >
> > >  The reference I used for this syntax is this pdf:
> > >
> > >
> http://www.bristol.ac.uk/cmm/learning/module-samples/5-concepts-sample.pdf
> > >
> > > On page 19, a similar analysis is described. It is said that by
> defining
> > > the random intercepts and slopes toghether, it is indirectly specified
> that
> > > we want the random intercepts and slopes to covary. Therefore, also the
> > > correlation between random slopes and random intercepts is estimated.
> > > Basically, exactly what I need for answering my research hyptohesis.
> > >
> > > However, when I look at the results:
> > >
> > >  summary(model3)
> > >
> > >  I am getting the following output:
> > >
> > > Linear mixed model fit by REML ['lmerMod']
> > > Formula: Y ~ SES + (1 + SES | schoolid)
> > >
> > > REML criterion at convergence: 8256.4
> > >
> > > Scaled residuals:
> > >     Min      1Q  Median      3Q     Max -3.1054 -0.6633 -0.0028
>  0.6810  3.5606
> > >
> > > Random effects:
> > >  Groups   Name        Variance  Std.Dev. Corr
> > >  schoolid (Intercept) 0.6427924 0.80174
> > >       SES             0.0009143 0.03024  1.00
> > >  Residual             0.3290902 0.57366
> > > Number of obs: 4376, groups: schoolid, 179
> > >
> > > Fixed effects:
> > >              Estimate Std. Error t value(Intercept) -0.036532
> 0.060582  -0.603
> > > SES          0.062491   0.009984   6.259
> > >
> > > Correlation of Fixed Effects:
> > >     (Intr)
> > > SES 0.226
> > >
> > >  As stated in the output, the correlation between the random slopes and
> > > random intercepts equals 1.00. I find this hard to believe. When I
> call in
> > > R:
> > >
> > > VarCorr(model3)$schoolid
> > >
> > >  I am getting the following output which gives the correlations and
> > > covariance matrix:
> > >
> > >                 (Intercept)          SES(Intercept)  0.64279243
> 0.0242429680
> > > SES          0.02424297 0.0009143255
> > >
> > > attr(,"stddev")(Intercept)         SES
> > >  0.80174337  0.03023782
> > >
> > > attr(,"correlation")
> > >         (Intercept) SES(Intercept)           1   1
> > > SES                   1   1
> > >
> > >  It seems as if the correlation between the slopes and intercepts was
> set
> > > to 1.00 by R. I did not see this in the output of anyone else when I
> was
> > > searching the internet on references on multilevel modelling.
> > >
> > > Does anybody know what can be the cause of this correlation? Can it be
> that
> > > the correlation is set to 1.00 because otherwise the model would not be
> > > identified? Or is it because the variance of the random slopes is so
> small
> > > (0.0009) that the correlation can not be estimated?
> > >
> > > I have tried to simulate data in order provide the code for a small
> > > reproducible dataset. I was however not yet able to reproduce this
> output
> > > by means of simulated data. As far as I have the code I will eidt my
> post
> > > and add the code.
> > >
> > >
> > > Many thanks for you time and effort,
> > > Kind regards,
> > >
> > > Inga
> > >
> >
> >   I was thinking about answering on StackOverflow, but since you've
> > asked here: yes, your second-to-last paragraph is exactly correct.  You
> > don't have enough information (or equivalently there is too much noise)
> > in your data to uniquely identify the full variance-covariance matrix,
> > so the results are 'singular'; that is, one of the underlying components
> > is identified as zero.  Common symptoms of this situation include either
> > variances equal to zero (or nearly equal, although I see in your case
> > that the variances while small are *not* exactly zero) or correlations
> > equal to +/- 1.
> >
> > http://rpubs.com/bbolker/4187 shows some simulated examples where the
> > estimated variance collapses to zero (despite the true, simulated model
> > having a non-zero among-group variance).  Presumably one could make up
> > similar examples (with randomized-block designs/random-slope models)
> > that would show analogous situations of correlations collapsing to +/- 1.
> >
> >  It's a little bit surprising that you have encountered this problem
> > since it is most commonly a problem with small numbers of levels in the
> > grouping variables, which doesn't appear to be the case here.
> >
> >   See http://glmm.wikidot.com/faq#singular_fits for more information ...
> >
> >   Ben Bolker
> >
> > _______________________________________________
> > R-sig-mixed-models@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]

