-
Notifications
You must be signed in to change notification settings - Fork 0
/
lab1.R
55 lines (43 loc) · 1.97 KB
/
lab1.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
toluca <- read.table("toluca.txt" , sep = "" , header = T)
lm_toluca <- lm(workHours ~ lotSize , data = toluca)
plot(toluca$lotSize , toluca$workHours , xlab = "Lot Size" ,
ylab = "Work Hours" , main = "Toluca Company: Refrigeration Manufacturer")
abline(lm_toluca)
summary(lm_toluca)
residuals(lm_toluca)
fitted(lm_toluca)
sqrt(sum((toluca$workHours - fitted.values(lm_toluca))^2)/(length(toluca$workHours)-2)) # sqrt(MSE)
sqrt(sum(residuals(lm_toluca)^2)/(length(toluca$workHours)-2)) # sqrt(MSE)
b1 = (sum((toluca$lotSize - mean(toluca$lotSize))*(toluca$workHours - mean(toluca$workHours))))/(sum((toluca$lotSize - mean(toluca$lotSize))^2))
MSE = sum(residuals(lm_toluca)^2)/(length(toluca$workHours)-2)
s2b1 = MSE/(sum((toluca$lotSize - mean(toluca$lotSize))^2))
(t_star = b1/sqrt(s2b1))
qt(1-.05/2, length(toluca$workHours-2))
ifelse(abs(t_star) <= qt(1-.05/2, length(toluca$workHours-2)), "Do not reject H0", "Reject H0")
confint(lm_toluca)
cbind(Estimate=coef(lm_toluca), confint(lm_toluca))
b0 = mean(toluca$workHours) - b1*mean(toluca$lotSize)
s2b0 = MSE*(1/length(toluca$lotSize) + mean(toluca$lotSize/(sum((toluca$lotSize - mean(toluca$lotSize))^2))))
(t_star = b0/sqrt(s2b0)) #something wrong here
qt(1-.05/2, length(toluca$lotSize-2))
predict(lm_toluca, data.frame(lotSize=33), interval="confidence")
predict(lm_toluca, data.frame(lotSize=c(33,99)), interval="confidence")
summary(lm_toluca)$fstatistic
anova(lm_toluca)
summary(lm_toluca)$r.squared
summary(lm_toluca)$coefficients
residplot <- function(fit, nbreaks=10) {
z <- rstudent(fit)
hist(z, breaks=nbreaks, freq=FALSE,
xlab="Studentized Residual",
main="Distribution of Errors")
rug(jitter(z), col="brown")
curve(dnorm(x, mean=mean(z), sd=sd(z)),
add=TRUE, col="blue", lwd=2)
lines(density(z)$x, density(z)$y,
col="red", lwd=2, lty=2)
legend("topright",
legend = c( "Normal Curve", "Kernel Density Curve"),
lty=1:2, col=c("blue","red"), cex=.7)
}
residplot(lm_toluca)