#############################
#############################
## Regression III          ##
## Lecture 7: Nonlinearity ##
## 7-02-08                 ##
#############################
#############################

## Page 4
## Box-Cox Transformation Example: Ornstein Data (1)

data(Ornstein)
mod1 <- lm(interlocks ~ assets + sector + nation, 
    data=Ornstein)


mod <- lm(interlocks ~ assets + sector + nation + box.cox.var(interlocks+1), 
    data=Ornstein)

y.trans <- Ornstein$interlocks^.31

mod <- lm(y.trans ~ assets + sector + nation , 
    data=Ornstein)

summary(mod)

## Page 5 
## Box-Cox Transformation Example: Ornstein Data (2)

av.plots(mod, "box.cox.var(interlocks + 1)", col="black", identify.points=F)

## Page 8
## Box-Tidwell Transformation Example: Prestige Data 


box.tidwell(prestige ~ income + education, ~poly(women, 2), data=Prestige)

## Page 9-10
## AV Plots for Box-Tidwell Constructed Variables
mod <- lm(prestige ~ income + education + women + I(women^2), data=Prestige)
mod1 <- lm(prestige ~ income + education + women + I(women^2) + I(income * log(income)) + 
    I(education *log(education)), data=Prestige)

av.plots(mod1, "I(income * log(income))", identify.points=F, col="black")
av.plots(mod1, "I(education * log(education))", identify.points=F, col="black")

## Page 12
## Income

data(Prestige)
inc.mod1 <- lm(prestige ~ log(income) + education + women, data=Prestige)

## Page 13
## Transform of Income

box.tidwell(prestige ~ income, ~ education + women, data=Prestige)    

## Page 14
## Transformation of Income
inc.mod2 <- lm(prestige ~ I(income^0.08073) + 
    education + women, data=Prestige)

## Page 15
## Diagnostics for Education
cr.plots(inc.mod2, "education")

## Page 16
## Polynomial in Education

ed.mod1 <- lm(prestige ~ I(income^0.08073) + 
    poly(education, 3) + women, data=Prestige)

## Page 17
## Women and Prestige

cr.plots(ed.mod1, "women")

## Page 18
## Polynomial for Women
wom.mod1 <- lm(prestige ~ I(income^0.08073) + 
    education + poly(women, 2), data=Prestige)

## Page 19
## Re-checking Income
box.tidwell(prestige ~ income, ~education + poly(women, 2), data=Prestige)

## Page 20
## Final Model

final.mod <- lm(prestige ~ log(income) + 
    education + poly(women, 2), data=Prestige)
 
## Page 21
## Using the Effects Library

library(effects)

plot(effect("income", final.mod))
plot(effect("women", final.mod))

## Page 27
## An Example
library(car)
data(Prestige)

mod <- lm(prestige ~ income, data=Prestige)
cr.plots(mod)

trans.mod <- lm(prestige ~ log(income), data=Prestige)
loess.mod <- loess(prestige ~ income, data=Prestige, span=.5, family="symmetric")



rss0 <- sum((Prestige$prestige - trans.mod$fitted)^(2))
rss1 <- sum((loess.mod$y - loess.mod$fitted)^2)
F0 <- ((rss0-rss1)/(loess.mod$trace.hat-length(trans.mod$coef)))/
    (rss1/(loess.mod$n-loess.mod$trace.hat))

1-pf(F0, (loess.mod$trace.hat-length(trans.mod$coef)), 
    (loess.mod$n-loess.mod$trace.hat))

## Page 28-30
## Plots of Predictions

inc.seq <- seq(from=min(Prestige$income), to=max(Prestige$income), length=1000)
para.pred <- predict(trans.mod, newdata=data.frame(income=inc.seq), se.fit=T)
nonpara.pred <- predict(loess.mod, newdata=data.frame(income=inc.seq), se=T)

pdf("paranonpara.pdf", height=8, width=8)
plot(para.pred$fit ~ inc.seq, type="l", xlab="Income", ylab="Predicted Prestige", 
    ylim=c(min(c(para.pred$fit, nonpara.pred$fit)), 
        max(c(para.pred$fit,nonpara.pred$fit))))
lines(nonpara.pred$fit ~ inc.seq, lty=2)
points(Prestige$prestige ~ Prestige$income, pch=".", cex=3)
legend("bottomright", c("Parametric Fit", "Nonparametric Fit"), lty=c(1,2), 
    inset=.01)
dev.off()

poly.x <- c(inc.seq, inc.seq[1000:1], inc.seq[1])
para.lower <- para.pred$fit - 2*para.pred$se.fit
para.upper <- para.pred$fit + 2*para.pred$se.fit
nonpara.lower <- nonpara.pred$fit - 2*nonpara.pred$se.fit
nonpara.upper <- nonpara.pred$fit + 2*nonpara.pred$se.fit

para.poly.y <- c(para.lower, para.upper[1000:1], para.lower[1])
nonpara.poly.y <- c(nonpara.lower, nonpara.upper[1000:1], nonpara.lower[1])

pdf("paranonparaci.pdf", height=8, width=8)
plot(para.pred$fit ~ inc.seq, type="n", xlab="Income", ylab="Predicted Prestige", 
    ylim=c(min(c(para.lower, nonpara.lower)), max(c(para.upper, nonpara.upper))))
polygon(poly.x, para.poly.y, col="gray65", border=NA)
polygon(poly.x, nonpara.poly.y, col=rgb(0,0,1,.25), border=NA)
lines(para.pred$fit ~ inc.seq, lty=1)
lines(nonpara.pred$fit ~ inc.seq, lty=2)
legend("bottomright", c("Parametric Fit", "Nonparametric Fit"), lty=c(1,2), 
    inset=.01)
dev.off()

## Page 35
## Regression Splines

x <- seq(1,100, by=1)

before <- function(x) ifelse (x<60, 60-x,0)
after <- function(x) ifelse (x<60,0, x-60)

X <- cbind(before(x), after(x))

y <- 1 + 1*before(x) + 1*after(x) + rnorm(100, 0, 5)

plot(y~x)

## Page 37
## Splines vs. Dummy Interactions

d <- as.numeric(x>60)
mod <- lm(y ~ x*d)

x.s <- seq(from=0, to=100, by=1)
d.s <- as.numeric(x.s>60)

pred <- predict(mod, newdata=data.frame(x=x.s, d=d.s))

plot(pred ~ x.s, type="l")

## Page 40
## Fixing the Discontinuity

x1p <- ifelse(x > 60, x-60, 0)

mod1 <- lm(y ~ x + x1p)
x1p.s <- ifelse(x.s > 60, x.s-60, 0)
pred1 <- predict(mod1, newdata=data.frame(x=x.s, x1p=x1p.s))

plot(pred1 ~ x.s, type="l")

## Page 46
## The Basis Function in R

rk <- function(x,z){
  ((z-0.5)^2-1/12)*((x-0.5)^2-1/12)/4-
  ((abs(x-z)-0.5)^4-(abs(x-z)-0.5)^2/2 + 7/240)/24
  }

#Model Matrix For Spline Regression
spl.X <- function(x, xk){
   # set up model matrix for cubic penalized regression spline
   q <- length(xk)+2 # number of parameters
   n <- length(x) #amount of data
   X <- matrix(1,n,q) #initialized model matrix with constant in first column
   X[,2] <- x #set second column to x
   X[,3:q] <- outer(x,xk,FUN=rk) # add in remaining to R(x,xk)
   X
   }

## Page 47-48
## Estimating the Model in $
library(foreign)
jacob <- read.dta("jacob.dta")
attach(jacob)

# Rescale Ind. Variable
perot <- perotvote - min(perotvote)
perot <- perot/max(perot)

#Now Choose Number of Knots - Here 4 Knots
xk <- 1:4/5
#Generate Model Matrix
X <- spl.X(perot,xk)
#Fit Regression Model n- Be sure to Remove the Constant
mod.1 <- lm(chal.vote~X-1)
#X values for prediction
xp <- 0:100/100
#Prediction Matrix
Xp <- spl.X(xp,xk)
#Figure 3.3
plot(perot, chal.vote, type="n", ylab="Challengers' Vote Share (%)", 
    xlab="Vote for Perot", bty="l")
points(perot, chal.vote, pch=".", cex=1.75)
lines(xp,Xp%*%coef(mod.1))

## Page 49 
## Other Basis Functions
library(splines)

## Page 50
## Comparison of Different Basis Functions
mod.ns <- lm(chal.vote ~ ns(perot, knots=xk))
mod.bs <- lm(chal.vote ~ bs(perot, knots=xk))

pred.ns <- predict(mod.ns, newdata=data.frame(perot=xp))
pred.bs <- predict(mod.bs, newdata=data.frame(perot=xp))

plot(perot, chal.vote, pch=".", cex=1.75, ylab="Challengers' Vote Share (%)", 
    xlab="Vote for Perot", 
bty="l")
lines(xp,Xp%*%coef(mod.1), lwd=3)
lines(pred.ns ~ xp, lty=2, col="red", lwd=3)
lines(pred.bs ~ xp, lty=3, col="blue", lwd=3)

legend("bottomright", c("Cubic Spline", "Natural Spline", "B-Spline"), lty=c(1,2,3), 
    lwd=c(3,3,3), col=c("black", "red", "blue"),inset=.01)

## Page 51
## Comparing Number of Knots

mod.bs2 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=4)[-c(1,4)]))
mod.bs3 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=5)[-c(1,5)]))
mod.bs4 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=6)[-c(1,6)]))
mod.bs5 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=7)[-c(1,7)]))
mod.bs6 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=8)[-c(1,8)]))
mod.bs7 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=9)[-c(1,9)]))
mod.bs8 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=10)[-c(1,10)]))
mod.bs9 <- lm(chal.vote ~ bs(perot, knots=seq(0,1, length=11)[-c(1,11)]))


preds <- list()
preds[[1]] <- predict(mod.bs2, newdata=data.frame(perot=xp))
preds[[2]] <- predict(mod.bs3, newdata=data.frame(perot=xp))
preds[[3]] <- predict(mod.bs4, newdata=data.frame(perot=xp))
preds[[4]] <- predict(mod.bs5, newdata=data.frame(perot=xp))
preds[[5]] <- predict(mod.bs6, newdata=data.frame(perot=xp))
preds[[6]] <- predict(mod.bs7, newdata=data.frame(perot=xp))
preds[[7]] <- predict(mod.bs8, newdata=data.frame(perot=xp))
preds[[8]] <- predict(mod.bs9, newdata=data.frame(perot=xp))

for(i in 1:8){
    pdf(paste("predbs", i, ".pdf", sep=""), height=6, width=6)
    plot(perot, chal.vote, pch=".", cex=1.75, ylab="Challengers' Vote Share (%)", 
        xlab="Vote for Perot")
    lines(xp, preds[[i]])
    dev.off()
}



## page 52-53
## How Many Knots

mods.aic <- c(
AIC(mod.bs2), 
AIC(mod.bs3), 
AIC(mod.bs4), 
AIC(mod.bs5), 
AIC(mod.bs6), 
AIC(mod.bs7), 
AIC(mod.bs8), 
AIC(mod.bs9))

mods.bic <- c(
BIC(mod.bs2), 
BIC(mod.bs3), 
BIC(mod.bs4), 
BIC(mod.bs5), 
BIC(mod.bs6), 
BIC(mod.bs7), 
BIC(mod.bs8), 
BIC(mod.bs9))
par(mar=c(5,4,4,4)
)
plot(2:9, mods.aic, xlab="# Knots", ylab="AIC", type="l")
par(new=T)
plot(2:9, mods.bic, axes=F, xlab="", ylab="", main="", type="l", lty=2)
axis(4)
mtext("BIC", 4, 3)