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