#  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')