x <- c(0,1,2,3,4,5) n <- c(70,70,70,70,70,70) z <- c( 0, 9,21,47,60,63) y <- n - z b <- cbind(y, z) fitx <- glm(b ~ x, family=binomial) print(fitx) plot(x,y/n,xlim=c(0,5),ylim=c(0,1),ylab="") beta0 <- fitx$coef[1] beta1 <- fitx$coef[2] fn <- function(x) exp(beta0+beta1*x)/(1+exp(beta0+beta1*x)) par(new=T) curve(fn,0,5,ylim=c(0,1),ylab="Propn surviving")