setwd("C:/Courses/Computer-Presentation")
rm(list=ls(all=TRUE))  # CLEAN UP THE CURRENT WORK DIRECTORY

options()$contrasts
# options(contrasts=c("contr.sum", "contr.poly"))




########################
# ONE-WAY ANOVA
########################

# THE DATA
productivity <- c(7.6, 8.2, 6.8, 5.8, 6.9, 6.6, 6.3, 7.7, 6.0, 
    6.7, 8.1, 9.4, 8.6, 7.8, 7.7, 8.9, 7.9, 8.5, 8.7, 7.1, 8.4, 
    8.5, 9.7, 10.1, 7.8, 9.6, 9.5)
RnD <- rep(c("low", "moderate", "high"), c(9, 12, 6))
dat <- data.frame(productivity=productivity, RnD=RnD)
dat 

# PLOT OF THE DATA
RnD0 <- rep(c(1,2,3), c(9, 12, 6))
par(bg="grey", mar=c(6, 8, 4, 8))
plot(x=RnD0, y=productivity,  type="p", yaxt="s", xaxt="n", ann=T, main="Plot of One-Way ANOVA Data",
    xlab="Expenditure RnD", ylab="prod improvement")
axis(1, at=c(1,2,3), labels = c("low", "moderate", "high"), tick =F)

# REGRESSION WITH DUMMY VARIABLES
fit <- lm(productivity~factor(RnD), x=T, data=dat)
summary(fit)
anova(fit)
fit$x  # THE DESIGN MATRIX
names(fit)



########################
# TWO-WAY ANOVA 
########################

growth <- c(1.4, 2.4, 2.2, 2.1, 1.7, .7, 1.1, 2.4, 1.8, 2.0, 2.5, .9, 1.3, .5)
gender <- rep(c("male", "female"), c(7, 7))
bone <- c(rep(c("severe", "moderate", "mild"), c(3, 2, 2)), 
    rep(c("severe", "moderate", "mild"), c(3, 3, 1)))
dat <- data.frame(growth=growth, gender=gender, bone=bone)
dat 

# INTERACTION PLOT
par(bg="white", mfrow=c(2, 1), mar=c(2, 2, 2, 2))
interaction.plot(gender, bone, growth, fixed=TRUE, col = 2:3, leg.bty = "o")
interaction.plot( bone, gender,growth, fixed=TRUE, col = 2:3, leg.bty = "o")

# REGRESSION WITH DUMMY VARIABLES
# ---------------------------------

# SEVERAL RELATED MODELS
fit1 <- lm(growth~factor(gender) + factor(bone) + factor(gender):factor(bone), x=T, data=dat)
fit1 <- lm(growth~factor(gender)*factor(bone), x=T, data=dat)
summary(fit1); anova(fit1)
fit1$x  # THE DESIGN MATRIX

fit2 <- lm(growth~factor(gender) + factor(bone), x=T, data=dat)
fit3 <- lm(growth~factor(gender), x=T, data=dat)
fit4 <- lm(growth~factor(bone), x=T, data=dat)

# TEST OF INTERACTION
anova(fit1, fit2)
# TEST OF MAIN EFFECT OF GENDER
anova(fit2, fit4)
# TEST OF MAIN EFFECT OF BONE
anova(fit2, fit3)




########################
# ANCOVA 
########################

yield <- scan()
144 140 171
98 99 122 
123 139 153 
147 113 115 
193 171 146 
136 113 145 
129 125 129 
141 109 134 
125 134 143 
136 142 133 

nitrogen <- scan()
42 54 77
19 33 57
31 56 74
44 42 53
64 69 67
40 43 68
34 52 60
44 39 61
31 52 67
36 59 57

regulator <- rep(c("A", "B", "C"), 10)
dat <- data.frame(yield, nitrogen, regulator)
dat


# PLOT OF ANCOVA Data
par(mfrow=c(1, 1), oma=rep(4, 4))
plot(range(nitrogen), range(yield), type="n", xlab="nitrogen", ylab="yield",
    main="Plot of Sugar Beets Data")
level <- sort(unique(regulator))
for (i in 1:length(level)){
    temp <- dat[dat$regulator==level[i],]
    points(temp$nitrogen, temp$yield, pch=level[i], col=i)
}

# REGRESSION WITH DUMMY VARIABLES
# --------------------------------
fit1 <- lm(yield ~ nitrogen + factor(regulator) + nitrogen:factor(regulator), data=dat, x=T)  # INTERACTION MODEL
fit1 <- lm(yield ~ nitrogen*factor(regulator), data=dat, x=T)
summary(fit1); anova(fit1)
fit1$x

fit2 <- lm(yield ~ nitrogen + factor(regulator), data=dat, x=T)  # ANCOVA MODEL
summary(fit2); anova(fit2)

fit3 <- lm(yield ~ factor(regulator), data=dat, x=T) # ONE-WAY ANOVA MODEL
summary(fit3); anova(fit3)

fit4 <- lm(yield ~ nitrogen, data=dat, x=T) 
summary(fit4); anova(fit4)

# ANOVA 
anova(fit3)
# ANCOVA
anova(fit2, fit4)
# CHECK THE ASSUMPTION OF NO INTERACTION
anova(fit1, fit2)

# COINCIDENCE
anova(fit1, fit4)
#