pdf("E:\\Rmisc\\graphs\\paddyfields_slr.pdf") bearcap <- read.table("http://www.stat.ufl.edu/~winner/data/paddyfields_slr.dat", header=F,col.names=c("depth","capacity")) attach(bearcap) # Makes data "available" for direct computations bearcap # Print the data (mean_X <- mean(depth)) # Obtain X-bar ... Parentheses has it printed out (mean_Y <- mean(capacity)) # Obtain Y-bar ... Parentheses has it printed out (n <- length(depth)) # Obtain the sample size, n (SS_XX <- (n-1)*var(depth)) # Obtain SS_XX (SS_YY <- (n-1)*var(capacity)) # Obtain SS_YY (SS_XY <- (n-1)*cov(depth,capacity)) # Obtain SS_XY (b1 <- SS_XY/SS_XX) # Obtain b1 (b0 <- mean_Y-b1*mean_X) # Obtain b0 yhat <- b0 + b1*depth # Obtain fitted values (one for each observation) e <- capacity - yhat # Obtain residuals (one for each observation) (SSE <- sum(e^2)) # Obtain the Error Sum of Squares (MSE <- SSE/(n-2)) # Obtain the Mean Square Error - estimate of sigma^2 (s_b1 <- sqrt(MSE/SS_XX)) # Obtain the standard error of b1 - s{b1} (s_b0 <- sqrt(MSE*((1/n)+((mean_X^2)/SS_XX)))) # Obtain s{b0} (t_b1 <- b1/s_b1) # Obtain t-statistic for H_0: beta_1 = 0 (t_b0 <- b0/s_b0) # Obtain t-statistic for H_0: beta_0 = 0 (p_b1 <- 2*(1-pt(abs(t_b1),n-2))) # Obtain 2-sided P-value for H_0: beta_1 = 0 (p_b0 <- 2*(1-pt(abs(t_b0),n-2))) # Obtain 2-sided P-value for H_0: beta_0 = 0 # Note: pt(a,df) = P(t(df) <= a) (ci_b1 <- c(b1-qt(.975,n-2)*s_b1,b1+qt(.975,n-2)*s_b1)) # Obtain 95% CI for beta_1 (ci_b0 <- c(b0-qt(.975,n-2)*s_b0,b0+qt(.975,n-2)*s_b0)) # Obtain 95% CI for beta_0 # Note: qt(a,df) = t(a;df) (SSTO <- SS_YY) # Obtain the Total Sum of Squares (SSR <- SSTO - SSE) # Obtain the Regression Sum of squares (F_stat <- (SSR/1)/MSE) # Compute the F-Stat for H_0: beta_1 = 0 (p_F <- 1-pf(F_stat,1,n-2)) # Obtain p-value (2-sided) for H_0: beta_1 = 0 # Note: pf(a,df1,df2) = P(F(df1,df2) <= a (R_square <- SSR/SSTO) # Obtain Coefficient of Determination - R^2 (R <- SS_XY/sqrt(SS_XX*SS_YY)) # Obtain Pearson Product Moment Coefficient of Correlation plot(depth,capacity) abline(a=b0,b=b1) ################################################################################################## ## Now running analysis with linear model (lm) function ######################################### ################################################################################################## bearcap.reg1 <- lm(capacity ~ depth) # Run the linear regression with Y=capacity, X=depth summary(bearcap.reg1) # Print out summary anova(bearcap.reg1) # Print out ANOVA table bearcap.regsum <- summary(bearcap.reg1) # Retrieve coefs, SE's from model (b0_lm <- coef(bearcap.regsum)[1,1]) # Retrieve b0 from the model (b1_lm <- coef(bearcap.regsum)[2,1]) # Retrieve b1 from the model (s_b0_lm <- coef(bearcap.regsum)[1,2]) # Retrieve s{b0} from the model (s_b1_lm <- coef(bearcap.regsum)[2,2]) # Retrieve s{b1} from the model plot(depth,capacity) abline(lm(capacity ~ depth)) dev.off()