Rense Nieuwenhuis
2007-Jun-10 14:35 UTC
[R] {nlme} Multilevel estimation heteroscedasticity
Dear All,
I'm trying to model heteroscedasticity using a multilevel model. To
do so, I make use of the nlme package and the weigths-parameter.
Let's say that I hypothesize that the exam score of students
(normexam) is influenced by their score on a standardized LR test
(standLRT). Students are of course nested in "schools". These
variables are contained in the Exam-data in the mlmRev package.
library(nlme)
library(mlmRev)
lme(fixed = normexam ~ standLRT,
data = Exam,
random = ~ 1 | school)
If I want to model only a few categories of variance, all works fine.
For instance, should I (for whatever reason) hypothesize that the
variance on the normexam-scores is larger in mixed schools than in
boys-schools, I'd use weights = varIdent(form = ~ 1 | type), leading to:
heteroscedastic <- lme(fixed = normexam ~ standLRT,
data = Exam,
weights = varIdent(form = ~ 1 | type),
random = ~ 1 | school)
This gives me nice and clear output, part of which is shown below:
Variance function:
Structure: Different standard deviations per stratum
Formula: ~normexam | type
Parameter estimates:
Mxd Sngl
1.000000 1.034607
Number of Observations: 4059
Number of Groups: 65
Though, should I hypothesize that the variance on the normexam-
variable is larger on schools that have a higher average score on
intake-exams (schavg), I run into troubles. I'd use weights = varIdent
(form = ~ 1 | schavg), leading to:
heteroscedastic <- lme(fixed = normexam ~ standLRT,
data = Exam,
weights = varIdent(form = ~ 1 | schavg),
random = ~ 1 | school)
This leads to estimation problems. R tells me:
Error in lme.formula(fixed = normexam ~ standLRT, data = Exam,
weights = varIdent(form = ~1 | :
nlminb problem, convergence error code = 1; message = iteration
limit reached without convergence (9)
Fiddling with maxiter and setting an unreasonable tolerance doesn't
help. I think the origin of this problem lies within the large number
of categories on "schavg" (65), that may make estimation troublesome.
This leads to my two questions:
- How to solve this estimation-problem?
- Is is possible that the varIdent (or more general: VarFunc) of lme
returns a single value, representing a coëfficiënt along which
variance is increasing / decreasing?
- In general: how can a variance-component / heteroscedasticity be
made dependent on some level-2 variable (school level in my examples) ?
Many thanks in advance,
Rense Nieuwenhuis
[[alternative HTML version deleted]]
Rense, how about weights = varPower(form = ~ schavg) or weights = varConstPower(form = ~ schavg) or even weights = varPower(form = ~ schavg | type) Yuo might find Pinheiro and Bates (2000) to be a valuable investment. I hope that this helps, Andrew On Sun, Jun 10, 2007 at 04:35:58PM +0200, Rense Nieuwenhuis wrote:> Dear All, > > I'm trying to model heteroscedasticity using a multilevel model. To > do so, I make use of the nlme package and the weigths-parameter. > Let's say that I hypothesize that the exam score of students > (normexam) is influenced by their score on a standardized LR test > (standLRT). Students are of course nested in "schools". These > variables are contained in the Exam-data in the mlmRev package. > > library(nlme) > library(mlmRev) > lme(fixed = normexam ~ standLRT, > data = Exam, > random = ~ 1 | school) > > > If I want to model only a few categories of variance, all works fine. > For instance, should I (for whatever reason) hypothesize that the > variance on the normexam-scores is larger in mixed schools than in > boys-schools, I'd use weights = varIdent(form = ~ 1 | type), leading to: > > heteroscedastic <- lme(fixed = normexam ~ standLRT, > data = Exam, > weights = varIdent(form = ~ 1 | type), > random = ~ 1 | school) > > This gives me nice and clear output, part of which is shown below: > Variance function: > Structure: Different standard deviations per stratum > Formula: ~normexam | type > Parameter estimates: > Mxd Sngl > 1.000000 1.034607 > Number of Observations: 4059 > Number of Groups: 65 > > > Though, should I hypothesize that the variance on the normexam- > variable is larger on schools that have a higher average score on > intake-exams (schavg), I run into troubles. I'd use weights = varIdent > (form = ~ 1 | schavg), leading to: > > heteroscedastic <- lme(fixed = normexam ~ standLRT, > data = Exam, > weights = varIdent(form = ~ 1 | schavg), > random = ~ 1 | school) > > This leads to estimation problems. R tells me: > Error in lme.formula(fixed = normexam ~ standLRT, data = Exam, > weights = varIdent(form = ~1 | : > nlminb problem, convergence error code = 1; message = iteration > limit reached without convergence (9) > > Fiddling with maxiter and setting an unreasonable tolerance doesn't > help. I think the origin of this problem lies within the large number > of categories on "schavg" (65), that may make estimation troublesome. > > This leads to my two questions: > - How to solve this estimation-problem? > - Is is possible that the varIdent (or more general: VarFunc) of lme > returns a single value, representing a co?ffici?nt along which > variance is increasing / decreasing? > > - In general: how can a variance-component / heteroscedasticity be > made dependent on some level-2 variable (school level in my examples) ? > > Many thanks in advance, > > Rense Nieuwenhuis > [[alternative HTML version deleted]] >> ______________________________________________ > R-help at stat.math.ethz.ch mailing list > 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.-- Andrew Robinson Department of Mathematics and Statistics Tel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/