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.
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-- Andrew Robinson Ph: 208 885 7115 Department of Forest Resources Fa: 208 885 6226 University of Idaho E : andrewr at uidaho.edu PO Box 441133 W : http://www.uidaho.edu/~andrewr Moscow ID 83843 Or: http://www.biometrics.uidaho.edu No statement above necessarily represents my employer's opinion.
> 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.You are aware that you need to use something like weigths= varPower (form= fitted (.)) and the plot residuals using e.g. scatter.smooth (fitted (model), resid (model, type= 'n')) Maybe the latter should also be ok with the default pearson residuals, but I am not sure. Possibly a look into the following would help? @Book{Pin:00a, author = {Pinheiro, Jose C and Bates, Douglas M}, title = {Mixed-Effects Models in {S} and {S}-{P}{L}{U}{S}}, publisher = {Springer}, year = {2000}, address = {New York} }> 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.Well, what error distribution do you have / do you expect? Regards, Lorenz - Lorenz Gygax, Dr. sc. nat. Centre for proper housing of ruminants and pigs Swiss Federal Veterinary Office
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
>>
>>
>
>
>
Christoph, ----- Original Message ----- From: Christoph Scherber <Christoph.Scherber at uni-jena.de> Date: Wednesday, January 26, 2005 7:57 pm Subject: Re: [R] lme and varFunc()> Dear all, > > I am expecting a Poisson error distribution in my lme with > weights=varFunc(). > > The "weigths= varPower (form= fitted (.))" doesn?t work due to > missing > values in the response: > > Problem in lme.formula(fixed = sqrt(nrmainaxes + 0...: Maximum > number of iterations reached without convergence. > Use traceback() to see the call stackAre you certain that this is caused by missing values in the response? It looks like something quite different to me. How did the response come to have missing values?> That?s why I?ve used one of my most important explanatory > variables as a variance covariate: > > weigths= varPower (form=~explanatory) > > With that, it worked out properly so far. > > What would your suggestion in such a case be?Does it satisfy your needs with regards to the model diagnostics? If so, then ask yourself: am I fitting a Poisson because it is an important and intrinsic part of the model that I need to construct, or am I fitting it in order to satisfy the model assumptions? If the latter then it seems that the Poisson might be unnecessary. I hope that this helps, Andrew