homework 2 mat613

data <- read.csv("/Users/caelynsobie/Downloads/MAT613/plastic-hardness.csv")

# a) construct a 98% confidence interval for the mean hardness with an elapsed time of 30 hours
# y_h^hat +- t(1-a/2,n-2)*SE(Y_h^hat)
y <- data$Y
x <- data$X
model <- lm(y~x)
av <- anova(model)
n <- length(x)
a = .02

x_mean <- mean(x)
y_mean <- mean(y)

b_1 <- (sum((x - x_mean)*(y - y_mean))) / (sum((x - x_mean)^2))
b_0 <- y_mean - (b_1 * (x_mean))
  
x_h <- 30
y_hat <- b_0 + (b_1*x_h)
t <- qt(1-a/2,n-2)

mse_ci <- av$"Mean Sq"[2]


sxx <- sum((x - x_mean)^2)

se <- sqrt(mse_ci*((1/n)+ ((x_h-x_mean)^2/(sxx)) ))

conf_u <- y_hat + (t*se)
conf_l <- y_hat - (t*se)
print(paste0("Confidence Interval: (", conf_l,", ", conf_u, ")"))

[1] "Confidence Interval: (227.456928039899, 231.805571960101)"

# b) construct a 98% prediction interval for the hardness of a newly molded test item with an elapsed time of 30 hours
# y_h^hat +- t(1-a/2,n-2)*SE(pred)

y_hat <- b_0 + (b_1*x_h)
t <- qt(1-a/2,n-2)

mse_pi <- av$"Mean Sq"[2]
sxx <- sum((x - x_mean)^2)
se <- sqrt(mse_pi*(1+(1/n)+ ((x_h-x_mean)^2/(sxx)) ))

pred_u <- y_hat + (t*se)
pred_l <- y_hat - (t*se)

print(paste0("Prediction Interval: (", pred_l,", ", pred_u, ")"))

[1] "Prediction Interval: (220.869489165647, 238.393010834353)"

#d) plot the boundary values of the 98% Working-Hoteling confidence band for the regression line

av <- anova(model)

y_hat <- b_0 + (b_1*x)
y_hat30 <- b_0 + (b_1*30)
mse_cb <- av$"Mean Sq"[2]
sxx <- sum((x - x_mean)^2)
se <- sqrt(mse_cb*((1/n)+ ((x-x_mean)^2/(sxx)) ))
se_30 <- sqrt(mse_cb*((1/n)+ ((30-x_mean)^2/(sxx)) ))
  
a = 0.02
n = length(x)
f <- qf(1-a, 2, n-2)

w <- sqrt(2*(f))

conf_band_u <- y_hat + (w*se)
conf_band_l <- y_hat - (w*se)
cb_30_u <- y_hat30 + (w*se_30)
cb_30_l <-  y_hat30 - (w*se_30)

print(paste0("Confidence band at X_h = 30: (", cb_30_l,", ", cb_30_u, ")"))

[1] "Confidence band at X_h = 30: (226.94905721092, 232.31344278908)"

plot(x,y, pch=4)

lines(x, conf_band_u, col = "steel blue")  # Upper band

lines(x, conf_band_l, col = "steel blue")  # Lower band
adjustcolor("blue", alpha.f=0.5) 

[1] "#0000FF80"

polygon(c(x, rev(x)), c(conf_band_u, rev(conf_band_l)),
        col = rgb(70/255, 130/255, 180/255, 0.4))

abline(model, col = rgb(95/255, 158/255, 72/255), lwd = 2)
legend("topleft", legend = c("Regression Line", "98% W-H Confidence Band"),
       col = c(rgb(95/255, 158/255, 72/255), "steel blue"), lty = c(1, 2))

#e) set up ANOVA table
SST <- sum((y - y_mean)^2)
SSR <- sum((y_hat - y_mean)^2)
SSE <- sum((y - y_hat)^2)

n = length(x)
df_r = 1
df_e = n-2
df_t = n-1

MSR_av <- SSR/df_r
MSE_av <- SSE/df_e

sov <- c("Regression","Error","Total")

anov <- data.frame(
  "Source of Variation" = sov,
  "SS" = c(SSR,SSE,SST),
  "df" = c(df_r, df_e, df_t),
  "MS" = c(MSR_av,MSE_av,NA)
)
anov

#f) use an F test to test whether or not there is a linear association between the hardness of the plastic and elapsed time
alpha = 0.01
n=length(x)
f_test <- qf(1-alpha, 1, n-2)
mse <- av$"Mean Sq"[2]
msr <- av$"Mean Sq"[1]
f_stat <- msr/mse

print(paste0("F critical-value: ", f_test))

[1] "F critical-value: 8.86159266517642"

print(paste0("F test statistic: ", f_stat))

[1] "F test statistic: 506.506231859314"

ifelse(f_stat > f_test, "Reject H_0","Fail to reject H_0")

[1] "Reject H_0"

Question 5

residuals <- residuals(model)
fit <- fitted(model)
plot(fit, residuals)

plot(x, residuals)

#c) q-q plot

qqnorm(residuals)
qqline(residuals,col="steel blue",lwd=2)

#c) Shapiro-Wilk test
shapiro.test(residuals)

    Shapiro-Wilk normality test

data:  residuals
W = 0.97348, p-value = 0.8914

#d) |residuals| vs Xi

abs_res <- abs(residuals)
plot(x, abs_res)

#d) residuals^2 vs Xi

sq_res <- residuals^2

plot(x, sq_res)

#e) Brown-Forsythe test
g <- rep(1, n)
g[x <= 24] <- 0

bftest(model,g)

       t.value   P.Value alpha df
[1,] 0.8557853 0.4065253  0.05 14

#f) runs test

mean_res <- mean(residuals)
bin_res <- ifelse(residuals > mean_res,1 ,0)
res_factor <- factor(bin_res)

runs.test(res_factor)

    Runs Test

data:  res_factor
Standard Normal = 1.1185, p-value = 0.2633
alternative hypothesis: two.sided