#
# 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.4.pdf",width=10.5,height=7.5,horiz=T)
par(mfrow=c(1,2))
#
lfit <- loess(prestige ~ income + education, span=0.5, degree=1)
summary(lfit)
#
# 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(lfit, 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="Local polynomial 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")
#
# Now do significance statistics
#
lfit.income <- loess(prestige ~ income, span=0.7, degree=1)
lfit.education <- loess(prestige ~ education, span=0.7, degree=1)
#
# Look at Analysis-of-Variance decomposition
# for both fits
#
# First do the incremental effect of Education
#
anova(lfit.income, lfit.education)
#
# Now do the incremental effect of Income
#
anova(lfit.education, lfit.income)