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

