[R] summary.lme: argument "adjustSigma"

[Spencer Graves]

<see inline>	

Christoph Buser wrote:
> 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'. 
> 
I would say that the documentation is ambiguous on the issue you raise. 
  Furthermore, I suspect that most people would likely expect the 
"sigma" component of summary(lme(..., method="ML")) to match that of 
lme(..., method="ML").

	  hope this helps,
	  Spencer Graves
> 
> 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
> 
> ______________________________________________
> 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