toluca <- read.table("http://www.stat.ufl.edu/~winner/sta4210/data/CH01TA01.txt",header=F,col.names=c("X","Y")) attach(toluca) toluca n <- length(Y) mean_x <- mean(X); mean_y <- mean(Y) var_x <- var(X); var_y <- var(Y); cov_xy <- cov(X,Y) SS_xx <- (n-1)*var_x SS_xy <- (n-1)*cov_xy SS_yy <- (n-1)*var_y b1 <- SS_xy/SS_xx b0 <- mean_y - b1*mean_x yhat <- b0 + b1*X e <- Y-yhat SSE <- sum(e^2) MSE <- SSE/(n-2) s <- sqrt(MSE) s_b1 <- s/sqrt(SS_xx) s_b0 <- s*sqrt((1/n)+(mean_x^2/SS_xx)) t_b1 <- b1/s_b1 t_b0 <- b0/s_b0 p_b1 <- 2*(1-pt(abs(t_b1),n-2)) p_b0 <- 2*(1-pt(abs(t_b0),n-2)) lb_b1_95 <- b1 - qt(.975,n-2)*s_b1; ub_b1_95 <- b1 + qt(.975,n-2)*s_b1; lb_b0_95 <- b0 - qt(.975,n-2)*s_b0; ub_b0_95 <- b0 + qt(.975,n-2)*s_b0; print(cbind(b0,s_b0,t_b0,p_b0,lb_b0_95,ub_b0_95)) print(cbind(b1,s_b1,t_b1,p_b1,lb_b1_95,ub_b1_95)) toluca.reg1 <- lm(Y ~ X) summary(toluca.reg1) anova(toluca.reg1) pdf("F:\\sta4210\\KLLN_TA2_1.pdf") plot(X,Y,xlim=c(0,125),ylim=c(0,600),main="Toluca Data - Confidence Bands and Loess") #abline(a=b0,b=b1) lines(lowess(X,Y)) Xh <- 0:125 lines(Xh,b0+b1*Xh - sqrt(2*qf(.95,2,n-2))*s*sqrt((1/n)+((Xh-mean_x)^2/SS_xx))) lines(Xh,b0+b1*Xh + sqrt(2*qf(.95,2,n-2))*s*sqrt((1/n)+((Xh-mean_x)^2/SS_xx))) dev.off() SSTO <- SS_yy; df_TO <- n-1 SSR <- SSTO - SSE; df_R <- 1; MSR <- SSR/df_R df_E <- n-2 F_stat <- MSR/MSE F_05 <- qf(.95,1,n-2) p_F <- 1-pf(F_stat,1,n-2) print(cbind(SSR,df_R,MSR,F_stat,F_05,p_F)) print(cbind(SSE,df_E,MSE)) print(cbind(SSTO,df_TO))