## Commands for chapter 8 lecture ## read in the data reg.data <- read.table("E:\Uvic folder\359\R code\c8\tannin.txt",header=T) attach(reg.data) names(reg.data) ## scatter plot plot(tannin,growth,pch=16) ################################################ ## rough guess at parameter values a<-12 b<-(growth[9]-growth[1])/(tannin[9]-tannin[1]) ## plot the guessed linear relationship lines(tannin,a+b*tannin,lty=2) ################################################# ## plot the errors simple.fit<-a+b*tannin for (i in 1:length(growth)) { lines(c(tannin[i],tannin[i]),c(simple.fit[i],growth[i])) } ################################################# ## compute least squares estimators ## change notation to match the notes y<-growth x<-tannin b.ls<-(sum(x*y) - sum(y)*sum(x)/length(y))/(sum(x^2) - sum(x)*sum(x)/length(x)) a.ls<-mean(y) - b.ls*mean(x) ## plot the least squares line plot(tannin,growth,pch=16) ## plot the guessed linear relationship lines(tannin,a+b*tannin,lty=2) ## plot the LS line lines(tannin,a.ls+b.ls*tannin,lty=2,col="red") ## can compute these using the lm() function lm(growth~tannin) ############################################3 ## error sum of squares SSE<-sum((y-a.ls-b.ls*x)^2) ## total sum of squares SSY<-sum(y^2) - (sum(y)^2)/length(y) ## regression SS SSR<-SSY-SSE ##########################################3 ## compute the p-value for F-test F.ratio<-(SSR/1)/(SSE/7) 1-pf(F.ratio,1,7) ## can produce the anova table in R model<-lm(growth~tannin) summary.aov(model) ################################ ## examine the standard errors and other output summary(model) ## want a 95% confidence interval for b lower.b<-b.ls-0.2186*qt(p=0.975,df=length(y)-2) upper.b<-b.ls+0.2186*qt(p=0.975,df=length(y)-2) CI.b<-c(lower.b,upper.b) ## want a 95% confidence interval for a lower.a<-a.ls-1.0408*qt(p=0.975,df=length(y)-2) upper.a<-a.ls+1.0408*qt(p=0.975,df=length(y)-2) CI.a<-c(lower.a,upper.a) ################################################# ## regression diagnostics par(mfrow=c(1,3)) plot(model,which=c(1,2,4)) ################################################ ## sensitivity analysis investigating impact of most influential point c(y[7],x[7]) # remove this point and re-fit the model model2<-update(model,subset=(tannin !=6)) summary(model) summary(model2) ################################################# ## examine data on the relationship between radioactive emissions (y) and time (x) ## this data is collected for a radioactive decay process detach(reg.data) rm(x,y) ## remove variables x and y from gloabal env par(mfrow=c(1,1)) curve <- read.table("E:\Uvic folder\359\R code\c8\decay.txt",header=T) attach(curve) names(curve) plot(x,y) abline(lm(y~x)) ##examine r^2 summary(lm(y~x)) ##########################################3 ## fit a quadratic x2<-x^2 quadratic<-lm(y~x+x2) plot(x,y) xv<-seq(0,30,.1) ## plot the line linear<-lm(y~x) yv.linear<-predict(linear,list(x=xv)) lines(xv,yv.linear) ## plot the quadratic model yv.quadratic<-106.3888-7.3448*xv+0.1506*xv^2 yv.quadratic<-predict(quadratic,list(x=xv,x2=xv^2)) lines(xv,yv.quadratic,col="red") ## testing summary(quadratic) ## or can use anova table anova(quadratic,linear) ########################################## ## compare a seqence of four nested models constant<-lm(y~1) linear<-lm(y~x) x3<-x^3 cubic<-lm(y~x+x2+x3) ## compare 4 nested models ## the F-statistic provides pairwise comparison of models ## numerator = difference in RSS/difference in DF ## denominator = RSS of largest model/DF anova(constant,linear,quadratic,cubic) ############################################# exponential<-lm(log(y)~x) summary(exponential) ## make predicitions on the log scale and exponentiate yv.exp<-exp(predict(exponential,list(x=xv))) lines(xv,yv.exp,col="blue") ## nice to add a legend