waldplot <- function(fit, var, steps=10, type=3, col=c("#FF4E37FF", "#008DFFFF"), ...) {
  require(MASS)
  require(splines)
  prof <- profile(fit, which(names(coef(fit))==var), maxsteps=steps)
  x <- prof[[var]]$par.vals[,var]
  y <- -prof[[var]]$z^2
  sp <- interpSpline(x,y)
  xx <- seq(min(x),max(x),len=101)
  par(mar=c(5,5,3,1))
  plot(xx, predict(sp,xx)$y, type="l", lwd=3, xlab=~beta, ylab=~"2{l"*(beta)-l(hat(beta))*"}", col=col[1], las=1, ...)
  if (type==1) return(invisible(NULL))
  lines(xx,-(xx-coef(fit)[var])^2/(vcov(fit)[var,var]),lwd=3, col=col[2])
  toplegend(c("Likelihood ratio", "Wald"), col=col, lwd=3, ncol=2)
  if (type==2) return(invisible(NULL))
  abline(h=-qchisq(.95,1), col="gray")
}

## Donner data
donner <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/donner.txt")
fit <- glm(Status~Age*Sex,data=donner,family=binomial)
waldplot(fit,"Age")

## Confidence intervals for parameters
confint(fit)          ## Likelihood ratio
confint.default(fit)  ## Wald

## LR confidence intervals for odds ratio: 10-year difference in age
Male <- 1*(donner$Sex=="Male")
Female <- 1-Male
donner$Status <- relevel(donner$Status, "Survived")
fit1 <- glm(Status~Age*Male, data=donner, family=binomial)
fit2 <- glm(Status~Age*Female, data=donner, family=binomial)
exp(10*confint(fit1, "Age"))
exp(10*confint(fit2, "Age"))

## Testing interaction
fit0 <- glm(Status~Sex+Age,data=donner,family=binomial)
anova(fit0, fit, test="Chisq")

## High agreement
x <- sort(rnorm(100))
y <- rep(0:1,50)
fit <- glm(y~x,family=binomial)
confint(fit)
confint.default(fit)
waldplot(fit,"x")

## Low agreement
x <- sort(rnorm(6))
y <- c(0,0,1,0,1,1)
fit <- glm(y~x,family=binomial)
waldplot(fit,"x")
confint(fit)
confint.default(fit)

## Complete separation
x <- sort(rnorm(10))
y <- c(0,0,0,0,0,1,1,1,1,1)
fit <- glm(y~x,family=binomial)

f <- function(bb) {
  ll <- numeric(length(bb))
  for (i in 1:length(bb)) {
    fit.i <- glm(y~offset(x*bb[i]),family=binomial)
    ll[i] <- 2*(logLik(fit.i))
  }
  ll
}
bb <- seq(0,100,len=51)
plot(bb, f(bb), type="l", lwd=3, xlab=~beta, ylab=expression(l(beta)))

## Note that one can still perform a hypothesis test:
fit0 <- glm(y~1, family=binomial)
anova(fit0, fit, test="Chisq")