> toluca <- read.table("http://www.stat.ufl.edu/~winner/sta4210/data/CH01TA01.txt",header=F,col.names=c("X","Y"))
> 
> attach(toluca)
The following object(s) are masked from 'toluca (position 5)':

    X, Y
The following object(s) are masked from 'toluca (position 6)':

    X, Y
The following object(s) are masked from 'toluca (position 7)':

    X, Y
> 
> 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) 
 lag Autocorrelation D-W Statistic p-value
   1       0.2593193       1.43179   0.132
 Alternative hypothesis: rho != 0
> 
> # 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)
[1] 0.9915055
> (critval <- 1.02-1/sqrt(10*length(Y)))
[1] 0.9567544
> 
> 
> # Shapiro-Wilk Normality test on Residuals
> 
> shapiro.test(Residuals)

        Shapiro-Wilk normality test

data:  Residuals 
W = 0.9789, p-value = 0.8626

> 
> # 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))))
[1] 1.316482
> (t_crit <- qt(.975,n_d1+n_d2-2))
[1] 2.068658
> (p_BF <- 2*(1-pt(abs(t_BF),n_d1+n_d2-2)))
[1] 0.2009812
> 
> 
> 
> 
> # 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
[1] 174292038
> 
> toluca.reg2 <- lm(E2 ~ X)                    # Fit regression of e^2 on X
> summary(toluca.reg2)                         

Call:
lm(formula = E2 ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-2760.1 -1765.4  -990.7   609.5  7726.3 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  3590.91    1442.14   2.490   0.0204 *
X             -19.97      19.12  -1.045   0.3070  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2690 on 23 degrees of freedom
Multiple R-squared: 0.0453,     Adjusted R-squared: 0.003796 
F-statistic: 1.091 on 1 and 23 DF,  p-value: 0.307 

> (SSE_E2 <- deviance(toluca.reg2))            # Obtain SSE from regression of e^2 on X
[1] 166395896
> (SSR_E2 <- SSE2E2 - SSE_E2)                  # Compute SSR*
[1] 7896142
> 
> (X2_BP <- (SSR_E2/2)/(sum(E2)/length(E2))^2)   # Compute Breusch-Pagan test statistic
[1] 0.8209192
> (X2_crit <- qchisq(.95,1))                     # Obtain critical value
[1] 3.841459
> (p_BP <- 1-pchisq(X2_BP,1))                    # Compute P-value
[1] 0.3649116
> 
> 
> # Using lmtest package
> 
> #install.packages("lmtest")   # You must have already Set CRAN Mirror under Packages Tab
> library(lmtest)
> 
> bptest(Y ~ X,studentize=FALSE)

        Breusch-Pagan test

data:  Y ~ X 
BP = 0.8209, df = 1, p-value = 0.3649