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

[Thierry Onkelinx]

Dear Andy,

Can you post a reproducible example?

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx using inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>


Op do 31 jan. 2019 om 21:47 schreef Andy Beet via R-sig-mixed-models <
r-sig-mixed-models using r-project.org>:

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

	[[alternative HTML version deleted]]