mixed lognormal hurdle model with multiple grouping factors

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 at 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 at r-project.org> ?? ?????? ??? ?????? Liming
Wang <lmwang at 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 at 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 at 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

on

behalf of Mollie Brooks <mollieebrooks at 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 at 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]]

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

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

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

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

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

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

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