Greetings, 1) Is there a nice way of extracting the variance estimates from an lme fit? They don't seem to be part of the lme object. 2) In a series of simulations, I am finding that with ML fitting one of my random effect variances is sometimes being estimated as essentially zero with massive CI instead of the finite value it should have, whilst using REML I get the expected value. I guess it is a numerical/optimisation problem but don't know enough about the lme fitting algorithm to know which tollerance/scale parameter to mess about with. Any suggestions where to start? Thanks, Steve. [[alternative HTML version deleted]]
There are two way to accomplish this in nlme. First try using the summary() command, which will produce all variance components and estimates for the fixed effects. Also, try the following to extract the point estimates and approximate CIs for the variance comonents.> intervals(model.lme, which="var")Harold -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch]On Behalf Of Steve Roberts Sent: Monday, April 05, 2004 3:32 PM To: r-help at stat.math.ethz.ch Cc: Steve Roberts Subject: [R] 2 lme questions Greetings, 1) Is there a nice way of extracting the variance estimates from an lme fit? They don't seem to be part of the lme object. 2) In a series of simulations, I am finding that with ML fitting one of my random effect variances is sometimes being estimated as essentially zero with massive CI instead of the finite value it should have, whilst using REML I get the expected value. I guess it is a numerical/optimisation problem but don't know enough about the lme fitting algorithm to know which tollerance/scale parameter to mess about with. Any suggestions where to start? Thanks, Steve. [[alternative HTML version deleted]] ______________________________________________ R-help at stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
For clarity, nlme produces only SDs while lme4 produces both the variance and sd. Harold -----Original Message----- From: Spencer Graves [mailto:spencer.graves@pdf.com] Sent: Mon 4/5/2004 6:54 PM To: Harold Doran Cc: Steve Roberts; r-help@stat.math.ethz.ch Subject: Re: [R] 2 lme questions The "print" method and the "summary" command display the STANDARD DEVIATIONS (not the variances) on the screen (or in a sink file). However, when I do attributes(lme(...)) and attributes(summary(lme(...))), I don't see anything I can use. Fortunately, the "interval" function produces a list, from which the variance estimates can be extracted. Consider the following example: DF <- data.frame(group=c(1,1,2,2), y=c(1, 2, 11, 12)) library(nlme) fit <- lme(y~1, random=~1|group, DF) Linear mixed-effects model fit by REML Data: DF Log-restricted-likelihood: -6.559401 Fixed: y ~ 1 (Intercept) 6.5 Random effects: Formula: ~1 | group (Intercept) Residual StdDev: 7.053597 0.7070954 Number of Observations: 4 Number of Groups: 2 > lme.int <- intervals(fit) > lme.int$reStruct^2 Error in lme.int$reStruct^2 : non-numeric argument to binary operator > lme.int$reStruct$group^2 lower est. upper sd((Intercept)) 3.06859 49.75323 806.6845 > lme.int$sigma^2 lower est. upper 0.0704255 0.4999839 3.5496217 attr(,"label") [1] "Within-group standard error:" There may be a better way; if there is, I hope someone will enlighten us all. If not, at least this works in R 1.8.1 hope this helps. spencer graves Harold Doran wrote: >There are two way to accomplish this in nlme. First try using the summary() command, which will produce all variance components and estimates for the fixed effects. Also, try the following to extract the point estimates and approximate CIs for the variance comonents. > > > >>intervals(model.lme, which="var") >> >> > >Harold > >-----Original Message----- >From: r-help-bounces@stat.math.ethz.ch >[mailto:r-help-bounces@stat.math.ethz.ch]On Behalf Of Steve Roberts >Sent: Monday, April 05, 2004 3:32 PM >To: r-help@stat.math.ethz.ch >Cc: Steve Roberts >Subject: [R] 2 lme questions > > >Greetings, > >1) Is there a nice way of extracting the variance estimates from an lme fit? They don't seem to be part of the lme object. > >2) In a series of simulations, I am finding that with ML fitting one of my random effect variances is sometimes being estimated as essentially zero with massive CI instead of the finite value it should have, whilst using REML I get the expected value. I guess it is a numerical/optimisation problem but don't know enough about the lme fitting algorithm to know which tollerance/scale parameter to mess about with. Any suggestions where to start? > >Thanks, > >Steve. > > [[alternative HTML version deleted]] > >______________________________________________ >R-help@stat.math.ethz.ch mailing list >https://www.stat.math.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://www.stat.math.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > > [[alternative HTML version deleted]]
Below are a couple functions written by Jose Pinheiro and posted to the
Nlme-help listserv (by someone else) a few years back. They will extract the
random effects matrix (for a given level) and the level-1 error.
I've used them periodically when I wanted to get at the variance components.
Hope that helps.
Dave
Dave Atkins
Assistant Professor in Clinical Psychology
Travis Research Institute
datkins at fuller.edu
# Message: 1
# Date: Thu, 22 Mar 2001 14:27:18 -0500 (EST)
# From: Kouros Owzar <owzar at email.unc.edu>
# To: Nlme-help at franz.stat.wisc.edu
# Subject: [Nlme-help] Useful Functions
#
# Dear All,
#
# Jose has been nice enough to write two useful functions
# which enable one to extract estimated variance matrices from
# lme and nlme objects.
#
# The details follow:
#
# Consider the model
#
# Y_i = g(X_i beta + Z_i b_i) + e_i
# nix1 nixp px1 nixq qX1 + n_i x1
# where
#
# Y_i : observations for the ith subject
# X_i : design matrix for the fixed effects for the ith subject
# beta : vector of fixed effects
# Z_i : design matrix for the random effects for the ith subject
# b_i : random effects vector for ith subject (~N(0,\Sigma_b)
# e_i : measurement errors for the ith subject (~N(0,\Sigma_e)
#
# The functions varRan and varWithin (shown below) extract the
# estimated matrices \hat(\Sigma_b) and \hat(\Sigma_e) of
# \Sigma_b and \Sigma_e respectively.
#
# NOTES and REMARKS:
#
# 1. I have NOT carried out extensive testing of these functions.
# So please use them with care. Some of my tests are shown
# below.
# 2. varWithin works even if one imposes structures on Sigma_e
# 3. varWithin also supports gls and gnls. The amount of testing
# done on these objects is considerably less than that done
# on lme and nlme objects.
# 4. I sincerely thank Jose once again for taking out time out of his busy
# schedule to write these functions.
#
#
# Sincerely,
#
# Kouros
#
############################### Begin Functions
#################################
varRan <- function(object, level = 1)
{
sigE <- object$sig^2
sigE*pdMatrix(object$modelStruct$reStruct)[[level]]
}
varWithin <- function(object, wch)
{
getMat <- function(wch, grps, ugrps, cS, corM, vF, std) {
wgrps <- grps[grps == ugrps[wch]]
nG <- length(wgrps)
if (!is.null(cS)) {
val <- corM[[wch]]
} else {
val <- diag(nG)
}
if (!is.null(vF)) {
std <- diag(std[[wch]])
val <- std %*% (val %*% std)
} else {
val <- std*std*val
}
val
}
grps <- getGroups(object)
if (is.null(grps)) { # gls/gnls case
grps <- rep(1, length(resid(object)))
ugrps <- "1"
} else {
ugrps <- unique(as.character(grps))
}
if (!is.null(vF <- object$modelStruct$varStruct)) {
## weights from variance function, if present
std <- split(object$sigma/varWeights(vF), grps)[ugrps]
} else {
std <- object$sigma
}
if (!is.null(cS <- object$modelStruct$corStruct)) {
corM <- corMatrix(cS)[ugrps]
} else {
corM <- NULL
}
if (!missing(wch)) { # particular group
return(getMat(wch, grps, ugrps, cS, corM, vF, std))
} else {
if (length(ugrps) == 1) { # single group
return(getMat(1, grps, ugrps, cS, corM, vF, std))
}
val <- vector("list", length(ugrps))
names(val) <- ugrps
for(i in 1:length(ugrps)) {
val[[i]] <- getMat(i, grps, ugrps, cS, corM, vF, std)
}
return(val)
}
}
############################## Some Tests
########################################
mod1<-lme(Orange)
mod2<-update(mod1,corr=corAR1())
OJ<-Orange[-1,]
mod3<-lme(OJ)
mod4<-update(mod3,corr=corARMA(p=1,q=1))
fm1 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal),data=Soybean)
fm2 <-update(fm1,corr = corAR1())
dim(varWithin(fm1))
sapply(varWithin(fm2),dim)
varWithin(fm2,1)
varWithin(fm2,2)
varWithin(fm2,9)
SB<-Soybean[-c(11,77),]
fm3 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal),data=SB)
fm4 <-update(fm3,corr = corAR1())
dim(varWithin(fm3))
sapply(varWithin(fm4),dim)
varWithin(fm4,1)
varWithin(fm4,2)
varWithin(fm4,9)
fm11 <- nlme(weight ~ SSlogis(Time, Asym, xmid,
scal),fixed=Asym+xmid+scal~1,sta
rt=c(19,55,8),data=Soybean)
fm22 <- nlme(weight ~ SSlogis(Time, Asym, xmid,
scal),fixed=Asym+xmid+scal~1,sta
rt=c(19,55,8),data=SB)
sapply(varWithin(fm11),dim)
sapply(varWithin(fm22),dim)
--__--__--
_______________________________________________
Nlme-help maillist - Nlme-help at nlme.stat.wisc.edu
http://nlme.stat.wisc.edu/mailman/listinfo/nlme-help
End of Nlme-help Digest