Hello all, This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?. Is this the expected behavior? And if so, are there any reasons for this choice? Reproducible example below. Thanks, Florin library(nlme) fit.reml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="REML") (se.reml = summary(fit.reml)$tTable[,2]) ## (Intercept) Time ## 0.019457141 0.005829144 (se.reml = sqrt(diag(summary(fit.reml)$varFix))) ## (Intercept) Time ## 0.019457141 0.005829144 fit.ml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="ML") (se.ml = summary(fit.ml)$tTable[,2]) # they match the REML SE?s ## (Intercept) Time ## 0.019457141 0.005829144 (se.ml = sqrt(diag(summary(fit.ml)$varFix))) # they do not match the tTable SE?s ## (Intercept) Time ## 0.019405324 0.005813621
lme varFix under ML fit does not match coefficients standard error
8 messages · Rolf Turner, Ben Bolker, Vaida, Florin +1 more
From ?summary.lme:
adjustSigma: an optional logical value. If ?TRUE? and the estimation
method used to obtain ?object? was maximum likelihood, the
residual standard error is multiplied by sqrt(nobs/(nobs -
npar)), converting it to a REML-like estimate. This argument
is only used when a single fitted object is passed to the
function. Default is ?TRUE?.
On 2024-02-18 2:34 p.m., Vaida, Florin via R-sig-mixed-models wrote:
Hello all, This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?. Is this the expected behavior? And if so, are there any reasons for this choice? Reproducible example below. Thanks, Florin library(nlme) fit.reml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="REML") (se.reml = summary(fit.reml)$tTable[,2]) ## (Intercept) Time ## 0.019457141 0.005829144 (se.reml = sqrt(diag(summary(fit.reml)$varFix))) ## (Intercept) Time ## 0.019457141 0.005829144 fit.ml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="ML") (se.ml = summary(fit.ml)$tTable[,2]) # they match the REML SE?s ## (Intercept) Time ## 0.019457141 0.005829144 (se.ml = sqrt(diag(summary(fit.ml)$varFix))) # they do not match the tTable SE?s ## (Intercept) Time ## 0.019405324 0.005813621 [[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
On Sun, 18 Feb 2024 19:34:58 +0000 "Vaida, Florin via R-sig-mixed-models"
<r-sig-mixed-models at r-project.org> wrote:
Hello all, This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?. Is this the expected behavior? And if so, are there any reasons for this choice? Reproducible example below.
<SNIP>
I am probably displaying my ignorance and na?vet? here, but isn't this
essentially the same thing as using \hat{sigma} = \sqrt(SSE/(n-p))
(REML) rather than \hat{sigma} = \sqrt(SSE/n) (ML) in "ordinary"
regression modelling?
cheers,
Rolf Turner
Honorary Research Fellow
Department of Statistics
University of Auckland
Stats. Dep't. (secretaries) phone:
+64-9-373-7599 ext. 89622
Home phone: +64-9-480-4619
Yes, as far as I know
On Sun, Feb 18, 2024, 4:15 PM Rolf Turner <rolfturner at posteo.net> wrote:
On Sun, 18 Feb 2024 19:34:58 +0000 "Vaida, Florin via R-sig-mixed-models" <r-sig-mixed-models at r-project.org> wrote:
Hello all, This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?. Is this the expected behavior? And if so, are there any reasons for this choice? Reproducible example below.
<SNIP>
I am probably displaying my ignorance and na?vet? here, but isn't this
essentially the same thing as using \hat{sigma} = \sqrt(SSE/(n-p))
(REML) rather than \hat{sigma} = \sqrt(SSE/n) (ML) in "ordinary"
regression modelling?
cheers,
Rolf Turner
--
Honorary Research Fellow
Department of Statistics
University of Auckland
Stats. Dep't. (secretaries) phone:
+64-9-373-7599 ext. 89622
Home phone: +64-9-480-4619
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
1 day later
Thanks Ben, interesting!
Any particular reason for this default choice?
It sounds like separating parameter estimation (ML) from SE of the parameters (REML), but presumably when someone chooses ML for estimation they assume this goes to how the SE's are computed too.
Florin
?On 2/18/24, 11:52?AM, "R-sig-mixed-models on behalf of Ben Bolker" <r-sig-mixed-models-bounces at r-project.org on behalf of bbolker at gmail.com> wrote:
From ?summary.lme:
adjustSigma: an optional logical value. If ?TRUE? and the estimation
method used to obtain ?object? was maximum likelihood, the
residual standard error is multiplied by sqrt(nobs/(nobs -
npar)), converting it to a REML-like estimate. This argument
is only used when a single fitted object is passed to the
function. Default is ?TRUE?.
On 2024-02-18 2:34 p.m., Vaida, Florin via R-sig-mixed-models wrote:
> Hello all,
>
> This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?.
> Is this the expected behavior? And if so, are there any reasons for this choice?
> Reproducible example below.
>
> Thanks,
> Florin
>
> library(nlme)
> fit.reml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="REML")
> (se.reml = summary(fit.reml)$tTable[,2])
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.reml = sqrt(diag(summary(fit.reml)$varFix)))
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> fit.ml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="ML")
> (se.ml = summary(fit.ml)$tTable[,2]) # they match the REML SE?s
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.ml = sqrt(diag(summary(fit.ml)$varFix))) # they do not match the tTable SE?s
> ## (Intercept) Time
> ## 0.019405324 0.005813621
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
In principle, the standard errors should be different with REML and ML. When I try different datasets, I see differences. For example,
fm_reml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject)
fm_ml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject,
method = "ML")
coef(summary(fm_reml))
coef(summary(fm_ml))
Best,
Dimitris
-----Original Message-----
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> On Behalf Of Vaida, Florin via R-sig-mixed-models
Sent: Tuesday, February 20, 2024 9:25 PM
To: Ben Bolker <bbolker at gmail.com>; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] lme varFix under ML fit does not match coefficients standard error
Waarschuwing: Deze e-mail is afkomstig van buiten de organisatie. Klik niet op links en open geen bijlagen, tenzij u de afzender herkent en weet dat de inhoud veilig is.
Caution: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.
Thanks Ben, interesting!
Any particular reason for this default choice?
It sounds like separating parameter estimation (ML) from SE of the parameters (REML), but presumably when someone chooses ML for estimation they assume this goes to how the SE's are computed too.
Florin
?On 2/18/24, 11:52?AM, "R-sig-mixed-models on behalf of Ben Bolker" <r-sig-mixed-models-bounces at r-project.org on behalf of bbolker at gmail.com> wrote:
From ?summary.lme:
adjustSigma: an optional logical value. If ?TRUE? and the estimation
method used to obtain ?object? was maximum likelihood, the
residual standard error is multiplied by sqrt(nobs/(nobs -
npar)), converting it to a REML-like estimate. This argument
is only used when a single fitted object is passed to the
function. Default is ?TRUE?.
On 2024-02-18 2:34 p.m., Vaida, Florin via R-sig-mixed-models wrote:
> Hello all,
>
> This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?.
> Is this the expected behavior? And if so, are there any reasons for this choice?
> Reproducible example below.
>
> Thanks,
> Florin
>
> library(nlme)
> fit.reml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="REML")
> (se.reml = summary(fit.reml)$tTable[,2])
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.reml = sqrt(diag(summary(fit.reml)$varFix)))
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> fit.ml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="ML")
> (se.ml = summary(fit.ml)$tTable[,2]) # they match the REML SE?s
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.ml = sqrt(diag(summary(fit.ml)$varFix))) # they do not match the tTable SE?s
> ## (Intercept) Time
> ## 0.019405324 0.005813621
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Hi Dimitris,
The point is that the ML standard errors reported by summary() do not match those computed under varFix.
The reason for this is that indicated by Ben, which is that by default the varFix SE's get multiplied by n/(n-p), to behave more like those from REML/unbiased variances.
fm_ml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject,
method = "ML")
coef(summary(fm_ml))
sqrt(diag(summary(fm_ml)$varFix)) # ML SE's different than those reported by coef(summary(fm_ml)) above
sqrt(diag(summary(fm_ml)$varFix))*sqrt(108/106) # same SE's as in coef(summary(fm_ml))
coef(summary(fm_ml, adjustSigma=FALSE)) # same SE's as given by varFix
Best,
Florin
?On 2/20/24, 12:40?PM, "Dimitris Rizopoulos" <d.rizopoulos at erasmusmc.nl> wrote:
In principle, the standard errors should be different with REML and ML. When I try different datasets, I see differences. For example,
fm_reml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject)
fm_ml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject,
method = "ML")
coef(summary(fm_reml))
coef(summary(fm_ml))
Best,
Dimitris
-----Original Message-----
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> On Behalf Of Vaida, Florin via R-sig-mixed-models
Sent: Tuesday, February 20, 2024 9:25 PM
To: Ben Bolker <bbolker at gmail.com>; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] lme varFix under ML fit does not match coefficients standard error
Waarschuwing: Deze e-mail is afkomstig van buiten de organisatie. Klik niet op links en open geen bijlagen, tenzij u de afzender herkent en weet dat de inhoud veilig is.
Caution: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.
Thanks Ben, interesting!
Any particular reason for this default choice?
It sounds like separating parameter estimation (ML) from SE of the parameters (REML), but presumably when someone chooses ML for estimation they assume this goes to how the SE's are computed too.
Florin
On 2/18/24, 11:52?AM, "R-sig-mixed-models on behalf of Ben Bolker" <r-sig-mixed-models-bounces at r-project.org on behalf of bbolker at gmail.com> wrote:
From ?summary.lme:
adjustSigma: an optional logical value. If ?TRUE? and the estimation
method used to obtain ?object? was maximum likelihood, the
residual standard error is multiplied by sqrt(nobs/(nobs -
npar)), converting it to a REML-like estimate. This argument
is only used when a single fitted object is passed to the
function. Default is ?TRUE?.
On 2024-02-18 2:34 p.m., Vaida, Florin via R-sig-mixed-models wrote:
> Hello all,
>
> This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?.
> Is this the expected behavior? And if so, are there any reasons for this choice?
> Reproducible example below.
>
> Thanks,
> Florin
>
> library(nlme)
> fit.reml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="REML")
> (se.reml = summary(fit.reml)$tTable[,2])
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.reml = sqrt(diag(summary(fit.reml)$varFix)))
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> fit.ml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="ML")
> (se.ml = summary(fit.ml)$tTable[,2]) # they match the REML SE?s
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.ml = sqrt(diag(summary(fit.ml)$varFix))) # they do not match the tTable SE?s
> ## (Intercept) Time
> ## 0.019405324 0.005813621
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!gR6OQQ6PtYOaMYJLAUkf5UjXFoNg68ARTl6TN_il5daK7QKo-cdFJqHS3Ax5SC3bsko0LV6QZTrCEqlsCAVJO4EGR74$
Note that the initial example provided gave a singular or close-to-singular fit, so we would expect very similar answers from REML and ML.
On 2024-02-20 3:40 p.m., Dimitris Rizopoulos wrote:
In principle, the standard errors should be different with REML and ML. When I try different datasets, I see differences. For example,
fm_reml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject)
fm_ml <- lme(distance ~ age, data = Orthodont, random = ~ age | Subject,
method = "ML")
coef(summary(fm_reml))
coef(summary(fm_ml))
Best,
Dimitris
-----Original Message-----
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> On Behalf Of Vaida, Florin via R-sig-mixed-models
Sent: Tuesday, February 20, 2024 9:25 PM
To: Ben Bolker <bbolker at gmail.com>; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] lme varFix under ML fit does not match coefficients standard error
Waarschuwing: Deze e-mail is afkomstig van buiten de organisatie. Klik niet op links en open geen bijlagen, tenzij u de afzender herkent en weet dat de inhoud veilig is.
Caution: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.
Thanks Ben, interesting!
Any particular reason for this default choice?
It sounds like separating parameter estimation (ML) from SE of the parameters (REML), but presumably when someone chooses ML for estimation they assume this goes to how the SE's are computed too.
Florin
?On 2/18/24, 11:52?AM, "R-sig-mixed-models on behalf of Ben Bolker" <r-sig-mixed-models-bounces at r-project.org on behalf of bbolker at gmail.com> wrote:
From ?summary.lme:
adjustSigma: an optional logical value. If ?TRUE? and the estimation
method used to obtain ?object? was maximum likelihood, the
residual standard error is multiplied by sqrt(nobs/(nobs -
npar)), converting it to a REML-like estimate. This argument
is only used when a single fitted object is passed to the
function. Default is ?TRUE?.
On 2024-02-18 2:34 p.m., Vaida, Florin via R-sig-mixed-models wrote:
> Hello all,
>
> This is probably known, but it?s news to me: the standard errors printed for the lme model fit run under method=?ML? are in fact those computed under method=?REML?.
> Is this the expected behavior? And if so, are there any reasons for this choice?
> Reproducible example below.
>
> Thanks,
> Florin
>
> library(nlme)
> fit.reml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="REML")
> (se.reml = summary(fit.reml)$tTable[,2])
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.reml = sqrt(diag(summary(fit.reml)$varFix)))
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> fit.ml = lme(log(conc) ~ Time, random=~1|Subject, data=Glucose, na.action=na.omit, method="ML")
> (se.ml = summary(fit.ml)$tTable[,2]) # they match the REML SE?s
> ## (Intercept) Time
> ## 0.019457141 0.005829144
> (se.ml = sqrt(diag(summary(fit.ml)$varFix))) # they do not match the tTable SE?s
> ## (Intercept) Time
> ## 0.019405324 0.005813621
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!maZY0riqf0o-YADm5wrg-eNFM4LHerS92ofHlC_Pv5pnBNSjJYPepDr4S5Y16jYfN9zHl3B6SPmJz71d$
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models