rm(list = ls()) # clear the memory thick = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/adaptedlayer_table_4_1.dat", header = F, quote= "") A = as.factor(thick[1:16,2]) B = as.factor(thick[1:16,3]) C = as.factor(thick[1:16,4]) D = as.factor(thick[1:16,5]) ymat = as.matrix(thick[,6:11]) y = c(ymat) data = data.frame(y, A, B, C, D) ybar = apply(ymat, 1, mean) s2 = apply(ymat, 1, var) logs2 = log(s2) plot.design(data) windows() par(mfrow=c(3,2)) interaction.plot(B,A,ybar) interaction.plot(C,A,ybar) interaction.plot(D,A,ybar) interaction.plot(C,B,ybar) interaction.plot(D,B,ybar) interaction.plot(D,C,ybar) windows() par(mfrow=c(1,2)) interaction.plot(D,C,ybar, sub = "synergistic") interaction.plot(C,D,ybar, sub = "antagonistic") fit = aov(y ~ (A + B + C + D)^4, data) summary(fit) # Of minor interest: tab = model.tables(fit, type = "effects") #tab xA = 2*as.numeric(data$A=="+")-1 xB = 2*as.numeric(data$B=="+")-1 xC = 2*as.numeric(data$C=="+")-1 xD = 2*as.numeric(data$D=="+")-1 #Example: 2*mean(y*xA) # Main effect of A fit2 = lm(y ~(xA + xB + xC + xD)^4) anova(fit2) #D, CD and possibly B are significant effects = 2*fit2$coef[-1] effects # Main effect of A and others # load the BsMD package library(BsMD) windows() par(mfrow=c(1,2)) DanielPlot(fit2, half = T) LenthPlot(fit2, alpha=.05) # Not clear if trimming is done. # alpha PSE ME SME # 0.0500000 0.0863125 0.2218733 0.4504348 # Simultaneous (1-alpha) CIs are thetahat +/- SME #LenthPlot(fit2, alpha=.01) # Not clear if trimming is done. # alpha PSE ME SME # 0.0100000 0.0863125 0.3480243 0.6466053 windows() par(mfrow = c(1,1)) LenthPlot(fit2, alpha=.05) library(car) leveneTest(y, group = as.factor(rep(1:16,6))) # Pooled estimate of variance: S=sqrt(mean(s2)) S # compare with 'fit' ####################################################### # Nominal-the-best analysis rm(list = ls()) # clear the memory # Compute log(s^2) for each of the 16 combinations of levels (6 replicates each); use the ORIGINAL data thick0 = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/originallayer_table_4_10.dat", header = F, quote= "") A = as.factor(thick0[1:16,1]) B = as.factor(thick0[1:16,2]) C = as.factor(thick0[1:16,3]) D = as.factor(thick0[1:16,4]) ymat = as.matrix(thick0[,5:10]) y = c(ymat) data0 = data.frame(y, A, B, C, D) logs2 = log(apply(ymat,1,var)) # View these as the data in one replicate of a 2^4 factorial xA = 2*as.numeric(data0$A=="+")-1 xB = 2*as.numeric(data0$B=="+")-1 xC = 2*as.numeric(data0$C=="+")-1 xD = 2*as.numeric(data0$D=="+")-1 ################### # Response is Y fit3 = lm(y ~(xA + xB + xC + xD)^4) summary(fit3) #effects = round(as.matrix(2*fit3$coef),3) #effects = as.matrix(effects[-1,]) #effects library(BsMD) windows() par(mfrow=c(1,2)) DanielPlot(fit3, half = T, sub = "location") fit3a = lm(y ~ xD) ###################### # Response is log(S^2) library(car) leveneTest(y, group = as.factor(rep(1:16,6))) # differences in variances highly significant xA = xA[1:16] xB = xB[1:16] xC = xC[1:16] xD = xD[1:16] fit4 = lm(logs2 ~(xA + xB + xC + xD)^4) summary(fit4) #vareffects = round(as.matrix(2*fit4$coef),3) #vareffects = as.matrix(vareffects[-1,]) #vareffects DanielPlot(fit4, half = T, sub = "dispersion") fit4a = lm(logs2 ~ xA) windows() par(mfrow=c(1,2)) LenthPlot(fit3, alpha=.05, sub = "location") LenthPlot(fit4, alpha=.05, sub = "dispersion") fit3a fit4a