Comparison between the ouputs of the lmer and the brm (from the brms package) functions
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 at phillipalday.com
<mailto:me at 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 at r-project.org
<mailto:R-sig-mixed-models at r-project.org> mailing list
_______________________________________________
R-sig-mixed-models at r-project.org
<mailto:R-sig-mixed-models at r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models