[R] hierarchical model with heteroscedastic variances

[Rafael Moral]

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]]