[R] summary.lme: argument "adjustSigma"

[Christoph Buser]

Dear Spencer 

Thank you for your answer. You are correct. The "adjustSigma"
argument is only used for the "ML" method. In the code there is: 

stdFixed <- sqrt(diag(as.matrix(object$varFix)))

if (object$method == "ML" && adjustSigma == TRUE) {
    stdFixed <- sqrt(object$dims$N/(object$dims$N - length(stdFixed))) * 
        stdFixed
}

But my question is still open:

Looking into the code above, "stdFixed" is adapted for the "ML"
method if "adjustSigma" is TRUE. Therefore the standard error
and the test results for the fixed effects is affected.

But I would expect that the residual standard error should be
adapted, too. But "object$sigma" of the lme object is unchanged
and therefore in the summary output, the estimate of the
residual standard error is identical, independent of the value
of "adjustSigma".

To my understanding this is a discrepancy to the help page that
says: 

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


Regards,

Christoph Buser

--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH Zurich	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------



Spencer Graves writes:
 > 	  It appears that the "adjustSigma" argument of 'summary.lme' does 
 > nothing with the default method, "REML".  To check, I tried the 
 > following modification of the 'summary.lme' example:
 > 
 > fm1 <- lme(distance ~ age, Orthodont, random = ~ age | Subject)
 > fm1sa <- summary(fm1)
 > fm1s <- summary(fm1, adjustSigma=FALSE)
 > fm1saa <- summary(fm1, adjustSigma=TRUE)
 > all.equal(fm1s, fm1sa)
 > TRUE
 > all.equal(fm1s, fm1saa)
 > TRUE
 > #################
 > 
 > 	  When I changed 'method' to "ML" in this example, the result suggested 
 > that the adjustment to sigma also affected the standard errors and t 
 > values for the fixed effects.  If the p-values had not been 0, they also 
 > would have been affected:
 > 
 > fm2 <- lme(distance ~ age, Orthodont, random = ~ age | Subject,
 >             method="ML")
 > (fm2sa <- summary(fm2))
 > Linear mixed-effects model fit by maximum likelihood
 >   Data: Orthodont
 >         AIC      BIC    logLik
 >    451.2116 467.3044 -219.6058
 > 
 > Random effects:
 >   Formula: ~age | Subject
 >   Structure: General positive-definite, Log-Cholesky parametrization
 >              StdDev    Corr
 > (Intercept) 2.1940998 (Intr)
 > age         0.2149245 -0.581
 > Residual    1.3100399
 > 
 > Fixed effects: distance ~ age
 >                  Value Std.Error DF   t-value p-value
 > (Intercept) 16.761111 0.7678975 80 21.827277       0
 > age          0.660185 0.0705779 80  9.353997       0
 >   Correlation:
 >      (Intr)
 > age -0.848
 > 
 > Standardized Within-Group Residuals:
 >           Min           Q1          Med           Q3          Max
 > -3.305969237 -0.487429631  0.007597973  0.482237063  3.922789795
 > 
 > Number of Observations: 108
 > Number of Groups: 27
 >  > (fm2s <- summary(fm2, adjustSigma=FALSE))
 > Linear mixed-effects model fit by maximum likelihood
 >   Data: Orthodont
 >         AIC      BIC    logLik
 >    451.2116 467.3044 -219.6058
 > 
 > Random effects:
 >   Formula: ~age | Subject
 >   Structure: General positive-definite, Log-Cholesky parametrization
 >              StdDev    Corr
 > (Intercept) 2.1940998 (Intr)
 > age         0.2149245 -0.581
 > Residual    1.3100399
 > 
 > Fixed effects: distance ~ age
 >                  Value Std.Error DF  t-value p-value
 > (Intercept) 16.761111 0.7607541 80 22.03223       0
 > age          0.660185 0.0699213 80  9.44183       0
 >   Correlation:
 >      (Intr)
 > age -0.848
 > 
 > Standardized Within-Group Residuals:
 >           Min           Q1          Med           Q3          Max
 > -3.305969237 -0.487429631  0.007597973  0.482237063  3.922789795
 > 
 > Number of Observations: 108
 > Number of Groups: 27
 >  >
 > 	  Does this answer the question?
 > 	  spencer graves
 > 	
 > Christoph Buser wrote:
 > 
 > > Dear R-list
 > > 
 > > I have a question concerning the argument "adjustSigma" in the
 > > function "lme" of the package "nlme". 
 > > 
 > > The help page says:
 > > 
 > > "the residual standard error is multiplied by sqrt(nobs/(nobs - 
 > > npar)), converting it to a REML-like estimate."
 > > 
 > > Having a look into the code I found:
 > > 
 > > stdFixed <- sqrt(diag(as.matrix(object$varFix)))
 > > 
 > > if (object$method == "ML" && adjustSigma == TRUE) {
 > >     stdFixed <- sqrt(object$dims$N/(object$dims$N - length(stdFixed))) * 
 > >         stdFixed
 > > }
 > > 
 > > tTable <- data.frame(fixed, stdFixed, object$fixDF[["X"]], 
 > >     fixed/stdFixed, fixed)
 > > 
 > > 
 > > To my understanding, only the standard error for the fixed
 > > coefficients is adapted and not the residual standard error. 
 > > 
 > > Therefore only the tTable of the output is affected by the
 > > argument "adjustSigma", but not the estimate for residual
 > > standard error (see the artificial example below). 
 > > 
 > > May someone explain to me if there is an error in my
 > > understanding of the help page and the R code? 
 > > Thank you very much.  
 > > 
 > > Best regards,
 > > 
 > > Christoph Buser
 > > 
 > > --------------------------------------------------------------
 > > Christoph Buser <buser at stat.math.ethz.ch>
 > > Seminar fuer Statistik, LEO C13
 > > ETH Zurich	8092 Zurich	 SWITZERLAND
 > > phone: x-41-44-632-4673		fax: 632-1228
 > > http://stat.ethz.ch/~buser/
 > > --------------------------------------------------------------
 > > 
 > > 
 > > Example
 > > -------
 > > 
 > > set.seed(1)
 > > dat <- data.frame(y = rnorm(16), fac1 = rep(1:4, each = 4),
 > >                   fac2 = rep(1:2,each = 8))
 > > 
 > > telme <- lme(y ~ fac1, data = dat, random = ~ 1 | fac2, method = "ML")
 > > summary(telme)
 > > summary(telme, adjustSigma = FALSE)
 > > 
 > > ______________________________________________
 > > R-help at stat.math.ethz.ch mailing list
 > > https://stat.ethz.ch/mailman/listinfo/r-help
 > > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html