[R-sig-ME] mixed lognormal hurdle model with multiple grouping factors
Chen, Gang (NIH/NIMH) [E]
g@ngchen @end|ng |rom m@||@n|h@gov
Tue Dec 24 11:53:43 CET 2019
Guillaume,
To make statistical inferences for each particular effect with the Bayesian approach, first you need to directly draw the posterior samples from the output of brm():
ps <- fixef(fit, summary = FALSE)
Based on the reference (or base) level of each factor, assemble the posterior samples for each effect of interest with the above output. Then, you can get various statistical measures (mean, mode, median, standard error, quantile intervals, etc.) from those posterior samples.
The population-level effects you previously obtained from summary(fit) can basically be retrieved from
fixef(fit, summary = TRUE)
Gang
--
Gang Chen, Ph. D.
SSCC/DIRP/NIMH
National Institutes of Health
Bethesda, MD 20892-1148
U.S.A.
On 12/24/19, 4:54 AM, "Guillaume Adeux" <guillaumesimon.a2 using gmail.com> wrote:
Thank you all for continuing to provide valuable information.
Indeed, I ruled out glmmadaptive because it does not handle crossed random
effects (such a pity in my case, as it is the only downside).
My dataset is not massive (480 observations), so computation time is not a
problem.
As I specified earlier, i successfully fitted this model with brms:
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))
but this bayesian framework is completely new to me and am having a hard
time finding equivalents to (i) test of main effects and (ii) comparison of
treatments means.
Any further advice is welcome. Thanks again for your help.
GA2
Le lun. 23 déc. 2019 à 20:58, Daniel Fulop <dfulop.ucd using gmail.com> a écrit :
> Another (Tweedie) option is using the cpglmm() function in the cplm
> package.
>
> On Mon, Dec 23, 2019 at 11:06 AM Liming Wang <lmwang using gmail.com> wrote:
>
> > Dimitris,
> >
> > Thanks for the correction and clarification! I apparently didn't note
> > that Guillaume has ruled out GLMMAdaptive and I didn't recognize that you
> > are the author of GLMMAdaptive.
> >
> > Best,
> >
> > Liming
> >
> > On Mon, Dec 23, 2019 at 10:47 AM D. Rizopoulos <
> d.rizopoulos using erasmusmc.nl>
> > wrote:
> >
> > > GLMMadaptive does support random effects for the zero part model.
> > >
> > > The reason why GLMMadaptive is slower is because it implements the
> > > adaptive Gaussian quadrature for approximating the log-likelihood. This
> > > approximation is in general considered more accurate (but slower) than
> > the
> > > Laplace approximation implemented in glmmTMB.
> > >
> > > Best,
> > > Dimitris
> > >
> > > —-
> > > Dimitris Rizopoulos
> > > Professor of Biostatistics
> > > Erasmus University Medical Center
> > > The Netherlands
> > >
> > > ------------------------------
> > > *Από:* Ο χρήστης R-sig-mixed-models <
> > > r-sig-mixed-models-bounces using r-project.org> εκ μέρους του χρήστη Liming
> > > Wang <lmwang using gmail.com>
> > > *Στάλθηκε:* Δευτέρα, Δεκεμβρίου 23, 2019 19:12
> > > *Προς:* Guillaume Adeux
> > > *Κοιν.:* R-mixed models mailing list
> > > *Θέμα:* Re: [R-sig-ME] mixed lognormal hurdle model with multiple
> > > grouping factors
> > >
> > > Hi Guillaume,
> > >
> > > The mixed_model in the GLMMadaptive package supports Two-Part Mixed
> > > Effects
> > > Model for Semi-Continuous Data (
> > >
> > >
> >
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdrizopoulos.github.io%2FGLMMadaptive%2Farticles%2FZeroInflated_and_TwoPart_Models.html%23zero-inflated-negative-binomial-mixed-effects-model&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536371685&sdata=cDReHBDB%2FLuHfeM06tARqolihNCpVErMgCcMpyr6BYc%3D&reserved=0
> > ).
> > >
> > > It doesn't seem to support mixed effect for the zero outcome component
> > > though. Another downside is that, since it is written in pure R, it is
> > > much
> > > slower than glmmTMB (for similar models).
> > >
> > > Liming
> > >
> > > On Fri, Dec 20, 2019 at 6:31 AM Guillaume Adeux <
> > > guillaumesimon.a2 using gmail.com>
> > > wrote:
> > >
> > > > 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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536371685&sdata=PL6VOZ8Q9M5DQa8xp%2FrKKe%2FqXO%2BzE5sfMP%2FvxZr2jHE%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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536371685&sdata=PL6VOZ8Q9M5DQa8xp%2FrKKe%2FqXO%2BzE5sfMP%2FvxZr2jHE%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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536381676&sdata=YRc7r3Ni5P345YiPoEuNhmhmte20pMrIb3ojmkmY6qA%3D&reserved=0
> > > > > >
> > > > > > [[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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536381676&sdata=YRc7r3Ni5P345YiPoEuNhmhmte20pMrIb3ojmkmY6qA%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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536381676&sdata=YRc7r3Ni5P345YiPoEuNhmhmte20pMrIb3ojmkmY6qA%3D&reserved=0
> > > > >
> > > >
> > > > [[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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536381676&sdata=YRc7r3Ni5P345YiPoEuNhmhmte20pMrIb3ojmkmY6qA%3D&reserved=0
> > > >
> > >
> > > [[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%7C94c75233c5084fb99c6908d787d3adf4%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C1%7C637127215536381676&sdata=YRc7r3Ni5P345YiPoEuNhmhmte20pMrIb3ojmkmY6qA%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
> >
> --
> --
> Daniel Fulop, Ph.D.
> Postdoctoral Scholar
> Dept. Plant Biology, UC Davis
> Maloof Lab, Rm. 2220
> Life Sciences Addition, One Shields Ave.
> Davis, CA 95616
>
> 510-253-7462
> dfulop using ucdavis.edu
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
_______________________________________________
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