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]]
______________________________________________
[email protected] 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.
