###################################################
### chunk number 1: fig2-1
###################################################
#line 67 "Lecture2.Rnw"
x<-sort(rnorm(50))
y<-10+3*x+rnorm(50,sd=2)
plot(x,y)
abline(v=mean(x))
lines(x,lm(y~x)$fit,lwd=2)
points(x[1:10],y[1:10],col=2,pch=8) ## marking the extreme observations
points(x[40:50],y[40:50],col=2,pch=8)


###################################################
### chunk number 2: fig2-2
###################################################
#line 76 "Lecture2.Rnw"
vv<-sum((x-mean(x))^2)
ki<-(x-mean(x))/vv ## computing the weights k
plot(x,ki)
points(x[1:10],ki[1:10],col=2,pch=8) ## marking the extreme observations
points(x[40:50],ki[40:50],col=2,pch=8)
abline(h=0,lty=2)


###################################################
### chunk number 3: l2-v1
###################################################
#line 115 "Lecture2.Rnw"
## here I create several data sets y and fit regression models to them
## I plot and compare with the true model colored red (col=2)
## Note - in all the figures below I use the same axes -4,4 to
## make direct comparisons between the figures possible
x<-sort(rnorm(50))
for (kk in (1:50)) {
 y1<-2+3*x+rnorm(50)
 if (kk==1) {plot(x,y1,xlim=c(-4,4),ylim=c(-12,12))}
 mm<-lm(y1~x)
 lines(c(-4,4),mm$coef[1]+mm$coef[2]*c(-4,4)) }
lines(seq(-4,4),2+3*seq(-4,4),lwd=2,col=2)


###################################################
### chunk number 4: l2-v1b
###################################################
#line 128 "Lecture2.Rnw"
x<-sort(rnorm(50))
for (kk in (1:50)) {
 y1<-2+3*x+rnorm(50,sd=4) ## increasing the noise level
 if (kk==1) {plot(x,y1,xlim=c(-4,4),ylim=c(-12,12))}
 mm<-lm(y1~x)
 lines(c(-4,4),mm$coef[1]+mm$coef[2]*c(-4,4)) }
lines(seq(-4,4),2+3*seq(-4,4),lwd=2,col=2)


###################################################
### chunk number 5: l2-v2
###################################################
#line 137 "Lecture2.Rnw"
x<-sort(rnorm(25)) ## sample size n=25
for (kk in (1:50)) {
 y1<-2+3*x+rnorm(25)
 if (kk==1) {plot(x,y1,xlim=c(-4,4),ylim=c(-12,12))}
 mm<-lm(y1~x)
 lines(c(-4,4),mm$coef[1]+mm$coef[2]*c(-4,4)) }
lines(seq(-4,4),2+3*seq(-4,4),lwd=2,col=2)


###################################################
### chunk number 6: l2-v2b
###################################################
#line 146 "Lecture2.Rnw"
x<-sort(rnorm(250)) ## sample size n=250
for (kk in (1:50)) {
 y1<-2+3*x+rnorm(250)
 if (kk==1) {plot(x,y1,xlim=c(-4,4),ylim=c(-12,12))}
 mm<-lm(y1~x)
 lines(c(-4,4),mm$coef[1]+mm$coef[2]*c(-4,4)) }
lines(seq(-4,4),2+3*seq(-4,4),lwd=2,col=2)


###################################################
### chunk number 7: l2-v3
###################################################
#line 155 "Lecture2.Rnw"
x<-sort(rnorm(50))
for (kk in (1:50)) {
 y1<-2+3*x+rnorm(50)
 if (kk==1) {plot(x,y1,xlim=c(-4,4),ylim=c(-12,12))}
 mm<-lm(y1~x)
 lines(c(-4,4),mm$coef[1]+mm$coef[2]*c(-4,4)) }
lines(seq(-4,4),2+3*seq(-4,4),lwd=2,col=2)


###################################################
### chunk number 8: l2-v3b
###################################################
#line 164 "Lecture2.Rnw"
x<-sort(rnorm(50,sd=5)) ## increasing the spread in x
for (kk in (1:50)) {
 y1<-2+3*x+rnorm(50)
 if (kk==1) {plot(x,y1,xlim=c(-4,4),ylim=c(-12,12))}
 mm<-lm(y1~x)
 lines(c(-4,4),mm$coef[1]+mm$coef[2]*c(-4,4)) }
lines(seq(-4,4),2+3*seq(-4,4),lwd=2,col=2)


###################################################
### chunk number 9: LSfit
###################################################
#line 196 "Lecture2.Rnw"
x<-sort(rnorm(25))
y<-2+3*x+rnorm(25)
plot(x,y)
mm<-lm(y~x)
lines(x,mm$fit,lwd=2)
for (zz in (1:length(x))) {
  lines(rep(x[zz],2),c(y[zz],mm$fit[zz]),lty=2) 
  points(x[zz],mm$fit[zz],pch=8,col=2)}


###################################################
### chunk number 10: resfita
###################################################
#line 217 "Lecture2.Rnw"
plot(x,mm$res)
id<-identify(x,mm$res,pos=T)


###################################################
### chunk number 11: resfit
###################################################
#line 221 "Lecture2.Rnw"
plot(x,mm$res)
if (max(id$ind)!=-Inf) {text(x[id$ind],mm$res[id$ind],id$ind,pos=id$pos) }
abline(h=0,lty=2)


###################################################
### chunk number 12: l2-lev1
###################################################
#line 262 "Lecture2.Rnw"
x<-sort(rnorm(50))
y<-10+3*x+rnorm(50,sd=2)
plot(x,y)
abline(v=mean(x))
mm<-lm(y~x)
lines(x,mm$fit,lwd=2)
points(x[1:10],y[1:10],col=2,pch=8)
points(x[40:50],y[40:50],col=2,pch=8)


###################################################
### chunk number 13: l2-lev2
###################################################
#line 272 "Lecture2.Rnw"
lmi<-lm.influence(mm) 
## the lm.influence function extracts both hat-values and other useful diagnostics
plot(x,lmi$hat,ylab="leverage")
points(x[1:10],lmi$hat[1:10],col=2,pch=8)
points(x[40:50],lmi$hat[40:50],col=2,pch=8)


###################################################
### chunk number 14: l2-diag1
###################################################
#line 313 "Lecture2.Rnw"
## here I use lmi$coef[,2] = the change in the slope parameter when observations are dropped
## lmi$coef[,1] is the change in the intercept
plot(lmi$coef[,2],ylab="Change in Slope")
abline(h=0)


###################################################
### chunk number 15: l2-diag2
###################################################
#line 319 "Lecture2.Rnw"
## lmi$sig is the change in the least squares criterion (divided by n-2) when observations are dropped
plot(lmi$sig,ylab="Change in Sum of Squares")


###################################################
### chunk number 16: ch1a
###################################################
#line 338 "Lecture2.Rnw"
TVdat<-read.table('TV.dat',sep='\t') ## reading the data into R
print(dim(TVdat)) ## printing the dimensions of the data
print(names(TVdat))  ## printing the names of the variables
print(row.names(TVdat)) ## printing the names of the countries included
plot(TVdat$ppD,TVdat$ppT,xlab="people per Dr",ylab="people per TV")
id<-identify(TVdat$ppD,TVdat$ppT,row.names(TVdat),pos=T)
## identify() allows you to mark observations in the data with a left mouse click
## end the identify session with a right mouse click


###################################################
### chunk number 17: ch1
###################################################
#line 348 "Lecture2.Rnw"
## this code is just a way for me to allow for output to be put in my handouts on the go
## the if-statement allows for the figure to be created even if there are no observations marked
plot(TVdat$ppD,TVdat$ppT,xlab="people per Dr",ylab="people per TV")
if (max(id$ind)!=-Inf) {
text(TVdat$ppD[id$ind],TVdat$ppT[id$ind],row.names(TVdat)[id$ind],pos=id$pos) }


###################################################
### chunk number 18: ch2a
###################################################
#line 366 "Lecture2.Rnw"
plot(log(TVdat$ppD),TVdat$ppT) ## taking logs of the x-variable
id<-identify(log(TVdat$ppD),TVdat$ppT,row.names(TVdat),pos=T)


###################################################
### chunk number 19: ch2
###################################################
#line 370 "Lecture2.Rnw"
plot(log(TVdat$ppD),TVdat$ppT,xlab="people per Dr",ylab="people per TV")
if (max(id$ind)!=-Inf) {text(log(TVdat$ppD)[id$ind],TVdat$ppT[id$ind],row.names(TVdat)[id$ind],pos=id$pos)}


###################################################
### chunk number 20: ch3a
###################################################
#line 386 "Lecture2.Rnw"
plot(log(TVdat$ppD),log(TVdat$ppT)) ## taking logs of both x and y
id<-identify(log(TVdat$ppD),log(TVdat$ppT),row.names(TVdat),pos=T)
mm<-lm(log(TVdat$ppT)~log(TVdat$ppD))
# missing values for ppT for two countries
lines(sort(log(TVdat$ppD)[is.na(TVdat$ppT)==F]),mm$fit[sort.list(log(TVdat$ppD)[is.na(TVdat$ppT)==F])])
## here I use is.na() to identify the observations with missing values, I only plot those that had
## is.na==F for FALSE


###################################################
### chunk number 21: ch3
###################################################
#line 395 "Lecture2.Rnw"
plot(log(TVdat$ppD),log(TVdat$ppT))
if (max(id$ind)!=-Inf) {text(log(TVdat$ppD)[id$ind],log(TVdat$ppT)[id$ind],row.names(TVdat)[id$ind],pos=id$pos)}
lines(sort(log(TVdat$ppD)[is.na(TVdat$ppT)==F]),mm$fit[sort.list(log(TVdat$ppD)[is.na(TVdat$ppT)==F])])


###################################################
### chunk number 22: ch4
###################################################
#line 411 "Lecture2.Rnw"
## here I create a regression summary in a LateX table format
## you can just use the command summary(mm) in R directly
library(xtable)
xtable(summary(mm), caption="Regression summary", label="tab:ch4")


###################################################
### chunk number 23: ch5a
###################################################
#line 420 "Lecture2.Rnw"
induse<-seq(1,dim(TVdat)[1])[is.na(TVdat$ppT)==F] ##id the data with non-missing y-values
plot(log(TVdat$ppD)[induse],mm$res)
abline(h=0)
id<-identify(log(TVdat$ppD)[induse],mm$res,row.names(TVdat)[induse],pos=T) ## mark large residuals


###################################################
### chunk number 24: ch5
###################################################
#line 426 "Lecture2.Rnw"
plot(log(TVdat$ppD)[induse],mm$res)
abline(h=0)
if (max(id$ind)!=-Inf) {text((log(TVdat$ppD)[induse])[id$ind],mm$res[id$ind],(row.names(TVdat)[induse])[id$ind],pos=id$pos)}


###################################################
### chunk number 25: ch6a
###################################################
#line 431 "Lecture2.Rnw"
lmi<-lm.influence(mm)
plot(log(TVdat$ppD)[induse],lmi$hat,ylab="leverage")
id<-identify(log(TVdat$ppD)[induse],lmi$hat, row.names(TVdat)[induse],pos=T)


###################################################
### chunk number 26: ch6
###################################################
#line 436 "Lecture2.Rnw"
plot(log(TVdat$ppD)[induse],lmi$hat,ylab="leverage")
if (max(id$ind)!=-Inf) {text((log(TVdat$ppD)[induse])[id$ind],lmi$hat[id$ind], (row.names(TVdat)[induse])[id$ind],pos=id$pos)}
p<-locator()


###################################################
### chunk number 27: ch7a
###################################################
#line 452 "Lecture2.Rnw"
plot(induse,lmi$coef[,2],ylab="Impact on Slope")
abline(h=0)
id<-identify(induse,lmi$coef[,2],label=row.names(TVdat)[induse],pos=T)


###################################################
### chunk number 28: ch7
###################################################
#line 457 "Lecture2.Rnw"
plot(induse,lmi$coef[,2],ylab="Impact on Slope")
abline(h=0)
if (max(id$ind)!=-Inf) {text(induse[id$ind],(lmi$coef[,2])[id$ind],label=(row.names(TVdat)[induse])[id$ind],pos=id$pos)}


###################################################
### chunk number 29: ch8a
###################################################
#line 462 "Lecture2.Rnw"
plot(induse,lmi$sig,ylab="Impact on Sum of Squares")
id<-identify(induse,lmi$sig,label=row.names(TVdat)[induse],pos=T)


###################################################
### chunk number 30: ch8
###################################################
#line 466 "Lecture2.Rnw"
plot(induse,lmi$sig,ylab="Impact on Sum of Squares")
if (max(id$ind)!=-Inf) {text(induse[id$ind],lmi$sig[id$ind],label=(row.names(TVdat)[induse])[id$ind],pos=id$pos)}