#### Part I: Demos of density estimation ####

# Load some data about cats from a standard library
library(MASS)
data(cats)
dim(cats)
summary(cats)  # Bwt = body weight (kilograms), Hwt = heart weight (grams)
# A scatterplot tells us about the joint distribution
plot(Hwt~Bwt,data=cats)
# Rug plots along the axes give us a sense of the marginal distributions
rug(cats$Bwt,side=1,col="grey")
rug(cats$Hwt,side=2,col="grey")

# Focus on heart weight
# Plot the histogram, with probability density estimate on vertical axis
plot(hist(cats$Hwt,n=41),freq=FALSE)
rug(cats$Hwt,side=1,col="grey",lwd=2)
# The values are all positive, with a mode and then a longish right tail,
# so try a gamma distribution

# Fit a gamma distribution by maximum likelihood
   # Function: negative mean log likelihood of a gamma distribution
     # Negative so that it works will with minimizers like optim()
     # Inputs: vector of parameters (shape and scale), vector of data values
     # Output: negative average log likelihood of that data with those parameters
gamma.negloglik <- function(params,data) {
    -mean(dgamma(data,shape=params[1],scale=params[2],log=TRUE))
}
# Find the parameter values by minimizing this function
mle <- optim(par=c(1,10),fn=gamma.negloglik,data=cats$Hwt)$par
# Add the pdf of the gamma with those values to the plot
curve(dgamma(x,shape=mle[1],scale=mle[2]),add=TRUE,lwd=2)
   # Exercise for the reader: make a new 3D plot which shows gamma.negloglik
   # as we sweep the shape and scale parameters, holding the data fixed to
   # cat$Hwt

# How about a non-parametric, kernel, density estimate?
# Rough bandwidth rule: in 1D, use 1.05*standard deviation/(n^0.2)
   # Derivation: this turns out to be the optimal bandwidth, if we presume
   # the true distribution is Gaussian.  The n^(1/5) rate is generic, though.
n <- length(cats$Hwt)  # So we don't have to keep writing that out
rough.bw <- 1.05*sd(cats$Hwt)/(n^(1/5))
# To illustrate the idea of adding a copy of the kernel from each data point
  # Pick a _few_ data points at random, and show the copy of the kernel
  # for each of them
representative.cats <- sample(1:n,size=5,replace=FALSE)
rug(cats$Hwt[representative.cats],col="blue",lwd=3)
for (i in representative.cats) {
    curve((1/n)*dnorm(x,mean=cats$Hwt[i],sd=rough.bw),add=TRUE,col="blue")
}
# Note: by dividing by 1/n, each one makes a small contribution.  Change to
# 1/5 above to see the shape of each one in a bit more detail.

# Now add the over-all density to the plot
lines(density(cats$Hwt,bw=rough.bw),col="blue",lwd=2)
  # Note: If we don't specify bandwith through bw, density() will pick it
  # automatically; see help(density) for details

# What about 2 D?
plot(Hwt~Bwt,data=cats)
rug(cats$Bwt,side=1,col="grey")
rug(cats$Hwt,side=2,col="grey")
# Nonparametric kernel density estimation in multiple dimensions is built
# in to the np package
library(np)
# Estimate the unconditional (joint) density, and time how long it takes
# to do so
   # Can slow down drastically with more data and/or more dimensions; look
   # at tricks for speeding up npreg
system.time(cats.dens <- npudens(~ Hwt+Bwt, data=cats))
# Plot in 3D, with rotation by default...
plot(cats.dens)
# Plot in 3D with a fixed perspective, and return the information used to
# make the plot, if we want to pick that apart later ourselves
plot(cats.dens,view="fixed",plot.behavior="plot-data")

# Conditional density
cats.hwt.on.bwt <- npcdens(Hwt ~ Bwt, data=cats)
plot(cats.hwt.on.bwt,view="fixed",phi=75)

# Conditioning on a categorical as well:
cats.hwt.on.bwt.and.sex <- npcdens(Hwt ~ Bwt + factor(Sex), data=cats)
# Default plot: sweep one variable at a time through its range, holding the
# other two fixed
plot(cats.hwt.on.bwt.and.sex)

# Following the notes, we can make side-by-side 3D plots for the each level
# of the categorical variable
# First, load the appropriate graphics library
library(lattice)
# Next, create a regular grid on which to evaluate the conditional pdf;
# we need all combinations of all 3 variables, so use expand.grid()
  # See recipe 12.10 in _The R Cookbook_ for a tutorial on this function
cats.grid <- expand.grid(Bwt=seq(from=min(cats$Bwt),to=max(cats$Bwt),length.out=50),
                         Hwt=seq(from=min(cats$Hwt),to=max(cats$Hwt),length.out=50),
                         Sex=c("F","M"))
# Now, call predict for the conditional pdf on this grid
fhat <- predict(cats.hwt.on.bwt.and.sex,newdata=cats.grid)
# Finally, make 3D wireframe plots where height is the value of the pdf,
# the two horizontal-plane variables are Hwt and Bwt, and we make another
# plot for each level of Sex
wireframe(fhat~cats.grid$Bwt*cats.grid$Hwt|cats.grid$Sex,
          scales=list(arrows=FALSE),xlab="Bwt",ylab="Hwt",zlab="pdf")


####### Part II: Some hints related to the exam ######

# Fit an additive model of heart weight on body weight and sex
library(mgcv)
cat.gam.1 <- gam(Hwt ~ s(Bwt)+factor(Sex),data=cats)
# Sex has two categories, F and M; we a single coefficient which tells us
# what to add in the marked category (here, M), over and above what we'd
# expect for the default category (here, F).
# Running summary gives us estimates for all the parametric terms in the GAM,
# starting with the intercept
summary(cat.gam.1)
  # Note: just as in the linear models world, we can't estimate an intercept
  # _and_ an effect for one level of the category _and_ an effect for the
  # other level of the category; either we estimate a global intercept and
  # a _constrast_ between categories (the usual way, which is what R does),
  # or we estimate a different intercept for each category
# If we just want the vector of parametric coefficient estimates:
summary(cat.gam.1)$p.coef
# By contrast, we've got a whole smooth function for the partial response to
# body weight (only a bit curved in this case):
plot(cat.gam.1,pages=1,residuals=TRUE)
# Coefficients are a bit treacherous for GAMs:
coefficients(cat.gam.1)
  # The coefficients vector begins with all the parametric coefficients, which
  # we can interpret in the old familiar way
  # As for the rest, remember that s() stands for "spline smoothing", and
  # splines are piecewise cubic polynomials.  We're getting the coefficients
  # of the smooth partial response function in a particular representation or
  # basis, along the lines described in Section 8.3 of the notes.  These are
  # NOT interpretable in anything like the same way.

# We can interact a smooth term with a categorical one by using "by=" inside
# the smoothing command.
  # See help(gam.models) for more on that, and many other options/issues for
  # specifying GAMs.

# This tries out including a smooth main effect for Bwt and an interaction
# with Sex, but it's over-ambitious
cat.gam.2 <- gam(Hwt ~ s(Bwt)
                 + factor(Sex)
                 + s(Bwt,by=factor(Sex)),
                     data=cats)
plot(cat.gam.2,pages=1,residuals=TRUE,scale=0,se=FALSE)
# Plotting, we see that one of the estimated smooth functions is pretty exactly
# flat.  (Why is it the one for the default category?)  We have the same problem
# that we can't include a main effect _and_ an interactive effect for _every_
# level.

# Remedy one: no main effect, just one smooth function per level.
cat.gam.3 <- gam(Hwt ~ s(Bwt,by=Sex)+factor(Sex),data=cats)
plot(cat.gam.3,pages=1)

# Remedy two: include a main effect, and a _contrast_ smooth function
  # Every cat gets the contribution from the main-effect smooth
  # Only male cats (in this case) get the extra contribution from the
  # contrast function
# This is actually easier to do if the category is numerically coded as
# 0 and 1, so add such a column
cats$SexNumber <- ifelse(cats$Sex=="M",1,0)
summary(cats)
cat.gam.4 <- gam(Hwt ~ s(Bwt) + SexNumber,
                   +s(Bwt,by=SexNumber),data=cats)
plot(cat.gam.4)
# The 2nd smooth function is now the _constrast_ in expected heart weight
# between male and female cats at a give body weight
  # Whereas in cat.gam.3, the 2nd plot was the _total_ effect of Bwt on
  # Hwt, but only for male cats

### Extracting the partial response functions
# Suppose we want to know what each term in the model contributes to
# the prediction; predict() will give this to us:
predict(cat.gam.4,type="terms")
  # We get an array where rows are data points, and columns are the additive
  # contribution of each term in the model
# By default, predict() gives predictions for the original data, but we
# can ask for predictions on newdata, and also ask what each term in the
# model contributes
predict(cat.gam.4,
        newdata=expand.grid(SexNumber=c(0,1),
                            Bwt=c(2,3,4)),
        type="terms")

# If we want to sweep just one variable at a time, expand.grid() is not
# needed
predict(cat.gam.1,
        newdata=data.frame(Sex="F",
                           Bwt=seq(from=2,to=4,length.out=10)),
        type="terms")

# If we are only interested in one of the partial response functions at
# a time, we can capture that by selecting just one column of the terms
# matrix
partial.response.on.grid <- predict(cat.gam.1,
  newdata=data.frame(Sex="F",
  Bwt=seq(from=2,to=4,length.out=10)),
  type="terms")[,2]
# alternately, we could select the column by the name of the term, here
# "s(Bwt)"

# We could now plot this, if we liked:
plot(seq(from=2,to=4,length.out=10),partial.response.on.grid, type="l")

# On the other hand, if we re-estimated the GAM on some new data, we would
# just have to change that part of line 198 above to get a new version
# of the partial response function evaluated on the same grid of points.
# So we can use what we learned about bootstrapping kernel regressions
# and splines to get confidence bands for the partial response functions...