Re: [R] lme and varFunc()
Christoph, - Original Message - From: Christoph Scherber <[EMAIL PROTECTED]> 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 stack Are 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 __ R-help@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
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 stack
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?
Regards
Christoph
[EMAIL PROTECTED] wrote:
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
__
R-help@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@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
Re: [R] lme and varFunc()
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" <[EMAIL PROTECTED]>
To:
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@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@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@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
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@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@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
RE: [R] lme and varFunc()
> 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
__
R-help@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
Re: [R] lme and varFunc()
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@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 : [EMAIL PROTECTED] PO Box 441133W : http://www.uidaho.edu/~andrewr Moscow ID 83843 Or: http://www.biometrics.uidaho.edu No statement above necessarily represents my employer's opinion. __ R-help@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
