Dear R users,
I am analyzing data using a gls model from the nlme package The model is: m2<-gls(C~LDH + SI3, method = "ML", weights = vpDist) ## I will later fit this model using REML C = total soil carbon, LDH = log of a distance that I am interested in, and SI3 = is a index of invasion severity, but these details are not important. I am trying to plot for the effects of LDH and SI3 individually on C while controlling for the other factor. The "effects" package achieves this through the code below: plot(effect('LDH', m2)) ## effects of LDH on C, while controlling for the effects of SI3 The graph that is produced shows predicted linear relationship and a 95% CI, but does not actual show data values for C. I tried to add the points using the code below, but the problem is that the points that I add to the effects graph are not on the correct scale. That is their values are not graphed correctly. y1<-c(min(C), max(C)) x1<-c(min(LDH), max(LDH)) plot.new() plot(effect('LDH', m2), xlim = x1, ylim = y1) ## called from the "effects" package points(LDH, C, ylim = y1, xlim = x1, pch = 16, col = "black") NOTE: I use the "plot.new()" because sometimes if I do not the points do not show up on the graph. So my question is, how do I place points over an effects graph in a way that the value (x,y coordinates) of the points are correct (i.e., on the same scale as the axes of the effects graph)? I hope this makes sense. I was not sure how to attach my data file to this email so that people can try this code out for themselves, but the made up example below generates the same issue. NOTE: I tried graphing two ways: one where I specify the scale of the points and the later where I do not. Neither way seems to work. X<-rnorm(100) Y<-rnorm(100) Z<-rnorm(100) m1<-gls(Y~X+Z) summary(m1) plot.new() y1<-c(min(Y), max(Y)) x1<-c(min(X), max(X)) plot(effect('X', m1), xlim = x1, ylim = y1) points(Y~X, ylim = y1, xlim = x1, pch = 16, col = "black") plot.new() y1<-c(min(Y), max(Y)) x1<-c(min(X), max(X)) plot(effect('X', m1), xlim = x1, ylim = y1) points(Y~X, pch = 16, col = "black") Thanks a bunch for any suggestions that anyone can offer. -- Basil Iannone University of Illinois at Chicago Department of Biological Sciences (MC 066) 845 W. Taylor St. Chicago, IL 60607-7060 Email: [hidden email] Phone: 312-355-3231 Fax: 312-413-2435 http://www2.uic.edu/~bianno2 [[alternative HTML version deleted]] _______________________________________________ R-sig-ecology mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-ecology |
Hi Basil,
I think you actually want par(new=TRUE) rather than plot.new() y1<-c(min(Y), max(Y)) x1<-c(min(X), max(X)) pp <- plot(effect('X', m1), xlim = x1, ylim = y1) par(new=TRUE) plot(Y~X, xlim = x1, ylim = y1, pch = 16, col = "black") but this then has another problem - effects alters the par(mai) setting by the look of it. is you inspect pp, str(pp), you might be able to spot it. Alternatively, create an effect object, pull the relevant parts and plot it yourself. effect('X', m1)$fit # gives the y coordinates for the think line effect('X', m1)$x # gives the x coordinates effect('X', m1)$se # gives the SE values effect('X', m1)$lower effect('X', m1)$upper # give the y values from the red lines y1<-c(min(Y), max(Y)) x1<-c(min(X), max(X)) plot(effect('X', m1)$fit ~ as.matrix(effect('X', m1)$x), xlim = x1, ylim = y1, ty="l") lines(effect('X', m1)$lower ~ as.matrix(effect('X', m1)$x), lty=2, col="red) # + the same for upper points(Y~X) The as.matrix seems to be needed to convert the x coord from a list. as.data.frame, as.vector... didnt work for some reason. Not the ideal solution, but its not overly difficult and you could make a function from it easily enough if you want. HTH Alan -------------------------------------------------- Email: [hidden email] Mobile: +41794385586 Skype: aghaynes On 28 June 2012 05:34, Basil Iannone <[hidden email]> wrote: > Dear R users, > > I am analyzing data using a gls model from the nlme package > > The model is: > > m2<-gls(C~LDH + SI3, method = "ML", weights = vpDist) ## I will later fit > this model using REML > > C = total soil carbon, LDH = log of a distance that I am interested in, and > SI3 = is a index of invasion severity, but these details are not important. > > I am trying to plot for the effects of LDH and SI3 individually on C while > controlling for the other factor. The "effects" package achieves this > through the code below: > > plot(effect('LDH', m2)) ## effects of LDH on C, while controlling for the > effects of SI3 > > The graph that is produced shows predicted linear relationship and a 95% > CI, but does not actual show data values for C. > > I tried to add the points using the code below, but the problem is that the > points that I add to the effects graph are not on the correct scale. That > is their values are not graphed correctly. > > y1<-c(min(C), max(C)) > x1<-c(min(LDH), max(LDH)) > plot.new() > plot(effect('LDH', m2), xlim = x1, ylim = y1) ## called from the "effects" > package > points(LDH, C, ylim = y1, xlim = x1, pch = 16, col = "black") > > NOTE: I use the "plot.new()" because sometimes if I do not the points do > not show up on the graph. > > So my question is, how do I place points over an effects graph in a way > that the value (x,y coordinates) of the points are correct (i.e., on the > same scale as the axes of the effects graph)? > > I hope this makes sense. I was not sure how to attach my data file to this > email so that people can try this code out for themselves, but the made up > example below generates the same issue. NOTE: I tried graphing two ways: > one where I specify the scale of the points and the later where I do not. > Neither way seems to work. > > X<-rnorm(100) > Y<-rnorm(100) > Z<-rnorm(100) > > m1<-gls(Y~X+Z) > summary(m1) > > plot.new() > y1<-c(min(Y), max(Y)) > x1<-c(min(X), max(X)) > plot(effect('X', m1), xlim = x1, ylim = y1) > points(Y~X, ylim = y1, xlim = x1, pch = 16, col = "black") > > plot.new() > y1<-c(min(Y), max(Y)) > x1<-c(min(X), max(X)) > plot(effect('X', m1), xlim = x1, ylim = y1) > points(Y~X, pch = 16, col = "black") > > Thanks a bunch for any suggestions that anyone can offer. > > -- > Basil Iannone > University of Illinois at Chicago > Department of Biological Sciences (MC 066) > 845 W. Taylor St. > Chicago, IL 60607-7060 > Email: [hidden email] > Phone: 312-355-3231 > Fax: 312-413-2435 > http://www2.uic.edu/~bianno2 > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-ecology mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-ecology > [[alternative HTML version deleted]] _______________________________________________ R-sig-ecology mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-ecology |
Free forum by Nabble | Edit this page |