#
# Load data on the Prestige of 102 Canadian Occupational
# categories studied in the great text: 
#
#  Fox, J. (1997) Applied Regression, Linear Models,
#  and Related Methods. Sage.
#
# We will focus on:
#
#     income: Average income of incumbents, dollars, in 1971.
#
#     prestige: Pineo-Porter prestige score for occupation, from a social
#          survey conducted in the mid-1960s.
#
#     education: Average education of incumbents, years, in 1971.
#
# It is available from the "car" package on CRAN.
#
library(car)
data(Prestige)
attach(Prestige)
#
# loess is in the "modreg" base package
#
library(modreg)
#
pdf(file = "FigLec8.5.pdf",width=10.5,height=7.5,horiz=T)
par(mfrow=c(1,2))
#
plot(x=income, y=prestige, xlab="Average Income", ylab="Prestige",
  main="Local polynomial regression fits to Prestige")
#
# Add a local polynomial regression line
# with 70% of the data as the kernel neighborhood range
#
# First create the fit
#
lfit <- loess(prestige ~ income, span=0.7, degree=1)
summary(lfit)
#
# Predict over a range of 200 points
#
xfit <- data.frame(income=seq(from=0,to=25000,length=200))
pfit <- predict(lfit, newdata=xfit)
#
lines(x=xfit$income, y=pfit, lwd=2)
#
# Add a smoothing spline fit 
# with 70% of the data as the kernel neighborhood range
#
sfit <- smooth.spline(x=income, y=prestige, df=3.85)
sfit
pfit <- predict(sfit, x=xfit$income)
#
lines(x=xfit$income, y=pfit$y, lwd=2, lty=2)
legend(x=c(8000,8000),y=c(20,20),lty=1:2,lwd=2,
  legend=c("loess linear 70% nghd", "smoothing spline d.f.=3.85"))
#
library(mgcv)
#
gam.fit <- gam(prestige ~ s(income) + s(education))
gam.fit
#
# Next predict over a grid of points
#
income.vals <- seq(from=0,to=25000,length=25)
education.vals <- seq(from=6,to=16,length=25)
#
# Form a grid of values
#
qgrid <- expand.grid(income=income.vals, education=education.vals)
#
# Predict on this grid 
#
pfit <- matrix(predict(gam.fit, newdata=qgrid), ncol=25, nrow=25)
#
# Create a perspective plot of the fit
#
persp(x=income.vals, y=education.vals, z=pfit,
      xlab="Average Income", ylab="Education", zlab="Prestige",
      expand=2/3, shade=0.5, theta=45,
      main="Additive nonparametric regression fits to Prestige \n based on Income and Education",
      sub="perspective plot")
##
## Add a local polynomial regression line
## with 50% of the data as the kernel neighborhood range
#contour(x=income.vals, y=education.vals, z=pfit,
#      xlab="Average Income", ylab="Education",
#      main="Local polynomial regression fits to Prestige \n based on Income and Education",
#      sub="contour plot")