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: <[email protected]>
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.
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
