problem placing points on an effects plot for gls model

[Iannone, Basil]

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: [phone]
Fax: [phone]
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

[Alan Haynes]

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: [phone]
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: [phone]
> Fax: [phone]
> 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