[R-sig-ME] mixed lognormal hurdle model with multiple grouping factors
Guillaume Adeux
gu|||@ume@|mon@@2 @end|ng |rom gm@||@com
Fri Dec 20 15:31:07 CET 2019
Yes, for now, the most seducing solution has been the following (with the
brms package):
fit=brm(bf(H_q~block+syst+(1|plot)+(1|year)+(1|plot:year),hu~block+syst+(1|plot)+(1|year)+(1|plot:year)),data=density,family=hurdle_lognormal(),control
= list(adapt_delta = 0.999))
Family: hurdle_lognormal
Links: mu = identity; sigma = identity; hu = logit
Formula: H_q ~ block + syst + (1 | plot) + (1 | year) + (1 | plot:year)
hu ~ block + syst + (1 | plot) + (1 | year) + (1 | plot:year)
Data: density (Number of observations: 480)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~plot (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.12 0.11 0.01 0.41 1.00 852 912
sd(hu_Intercept) 1.21 1.01 0.05 3.86 1.01 637 1011
~plot:year (Number of levels: 60)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.13 0.06 0.01 0.23 1.00 737 546
sd(hu_Intercept) 2.10 0.46 1.34 3.15 1.00 1427 2435
~year (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.06 0.07 0.00 0.21 1.00 1497 1843
sd(hu_Intercept) 0.83 0.71 0.03 2.55 1.01 935 1469
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.23 0.07 -0.38 -0.08 1.00 2159 1513
hu_Intercept -1.04 0.66 -2.26 0.43 1.00 1832 2061
block1 -0.06 0.06 -0.17 0.07 1.00 2394 1933
syst1 -0.32 0.17 -0.64 0.03 1.00 1930 1487
syst2 -0.11 0.12 -0.35 0.13 1.00 2023 1617
syst3 0.18 0.11 -0.06 0.41 1.00 1724 1267
syst4 0.08 0.12 -0.16 0.32 1.00 2222 1691
hu_block1 0.36 0.60 -0.81 1.63 1.00 1658 1684
hu_syst1 4.19 1.34 1.61 6.91 1.00 1236 1233
hu_syst2 -0.09 1.17 -2.40 2.38 1.00 1266 1419
hu_syst3 -1.98 1.30 -4.59 0.54 1.01 1183 1235
hu_syst4 -0.52 1.23 -2.95 1.88 1.00 1170 1108
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.48 0.02 0.44 0.52 1.00 2632 3025
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
I am currently struggling to obtain "syst" means averaged over blocks (i.e.
what would be done traditionally with emmeans(fit,~syst)) and to compare
them pairwise.
As of now, I am only capable of retrieving a mean for each combination of
"syst" and "block" :
newdata =
expand.grid(syst=levels(density$syst),block=levels(density$block))
means=fitted(fit,newdata,re_formula=NA,summary=TRUE)
colnames(means)=c("fit","se","lwr","upr")
df_plot=cbind(newdata,means)
I am sure this must be simple, particularly with the hypothesis() function
but I am not sure how to include the "zero-non zero part" (i.e. hu_...).
If anyone can foward some information on this, I would be eternally
grateful (or even more should I say).
Sincerely,
Guillaume ADEUX
Le ven. 20 déc. 2019 à 15:12, Ben Bolker <bbolker using gmail.com> a écrit :
>
> Good point.
> This might be manageable in MCMCglmm or brms (or JAGS) ...
>
> On 2019-12-20 2:58 a.m., D. Rizopoulos wrote:
> > For a hurdle model for repeated measurements data, the dichotomous
> outcome I(diversity > 0) is also a repeated measurements outcome. Hence, in
> the logistic regression model for this dichotomous outcome you will need to
> include random effects to account for the correlations. And it is logical
> to assume that the random effects from this logistic regression model will
> be correlated with the random effects of the linear mixed model for only
> the positive responses.
> >
> > In this case the likelihood of the two parts does not split in two
> functionally independent parts that can be separately maximized. If this is
> indeed the case, then fitting the two parts separately may cause bias and
> loss of efficiency.
> >
> > Best,
> > Dimitris
> >
> > ��
> > Dimitris Rizopoulos
> > Professor of Biostatistics
> > Erasmus University Medical Center
> > The Netherlands
> >
> > ________________________________
> > From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on
> behalf of Mollie Brooks <mollieebrooks using gmail.com>
> > Sent: Wednesday, December 18, 2019 3:08 PM
> > To: Help Mixed Models; Guillaume Adeux
> > Subject: Re: [R-sig-ME] mixed lognormal hurdle model with multiple
> grouping factors
> >
> > Hi Guillaume,
> >
> > I don�t think the hurdle lognormal can be fit in a single function call
> to glmmTMB since the model for the non-zero response requires
> log-transforming the response. Other types of hurdle models could be fit in
> glmmTMB using the zero-inflation model.
> >
> > I don�t think you gain much information in hurdle models by modeling the
> two parts (zeros and non-zeros) in one function call. The only potential
> benefit to fitting a hurdle in a single function call is that you get
> likelihood and AIC for the entire data set, but I don�t know if those are
> produced by brms.
> >
> > You could just fit a binomial model for the zero-non-zero process (i.e.
> monoculture) like
> > mod_binom=lmer((diversity>0)~block+syst+(1|plot)+(1|year)+(1|plot:year),
> data=density, family=binomial)
> >
> > and then fit a model to the log of the positive data
> > mod_gaus=lmer(log(diversity)~block+syst+(1|plot)+(1|year)+(1|plot:year),
> data=subset(density, diversity>0))
> >
> > Or, given that the outcome is non-negative and continuous, it might make
> sense to try a Tweedie distribution, but I�m not sure I�ve seen this
> applied to diversity indices in the literature. Has anyone else seen this
> done?
> > mod_twe = glmmTMB(diversity~block+syst+(1|plot)+(1|year)+(1|plot:year),
> data=density, family=tweedie)
> >
> > Cheers,
> > Mollie
> >
> >> On 18Dec 2019, at 14:45, Cesko Voeten <c.c.voeten using hum.leidenuniv.nl>
> wrote:
> >>
> >> Hi Guillaume,
> >>
> >> If you're not afraid to go Bayesian, brms can do it. Alternatively, you
> may be able to use glmmTMB and treat the hurdle part as zero inflation, but
> this is conceptually not the same thing as a hurdle model so you would need
> to judge whether that would make sense at all for your application.
> >>
> >> HTH,
> >> Cesko
> >>
> >> Op 18-12-2019 om 13:48 schreef Guillaume Adeux:
> >>> Hi everyone,
> >>> I am looking for a package which can handle "hurdle.lognormal"
> distribution
> >>> family and multiple grouping factors.
> >>> GLMMadaptive seemed as the way to go but unfortunately, to the best of
> my
> >>> knowledge, it does not handle multiple grouping factors (random
> effects).
> >>> You may ask why? I am analyzing plant diversity and one of the
> treatments
> >>> led to plots which were dominated by one species. Hence, certain
> diversity
> >>> indices are estimated as zero in these plots, and produces a mass at
> zero.
> >>> All other values are positive and continuous.
> >>> Anyone have an idea of a package/function which can handle this? Or any
> >>> alternative approach?
> >>> In lmer syntax, the model is the following:
> >>>
> mod=lmer(diversity~block+syst+(1|plot)+(1|year)+(1|plot:year),data=density,REML=F)
> >>> Thank you for your time and help.
> >>> Sincerely,
> >>> Guillaume ADEUX
> >>> [[alternative HTML version deleted]]
> >>> _______________________________________________
> >>> R-sig-mixed-models using r-project.org mailing list
> >>>
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C0f8f22d15f244f0190f208d783c3ca4b%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637122749245106064&sdata=XWAMTSylM6j81Lf%2BLb8qikDEho3QIo8l%2Bh%2BnxSpwjHM%3D&reserved=0
> >>>
> >>
> >> _______________________________________________
> >> R-sig-mixed-models using r-project.org mailing list
> >>
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C0f8f22d15f244f0190f208d783c3ca4b%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637122749245106064&sdata=XWAMTSylM6j81Lf%2BLb8qikDEho3QIo8l%2Bh%2BnxSpwjHM%3D&reserved=0
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> >
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C0f8f22d15f244f0190f208d783c3ca4b%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637122749245106064&sdata=XWAMTSylM6j81Lf%2BLb8qikDEho3QIo8l%2Bh%2BnxSpwjHM%3D&reserved=0
> >
> > [[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]]
More information about the R-sig-mixed-models
mailing list