# This Splus code fits a logisitic regression model # to the Illinois turtle data (see example 3.4 in # Chapter 3. # First enter the data temp<-c(27.2, 27.7, 28.3, 28.4, 29.9) nmale <- c(2, 17, 26, 19, 27) nfemale <- c(25, 7, 4, 8, 1) ntur <- nmale+nfemale pmale <- nmale/ntur # Fit a logisitic regression model using # maximum likelihood estimation. Here the # response is the proportion of males and # the weight is the number of turtle eggs # exposed to a certain temperature. tur.logit <- glm(pmale ~ temp, family=binomial, weights=ntur) summary(tur.logit) anova(tur.logit, test="Chisq") # Make a plot of the data along with the # estimated curve. # The par( ) function is used to specify # features of the plot such as the shape # and size of plotting symbols and the size # of the plot. The following options are # specified # fin=c(w,h) requests a plot that us w inches # wide and h inches high # pch=18 requests a filled rectangle as a # plotting symbol # mkh=b requests plotting symbols that # that are b inched high # cex=1.2 increases size of printed letters # lwd=3 controls widths of lines and curves par(fin=c(7,7), pch=18, mkh=.15, mex=1.5, cex=1.2, lwd=3) # Now plot the data points plot(temp, pmale, xlab="Temperature (C)", ylab="Proportion Male", main="Proportion of Male Turtles vs. Incubation Temperature") # Now compute points on the estimated curve tempnew <- seq(27.0,30.0,0.1) logitnew <- tur.logit\$coef[1]+tur.logit\$coef[2]*tempnew pihatnew <- exp(logitnew)/(1+exp(logitnew)) # Add the curve to the previous plot lines(tempnew, pihatnew) # Compute confidence interals for parameters using # standard large sample theory. # First add the MASS library to obtain the variance # covariance matrix. library(MASS) # Compute the lower and upper limits of the CI. est.b <- function(obj, level=0.95) { b <- coef(obj) stderr <- sqrt(diag(vcov(obj))) low <- b - qnorm(1-(1-level)/2)*stderr high <- b+ qnorm(1-(1-level)/2)*stderr a <- cbind(b, stderr, low, high) a } est.b(tur.logit) # create a function to compute estimates of # the probability of success, standard errors, # and confidence limits at specific values # of the explanatory variable est.prob <- function (x0, obj, level=.95) { b <- coef(obj) x <- cbind(rep(1,length(x0)),x0) logitnew <- x%*%as.matrix(b) pihatnew <- exp(logitnew)/(1+exp(logitnew)) stderr <- sqrt(diag(x%*%vcov(obj)%*%t(x))) low <- logitnew - qnorm(1-(1-level)/2)*stderr low <- exp(low)/(1+exp(low)) high <- logitnew + qnorm(1-(1-level)/2)*stderr high <- exp(high)/(1+exp(high)) stdp <- stderr*pihatnew*(1-pihatnew) a <- cbind(x0, pihatnew, stdp, low, high) dimnames(a) <- list(NULL, c("x", "Prob", "Stderr", "low", "high")) a } newtemp <- c(27.0, 27.5, 28.0, 28.5, 29.0, 29.5, 30.0) est.prob(newtemp, tur.logit) # Confidence intervals for the value of the explanatory # variable at which the probability of success is p est.temp <- function(obj, p=0.5, level=.95){ eta <- family(obj)\$link(p) b <- coef(obj) x.p <- (eta-b[1])/b[2] pd <- -cbind(1, x.p)/b[2] stderr <- sqrt(diag((pd%*%vcov(obj)%*%t(pd)))) low <- x.p - qnorm(1-(1-level)/2)*stderr high <- x.p + qnorm(1-(1-level)/2)*stderr a <- cbind(p,x.p, stderr, low, high) dimnames(a) <- list(NULL, c("Prob","x","Stderr","low","high")) a } est.temp(tur.logit,,p=c(.10,.25,.50,.75,.90)) # Estimate a logistic curve for a second sample and # display it on the same plot temp2<-c(27.2, 27.7, 28.3, 28.4, 29.9) nmale2 <- c(1, 6, 14, 19, 26) nfemale2 <- c(29, 22, 14, 9, 2) ntur2 <- nmale2+nfemale2 pmale2 <- nmale2/ntur2 tur.logit2 <- glm(pmale2 ~ temp2, family=binomial, weights=ntur2) summary(tur.logit2) tempnew <- seq(27.0,30.0,0.1) logitnew2 <- tur.logit2\$coef[1]+tur.logit2\$coef[2]*tempnew pihatnew2 <- exp(logitnew2)/(1+exp(logitnew2)) # Add the curve to the previous plot lines(tempnew, pihatnew2, lty=3) # Create a legend. Note that the first two arguments legend # in the legend function are the (x,y) values for the position # of the upper left corner of the legend on the graph. legend(28.4,0.3, legend=c('Illinois turtles', 'Other turtles'), lty=c(1,3),bty='n')