[R-sig-ME] Extracting variance components when using weights?
Manuel Morales
Manuel.A.Morales at williams.edu
Fri Feb 22 16:14:17 CET 2008
I have a question about how to extract variance components from an lme
object when using weights. In particular, it appears that the w/i group
SD that is returned is estimated for group 1 only. Does that mean I need
to calculate the mean of the sds from all groups to estimate the overall
within-group variance component?
An example may clarify my question:
library(nlme)
### Create data
mns <- rnorm(10, 1500, 100)
sds <- rlnorm(10, log(300), .5)
obs <- mapply(function(x, y) rnorm(20, x, y), mns, sds)
ID <- rep(1:10, each=20)
resp <- c(as.vector(obs))
### Fit models
m1 <- lme(resp~1, rand=~1|ID)
m2 <- lme(resp~1, rand=~1|ID, weights=varIdent(form=~1|ID))
anova(m1,m2) # Model 2 much better than model 1
### Extract w/i grp variance
as.numeric(VarCorr(m1)[2,2]) # w/i grp sd for m1
(s1 <- as.numeric(VarCorr(m2)[2,2])) # w/i grp sd for grp 1 of m2???
var.wts <- varWeights(m2$modelStruct$varStruct) # var weights for m2
(m2.wi.sds <- s1*1/unique(var.wts)) # w/i grp sds for all grps of m2???
(mean(m2.wi.sds)) # w/i grp sd for m2???
--
http://mutualism.williams.edu
More information about the R-sig-mixed-models
mailing list