[R] hierarchical model with heteroscedastic variances

[John McKown]

Please repost your message in plain text. On my gmail account, it was
totally unreadable due to the fact that the forum software strips out all
the HTML "stuff".

On Mon, Dec 1, 2014 at 7:51 PM, Rafael Moral <rafa_moral2004 at yahoo.com.br>
wrote:

> 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 at 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.
>



-- 
The temperature of the aqueous content of an unremittingly ogled
culinary vessel will not achieve 100 degrees on the Celsius scale.

Maranatha! <><
John McKown

	[[alternative HTML version deleted]]