[R-sig-ME] lme varFix under ML fit does not match coefficients standard error

[Vaida, Florin]

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 using 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 using 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 using gmail.com>; r-sig-mixed-models using 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 using r-project.org on behalf of bbolker using 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 using 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 using 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 using r-project.org mailing list
    https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!LLK065n_VXAQ!gR6OQQ6PtYOaMYJLAUkf5UjXFoNg68ARTl6TN_il5daK7QKo-cdFJqHS3Ax5SC3bsko0LV6QZTrCEqlsCAVJO4EGR74$