################################################### ### chunk number 1: l3-fitlev ################################################### #line 42 "Lecture3.Rnw" x<-rnorm(50) mma<-rep(0,2) for (kk in (1:50)) { y<-3+2*x+rnorm(50) if (kk==1 ) {plot(x, y)} mm<-lm(y~x) lines(x,mm$fit,col=3,lty=2) mma<-mma+(1/50)*mm$coef } lines(x,3+2*x,lwd=2) lines(x,mma[1]+mma[2]*x,lwd=2,col=2) ################################################### ### chunk number 2: l3-stdres ################################################### #line 92 "Lecture3.Rnw" x<-rnorm(100) x<-sort(x) y<-3+2*x+rnorm(100) mm<-lm(y~x) plot(x,mm$res) abline(h=0,lty=2) points(x,rstandard(mm),pch=2,col=2) ################################################### ### chunk number 3: l3-varcomp ################################################### #line 132 "Lecture3.Rnw" x<-rnorm(1000) x<-sort(x) y<-3+2*x+rnorm(1000) mm<-lm(y~x) plot(x,y) lines(x,mm$fit,lwd=2) ## here I'm marking a vertical slot between the 60th and 80th percentile of x abline(v=quantile(x,.6),lwd=2,lty=2) abline(v=quantile(x,.8),lwd=2,lty=2) ## here I mark the 2.5 and 97.5 percentiles of all y (col=4, blue) abline(h=quantile(y,.025),lwd=2,col=4) abline(h=quantile(y,.975),lwd=2,col=4) ysub<-y[xquantile(x,.6)] xsub<-x[xquantile(x,.6)] xm<-mean(xsub) points(xsub,ysub,col=2,pch=8) ## here I mark the percentiles of the y's in the vertical slice abline(h=quantile(ysub,.025),lwd=2,col=2) abline(h=quantile(ysub,.975),lwd=2,col=2) text(quantile(x,.05),quantile(y,.99),"Spread in y",col="blue") text(xm,min(y),"x*") text(quantile(x,.025),quantile(ysub,.85),"Spread in y given x near x*",col="red") ################################################### ### chunk number 4: l3-varcompb ################################################### #line 167 "Lecture3.Rnw" x<-rnorm(1000) x<-sort(x) y<-3+.5*x+rnorm(1000) mm<-lm(y~x) plot(x,y) lines(x,mm$fit,lwd=2) abline(v=quantile(x,.6),lwd=2,lty=2) abline(v=quantile(x,.8),lwd=2,lty=2) abline(h=quantile(y,.025),lwd=2,col=4) abline(h=quantile(y,.975),lwd=2,col=4) ysub<-y[xquantile(x,.6)] xsub<-x[xquantile(x,.6)] xm<-mean(xsub) points(xsub,ysub,col=2,pch=8) abline(h=quantile(ysub,.025),lwd=2,col=2) abline(h=quantile(ysub,.975),lwd=2,col=2) text(quantile(x,.025),quantile(y,.99),"Spread in y",col="blue") text(xm,min(y),"x*") text(quantile(x,.1),quantile(ysub,.85),"Spread in y given x near x*",col="red") ################################################### ### chunk number 5: l3-r2a ################################################### #line 209 "Lecture3.Rnw" x<-rnorm(1000) x<-sort(x) y<-3+.5*x+rnorm(1000) mm<-lm(y~x) plot(x,y) lines(x,mm$fit,lwd=2) rsqa<-round(summary(mm)$r.sq,2) text(quantile(x,.1),quantile(y,.99),paste("R-sq =",rsqa),col="red") ################################################### ### chunk number 6: l3-r2b ################################################### #line 219 "Lecture3.Rnw" x<-rnorm(1000) x<-sort(x) y<-2+5*x+rnorm(1000) mm<-lm(y~x) plot(x,y) lines(x,mm$fit,lwd=2) rsqb<-round(summary(mm)$r.sq,2) text(quantile(x,.1),quantile(y,.99),paste("R-sq =",rsqb),col="red") ################################################### ### chunk number 7: l3-r2c ################################################### #line 229 "Lecture3.Rnw" x<-rnorm(1000) x<-sort(x) y<-2+.5*x+rnorm(1000,sd=.2) mm<-lm(y~x) rsqc<-round(summary(mm)$r.sq,2) plot(x,y) lines(x,mm$fit,lwd=2) text(quantile(x,.1),quantile(y,.99),paste("R-sq =",rsqc),col="red") ################################################### ### chunk number 8: ch1a ################################################### #line 257 "Lecture3.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 9: ch1 ################################################### #line 267 "Lecture3.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 10: ch2a ################################################### #line 285 "Lecture3.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 11: ch2 ################################################### #line 289 "Lecture3.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 12: ch3a ################################################### #line 305 "Lecture3.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 13: ch3 ################################################### #line 314 "Lecture3.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 14: ch4 ################################################### #line 330 "Lecture3.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 15: ch5a ################################################### #line 339 "Lecture3.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 16: ch5 ################################################### #line 345 "Lecture3.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 17: ch6a ################################################### #line 350 "Lecture3.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 18: ch6 ################################################### #line 355 "Lecture3.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 19: ch7a ################################################### #line 371 "Lecture3.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 20: ch7 ################################################### #line 376 "Lecture3.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 21: ch8a ################################################### #line 381 "Lecture3.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 22: ch8 ################################################### #line 385 "Lecture3.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)} ################################################### ### chunk number 23: ch9 ################################################### #line 399 "Lecture3.Rnw" print(summary(mm))