[R-sig-ME] glmmTMB and ar1

[Ben Bolker]

  Without looking too closely at this, I think you also need to include
(1|Common.Name) in the model; otherwise the assumption is that there is
no within-group variance?

This formula is taken from
http://bbolker.github.io/mixedmodels-misc/notes/corr_braindump.html

glmmTMB_simple_fit <- glmmTMB(y~1 + (1|f) + ar1(tt-1|f),
data=d,family=gaussian)

(you used +0 rather than -1 to suppress the intercept in the ar1() term,
and you have a non-trivial fixed-effect term, but otherwise these are
similar).

  Suggestions for doc improvements/pull requests welcome ...

  Ben Bolker


On 17-11-22 05:26 AM, Jarrod Hadfield wrote:
> Hi,
> 
> It is *really* great that glmmTMB allows ar1 models. However, I'm having
> some trouble understanding the output and reconciling the estimates with
> asreml.
> 
> The data consist of the number of birds censused each year for 34 years.
> In 13 years the birds were censused twice.
> 
> The model I would like to fit has year as a continuous fixed effect, and
> then an ar1 process across years. The residual variance should pick up
> the within-year variance.
> 
>  m1.glmmTMB<-glmmTMB(log(pop)~year+ar1(year.factor+0|Common.Name),
> data=shag_data)
> 
> However, this gives one fewer parameters than I was expecting:
> 
>  summary( m1.glmmTMB)
>  Family: gaussian  ( identity )
> Formula:          log(pop) ~ year + ar1(year.factor + 0 | Common.Name)
> Data: shag_data
> 
>      AIC      BIC   logLik deviance df.resid
>     12.5     21.7     -1.2      2.5       42
> 
> Random effects:
> 
> Conditional model:
>  Groups      Name            Variance Std.Dev. Corr
>  Common.Name year.factor1973 0.214434 0.46307   (ar1)
>  Residual                    0.004099 0.06403
> Number of obs: 47, groups:  Common.Name, 1
> 
> Dispersion estimate for gaussian family (sigma^2): 0.0041
> 
> Conditional model:
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept) 55.73041   29.70301   1.876   0.0606 .
> year        -0.02464    0.01493  -1.650   0.0989 .
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> Is 0.214434 the process variance for the ar1?  But then where is the
> autocorrelation parameter?
> 
> What I hoped was the equivalent model in asreml gives different answers
> 
> m1.asreml<-asreml(log(pop)~year, random=~ar1v(year.factor), data=shag_data)
> 
> The estimate are 0.74 (autocorrelation), 0.30 (process variance) and
> 0.0041 (the residual variance). Asreml uses REML not ML so this might
> explain some of the discrepancy but I'd be surprised if it explained all.
> 
> Cheers,
> 
> Jarrod
> 
> 
> 
> 
>