# [R] 2 lme questions

Dave Atkins datkins at u.washington.edu
Tue Apr 6 15:11:33 CEST 2004

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