[R] {nlme} Multilevel estimation heteroscedasticity

[Andrew Robinson]

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/