### Class examples from Week 9: Regresion
### Ken Benoit Mar 2010


## Table/Figure 11.1 example
par(mar=c(4,4,1,1))
x <- c(0,3,1,0,6,5,3,4,10,8)
y <- c(12,13,15,19,26,27,29,31,40,48)
plot(x, y, xlab="Number of prior convictions (X)",
     ylab="Sentence length (Y)", pch=19)
abline(h=c(10,20,30,40), col="grey70")

## calculations from Table 11.2
(tab11.2 <- data.frame(x, y,
                       xdev=(x-mean(x)),
                       ydev=(y-mean(y)),
                       xdevydev=((x-mean(x))*(y-mean(y))),
                       xdev2=(x-mean(x))^2,
                       ydev2=(y-mean(y))^2))
(SP <- sum(tab11.2$xdevydev))
(SSx <- sum(tab11.2$xdev2))
(SSy <- sum(tab11.2$ydev2))
(b <- SP / SSx)
(a <- mean(y) - b*mean(x))

## Figure 11.2
par(mar=c(4,4,1,1))
plot(x, y, xlab="Number of prior convictions (X)",
     ylab="Sentence length (Y)", pch=19)
abline(lm(y~x))
text(8,45,"Regression line")
text(8,43.5,"Yhat = 14 + 3X")
abline(v=mean(x), h=mean(y), lty="dashed")
points(0,a,pch=0)
points(mean(x),mean(y),pch=0)
text(.8,13,"Y intercept")

## Example pp267-268
x <- c(49,50,53,55,60,55,60,50)
y <- c(81,88,87,99,91,89,95,90)
N <- length(x)
(SSx <- sum(x^2) - N*mean(x)^2)
(SSy <- sum(y^2) - N*mean(y)^2)
(SP <- sum(x*y) - N*mean(x)*mean(y))
(b <- SP/SSx)
(a <- mean(y) - b*mean(x))
(r <- SP / sqrt(SSx*SSy))
(r2 <- r^2)

# plot the regression model
plot(x, y, pch=19)
abline(lm(y~x))

# regression models in R
regmodel <- lm(y ~ x)
summary(regmodel)
coef(summary(regmodel))


## Illustration using Anscombe dataset
data(anscombe)
attach(anscombe)
round(coef(lm(y1~x1)), 2)           # same b0, b1
round(coef(lm(y2~x2)), 2)
round(coef(lm(y3~x3)), 2)
round(coef(lm(y4~x4)), 2)
round(summary(lm(y1~x1))$r.squared, 2)   # same R^2
round(summary(lm(y2~x2))$r.squared, 2)
round(summary(lm(y3~x3))$r.squared, 2)
round(summary(lm(y4~x4))$r.squared, 2)
# plot the four x-y pairs
par(mfrow=c(2,2), mar=c(4, 4, 1, 1)+0.1) # 4 plots in one graphic window
plot(x1,y1)
abline(lm(y1~x1), col="red", lty="dashed")
plot(x2,y2)
abline(lm(y2~x2), col="red", lty="dashed")
plot(x3,y3)
abline(lm(y3~x3), col="red", lty="dashed")
plot(x4,y4)
abline(lm(y4~x4), col="red", lty="dashed")
detach(anscombe)


## Multiple regression Dail spending example
require(foreign)
dail <- read.dta("dailcorrected.dta")
summary(lm(votes1st ~ spend_total + incumb + electorate + minister, data=dail))