[R-sig-ME] mixed lognormal hurdle model with multiple grouping factors

[Daniel Fulop]

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]]