### Example 0
### Motivating Example
### "Regression with Factors"

y<-c(1.2, 2.3, 1.1, 4.3, 2.8, 2.4, 3.7, 2.7, 4.0)
x<-c(rep("A",5),rep("B",4))

ex0.lm<-lm(y~x)
summary(ex0.lm)
anova(ex0.lm)

plot(y~x)
boxplot(y~x)

# Dummy variable
# Convert to 0 and 1

x
X<-matrix(c(rep(1,9),0,0,0,0,0,1,1,1,1),ncol=2)
X
bhat<-solve(t(X)%*%X)%*%t(X)%*%y
bhat


### Example 1

a <- c(1,1,1,1,2,2,3,3,3,3,3)
y <- c(3,4,5,5,3,2,9,12,5,8,5)

aov(y~a)
aov(y~factor(a))
summary(aov(y~factor(a)))
lm(y~factor(a))
anova(lm(y~factor(a)))
summary(lm(y~factor(a)))

boxplot(y~factor(a))


### Example 2

?lm

ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
anova(lm.D9)
summary(lm.D9)
lm.D90 <- lm(weight ~ group - 1) # omitting intercept
summary(lm.D90)

t.test(ctl,trt,var.equal=T)

X<-cbind(rep(1,20),rep(c(1,0),each=10),
         rep(c(0,1),each=10))
X1<-cbind(rep(1,20),rep(c(0,1),each=10))
X2<-cbind(rep(c(1,0),each=10),
          rep(c(0,1),each=10))
y<-weight
mu.hat<-solve(t(X2)%*%X2)%*%t(X2)%*%y
mu.hat

boxplot(weight ~ group)


### Example 3

N <- c(0,1,0,1,1,1,0,0,0,1,1,0,1,1,0,0,1,0,1,0,1,1,0,0)
K <- c(1,0,0,1,0,1,1,0,0,1,0,1,0,1,1,0,0,0,1,1,1,0,1,0)
yield <- c(49.5,62.8,46.8,57.0,59.8,58.5,55.5,56.0,62.8,
           55.8,69.5,55.0,62.0,48.8,45.5,44.2,52.0,51.5,
           49.8,48.8,57.2,59.0,53.2,56.0)
table(N,K)
anova(lm(yield~factor(N)+factor(K)))
pf(6.7157, 1, 21, lower.tail=F)
pf(3.3778, 1, 21, lower.tail=F)
anova(lm(yield~factor(N)+factor(K)+factor(N):factor(K)))
anova(lm(yield~factor(N)*factor(K)))
interaction.plot(N,K,yield)
interaction.plot(K,N,yield)

P <- c(1,1,0,0,0,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,0,1,1,0)
anova(lm(yield~factor(N)*factor(K)*factor(P)))


### Example 4

warpbreaks
?warpbreaks

wbt.lm<-lm(breaks~tension, data=warpbreaks)
summary(wbt.lm)
anova(wbt.lm)

boxplot(breaks~tension, data=warpbreaks)

x1<-rep(0,nrow(warpbreaks))
x1[warpbreaks$tension=="M"]<-1
x2<-rep(0,nrow(warpbreaks))
x2[warpbreaks$tension=="H"]<-1
x1
x2
X<-matrix(c(rep(1,nrow(warpbreaks)),x1,x2),ncol=3)
X
y<-warpbreaks$breaks

bhat<-solve(t(X)%*%X)%*%t(X)%*%y
bhat

wbw.lm<-lm(breaks~wool, data=warpbreaks)
summary(wbw.lm)
anova(wbw.lm)

boxplot(breaks~wool, data=warpbreaks)

wb2.lm<-lm(breaks~wool+tension, data=warpbreaks)
summary(wb2.lm)
anova(wb2.lm)

wb2i.lm<-lm(breaks~wool*tension, data=warpbreaks)
summary(wb2i.lm)
anova(wb2i.lm)

interaction.plot(warpbreaks$tension, warpbreaks$wool, warpbreaks$breaks)
interaction.plot(warpbreaks$wool, warpbreaks$tension, warpbreaks$breaks)

### Example 5

setwd("~/Class/MATH5910/Fall2024/R")

anova1<-read.table("anova1.txt", header=T)

summary(anova1)
table(anova1$A)
table(anova1$B)
table(anova1$C)

anova(lm(Resp~A*B*C, data=anova1))
anova(lm(Resp~A*B*C-A:B:C, data=anova1))
anova(lm(Resp~A*B*C-A:B:C-B:C, data=anova1))
anova(lm(Resp~A*B*C-A:B:C-B:C-A:B, data=anova1))

anova(lm(Resp~A+B+C, data=anova1))
anova(lm(Resp~A+B+C+A:B, data=anova1))
anova(lm(Resp~A+B+C+B:C, data=anova1))

anova(lm(Resp~A+B+C+A:C, data=anova1))

boxplot(Resp~A, data=anova1)
boxplot(Resp~B, data=anova1)
boxplot(Resp~C, data=anova1)

interaction.plot(anova1$A, anova1$C, anova1$Resp)