[R-sig-ME] How to specify correlation structure along with modeling heteroskedasticity in glmmTMB?

[Ben Bolker]

On 2023-10-08 11:32 a.m., Santosh Srinivas wrote:
> Hello list members,
> 
> For one of our papers, we are trying to understand if there is an equivalent of the following nlme::lme() specification in glmmTMB:
> m = lme(y ~ gender + time, data = df, random= ~1+time|id, correlation = corAR1(), weights = varExp(form = ~time))
> 
> In our case, time = 0, 1, 2, ....8, where year 2008 was converted to 0, 2009 to 1, and so on for all individuals. Not all individuals have data on all timepoints. They also may start at different timepoints.
> 
> ------------------------------------
> id --- y --- time --- gender
> ------------------------------------
> ab --- 0.2 --- 1 --- 2
> ab --- 0.3 --- 2 --- 2
> ac --- 0.5 --- 0 --- 1
> ac --- 0.2 --- 1 --- 1
> ac --- 0.7 --- 3 --- 1
> bd --- 0.3 --- 2 --- 2
> bd --- 0.2 --- 5 --- 2
> bd --- 0.3 --- 8 --- 2
> :
> :
> 
> In essence, we are wondering if it is possible to include a) first-order autoregressive structure, along with b) modeling heteroscedasticity as an exponent of the covariate instead of as a power function of the covariate in glmmTMB.
> 
>  From the documentation, it appears that the below specification achieves the purpose:
> m_alt = glmmTMB(y ~ gender + time + ar1(0 + factor(time)|id) + (1|id), data= df, dispformula = ~factor(time))
> 
> Not sure if our understanding of the model specification is right.
> 
> 
>    1.  Is ar1() as specified here equivalent to corAr1() above?
>    2.  Is dispformula = ~factor(time) equivalent to weights = varExp(form = ~time)?
>    3.  We are also unsure if we are messing up with factorizing time.


   Almost.

   1. Yes.
   2. I don't think so. I think it should be dispformula = ~ time; the 
dispersion model uses a log link, so this should correspond to var = 
exp(a + b*time) = v0*exp(b*time), which should be the same as varExp
   3. You should be using time as a factor in the ar1() specification 
but numeric in the dispformula specification.

   When in doubt, I'd recommend fitting both packages/formulations to a 
small made-up data set to see if they give answers that are identical 
(up to rounding errors/variation)

> 
> Request your expertise and help.
> 
> Thank you.
> sbs
> 
> 
> 
> 
> 
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

-- 
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
(Acting) Graduate chair, Mathematics & Statistics
 > E-mail is sent at my convenience; I don't expect replies outside of 
working hours.