Dear All, With kind help from several friends on the list, I am getting close. Now here are something interesting I just realized: for random effects, lmer reports standard deviation instead of standard error! Is there a hidden option that tells lmer to report standard error of random effects, like most other multilevel or mixed modeling software, so that we can say something like "randome effect for xxx is significant, while randome effect for xxx is not significant"? Thanks! Best, Shige
Doran, Harold
2005-Aug-17 11:16 UTC
[R] How to assess significance of random effect in lme4
You can extract the posterior variance of the random effect from the bVar slot of the fitted lmer model. It is not a hidden option, but a part of the fitted model. It just doesn't show up when you use summary(). Look at the structure of your object to see what is available using str(). However, your comment below seems to imply that it is incorrect for lmer to report SDs instead of the standard error, which is not true. That is a quantity of direct interest. Other multilevel programs report the same exact statistics (for the most part). For instance, HLM reports the variances as well. If you want the posterior variance of an HLM model you need to extract it. -----Original Message----- From: r-help-bounces@stat.math.ethz.ch on behalf of Shige Song Sent: Wed 8/17/2005 6:30 AM To: r-help@stat.math.ethz.ch Cc: Subject: [R] How to assess significance of random effect in lme4 Dear All, With kind help from several friends on the list, I am getting close. Now here are something interesting I just realized: for random effects, lmer reports standard deviation instead of standard error! Is there a hidden option that tells lmer to report standard error of random effects, like most other multilevel or mixed modeling software, so that we can say something like "randome effect for xxx is significant, while randome effect for xxx is not significant"? Thanks! Best, Shige ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html [[alternative HTML version deleted]]
Doran, Harold
2005-Aug-17 13:26 UTC
[R] How to assess significance of random effect in lme4
These are the posterior variances of the random effects (I think more properly termed "empirical" posteriors). Your model apparently includes three levels of random variation (commu, bcohort, residual). The first are the variances associated with your commu random effect and the second are the variances associated with the bcohort random effect. Accessing either one would require fm at bVar$commu or fm at bVar$bcohort Obviously, replace "fm" with the name of your fitted model. -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Shige Song Sent: Wednesday, August 17, 2005 7:50 AM To: r-help at stat.math.ethz.ch Subject: Re: [R] How to assess significance of random effect in lme4 Hi Harold, Thanks for the reply. I looked at my outputs using str() as you suggested, here is the part you mentioned: ..@ bVar :List of 2 .. ..$ commu : num [1, 1, 1:29] 5e-10 5e-10 5e-10 5e-10 5e-10 ... .. ..$ bcohort: num [1, 1, 1:6] 1.05e-05 7.45e-06 6.53e-06 8.25e-06 7.11e-06 ... where commu and bcohort are the two second-level units. Are these standard errors? Why the second vector contains a series of different numbers? Thanks! Shige On 8/17/05, Doran, Harold <HDoran at air.org> wrote:> > > You can extract the posterior variance of the random effect from the > bVar slot of the fitted lmer model. It is not a hidden option, but a > part of the fitted model. It just doesn't show up when you usesummary().> > Look at the structure of your object to see what is available usingstr().> > However, your comment below seems to imply that it is incorrect for > lmer to report SDs instead of the standard error, which is not true. > That is a quantity of direct interest. > > Other multilevel programs report the same exact statistics (for the > most part). For instance, HLM reports the variances as well. If you > want the posterior variance of an HLM model you need to extract it. > > > > -----Original Message----- > From: r-help-bounces at stat.math.ethz.ch on behalf of > Shige Song > Sent: Wed 8/17/2005 6:30 AM > To: r-help at stat.math.ethz.ch > Cc: > Subject: [R] How to assess significance of random effect inlme4> > Dear All, > > With kind help from several friends on the list, I am getting close. > Now here are something interesting I just realized: for random > effects, lmer reports standard deviation instead of standard error! Is> there a hidden option that tells lmer to report standard error of > random effects, like most other multilevel or mixed modeling software,> so that we can say something like "randome effect for xxx is > significant, while randome effect for xxx is not significant"? Thanks! > > Best, > Shige > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > > > > >______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Doran, Harold
2005-Aug-18 10:17 UTC
[R] How to assess significance of random effect in lme4
Yes, it is a different issue. ranef() extracts the empirical Bayes estimates,
which are the empirical posterior modes. The bVar slot holds the corresponding
posterior variances of these modes.
Technically, (according to D. Bates) the values in the bVar slot are the the
diagonal elements of
(Z'Z+\Omega)^{-1}.
The original post was asking how to test and compare a specific random effect,
not a general assessment of how much information is provided by the data via
LRT.
Shige asked how to test whether a specific EB estimate is different than some
other value.
LRT doesn't answer this question, but the values in the bVar slot do.
-----Original Message-----
From: Spencer Graves [mailto:spencer.graves@pdf.com]
Sent: Wed 8/17/2005 10:08 PM
To: Doran, Harold
Cc: Shige Song; r-help@stat.math.ethz.ch
Subject: Re: [R] How to assess significance of random effect in lme4
Is there some reason you are NOT using "anova", as in
"Examples"
section of "?lmer"?
Permit me to summarize what I know about this, and I'll be pleased if
someone else who thinks they know different would kindly enlighten me
and others who might otherwise be misled if anything I say is
inconsistent with the best literature available at the moment:
1. Doug Bates in his PhD dissertation and later in his book with Don
Watts (1988) Nonlinear Regression Analysis and Its Applications (Wiley)
split approximation errors in nonlinear least squares into "intrinsic
curvature" and "parameter effects curvature". He quantified
these two
problems in the context of roughly three dozen published examples, if my
memory is correct, and found that in not quite all cases, the parameter
effects were at least an order of magnitude greater than the intrinsic
curvature.
2. In nonnormal situations, maximum likelihood is subject to more
approximation error -- intrinsic curvature -- than "simple" nonlinear
least squares. However, I would expect this comparison to still be
fairly accurate, even if the differences may not be quite as stark.
3. The traditional use of "standard errors" to judge statistical
significance is subject to both intrinsic and parameter effects errors,
while likelihood ratio procedures such as anova are subject only to the
intrinsic curvature (assuming there are no substantive problems with
nonconvergence). Consequently, to judge statistical significance of an
effect, anova is usually substantially better than the so-called Wald
procedure using approximate standard errors, and is almost never worse.
If anyone knows of a case where this is NOT true, I'd like to know.
4. With parameters at a boundary as with variance components, the
best procedure seems to double the p-value from a nested anova (unless
the reported p-value is already large). This is because the
2*log(likelihood ratio) in such cases is roughly a 50-50 mixture of 0
and chi-square(1) [if testing only 1 variance component parameter].
This is supported by a substantial amount of research, including
simulations discussed in a chapter in Pinheiro and Bates (2000)
Mixed-Effects Models in S and S-Plus (Springer). The may be more
accurate procedures available in the literature, but none so simple as
this as far as I know.
Comments?
spencer graves
p.s. It looks like fm@bVars is a list containing vectors of length 29
and 6 in your example. I don't know what they are, but I don't see how
they can be standard errors in the usual sense.
Doran, Harold wrote:
> These are the posterior variances of the random effects (I think more
> properly termed "empirical" posteriors). Your model apparently
includes
> three levels of random variation (commu, bcohort, residual). The first
> are the variances associated with your commu random effect and the
> second are the variances associated with the bcohort random effect.
>
> Accessing either one would require
>
> fm@bVar$commu or fm@bVar$bcohort
>
> Obviously, replace "fm" with the name of your fitted model.
>
> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Shige Song
> Sent: Wednesday, August 17, 2005 7:50 AM
> To: r-help@stat.math.ethz.ch
> Subject: Re: [R] How to assess significance of random effect in lme4
>
> Hi Harold,
>
> Thanks for the reply. I looked at my outputs using str() as you
> suggested, here is the part you mentioned:
>
> ..@ bVar :List of 2
> .. ..$ commu : num [1, 1, 1:29] 5e-10 5e-10 5e-10 5e-10 5e-10 ...
> .. ..$ bcohort: num [1, 1, 1:6] 1.05e-05 7.45e-06 6.53e-06 8.25e-06
> 7.11e-06 ...
>
> where commu and bcohort are the two second-level units. Are these
> standard errors? Why the second vector contains a series of different
> numbers?
>
> Thanks!
>
> Shige
>
> On 8/17/05, Doran, Harold <HDoran@air.org> wrote:
>
>>
>>
>>You can extract the posterior variance of the random effect from the
>>bVar slot of the fitted lmer model. It is not a hidden option, but a
>>part of the fitted model. It just doesn't show up when you use
>
> summary().
>
>>
>> Look at the structure of your object to see what is available using
>
> str().
>
>>
>> However, your comment below seems to imply that it is incorrect for
>>lmer to report SDs instead of the standard error, which is not true.
>>That is a quantity of direct interest.
>>
>> Other multilevel programs report the same exact statistics (for the
>>most part). For instance, HLM reports the variances as well. If you
>>want the posterior variance of an HLM model you need to extract it.
>>
>>
>>
>> -----Original Message-----
>> From: r-help-bounces@stat.math.ethz.ch on behalf of
>>Shige Song
>> Sent: Wed 8/17/2005 6:30 AM
>> To: r-help@stat.math.ethz.ch
>> Cc:
>> Subject: [R] How to assess significance of random effect in
>
> lme4
>
>>
>> Dear All,
>>
>> With kind help from several friends on the list, I am getting close.
>> Now here are something interesting I just realized: for random
>>effects, lmer reports standard deviation instead of standard error! Is
>
>
>>there a hidden option that tells lmer to report standard error of
>>random effects, like most other multilevel or mixed modeling software,
>
>
>>so that we can say something like "randome effect for xxx is
>>significant, while randome effect for xxx is not significant"?
Thanks!
>>
>> Best,
>> Shige
>>
>> ______________________________________________
>> R-help@stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide!
>>http://www.R-project.org/posting-guide.html
>>
>>
>>
>>
>>
>
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA
spencer.graves@pdf.com
www.pdf.com <http://www.pdf.com>
Tel: 408-938-4420
Fax: 408-280-7915
[[alternative HTML version deleted]]
Douglas Bates
2005-Aug-19 16:44 UTC
[R] How to assess significance of random effect in lme4
On 8/17/05, Shige Song <shigesong at gmail.com> wrote:> Dear All, > > With kind help from several friends on the list, I am getting close. > Now here are something interesting I just realized: for random > effects, lmer reports standard deviation instead of standard error! Is > there a hidden option that tells lmer to report standard error of > random effects, like most other multilevel or mixed modeling software, > so that we can say something like "randome effect for xxx is > significant, while randome effect for xxx is not significant"? Thanks!Sorry to be entering this discussion late. I know there have been several responses and a considerable amount of discussion following the original question but I think it would be worthwhile returning to it. Shige asks why standard errors of the estimates of the variance components are not reported. All other known software packages for fitting mixed-effects models report the estimates of variances and covariances of the random effects and the standard errors of these estimates. However, the summary for a fitted lmer model or a fitted lme model does not. This is intentional. It is not an oversight. The reason these are not reported is because they are not useful and, in fact, are quite misleading. It is useful to have a standard error of an estimate when the distribution of the estimator is approximately symmetric. Then the standard error can be used to create an approximate confidence interval or perform a hypothesis test. However, the distribution of the estimate of a variance generally looks like a scaled chi-squared. Try the following after installing version 0.98-4 or later of the Matrix package> library(Matrix) > Sm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) > Ss1 <- mcmcsamp(Sm1, 10000, trans = FALSE) > str(Ss1)mcmc [1:10000, 1:6] 238 252 253 251 245 ... - attr(*, "class")= chr "mcmc" - attr(*, "mcpar")= int [1:3] 1 10000 1 - attr(*, "dimnames")=List of 2 ..$ : NULL ..$ : chr [1:6] "(Intercept)" "Days" "sigma^2" "Sbjc.(In)" ... This fits an lmer model to a longitudinal data set then creates a Markov Chain Monte Carlo sample from the parameters of the model. The default is to create an MCMC sample of transformed parameters that correspond to the logarithm of any variances and Fisher's 'z' transformation of the correlations. The "trans = FALSE" optional argument allows for the sample to be on the scale of the variances and the covariances. Check the density plots or the normal probability plots of the variance components and the covariance. The distribution is not at all normal or even symmetric. A symmetric confidence interval based upon a standard error or, even worse, a hypothesis test based on the estimate of a variance component and its standard error are clearly inappropriate. BTW, that mcmcsamp function is frighteningly fast when applied to a linear mixed model. It takes less than 2 seconds to generate a sample of size 10000 on the machine that I use.