[R] ggplot2: stat_smooth for family=binomial with cbind(Y, N) formula

ONKELINX, Thierry Thierry.ONKELINX at inbo.be
Tue Dec 17 17:26:38 CET 2013


Dear Michael,

Calculate the propotions. Then it is easy to use the weight option of glm

data("SpaceShuttle", package="vcd")
SpaceShuttle$trials <- 6

fm <- glm(cbind(nFailures, 6 - nFailures) ~ Temperature, data = SpaceShuttle, family = binomial)
fm2 <- glm(nFailures/trials ~ Temperature, data = SpaceShuttle, family = binomial, weight = trials)
all.equal(coef(fm), coef(fm2))

ggplot(SpaceShuttle, aes(x = Temperature, y = nFailures / trials)) + geom_point() + geom_smooth(method = "glm", family = binomial, aes(weight = trials))

Best regards,

Thierry

-----Oorspronkelijk bericht-----
Van: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org] Namens Michael Friendly
Verzonden: dinsdag 17 december 2013 14:58
Aan: R-help
Onderwerp: [R] ggplot2: stat_smooth for family=binomial with cbind(Y, N) formula

With ggplot2, I can plot the glm stat_smooth for binomial data when the response is binary or a two-level factor as follows:

data("Donner", package="vcdExtra")
ggplot(Donner, aes(age, survived)) +
geom_point(position = position_jitter(height = 0.02, width = 0)) + stat_smooth(method = "glm", family = binomial, formula = y ~ x, alpha = 0.2, size=2)

But how can I specify the formula for stat_smooth when the response is cbind(successes, failures)?
The equivalent with plot (minus the confidence band) for the example I want is:

data("SpaceShuttle", package="vcd")

 > head(SpaceShuttle, 5)
   FlightNumber Temperature Pressure Fail nFailures Damage
1            1          66       50   no         0      0
2            2          70       50  yes         1      4
3            3          69       50   no         0      0
4            4          80       50 <NA>        NA     NA
5            5          68       50   no         0      0
 >

plot(nFailures/6 ~ Temperature, data = SpaceShuttle,
      xlim = c(30, 81), ylim = c(0,1),
      main = "NASA Space Shuttle O-Ring Failures",
      ylab = "Estimated failure probability",
      xlab = "Temperature (degrees F)",
      pch = 19, col = "blue", cex=1.2)
fm <- glm(cbind(nFailures, 6 - nFailures) ~ Temperature,
           data = SpaceShuttle,
           family = binomial)
pred <- predict(fm, data.frame(Temperature = 30 : 81), se=TRUE)
lines(30 : 81,
       predict(fm, data.frame(Temperature = 30 : 81), type = "response"),
       lwd = 3)

--
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

______________________________________________
R-help op r-project.org 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.
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.



More information about the R-help mailing list