toluca <- read.table("http://www.stat.ufl.edu/~winner/sta4210/data/CH01TA01.txt",header=F,col.names=c("X","Y")) attach(toluca) pdf("F:\\Rmisc\\graphs\\KLLN_TA3_toluca.pdf") freq.X <- table(X) barplot(freq.X,main="Barplot of X Levels") plot(X,type="o",main="X vs Run Number") boxplot(X, main="Boxplot of X") toluca.reg1 <- lm(Y~X) Residuals <- resid(toluca.reg1) plot(X,Residuals,xlab="Lot Size",ylab="Residuals", main="Residuals vs X") boxplot(Residuals,main="Residual boxplot") qqnorm(Residuals,main="Normal Probability Plot of Residuals") qqline(Residuals) plot(toluca.reg1) # Durbin-Watson Test for Autocorrelated Residuals #install.packages("car") # You must have already Set CRAN Mirror under Packages Tab library(car) durbinWatsonTest(toluca.reg1) # Correlation test for Normality of Residuals rank_e <- rank(Residuals) exp_e <- qnorm((rank_e-0.375)/(length(Y)+0.25)) cor(Residuals,exp_e) (critval <- 1.02-1/sqrt(10*length(Y))) # Shapiro-Wilk Normality test on Residuals shapiro.test(Residuals) # Conduct Brown-Forsythe Test of Homogeneity of Variance group <- numeric(length(Y)) for (i in 1:length(Y)) { if (X[i] <= 70) group[i]=1 else group[i]=2 } d1 <- abs(Residuals[group==1] - median(Residuals[group==1])) d2 <- abs(Residuals[group==2] - median(Residuals[group==2])) n_d1 <- length(d1); n_d2 <- length(d2) mean_d1 <- mean(d1); mean_d2 <- mean(d2) var_d1 <- var(d1); var_d2 <- var(d2) s_BF <- sqrt(((n_d1-1)*var_d1 + (n_d2-1)*var_d2)/(n_d1+n_d2-2)) (t_BF <- (mean_d1 - mean_d2)/(s_BF*sqrt((1/n_d1)+(1/n_d2)))) (t_crit <- qt(.975,n_d1+n_d2-2)) (p_BF <- 2*(1-pt(abs(t_BF),n_d1+n_d2-2))) # Conduct Breusch-Pagan Test of Homogeneity of Variance # Brute Force Approach E2 <- Residuals^2 # Compute e^2 for each observation (SSE2E2 <- (length(E2)-1)*var(E2)) # Compute SSTO for e^2 values toluca.reg2 <- lm(E2 ~ X) # Fit regression of e^2 on X summary(toluca.reg2) (SSE_E2 <- deviance(toluca.reg2)) # Obtain SSE from regression of e^2 on X (SSR_E2 <- SSE2E2 - SSE_E2) # Compute SSR* (X2_BP <- (SSR_E2/2)/(sum(E2)/length(E2))^2) # Compute Breusch-Pagan test statistic (X2_crit <- qchisq(.95,1)) # Obtain critical value (p_BP <- 1-pchisq(X2_BP,1)) # Compute P-value # Using lmtest package #install.packages("lmtest") # You must have already Set CRAN Mirror under Packages Tab library(lmtest) bptest(Y ~ X,studentize=FALSE) dev.off()