[R] lme and varFunc()

[Dimitris Rizopoulos]

Hi Chris,

You could perform a graphical check before deciding which variance 
function is reasonable to use. For example, in your case maybe 
something like:

plot(model1, resid(., type="p")~Block)

would have shown that the variability depends on `Block' (note: 
`Block' sounds like a categorical variable, if so probably you could 
also consider `varIdent()')

I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat
     http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Christoph Scherber" <Christoph.Scherber at uni-jena.de>
To: <r-help at stat.math.ethz.ch>
Sent: Tuesday, January 25, 2005 9:57 AM
Subject: Re: [R] lme and varFunc()


> Dear all,
>
> Regarding the lme with varFunc() question I posted a few days ago: I 
> have used the following two approaches:
>
> model1<-lme(response~Covariate+Block+TreatmentA+TreatmentB,random=~1|Plot/Subplot,method="ML")
> model2a<-update(model1,weights=varPower(form=~ fitted(.)))
> model2b<-update(model1,weights=varPower(form=~block))
>
> While model2a produces an error
> "Problem in .C("mixed_loglik",: subroutine mixed_loglik: Missing 
> values in argument 1 Use traceback() to see the call stack"
>
> Model 2b seems to work fine, now.
>
> I´m not sure why model2a doesn´t work, but using an important 
> explanatory variable (block) as a variance covariate seems to do a 
> better job (although I don´t really understand why)
>
> Does anyone have an explanation for this?
>
> Regards,
> Chris.
>
>
>
>
>
> Andrew Robinson wrote:
>
>>Dear Christoph,
>>
>>what command are you using to plot the residuals?  If you use the
>>default residuals it will not reflect the variance model.  If you 
>>use
>>the argument
>>
>>type="p"
>>
>>then you get the Pearson residuals, which will reflect the weights
>>model.  Try something like this:
>>
>>plot(model, resid(., type = "p") ~ fitted(.), abline = 0)
>>
>>I hope that this helps,
>>
>>Andrew
>>
>>On Mon, Jan 24, 2005 at 02:28:44PM +0100, Christoph Scherber wrote:
>>
>>>Dear R users,
>>>
>>>I am currently analyzing a dataset using lme(). The model I use has 
>>>the following structure:
>>>
>>>model<-lme(response~Covariate+TreatmentA+TreatmentB,random=~1|Block/Plot,method="ML")
>>>
>>>When I plot the residuals against the fitted values, I see a clear 
>>>positive trend (meaning that the variance increases with the mean).
>>>
>>>I tried to solve this issue using weights=varPower(), but it 
>>>doesn?t change the residual plot at all.
>>>
>>>How would you implement such a positive trend in the variance? I?ve 
>>>tried glmmPQL (which works great with poisson errors), but using 
>>>glmmPQL I can?t do model simplification.
>>>
>>>Many thanks for your help!
>>>
>>>Regards
>>>Chris.
>>>
>>>______________________________________________
>>>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
>>>
>>
>>
>
> ______________________________________________
> 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
>