Dear useRs,I have been wondering whether it would be possible to fit a linear
mixed model including heteroscedastic variances for a 2-level hierarchical
study with subsampling.Suppose that I had three levels for the first
hierarchical level and 4 for the second, with 2 subsamples, e.g.
set.seed(2014)y <- rnorm(24, 20, 2)level1 <- gl(3, 8)level2 <- gl(4, 2)my.data
<- data.frame(y, level1, level2)
Then, I can easily fit a model with nested random effects, i.e.
y_{ijk} = \mu + \alpha_i + \beta_{ij} + \epsilon_{ijk}
with \alpha_i ~ N(0, \sigma^2_1) the random effect for level 1, \beta_{ij} ~
N(0, \sigma^2_2) the random effect for level 2 and \epsilon_{ijk} ~ N(0,
sigma^2) the error,
using the following code
require(nlme)fit1 <- lme(y ~ 1, random=~1|level1/level2, my.data)
# which is equivalent tofit1 <- lme(y ~ 1, random=list(level1=pdDiag(~1),
level2=pdDiag(~1)), my.data)
However, I would like to have different variances per each "level1" level
(model i) and then a different model considering different variances per each
"level1" and per each "level2" level within "level1" (model ii).So we would have
Model (i):\alpha_i ~ N(0, \sigma^2_1_i), \beta_{ij} ~ N(0, \sigma^2_2) and
\epsilon_{ijk} ~ N(0, sigma^2),
that is, different sigma^2_1 per "level1" level
Model (ii):\alpha_i ~ N(0, \sigma^2_1_i), \beta_{ij} ~ N(0, \sigma^2_2_{ij})
and \epsilon_{ijk} ~ N(0, sigma^2),
that is, different sigma^2_1 per "level1" level and different sigma^2_2 per
"level1:level2" level combination
I tried the following for model (i):
fit2 <- lme(y ~ 1, random=list(level1=pdDiag(~level1)), my.data)
and this for model (ii):
fit3 <- lme(y ~ 1, random=list(level1=pdDiag(~level1),
level2=pdDiag(~level1:level2)), my.data)
This second fit also gives a warning, and I don't believe that these are doing
what I intend to.I also tried using the weights argument, i.e., for model (i)
fit2.2 <- lme(y ~ 1, random=~1|level1/level2, weights=varIdent(form=~1|level1),
my.data)
Perhaps this is closer to what I'm trying to do, but when it comes to model
(ii) I can't properly work the syntax, as the "/" creates an error, so I tried
fit2.3 <- lme(y ~ 1, random=~1|level1/level2,
weights=varIdent(form=~1|level1*level2), my.data)
but I'm still doubtful.
Any suggestions on how I might be able to fit these models?
Best wishes,Rafael.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
