## Epitaxial layer growth data - location/dispersion modelling
rm(list = ls())  # clear the memory
id1 = c(rep(1,8),rep(7,8))
thick = read.fwf("http://www.stat.ualberta.ca/~wiens/stat568/datasets/LayerGrowth_table_11_2.dat",widths = id1, skip = 1) 
thick

A = as.factor(thick[1:16,1])
B = as.factor(thick[1:16,2])
C = as.factor(thick[1:16,3])
D = as.factor(thick[1:16,4])
E = as.factor(thick[1:16,5])
F = as.factor(thick[1:16,6])
G = as.factor(thick[1:16,7])
H = as.factor(thick[1:16,8])

ymat = as.matrix(thick[,9:16])
y = c(ymat)

data = data.frame(y, A, B, C, D, E, F, G, H)
ybar = apply(ymat, 1, mean)
s2 = apply(ymat, 1, var)
logs2 = log(s2)


############ Location model:

par(mfrow=c(2,1))

xA =  2*as.numeric(A=="+")-1
xB =  2*as.numeric(B=="+")-1
xC =  2*as.numeric(C=="+")-1
xD =  2*as.numeric(D=="+")-1
xE =  2*as.numeric(E=="+")-1
xF =  2*as.numeric(F=="+")-1
xG =  2*as.numeric(G=="+")-1
xH =  2*as.numeric(H=="+")-1

xA = rep(xA,8)
xB = rep(xB,8)
xC = rep(xC,8)
xD = rep(xD,8)
xE = rep(xE,8)
xF = rep(xF,8)
xG = rep(xG,8)
xH = rep(xH,8)

fit1 = lm(y ~(xA + xB + xC + xD + xE + xF + xG + xH)^2)
summary(fit1)
# Fits only the main effects and 2 fi's with A
# same as aov(y ~(xA + xB + xC + xD + xE + xF + xG + xH)^3)
coefs = na.omit(fit1$coef)
effects.location = 2*coefs[-1]

len = length(effects.location)
quants = qnorm(.5*(1 + 1:len/(1+len)))
plot(quants, sort(abs(effects.location)), xlab = "HalfNormal quantiles", ylab = "location effects", xlim = c(0,2), ylim=c(0,1))

# Which are the significant ones?
qwe = sort(abs(effects.location))
qwe.location = qwe[-c(1:(len-3))]
qwe.location
#        xC        xH        xD 
# 0.1142672 0.1734266 0.8039047 
text(x=quants[(len-2):len], y = c(qwe.location[c(1:(len-3))]), labels = c("C", "H", "D"), pos = 4)

# Fit xD only:
fit.location = lm(y ~ xD)
summary(fit.location)


############ Dispersion model:

xA =  2*as.numeric(A=="+")-1
xB =  2*as.numeric(B=="+")-1
xC =  2*as.numeric(C=="+")-1
xD =  2*as.numeric(D=="+")-1
xE =  2*as.numeric(E=="+")-1
xF =  2*as.numeric(F=="+")-1
xG =  2*as.numeric(G=="+")-1
xH =  2*as.numeric(H=="+")-1

fit2 = lm(logs2 ~(xA + xB + xC + xD + xE + xF + xG + xH)^2)
summary(fit2)
# Fits only the main effects and 2 fi's with A
# same as aov(logs2 ~(xA + xB + xC + xD + xE + xF + xG + xH)^3)
coefs = na.omit(fit2$coef)
effects.dispersion = 2*coefs[-1]

len = length(effects.dispersion)
quants = qnorm(.5*(1 + 1:len/(1+len)))
plot(quants, sort(abs(effects.dispersion)), xlab = "HalfNormal quantiles", ylab = "dispersion effects", xlim = c(0,2), ylim = c(0, 2.5))


# Which are the significant ones?
qwe = sort(abs(effects.dispersion))
qwe.dispersion = qwe[-c(1:(len-3))]
qwe.dispersion
#       xD       xA       xH 
# 0.848132 1.233918 1.958919 
text(x=quants[(len-2):len], y = c(qwe.dispersion[c(1:(len-3))]), labels = c("D", "A", "H"), pos = 4)

# Fit xA and xH only:
fit.dispersion = lm(logs2 ~ xA + xH)
summary(fit.dispersion)


#######################################################################################
#######################################################################################
## Epitaxial layer growth data - response modelling
rm(list = ls())  # clear the memory
id1 = c(rep(1,8),rep(7,8))
thick = read.fwf("http://www.stat.ualberta.ca/~wiens/stat568/datasets/LayerGrowth_table_11_2.dat",widths = id1, skip = 1) 
thick

A = as.factor(thick[1:16,1])
B = as.factor(thick[1:16,2])
C = as.factor(thick[1:16,3])
D = as.factor(thick[1:16,4])
E = as.factor(thick[1:16,5])
F = as.factor(thick[1:16,6])
G = as.factor(thick[1:16,7])
H = as.factor(thick[1:16,8])


xA =  2*as.numeric(A=="+")-1
xB =  2*as.numeric(B=="+")-1
xC =  2*as.numeric(C=="+")-1
xD =  2*as.numeric(D=="+")-1
xE =  2*as.numeric(E=="+")-1
xF =  2*as.numeric(F=="+")-1
xG =  2*as.numeric(G=="+")-1
xH =  2*as.numeric(H=="+")-1

## Incorporate the values of L and M

dataL0M1 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(-1,16), rep(1,16), thick[,9])
dataL0M2 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(-1,16), rep(2,16), thick[,10])
dataL0M3 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(-1,16), rep(3,16), thick[,11])
dataL0M4 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(-1,16), rep(4,16), thick[,12])
dataL1M1 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(1,16), rep(1,16), thick[,13])
dataL1M2 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(1,16), rep(2,16), thick[,14])
dataL1M3 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(1,16), rep(3,16), thick[,15])
dataL1M4 = cbind(xA,xB,xC,xD,xE,xF,xG,xH,rep(1,16), rep(4,16), thick[,16])

data = rbind(dataL0M1,dataL0M2,dataL0M3,dataL0M4,dataL1M1,dataL1M2,dataL1M3,dataL1M4)
data  # All of the data

# Form the three orthogonal contrasts for M

xMl = (data[,10]==1 | data[,10]==2) - (data[,10]==3 | data[,10]==4)
xMq = (data[,10]==1 | data[,10]==4) - (data[,10]==2 | data[,10]==3)
xMc = (data[,10]==1 | data[,10]==3) - (data[,10]==2 | data[,10]==4)

# Re-label, using all 16 x 8 runs

xA = data[,1]
xB = data[,2]
xC = data[,3]
xD = data[,4]
xE = data[,5]
xF = data[,6]
xG = data[,7]
xH = data[,8]
xL = data[,9]
y = data[,11]

fit1 = lm(y ~(xA + xB + xC + xD + xE + xF + xG + xH + xL + xMl + xMq + xMc)^3)

# same as aov(y ~(xA + xB + xC + xD + xE + xF + xG + xH + xL + xMl + xMq + xMc)^3)

summary(fit1)
coefs = na.omit(fit1$coef)
effects = 2*coefs[-1]


par(mfrow=c(1,1))
len = length(effects)
quants = qnorm(.5*(1 + 1:len/(1+len)))
plot(quants, sort(abs(effects)), xlab = "HalfNormal quantiles", ylab = "effects")

# Which are the significant ones?
qwe = sort(abs(effects))
qwe = qwe[-c(1:(len-7))]
qwe
# xA:xH:xMq    xC:xMl        xH       xMl     xH:xL        xL        xD 
# 0.1633234 0.1660766 0.1734266 0.1803953 0.4776078 0.6591359 0.8039047 

# Fit the significant effects:

fit2 = lm(y ~ xD + xH + xL + xMl + xH:xL + xC:xMl + xA:xH:xMq)
summary(fit2)

windows()
par(mfrow=c(2,2))
interaction.plot(xL,xH,y)
xM = data[,10]
interaction.plot(xM,xC,y)
interaction.plot(xM[xH==1],xA[xH==1],y[xH==1])
interaction.plot(xM[xH==-1],xA[xH==-1],y[xH==-1])


########################################################################################
########################################################################################