################################################### ### chunk number 1: fig1plot ################################################### #line 68 "Lecture1.Rnw" n<-50 ## number of observations noisesd<-2 ## noise level x<-sort(rnorm(n)) ## the independent variable y<-2*x-.5*x^2+rnorm(n,sd=noisesd) ## generating data plot(x, y) ## scatter plot lines(x,2*x-.5*x^2,lwd=2) ## plotting true model mm<-lm(y~x) ## linear model lines(x,mm$fit,lwd=2,col=3) mm<-lm(y~x+I(x^2)) ## quadratic model lines(x,mm$fit,lty=2,lwd=2,col=2) mm3<-loess(y~x,degree=0,span=.1) ## local average model lines(x,mm3$fit,lty=5,lwd=3,col=4) ################################################### ### chunk number 2: fig2plot ################################################### #line 82 "Lecture1.Rnw" ## use same x y2<-2*x-.5*x^2+rnorm(n,sd=noisesd) plot(x, y2,lwd=2) lines(x,2*x-.5*x^2) mm<-lm(y2~x) lines(x,mm$fit,lwd=2,col=3) mm<-lm(y2~x+I(x^2)) lines(x,mm$fit,lty=2,lwd=2,col=2) mm3<-loess(y2~x,degree=0,span=.1) lines(x,mm3$fit,lty=5,lwd=3,col=4) ################################################### ### chunk number 3: fig3plot ################################################### #line 94 "Lecture1.Rnw" y3<-2*x-.5*x^2+rnorm(n,sd=noisesd) plot(x, y3) lines(x,2*x-.5*x^2,lwd=2) mm<-lm(y3~x) lines(x,mm$fit,lwd=2,col=3) mm<-lm(y3~x+I(x^2)) lines(x,mm$fit,lty=2,lwd=2,col=2) mm3<-loess(y3~x,degree=0,span=.1) lines(x,mm3$fit,lty=5,lwd=3,col=4) ################################################### ### chunk number 4: fig4plot ################################################### #line 121 "Lecture1.Rnw" K<-50 ## how many data sets to generate mma<-rep(0,2) ## creating a vector to store the model coefficients in ## The for-loop runs everything inside the curly brackets K times for (kk in (1:K)) { y<-2*x-.5*x^2+rnorm(n,sd=noisesd) if (kk==1 ) {plot(x, y)} ## plot the data once to get the axes fixed mm<-lm(y~x) lines(x,mm$fit,col=3,lty=2) ## plotting the different fites in col=3 (green). mma<-mma+(1/K)*mm$coef } ## storing the coefficients lines(x,2*x-.5*x^2,lwd=2) ## true model lines(x,mma[1]+mma[2]*x,lwd=2,col=2) ## average model ################################################### ### chunk number 5: fig5plot ################################################### #line 134 "Lecture1.Rnw" mma<-rep(0,3) for (kk in (1:K)) { y2<-2*x-.5*x^2+rnorm(n,sd=noisesd) if (kk==1) {plot(x, y2,lwd=2)} mm<-lm(y2~x+I(x^2)) lines(x,mm$fit,col=3,lty=2) mma<-mma+(1/K)*mm$coef} lines(x,2*x-.5*x^2,lwd=2) lines(x,mma[1]+mma[2]*x+mma[3]*x^2,lwd=2,col=2) ################################################### ### chunk number 6: fig6plot ################################################### #line 145 "Lecture1.Rnw" mma<-rep(0,length(x)) for (kk in (1:K)) { y3<-2*x-.5*x^2+rnorm(n,sd=noisesd) if (kk==1) {plot(x, y3) } mm3<-loess(y3~x,degree=0,span=.1) lines(x,mm3$fit,col=3,lty=2) mma<-mma+(1/K)*mm3$fit } lines(x,2*x-.5*x^2,lwd=2) lines(x,mma,lwd=2,col=2) ################################################### ### chunk number 7: LS ################################################### #line 196 "Lecture1.Rnw" x<-sort(rnorm(25)) y<-2+3*x+rnorm(25) plot(x,y) mm<-lm(y~x) mm2<-3+2*x lines(x,mm$fit,lwd=2,col=2) lines(x,mm2,col=4) for (zz in (1:length(x))) { lines(rep(x[zz],2),c(y[zz],mm$fit[zz])) ## plotting the vertical distance to the model = deviations lines(rep(x[zz],2),c(y[zz],mm2[zz]),lty=2) } ################################################### ### chunk number 8: B1a ################################################### #line 223 "Lecture1.Rnw" x<-sort(rnorm(50)) y1<-2+3*x+rnorm(50) plot(x,y1) lines(x,lm(y1~x)$fit) lines(x,2+3*x,lwd=2,col=2) ################################################### ### chunk number 9: B1b ################################################### #line 230 "Lecture1.Rnw" y2<-y1 y2[26:50]<-rnorm(25,sd=.5) ## creating a subgroup in the data plot(x,y2) lines(x,lm(y2~x)$fit) lines(x[1:25],2+3*x[1:25],col=2,lwd=2) lines(x[26:50],rep(0,25),col=2,lwd=2) ################################################### ### chunk number 10: B1c ################################################### #line 238 "Lecture1.Rnw" y3<-1+2*x-x^2+rnorm(50) ## true model is quadratic plot(x,y3) lines(x,lm(y3~x)$fit) lines(x,1+2*x-x^2,col=2,lwd=2) ################################################### ### chunk number 11: b2a ################################################### #line 265 "Lecture1.Rnw" x<-sort(rnorm(50)) y1<-10+3*x+rnorm(50) plot(x,y1) lines(x,lm(y1~x)$fit) lines(x,10+3*x,col=2,lwd=2) ################################################### ### chunk number 12: b2b ################################################### #line 272 "Lecture1.Rnw" y2<-10+3*x+(rchisq(50,1)) plot(x,y2) lines(x,lm(y2~x)$fit) lines(x,10+3*x,col=2,lwd=2) ################################################### ### chunk number 13: b2c ################################################### #line 278 "Lecture1.Rnw" #y2<-10+3*x+(rchisq(50,1)) plot(x,y2) lines(x,lm(y2~x)$fit) lines(x,10+3*x,col=2,lwd=2) lines(x,exp(lm(log(y2)~x)$fit),lty=2) ################################################### ### chunk number 14: b2d ################################################### #line 285 "Lecture1.Rnw" #y2<-10+3*x+(rchisq(50,1)) plot(x,log(y2)) #lines(x,log(10+3*x),col=2,lwd=2) lines(x,(lm(log(y2)~x)$fit)) ################################################### ### chunk number 15: b4a ################################################### #line 321 "Lecture1.Rnw" x<-sort(rnorm(100)) for (kk in (1:25)) { y1<-10+3*x+rnorm(100,sd=2) if (kk==1) {plot(x,y1)} lines(x,lm(y1~x)$fit) } lines(x,10+3*x,col=2,lwd=2) ################################################### ### chunk number 16: b4b ################################################### #line 329 "Lecture1.Rnw" for (kk in (1:25)) { y2<-10+3*x+rnorm(100,sd=seq(.1,6,length.out=100)) ## creating a data set with increasing error variance if (kk==1) {plot(x,y2)} lines(x,lm(y2~x)$fit) } lines(x,10+3*x,col=2,lwd=2) ################################################### ### chunk number 17: b4c ################################################### #line 336 "Lecture1.Rnw" y2<-10+3*x+rnorm(100,sd=seq(.1,6,length.out=100)) plot(x,(y2)) ################################################### ### chunk number 18: b4d ################################################### #line 340 "Lecture1.Rnw" y2<-10+3*x+rnorm(100,sd=seq(.1,6,length.out=100)) ## taking logs to try to even out the variance plot(x,log(y2)) ################################################### ### chunk number 19: b5a ################################################### #line 369 "Lecture1.Rnw" x<-sort(rnorm(50)) y1<-10+3*x+rnorm(50,sd=2) plot(x,y1) lines(x,lm(y1~x)$fit) lines(x,10+3*x,col=2,lwd=2) ################################################### ### chunk number 20: b5b ################################################### #line 376 "Lecture1.Rnw" xnew<-x xnew[50]<-xnew[50]+10 ## creating an extreme x-value y2<-10+3*xnew+rnorm(50,sd=2) plot(xnew,y2) lines(xnew,lm(y2~xnew)$fit) lines(xnew,10+3*xnew,col=2,lwd=2) points(xnew[50],y2[50],pch=8,col=2) ################################################### ### chunk number 21: b5c ################################################### #line 385 "Lecture1.Rnw" y2<-10+3*x+rnorm(50,sd=2) y2[50]<-y2[50]-20 ## creating an extreme y-value for an extreme x plot(x,y2) lines(x,lm(y2~x)$fit) lines(x,10+3*x,col=2,lwd=2) points(x[50],y2[50],pch=8,col=2) ################################################### ### chunk number 22: b5d ################################################### #line 393 "Lecture1.Rnw" y2<-10+3*x+rnorm(50,sd=2) y2[20]<-y2[20]-15 ## creating an extreme y-value for an x near mean(x) plot(x,y2) lines(x,lm(y2~x)$fit) lines(x,10+3*x,col=2,lwd=2) points(x[20],y2[20],pch=8,col=2) ################################################### ### chunk number 23: b6a ################################################### #line 436 "Lecture1.Rnw" x<-sort(rnorm(25)) x2<-c(x,sort(rnorm(25))+5) y1<-10+3*x+rnorm(25,sd=2) y2<-3+x2+rnorm(25) ## creating 2 groups of responses plot(c(x,x2),c(y1,y2)) ################################################### ### chunk number 24: b6b ################################################### #line 443 "Lecture1.Rnw" x<-sort(rexp(50,.1)) ## uneven spread in x y2<-10+3*x+rnorm(50,sd=5) plot(x,y2) ################################################### ### chunk number 25: b6c ################################################### #line 448 "Lecture1.Rnw" x<-rnorm(50) y2<-10+3*x+rnorm(50,sd=2) plot(x,y2,xlim=c(-5,10),ylim=c(-5,40)) lines(x,lm(y2~x)$fit) lines(seq(-5,10),10+3*seq(-5,10),col=2,lwd=2) ## looking at the model outside observed range of x ################################################### ### chunk number 26: ch1a ################################################### #line 483 "Lecture1.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 27: ch1 ################################################### #line 493 "Lecture1.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 28: ch2a ################################################### #line 511 "Lecture1.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 29: ch2 ################################################### #line 515 "Lecture1.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 30: ch3a ################################################### #line 531 "Lecture1.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 31: ch3 ################################################### #line 540 "Lecture1.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 32: ch4 ################################################### #line 556 "Lecture1.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")