Kiyoshi Sasaki
2013-Sep-22 03:12 UTC
[R] Conditioning plots (wth coplot function) with logistic regression curves
I have been trying to produce a conditional plot using coplot function
(http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/coplot.html) for a
binary response ("Presence" in my case) variable and one continuous
variable ("Overstory") given a specific levels of the other continuous
variable ("Ivy"). But, my codes produces an overlapping graph. Also, I
want to use three equal intervals for "Ivy" (i.e.,33.3 each), but I
could not figure out how. Here is my data and codes I used:
dat <- structure(list(Presence = c(0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L,
0L, 1L, 1L, 1L, 0L, 0L), Cover.size = c(7400, 2730, 2782, 1408,
5400, 4170, 630, 1120, 720, 1080, 1197, 2160, 1638, 2850, 2420,
90, 1017, 1260, 700, 1540, 100, 2720, 792, 360, 3048, 2620, 6253,
6253, 740, 2520, 3520), Ivy = c(10L, 90L, 0L, 40L, 100L, 100L,
100L, 90L, 50L, 95L, 95L, 95L, 0L, 85L, 10L, 5L, 5L, 80L, 70L,
70L, 0L, 0L, 0L, 0L, 0L, 90L, 0L, 0L, 6L, 50L, 15L), Overstory = c(6L,
4L, 4L, 4L, 3L, 6L, 5L, 4L, 5L, 8L, 8L, 8L, 8L, 8L, 6L, 6L, 9L,
12L, 12L, 7L, 12L, 16L, 16L, 16L, 16L, 13L, 14L, 14L, 14L, 12L,
12L), Moist = c(13.65, 15.3, 15.95, 17.2, 22.95, 18.25, 19.3,
18.75, 17, 21.3, 17.7, 28, 20.85, 19.85, 24.75, 22.2, 26.6, 21.15,
22, 15, 22.5, 30, 14.85, 27, 26.85, 19.55, 27.2, 27.2, 21.55,
18, 19.9), Leaf = c(95L, 95L, 95L, 90L, 80L, 90L, 95L, 95L, 90L,
80L, 80L, 80L, 95L, 90L, 95L, 95L, 95L, 100L, 100L, 95L, 5L,
80L, 80L, 80L, 80L, 90L, 95L, 95L, 95L, 85L, 90L), Prey = c(0.047,
0.051, 0.058, 0.057, 0.049, 0.046, 0.155, 0.108, 0.087, 0.093,
0.133, 0.095, 0.081, 0.104, 0.11, 0.086, 0.124, 0.134, 0.12,
0.139, 0.026, 0.045, 0.096, 0.035, 0.062, 0.065, 0.068, 0.068,
0.087, 0.052, 0.082), Abundance = c(0.141542817, 0.092071611,
0.077821012, 0.09906595, 0.113821138, 0.046202181, 0.641167367,
1.131059246, 0.346020761, 0.647748304, 1.424230915, 2.098280764,
0.23395382, 5.492325363, 1.713972302, 2.584891263, 1.049433858,
0.81246522, 1.895848332, 2.8125, 0.089866484, 0.15392509, 0.561643836,
0.400477763, 0.173869846, 1.158351129, 0.453206436, 0.50608052,
0.347007064, 0.461928178, 0.264333242)), .Names = c("Presence",
"Cover.size", "Ivy", "Overstory",
"Moist", "Leaf", "Prey", "Abundance"
), row.names = 106:136, class = "data.frame")
xrange <- seq(from=min(dat$Overstory), to=max(dat$Overstory), length=100)
panelFunc <- function(x,y,data,...){
plot(x,y,...)
fit <- glm(y ~ x, data = dat, family = binomial(link = logit), na.action =
na.exclude)
lines(xrange, inv.logit(fit$coef[1]+ fit$coef[2]*xrange))
}
coplot(dat$Presence ~ dat$Overstory | dat$Ivy, number=3, panel=panelFunc,
rows=1)
Also, if you have a better alternative code, please share with me.
Thank you in advance for your time and help!
Sincerely,
Kiyoshi
[[alternative HTML version deleted]]
Michael Friendly
2013-Sep-22 23:46 UTC
[R] Conditioning plots (wth coplot function) with logistic regression curves
On 9/21/2013 11:12 PM, Kiyoshi Sasaki wrote:> I have been trying to produce a conditional plot using coplot function (http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/coplot.html) for a binary response ("Presence" in my case) variable and one continuous variable ("Overstory") given a specific levels of the other continuous variable ("Ivy"). But, my codes produces an overlapping graph. Also, I want to use three equal intervals for "Ivy" (i.e.,33.3 each), but I could not figure out how. Here is my data and codes I used: >If you feel it hurts because you are banging your head into a wall trying to get coplot to do this, the answer is: Don't do that! Instead, you might consider using ggplot2, which handles this case nicely, as far as I can tell from your description. But first a due diligence caveat: Say you come to me for consulting on this little plotting question. I look at your data frame, dat, and I see there are a number of other variables that might explain Presence, so maybe the marginal plot that ignores them could be misleading, e.g., any of Moist, Leaf, Prey, ... could moderate the relation between overstory and presence, but you won't see that in a marginal plot. Here are a couple of quick ggplot examples, plotting classes of Ivy in the same plot frame, and on the probability scale. library(ggplot2) ggplot(dat, aes(Overstory, Presence), color=Ivy>50 ) + stat_smooth(method="glm", family=binomial, formula= y~x, alpha=0.3, aes(fill=Ivy>50)) dat$Ivy3 <- factor(cut(dat$Ivy,3)) plt <- ggplot(dat, aes(Overstory, Presence), color=Ivy3 ) + stat_smooth(method="glm", family=binomial, formula= y~x, alpha=0.3, aes(fill=Ivy3)) plt If you want separate panels, try something like plt + facet_grid(. ~ Ivy3) For plots on the logit scale, try something like logit <- function(p) log(p)/log(1-p), then plt + coord_trans(y="logit") -- Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. & Chair, Quantitative Methods York University Voice: 416 736-2100 x66249 Fax: 416 736-5814 4700 Keele Street Web: http://www.datavis.ca Toronto, ONT M3J 1P3 CANADA