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 <[EMAIL PROTECTED]>
> 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@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@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
