[R] 2 lme questions

[Dave Atkins]

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