[R-sig-ME] nlme::gls potential bug

[Ben Bolker]

  How different are the answers? Is it possibly an issue of dividing
by n vs. (n-1) ?

On Thu, Jan 31, 2019 at 3:47 PM Andy Beet via R-sig-mixed-models
<r-sig-mixed-models using r-project.org> wrote:
>
> Hi there,
>
> I was referred to this group. I hope someone can help out.
>
> I have been using the nlme::gls package in R to fit a pretty
> simple model (linear with AR error)
>
> y(t) = beta*x(t) + e(t)              where e(t) ~ rho*e(t-1) + Z(t)
>       and Z(t)~ N(0,sig^2)
>
> I call the R routine
>
> glsObj <- nlme::gls(y ~ x -1, data=data, correlation =
> nlme::corAR1(form= ~x), method="ML")
>
> All seems fine.
>
>
> In addition, I have also coded the likelihood myself and maximized it
> for beta, rho and sigma.
>
> I get the exact same estimates of beta and rho, (as nlme::gls) but the
> estimate of sigma is not the same and i can not figure out why.
>
> The maximum likelihood estimator for sigma under this model is
>
> sig^2 = (( 1-rho^2)u(1)^2 + sum((u(t)- rho*u(t-1))^2)/n
>
> where the sum is t=2,...,n and
>
> u(t) = y(t) - X(t)*beta
>
>
> I have read the mixed-effects models in S and S-Plus book (nlme::gls
> code is based directly on this) and this problem is specified on page
> 204 eq (5.5). I have also calculated sigma based on (5.7) -after the
> transformation documented (5.2) -and i do not get the same value as
> either the package or my implementation.
>
> Any advice would be most welcomed. Is there a bug in the estimation of
> sigma in this package?
>
> Thanks
>
> Andy
>
> --
> Andy Beet
> Ecosystem Dynamics & Assessment Branch
> Northeast Fisheries Science Center
> NOAA Fisheries Service
> 166 Water Street
> Woods Hole, MA 02543
> tel: 508-495-2073
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models