[R-sig-ME] Comparison between the ouputs of the lmer and the brm (from the brms package) functions

[Thierry Onkelinx]

Dear Vincent,

Also, notice that you have a tiny data set. Hence the prior distribution
will dominate the posterior distribution.
You have too few levels for a reasonable estimate of the random effect
variance.

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx using inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>


Op vr 16 apr. 2021 om 12:50 schreef Phillip Alday <me using phillipalday.com>:

> 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 mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]