pdf("rpd4_14r.pdf")

XREG <- matrix(c(rep(1,20),rep(0.02,5),rep(0.07,5),rep(0.11,5),rep(0.15,5)),byrow=F,ncol=2)
XAOV <- matrix(c(rep(c(1,0,0,0),5),rep(c(0,1,0,0),5),rep(c(0,0,1,0),5),rep(c(0,0,0,1),5)),
    byrow=T,ncol=4)

Y <- matrix(c(238.3,270.7,210.0,248.7,242.4,235.1,228.9,236.2,255.0,228.9,
              236.2,208.0,243.5,233.0,233.0,178.7,186.0,206.9,215.3,219.5),ncol=1)

BETAREG <- solve(t(XREG)%*%XREG) %*% t(XREG) %*% Y
BETAAOV <- solve(t(XAOV)%*%XAOV) %*% t(XAOV) %*% Y

PREG <- XREG %*% solve(t(XREG)%*%XREG) %*% t(XREG)
PAOV <- XAOV %*% solve(t(XAOV)%*%XAOV) %*% t(XAOV)

ix <- diag(1,ncol=nrow(Y),nrow=nrow(Y))
jx <- matrix(rep(1,nrow(Y)*nrow(Y)),ncol=nrow(Y))

SSTOTC <- t(Y) %*% (ix - (1/nrow(Y))*jx) %*% Y
DFTOTC <- nrow(Y)-1

SSMODREG <- t(Y) %*% (PREG - (1/nrow(Y))*jx) %*% Y
DFMODREG <- ncol(XREG)-1
MSMODREG <- SSMODREG/DFMODREG

SSMODAOV <- t(Y) %*% (PAOV - (1/nrow(Y))*jx) %*% Y
DFMODAOV <- ncol(XAOV)-1
MSMODAOV <- SSMODAOV/DFMODAOV

YHATREG <- XREG %*% BETAREG
YHATAOV <- XAOV %*% BETAAOV

SSRESREG <- t(Y) %*% (ix - PREG) %*% Y
DFRESREG <- nrow(Y)-ncol(XREG)
MSRESREG <- SSRESREG/DFRESREG

SSPE <- t(Y-YHATAOV) %*% (Y-YHATAOV)
DFPE <- nrow(Y)-ncol(XAOV)
MSPE <- SSPE/DFPE

SSLF <- t(YHATREG-YHATAOV) %*% (YHATREG-YHATAOV)
DFLF <- ncol(XAOV)-ncol(XREG)
MSLF <- SSLF/DFLF

print(cbind(c("Total(Corr)"),DFTOTC,SSTOTC))
print(cbind(c("Treatments"),DFMODAOV,SSMODAOV,MSMODAOV))
print(cbind(c("Regression"),DFMODREG,SSMODREG,MSMODREG))
print(cbind(c("Lack of Fit"),DFLF,SSLF,MSLF))
print(cbind(c("Pure Error"),DFPE,SSPE,MSPE))

F_LOF <- MSLF/MSPE
FCRIT_LOF <- qf(.95,DFLF,DFPE)
PROBF_LOF <- 1-pf(F_LOF,DFLF,DFPE)

F_LOF
FCRIT_LOF
PROBF_LOF

F_B1A <- MSMODREG/MSPE
FCRIT_B1A <- qf(.95,DFMODREG,DFPE)
PROBF_B1A <- 1-pf(F_B1A,DFMODREG,DFPE)

F_B1A
FCRIT_B1A
PROBF_B1A

F_B1B <- MSMODREG/MSRESREG
FCRIT_B1B <- qf(.95,DFMODREG,DFRESREG)
PROBF_B1B <- 1-pf(F_B1B,DFMODREG,DFRESREG)

F_B1B
FCRIT_B1B
PROBF_B1B

dev.off()