All, I'm trying to create an augPred plot in the nlme library, similar to the plot on p.43 of Pinheiro & Bates (Mixed Effects Models in S and S-Plus) for their Pixel data. My data structure is the same as the example but I still get the error msg below.> comp.adj.UKV <- groupedData(adj.UKV ~ Time | Patient_no/Lisinopril,data = comp.adj.UKV.frm, order.groups = F)> fm1comp = lme(adj.UKV ~ Time + Time.sq, data = comp.adj.UKV, random list(Patient_no = ~ 1) )> plot(augPred(fm1comp, level= 1))Error in model.frame(formula, rownames, variables, varnames, extras, extranames, : variable lengths differ I've checked all the variale lengths, and have also made sure that factors are correctly defined as factors. Is there anything special I need to be doing for augPred to work correctly? I checked the help but didn't find much. cheers, dave David Afshartous, PhD University of Miami School of Business Rm KE-408 Coral Gables, FL 33124 [[alternative HTML version deleted]]
All, I've solved part of the problem below by making sure that the formula in the grouped data object is the same as the formula specified within lme (this isn't the case in the cited example from Pinheiro & Bates). However, augPred seems to plot only a linear model instead of the polynomial model. Does anyone know how to make sure that augPred plots the same model as that specified in the model (as below)? cheers, dave ________________________________ From: Afshartous, David Sent: Thursday, September 07, 2006 6:18 PM To: r-help@stat.math.ethz.ch Cc: Afshartous, David Subject: augPred plot in nlme library All, I'm trying to create an augPred plot in the nlme library, similar to the plot on p.43 of Pinheiro & Bates (Mixed Effects Models in S and S-Plus) for their Pixel data. My data structure is the same as the example but I still get the error msg below.> comp.adj.UKV <- groupedData(adj.UKV ~ Time | Patient_no/Lisinopril,data = comp.adj.UKV.frm, order.groups = F)> fm1comp = lme(adj.UKV ~ Time + Time.sq, data = comp.adj.UKV, random list(Patient_no = ~ 1) )> plot(augPred(fm1comp, level= 1))Error in model.frame(formula, rownames, variables, varnames, extras, extranames, : variable lengths differ I've checked all the variale lengths, and have also made sure that factors are correctly defined as factors. Is there anything special I need to be doing for augPred to work correctly? I checked the help but didn't find much. cheers, dave David Afshartous, PhD University of Miami School of Business Rm KE-408 Coral Gables, FL 33124 [[alternative HTML version deleted]]
Deepayan, Thanks for your suggestion. Here are more details: I have a grouped data object for repeated measures data just like the Pixel grouped data object on p.42 of Pinheiro and Bates (2000). comp.adj.UKV.3 <- groupedData(adj.UKV ~ Time | Patient_no/Lisinopril, data = comp.adj.UKV.frm, order.groups = F #labels = list(x = "Hour", y = "adj.UKV") ) i.e., the response is continuous, Time is not treated as a factor, and there exists two factors, one nested within the other (Lisinopril nested witin patient, similar to Side within Dog on p.42). I also fit a model very similar to their model: fm1comp = lme(adj.UKV ~ Time + Time.sq, data = comp.adj.UKV.3, random list(Patient_no = ~ 1 , Lisinopril = ~ 1) ) However, the command below does not produce the fitted curves from this model, but rather it seems to be the fitted curves from a linear model. plot(augPred(fm3comp)) Possibly augPred behaves differently in R than in S, but reading the R help and trying various other approaches has not solved this. Thanks! Dave -----Original Message----- From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] Sent: Friday, September 08, 2006 4:37 PM To: Afshartous, David Cc: r-help at stat.math.ethz.ch Subject: Re: [R] augPred plot in nlme library On 9/8/06, Afshartous, David <afshart at exchange.sba.miami.edu> wrote:> All, > > I've solved part of the problem below by making sure that the formula > in the grouped data object is the same as the formula specified within> lme (this isn't the case in the cited example from Pinheiro & Bates). > > However, augPred seems to plot only a linear model instead of the > polynomial model. Does anyone know how to make sure that augPred > plots the same model as that specified in the model (as below)?You are unlikely to get any helpful answers unless you give us more information, as every r-help message asks you to do:> 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 > and provide commented, minimal, self-contained, reproducible code.-Deepayan
Hi Andrew, Thanks for your email. I assume you mean age^2 instead of age.2 for fm2a, and for fm2b, I get the following error:> fm2b <- lme(distance ~ age + I(age^2), data = Orthodont)Error in lme.formula(distance ~ age + I(age^2), data = Orthodont) : iteration limit reached without convergence (9) do you get his error as well? Finally, the Pixel example on p.42 of Pinheiro & Bates gets the quadratic plot w/o using I() as you do below; is this due to a difference between S and R? thanks! dave ps - sorry for not making the data available; if anyone is interested please let me know and I'll send it directly. -----Original Message----- From: Andrew Robinson [mailto:A.Robinson at ms.unimelb.edu.au] Sent: Friday, September 08, 2006 5:46 PM To: Afshartous, David Cc: Deepayan Sarkar; r-help at stat.math.ethz.ch Subject: Re: [R] augPred plot in nlme library Hi David, this is the sort of thing that Deepayan meant. Make a dataset available to us, or use one that will be installed by default on R. eg require(nlme) fm1 <- lme(distance ~ age, data = Orthodont) plot(augPred(fm1)) # All linear fm2a <- lme(distance ~ age + age.2, data = Orthodont) plot(augPred(fm2a)) # Still linear fm2b <- lme(distance ~ age + I(age^2), data = Orthodont) plot(augPred(fm2b)) # Quadratic! I hope that this helps you resolve the problem. Andrew On Fri, Sep 08, 2006 at 05:18:13PM -0400, Afshartous, David wrote:> > Deepayan, > > Thanks for your suggestion. Here are more details: > > I have a grouped data object for repeated measures data just like the > Pixel grouped data object on p.42 of Pinheiro and Bates (2000). > > comp.adj.UKV.3 <- groupedData(adj.UKV ~ Time | Patient_no/Lisinopril, > data = comp.adj.UKV.frm, order.groups = F > #labels = list(x = "Hour", y = "adj.UKV") > ) > > i.e., the response is continuous, Time is not treated as a factor, and> there exists two factors, one nested within the other (Lisinopril > nested > > witin patient, similar to Side within Dog on p.42). > > I also fit a model very similar to their model: > > fm1comp = lme(adj.UKV ~ Time + Time.sq, data = comp.adj.UKV.3, random > = list(Patient_no = ~ 1 , Lisinopril = ~ 1) ) > > > However, the command below does not produce the fitted curves from > this model, but rather it seems to be the fitted curves from a linear > model. > > plot(augPred(fm3comp)) > > Possibly augPred behaves differently in R than in S, but reading the R> help and trying various other approaches has not solved this. > > Thanks! > Dave > > > > > -----Original Message----- > From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] > Sent: Friday, September 08, 2006 4:37 PM > To: Afshartous, David > Cc: r-help at stat.math.ethz.ch > Subject: Re: [R] augPred plot in nlme library > > On 9/8/06, Afshartous, David <afshart at exchange.sba.miami.edu> wrote: > > All, > > > > I've solved part of the problem below by making sure that the > > formula in the grouped data object is the same as the formula > > specified within > > > lme (this isn't the case in the cited example from Pinheiro &Bates).> > > > However, augPred seems to plot only a linear model instead of the > > polynomial model. Does anyone know how to make sure that augPred > > plots the same model as that specified in the model (as below)? > > You are unlikely to get any helpful answers unless you give us more > information, as every r-help message asks you to do: > > > 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 > > and provide commented, minimal, self-contained, reproducible code. > > -Deepayan > > ______________________________________________ > 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 > and provide commented, minimal, self-contained, reproducible code.-- Andrew Robinson Department of Mathematics and Statistics Tel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: a.robinson at ms.unimelb.edu.au http://www.ms.unimelb.edu.au
Thanks Deepayan and Andrew. msMaxIter solved the convergence problem and plot(augPred) works with my data when I employ I() in the function call. One other strange thing I noticed is that when I take logs of dependent variable in the function call, the plot of augPred doesn't graph any prediction line at all. contr=nlmeControl(msMaxIter = 500) fm2c <- lme(log(distance) ~ age + I(age^2), data = Orthodont, control=contr) plot(augPred(fm2c)) However, this is fixed by hard coding the dependent variable: log.dist = log(distance) fm2c <- lme(log.dist ~ age + I(age^2), data = Orthodont, control=contr) plot(augPred(fm2c)) -----Original Message----- From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] Sent: Saturday, September 09, 2006 11:46 AM To: Afshartous, David Cc: Andrew Robinson; r-help at stat.math.ethz.ch Subject: Re: [R] augPred plot in nlme library On 9/9/06, Afshartous, David <afshart at exchange.sba.miami.edu> wrote:> Hi Andrew, > > Thanks for your email. I assume you mean age^2 instead of age.2 for > fm2a, and for fm2b, I get the following error: > > > fm2b <- lme(distance ~ age + I(age^2), data = Orthodont) > Error in lme.formula(distance ~ age + I(age^2), data = Orthodont) : > iteration limit reached without convergence (9) > > do you get his error as well?For me, adding 'control = list(msMaxIter = 500)' worked. I'm writing from memory, so the name may not be exactly right, see ?nlmeControl.> Finally, the Pixel example on p.42 of Pinheiro & Bates gets the > quadratic plot w/o using I() as you do below; is this due to a > difference between S and R?Yes.> > thanks! > dave > > ps - sorry for not making the data available; if anyone is interested > please let me know and I'll send it directly. > > > > > -----Original Message----- > From: Andrew Robinson [mailto:A.Robinson at ms.unimelb.edu.au] > Sent: Friday, September 08, 2006 5:46 PM > To: Afshartous, David > Cc: Deepayan Sarkar; r-help at stat.math.ethz.ch > Subject: Re: [R] augPred plot in nlme library > > Hi David, > > this is the sort of thing that Deepayan meant. Make a dataset > available to us, or use one that will be installed by default on R. > > eg > > require(nlme) > fm1 <- lme(distance ~ age, data = Orthodont) > plot(augPred(fm1)) > > # All linear > > fm2a <- lme(distance ~ age + age.2, data = Orthodont) > plot(augPred(fm2a)) > > # Still linear > > fm2b <- lme(distance ~ age + I(age^2), data = Orthodont) > plot(augPred(fm2b)) > > # Quadratic! > > I hope that this helps you resolve the problem. > > Andrew > > > > On Fri, Sep 08, 2006 at 05:18:13PM -0400, Afshartous, David wrote: > > > > Deepayan, > > > > Thanks for your suggestion. Here are more details: > > > > I have a grouped data object for repeated measures data just like > > the Pixel grouped data object on p.42 of Pinheiro and Bates (2000). > > > > comp.adj.UKV.3 <- groupedData(adj.UKV ~ Time |Patient_no/Lisinopril,> > data = comp.adj.UKV.frm, order.groups = F > > #labels = list(x = "Hour", y = "adj.UKV") > > ) > > > > i.e., the response is continuous, Time is not treated as a factor, > > and > > > there exists two factors, one nested within the other (Lisinopril > > nested > > > > witin patient, similar to Side within Dog on p.42). > > > > I also fit a model very similar to their model: > > > > fm1comp = lme(adj.UKV ~ Time + Time.sq, data = comp.adj.UKV.3, > > random = list(Patient_no = ~ 1 , Lisinopril = ~ 1) ) > > > > > > However, the command below does not produce the fitted curves from > > this model, but rather it seems to be the fitted curves from a > > linear model. > > > > plot(augPred(fm3comp)) > > > > Possibly augPred behaves differently in R than in S, but reading the> > R > > > help and trying various other approaches has not solved this. > > > > Thanks! > > Dave > > > > > > > > > > -----Original Message----- > > From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] > > Sent: Friday, September 08, 2006 4:37 PM > > To: Afshartous, David > > Cc: r-help at stat.math.ethz.ch > > Subject: Re: [R] augPred plot in nlme library > > > > On 9/8/06, Afshartous, David <afshart at exchange.sba.miami.edu> wrote: > > > All, > > > > > > I've solved part of the problem below by making sure that the > > > formula in the grouped data object is the same as the formula > > > specified within > > > > > lme (this isn't the case in the cited example from Pinheiro & > Bates). > > > > > > However, augPred seems to plot only a linear model instead of the > > > polynomial model. Does anyone know how to make sure that augPred > > > plots the same model as that specified in the model (as below)? > > > > You are unlikely to get any helpful answers unless you give us more > > information, as every r-help message asks you to do: > > > > > 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 > > > and provide commented, minimal, self-contained, reproducible code. > > > > -Deepayan
Hi Peter, Thanks for the email. No, plot(augPred(fm2c)) still graphs the original points but doesn't graph the prediction curve. However, this is solved by hardcording the dependent variable as below. yes, I attached Orthodont. cheers, dave -----Original Message----- From: Petr Pikal [mailto:petr.pikal at precheza.cz] Sent: Monday, September 11, 2006 2:52 AM To: Afshartous, David; r-help at stat.math.ethz.ch Subject: Re: [R] augPred plot in nlme library Hi please try not to hide an information> plot(augPred(fm2c))Error in log(distance) : object "distance" not found Is it what you have got and what you mean by "does not graph any prediction at all"? If not did you attached Orthodont before? I suppose plot.augPred probably expects the same name for original and and fitted data. You can go through source code to see what happens by nlme:::plot.augPred HTH Petr On 9 Sep 2006 at 13:39, Afshartous, David wrote: Date sent: Sat, 9 Sep 2006 13:39:04 -0400 From: "Afshartous, David" <afshart at exchange.sba.miami.edu> To: "Deepayan Sarkar" <deepayan.sarkar at gmail.com> Copies to: r-help at stat.math.ethz.ch Subject: Re: [R] augPred plot in nlme library> > Thanks Deepayan and Andrew. > > msMaxIter solved the convergence problem and plot(augPred) works with > my data when I employ I() in the function call. > > One other strange thing I noticed is that when I take logs of > dependent variable in the function call, the plot of augPred > doesn't graph any prediction line at all. > > contr=nlmeControl(msMaxIter = 500) > fm2c <- lme(log(distance) ~ age + I(age^2), data = Orthodont, > control=contr) > plot(augPred(fm2c)) > > However, this is fixed by hard coding the dependent variable: > > log.dist = log(distance) > fm2c <- lme(log.dist ~ age + I(age^2), data = Orthodont, > control=contr) plot(augPred(fm2c)) > > > > -----Original Message----- > From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] > Sent: Saturday, September 09, 2006 11:46 AM > To: Afshartous, David > Cc: Andrew Robinson; r-help at stat.math.ethz.ch > Subject: Re: [R] augPred plot in nlme library > > On 9/9/06, Afshartous, David <afshart at exchange.sba.miami.edu> wrote: > > Hi Andrew, > > Thanks for your email. I assume you mean age^2 instead > of age.2 for > fm2a, and for fm2b, I get the following error: > > > > fm2b <- lme(distance ~ age + I(age^2), data = Orthodont) > Error in > lme.formula(distance ~ age + I(age^2), data = Orthodont) : > > iteration limit reached without convergence (9) > > do you get his > error as well? > > For me, adding 'control = list(msMaxIter = 500)' worked. I'm writing > from memory, so the name may not be exactly right, see ?nlmeControl. > > > Finally, the Pixel example on p.42 of Pinheiro & Bates gets the > > quadratic plot w/o using I() as you do below; is this due to a > > difference between S and R? > > Yes. > > > > > thanks! > > dave > > > > ps - sorry for not making the data available; if anyone is > > interested please let me know and I'll send it directly. > > > > > > > > > > -----Original Message----- > > From: Andrew Robinson [mailto:A.Robinson at ms.unimelb.edu.au] > > Sent: Friday, September 08, 2006 5:46 PM > > To: Afshartous, David > > Cc: Deepayan Sarkar; r-help at stat.math.ethz.ch > > Subject: Re: [R] augPred plot in nlme library > > > > Hi David, > > > > this is the sort of thing that Deepayan meant. Make a dataset > > available to us, or use one that will be installed by default on R. > > > > eg > > > > require(nlme) > > fm1 <- lme(distance ~ age, data = Orthodont) > > plot(augPred(fm1)) > > > > # All linear > > > > fm2a <- lme(distance ~ age + age.2, data = Orthodont) > > plot(augPred(fm2a)) > > > > # Still linear > > > > fm2b <- lme(distance ~ age + I(age^2), data = Orthodont) > > plot(augPred(fm2b)) > > > > # Quadratic! > > > > I hope that this helps you resolve the problem. > > > > Andrew > > > > > > > > On Fri, Sep 08, 2006 at 05:18:13PM -0400, Afshartous, David wrote: > > > > > > Deepayan, > > > > > > Thanks for your suggestion. Here are more details: > > > > > > I have a grouped data object for repeated measures data just like > > > the Pixel grouped data object on p.42 of Pinheiro and Bates > > > (2000). > > > > > > comp.adj.UKV.3 <- groupedData(adj.UKV ~ Time | > Patient_no/Lisinopril, > > > data = comp.adj.UKV.frm, order.groups = F > > > #labels = list(x = "Hour", y = "adj.UKV") > > > ) > > > > > > i.e., the response is continuous, Time is not treated as a factor, > > > and > > > > > there exists two factors, one nested within the other (Lisinopril > > > nested > > > > > > witin patient, similar to Side within Dog on p.42). > > > > > > I also fit a model very similar to their model: > > > > > > fm1comp = lme(adj.UKV ~ Time + Time.sq, data = comp.adj.UKV.3, > > > random = list(Patient_no = ~ 1 , Lisinopril = ~ 1) ) > > > > > > > > > However, the command below does not produce the fitted curves from > > > this model, but rather it seems to be the fitted curves from a > > > linear model. > > > > > > plot(augPred(fm3comp)) > > > > > > Possibly augPred behaves differently in R than in S, but reading > > > the > > > > R > > > > > help and trying various other approaches has not solved this. > > > > > > Thanks! > > > Dave > > > > > > > > > > > > > > > -----Original Message----- > > > From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] > > > Sent: Friday, September 08, 2006 4:37 PM > > > To: Afshartous, David > > > Cc: r-help at stat.math.ethz.ch > > > Subject: Re: [R] augPred plot in nlme library > > > > > > On 9/8/06, Afshartous, David <afshart at exchange.sba.miami.edu> > > > wrote: > > > > All, > > > > > > > > I've solved part of the problem below by making sure that the > > > > formula in the grouped data object is the same as the formula > > > > specified within > > > > > > > lme (this isn't the case in the cited example from Pinheiro & > > Bates). > > > > > > > > However, augPred seems to plot only a linear model instead of > > > > the polynomial model. Does anyone know how to make sure that > > > > augPred plots the same model as that specified in the model (as > > > > below)? > > > > > > You are unlikely to get any helpful answers unless you give us > > > more information, as every r-help message asks you to do: > > > > > > > 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 > > > > and provide commented, minimal, self-contained, reproducible > > > > code. > > > > > > -Deepayan > > ______________________________________________ > 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 and provide commented, > minimal, self-contained, reproducible code.Petr Pikal petr.pikal at precheza.cz