[R] lme and varFunc()

[Christoph Scherber]

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