library(faraway)
# Snoqualmie Falls precipitation example
  # Data consists of precipitation, in 1/100 of an inch, for each day for
  # each year from 1948 through 1983, at Snoqualmie Falls, Wash.
# load:
snoqualmie <- read.csv("snoqualmie.csv",header=FALSE)
# Let's get rid of the year breaks
snoqualmie <- unlist(snoqualmie)
# Most years have NAs for the 366th day, let's remove that
snoqualmie <- na.omit(snoqualmie)

# Look at the histogram,
plot(hist(snoqualmie,n=50,probability=TRUE),xlab="Precipitation (1/100 inch)")
rug(snoqualmie,col="grey")
# and the scatter-plot of successive values
plot(snoqualmie[-length(snoqualmie)],snoqualmie[-1],
     xlab="Precipitation today (1/100 inch)",
     ylab="Precipitation tomorrow (1/100 inch)",cex=0.1)
rug(snoqualmie[-length(snoqualmie)],side=1,col="grey")
rug(snoqualmie[-1],side=2,col="grey")
# both show big spikes at zero precipitation.  We'll try predicting whether
# there will be any precipitation tomorrow, from how much there was today.

# Re-shape so we have pairs of sequential days
vector.to.pairs <- function(v) {
  v <- as.numeric(v)
  n <- length(v)
  return(cbind(v[-1],v[-n]))
}

snoq.pairs <- vector.to.pairs(snoqualmie)
colnames(snoq.pairs) <- c("tomorrow","today")
snoq <- as.data.frame(snoq.pairs)

# Fit the logistic regression model
snoq.logistic <- glm((tomorrow > 0) ~ today, data=snoq, family=binomial)
# Look at coeffcients etc. --- positive coefficient on today's precipitation,
# as we'd expected from the scatterplot and background knowledge
summary(snoq.logistic)

# Plot data, plus predicted probabilities with confidence interval on predicted
# probability
plot((tomorrow>0)~today,data=snoq,xlab="Precipitation today (1/100 inch)",
     ylab="Positive precipitation tomorrow?")
rug(snoq$today,col="grey")

# Fake data to give visually-nice, evenly-spaced points for plotting
data.plot <- data.frame(today=(0:500))
# Work out confidence interval on the logit scale (where things are roughly
# Gaussian), then transform to probability scale
logistic.predictions <- predict(snoq.logistic,newdata=data.plot,se.fit=TRUE)
lines(0:500,ilogit(logistic.predictions$fit))
lines(0:500,ilogit(logistic.predictions$fit+1.96*logistic.predictions$se.fit),
      lty=2)
lines(0:500,ilogit(logistic.predictions$fit-1.96*logistic.predictions$se.fit),
      lty=2)

# Simple spline smoothing for comparison
snoq.spline <- smooth.spline(x=snoq$today,y=(snoq$tomorrow>0))
# Add the spline to the plot
lines(snoq.spline,col="red")
  # Pretty good match to the logistic regression up to about 1.3 inches

# Try a GAM, which moves the smoother _inside_ the likelihood model
library(mgcv)
  # Note: you may need to un-load the gam library with detach()
snoq.gam <- gam((tomorrow>0)~s(today),data=snoq,family=binomial)
# Add the GAM to the plot (again using the logit-scale trick)
gam.predictions <- predict.gam(snoq.gam,newdata=data.plot,se.fit=TRUE)
lines(0:500,ilogit(gam.predictions$fit),col="blue")
lines(0:500,ilogit(gam.predictions$fit+1.96*gam.predictions$se.fit),col="blue",
      lty=2)
lines(0:500,ilogit(gam.predictions$fit-1.96*gam.predictions$se.fit),col="blue",
      lty=2)
  # Note: similar to the simple spline smoothing, but not quite identical!
    # Also note: huge error bands at the right, where there's very little data

# Width of GLM confidence bands:
plot(0:500,ilogit(logistic.predictions$fit)
           -ilogit(logistic.predictions$fit-1.96*logistic.predictions$se.fit),
    type="l",col="blue",xlab="Precipitation today (1/100 inch)",
    main="Difference in probability between prediction\n
and confidence limit for prediction",
    ylab = expression(paste(Delta,"probability")))
lines(0:500,ilogit(logistic.predictions$fit+1.96*logistic.predictions$se.fit)
            -ilogit(logistic.predictions$fit))
  # Notice that these are asymmetric, and actually narrower at really large
  # x than at moderate x; why?


# Significant mis-specification for the GLM?
# Observed difference in deviance
diff.dev.obs <- snoq.logistic$deviance - snoq.gam$deviance

# Simulate from the fitted logistic regression model for Snoqualmie
# Presumes: fitted values of the model are probabilities.
snoq.sim <- function(model=snoq.logistic) {
  fitted.probs <- fitted(model)
  n <- length(fitted.probs)
  new.binary <- rbinom(n,size=1,prob=fitted.probs)
  return(new.binary)
}

# Simulate from fitted logistic regression, re-fit logistic regression and
# GAM, calculate difference in deviances
diff.dev <- function(model=snoq.logistic,x=snoq[,2]) {
  y.new <- snoq.sim(model)
  GLM.dev <- glm(y.new ~ x,family=binomial)$deviance
  GAM.dev <- gam(y.new ~ s(x),family=binomial)$deviance
  return(GLM.dev-GAM.dev)
}

# The following takes ~1500 seconds on my computer; reduce number of replicates
# to taste
null.dist.of.diff.dev <- replicate(1000,diff.dev())
p.value <- (1+sum(null.dist.of.diff.dev > diff.dev.obs))/(1+length(null.dist.of.diff.dev))
  # p-value is < 10^-3; the observed difference is much bigger than what's
  # generated by simulating from the logistic model

# What went wrong with the logistic model?  Inspecting the plots, the spline
# and the GAM both predict lower probabilities at X=0 than the GLM does, and
# their predicted probabilities increase very steeply there.
# GLM's predicted probability of precipitation following a dry day:
ilogit(snoq.logistic$coefficients[1])
  # which is 0.502
# Actual frequency
mean(snoq$tomorrow[snoq$today==0]>0)
  # which is 0.470
# Since there are a _lot_ of dry days (also about 47%), a couple of percentage
# points difference in probability here has a big contribution in the
# likelihood.
# Therefore: Try telling the GLM that no precipitation today is a special value

# Expand the data with an extra column, indicating whether today was dry or
# not
snoq2 <- data.frame(snoq,dry=ifelse(snoq$today==0,1,0))
# Fit a new logistic model with that dummy variable
snoq2.logistic <- glm((tomorrow > 0) ~ today + dry,data=snoq2,family=binomial)
# GAM ditto
snoq2.gam <- gam((tomorrow > 0) ~ s(today) + dry,data=snoq2,family=binomial)

# Plot the new models
plot((tomorrow>0)~today,data=snoq,xlab="Precipitation today (1/100 inch)",
     ylab="Positive precipitation tomorrow?")
rug(snoq$today,side=1,col="grey")

# Need to add a "dry" column to the fake data for plotting
data.plot=data.frame(data.plot,dry=ifelse(data.plot$today==0,1,0))
logistic.predictions2 <- predict(snoq2.logistic,newdata=data.plot,se.fit=TRUE)
lines(0:500,ilogit(logistic.predictions2$fit))
lines(0:500,ilogit(logistic.predictions2$fit+1.96*logistic.predictions2$se.fit),
      lty=2)
lines(0:500,ilogit(logistic.predictions2$fit-1.96*logistic.predictions2$se.fit),
      lty=2)
gam.predictions2 <- predict.gam(snoq2.gam,newdata=data.plot,se.fit=TRUE)
lines(0:500,ilogit(gam.predictions2$fit),col="blue")
lines(0:500,ilogit(gam.predictions2$fit+1.96*gam.predictions2$se.fit),
      col="blue",lty=2)
lines(0:500,ilogit(gam.predictions2$fit-1.96*gam.predictions2$se.fit),
      col="blue",lty=2)
  # Notice that the GLM and the GAM now coincide almost exactly in their
  # predictions.

# Look at snoq2.gam to make sure I didn't accidentally forget to turn on
# smoothing!
summary(snoq2.gam)
plot(snoq2.gam)
  # The "today" term is in fact a non-parametric smooth with a little bit of
  # curvature

# Difference still significant?
diff.dev.obs2 <- snoq2.logistic$deviance - snoq2.gam$dev
diff.dev.2 <- function(model=snoq2.logistic,x=snoq2[,-1]) {
  y.new <- snoq.sim(model)
  sim.data <- data.frame(x,tomorrow=y.new)
  GLM.dev <- glm(tomorrow ~ today + dry,data=sim.data,family=binomial)$deviance
  GAM.dev <- gam(tomorrow ~ s(today)+dry,data=sim.data,family=binomial)$deviance
  return(GLM.dev-GAM.dev)
}
null.dist.of.diff.dev2 <- replicate(100,diff.dev.2())
p.value <- (1+sum(null.dist.of.diff.dev2 > diff.dev.obs2))/(1+length(null.dist.of.diff.dev2))
# As close to 1 as computation time allows

# So let's look at calibration
# First, what's the predicted probability of precipitation after dry days?
predict(snoq2.logistic,newdata=data.frame(today=0,dry=1),type="response")
  # Matches the actual frequency to seven decimal places --- it should, we've
  # introduced a coefficient to do just this!

# Look at the histogram of predicted probabilities
hist(fitted(snoq2.logistic))
  # Big gap between 0.47 and 0.55
# Are we calibrated between say 0.55 and 0.56?
mean(snoq$tomorrow[(fitted(snoq2.logistic) >= 0.55)
                   & (fitted(snoq2.logistic) < 0.56)] > 0)
  # Gives 0.547, pretty reasonable
# Write a function to do this less manually
frequency.vs.probability <- function(p.lower,p.upper=p.lower+0.01,
                                     model=snoq2.logistic,
                                     data=(snoq$tomorrow>0)) {
  # Presume that fitted values are probabilities
  fitted.probs <- fitted(model)
  # Which rows of the data have fitted probabilities in our range of interest?
  indices <- (fitted.probs >= p.lower) & (fitted.probs < p.upper)
  # For plotting purposes, what's our average predicted probability over this
  # range?
  ave.prob <- mean(fitted.probs[indices])
  # How big a standard error should we see for binomial data with this predicted
  # probability?
    # Rough calculation, could improve by taking account of variation in
    # predicted probabilities, not reliable if difference between p.lower
    # and p.upper was substantial
  se <- sqrt(ave.prob*(1-ave.prob)/sum(indices))
  # How often does the event actually happen?
    # Presumes data is either 0/1 or TRUE/FALSE valued
  frequency <- mean(data[indices])
  out <- list(frequency=frequency,ave.prob=ave.prob,se=se)
  return(out)
}

# Look at the whole range of probabilities between 0.55 and 0.74 (because
# that's most of our predictions, and there's a gap about 0.75 which is
# annoying to deal with)
f.vs.p <- sapply((55:74)/100,frequency.vs.probability)
# Re-shape the data to remove some R weirdness
f.vs.p <- data.frame(frequency=unlist(f.vs.p["frequency",]),
                     ave.prob=unlist(f.vs.p["ave.prob",]),
                     se=unlist(f.vs.p["se",]))

# Calibration plot
plot(f.vs.p$ave.prob,f.vs.p$frequency,xlim=c(0,1),ylim=c(0,1),
     xlab="Predicted probabilities",ylab="Observed frequencies")
rug(fitted(snoq2.logistic),col="grey")
abline(0,1,col="grey")
segments(x0=f.vs.p$ave.prob,y0=f.vs.p$ave.prob-1.96*f.vs.p$se,
         y1=f.vs.p$ave.prob+1.96*f.vs.p$se)
# Vertical lines are (approximate) 95% sampling intervals for the frequency,
# given the predicted probability and the number of observations