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 [EMAIL PROTECTED]
# Message: 1 # Date: Thu, 22 Mar 2001 14:27:18 -0500 (EST) # From: Kouros Owzar <[EMAIL PROTECTED]> # To: [EMAIL PROTECTED] # 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 - [EMAIL PROTECTED] http://nlme.stat.wisc.edu/mailman/listinfo/nlme-help
End of Nlme-help Digest
______________________________________________ [EMAIL PROTECTED] 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
