[R-sig-ME] Comparison between the ouputs of the lmer and the brm (from the brms package) functions
Phillip Alday
me @end|ng |rom ph||||p@|d@y@com
Fri Apr 16 13:00:58 CEST 2021
sample_prior generates extra return values, it doesn't block sampling
from the posterior.
Thierry is right though -- you have < 100 observation and only 4 levels
of your grouping variable. Recall that random effects are *variance*
components -- estimating the variance with only 4 observations won't
give you great estimates. The maximum-likelihood estimate is well
defined, but the likelihood surface is probably quite flat and so you
can have a "highest point" on it, but there will be a lot of points that
are as nearly as high. And that's what the Bayesian analysis catches
better because it's not using a point estimate of the variance components.
Phillip
On 16/4/21 12:57 pm, Torsten Hauffe wrote:
> What is happening when you changeyour sample_prior = "yes" in brm() back
> to the default "no"?
> (Disclaimer: From the brm help, I only understand the settings "no" and
> "only" sampling the prior).
>
> Cheers!
>
> On Fri, Apr 16, 2021 at 12:50 PM Phillip Alday <me using phillipalday.com
> <mailto:me using phillipalday.com>> wrote:
>
> Note that brms does some special reparameterization magic with the
> intercept and uses a weakly informative prior for both the
> reparameterized intercept and the random effects.
>
> Also, if you look, you'll see that youR ESS is much lower for the
> intercept and the RE than for the other parameters. This can be
> indicative of the model having trouble exploring how those parameters
> relate and thus still having a large amount of uncertainty.
>
> In other words, the uncertainty in the random-effect of Intercept
> results in increased uncertainty for the fixed-effect of Intercept.
>
> Phillip
>
> On 16/4/21 11:54 am, Vincent Bremhorst wrote:
> > Dear all,
> >
> > I fitted the same mixed model using two different functions : lmer
> and brm.
> >
> > The estimation of the standard deviation of the random effect and
> the estimation of the standard errors of the intercept differ. Both
> estimates are higher with the Bayesian procedure.
> > Since I use non-informative prior in the brm specification, I
> would expect similar results.
> > The other estimates are similar for both procedures.
> >
> > Do you have any idea what's happen here?
> > Thanks for your help,
> > Vincent Bremhorst.
> >
> > lmer (model assumptions are met):
> >
> > res <- lmer(dTmeanoff ~ habitat + (1|week), data=trh)
> >
> >
> > Linear mixed model fit by REML. t-tests use Satterthwaite's method
> ['lmerModLmerTest']
> >
> > Formula: dTmeanoff ~ habitat + (1 | week)
> >
> > Data: trh
> >
> >
> >
> > REML criterion at convergence: 312.5
> >
> >
> >
> > Scaled residuals:
> >
> > Min 1Q Median 3Q Max
> >
> > -3.5949 -0.5149 -0.0181 0.4792 2.2995
> >
> >
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev.
> >
> > week (Intercept) 0.06142 0.2478
> >
> > Residual 1.93221 1.3900
> >
> > Number of obs: 89, groups: week, 4
> >
> >
> >
> > Fixed effects:
> >
> > Estimate Std. Error df t value Pr(>|t|)
> >
> > (Intercept) 0.6200 0.2788 12.5997 2.224 0.0451 *
> >
> > habitatu -0.8101 0.3503 83.0542 -2.312 0.0232 *
> >
> > habitatw -1.6366 0.3698 83.2151 -4.425 2.9e-05 ***
> >
> > ---
> >
> > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >
> >
> >
> > Correlation of Fixed Effects:
> >
> > (Intr) habitt
> >
> > habitatu -0.638
> >
> > habitatw -0.605 0.481
> >
> > brm (convergence of the posterior chains ok)
> >
> >
> >> priors<- c(prior(normal(0,100), class="b"),
> >
> > + prior(normal(0, 100), class="Intercept"),
> >
> > + prior(exponential(0.1), class="sigma"),
> >
> > + prior(exponential(0.1), class="sd", group= "week")
> >
> > +
> >
> > + )
> >
> >
> >
> >> fit1<- brm(dTmeanoff~habitat+(1|week), data=trh,
> >
> > + prior=priors,
> >
> > + iter=4000, warmup=2000,chains=2,
> >
> > + family = gaussian(),#no logit function applied
> >
> > + file = "output.rds",
> >
> > + sample_prior = "yes",
> >
> > + control = list(adapt_delta = .9))
> >
> >
> >
> > Family: gaussian
> >
> > Links: mu = identity; sigma = identity
> >
> > Formula: dTmeanoff ~ habitat + (1 | week)
> >
> > Data: trh (Number of observations: 89)
> >
> > Samples: 2 chains, each with iter = 4000; warmup = 2000; thin = 1;
> >
> > total post-warmup samples = 4000
> >
> >
> >
> > Group-Level Effects:
> >
> > ~week (Number of levels: 4)
> >
> > Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
> Tail_ESS
> >
> > sd(Intercept) 0.50 0.48 0.02 1.82 1.00 900
> 1048
> >
> >
> >
> > Population-Level Effects:
> >
> > Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
> >
> > Intercept 0.63 0.36 -0.11 1.40 1.00 1159 901
> >
> > habitatu -0.82 0.35 -1.51 -0.13 1.00 2235 2581
> >
> > habitatw -1.65 0.37 -2.37 -0.89 1.00 2091 2500
> >
> >
> >
> > Family Specific Parameters:
> >
> > Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
> >
> > sigma 1.41 0.11 1.21 1.65 1.00 2863 2199
> >
> >
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org
> <mailto:R-sig-mixed-models using r-project.org> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> _______________________________________________
> R-sig-mixed-models using r-project.org
> <mailto:R-sig-mixed-models using r-project.org> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list