#APTS April 2016 - helicoper experiment for the design lecture library(FrF2) library(gplots) #some useful functions randomise <- function(n) { r<-sample(n) plot(c(0,n+1),c(0,1),type="n",xlab="",ylab="",xaxt="n",yaxt="n",bty="n") for (i in 1:n) text(i,0.5,r[i],cex=3) } contr.twolevel<-function() { contr<-c(-1,1) cbind(contr,deparse.level=0) } Amatrix<-function(X1,X2) { X1<-as.matrix(X1) X2<-as.matrix(X2) A<-solve(t(X1)%*%X1)%*%t(X1)%*%X2 A } #run the experiment and collect the data randomise(16) setwd("helicopters") times <- read.table("helicopter16.csv",sep=",",header=T,row.names=1) times$W <- factor(times$W) times$B <- factor(times$B) times$M <- factor(times$M) times$C <- factor(times$C) #code factor levels -1,1 contrasts(times$W)<-contr.twolevel() contrasts(times$B)<-contr.twolevel() contrasts(times$M)<-contr.twolevel() contrasts(times$C)<-contr.twolevel() #AoV and normal plot timesFULL.aov <- aov(Y~W*B*M*C,data=times) coeff.FULL<-coef(timesFULL.aov) X.FULL<-model.matrix(timesFULL.aov) XtX.FULL<-t(X.FULL)%*%X.FULL #par(mfrow=c(1,1)) windows(title="QQplot Full") #quartz(title="QQplot Full") qqnorm.aov(timesFULL.aov,label=T) #main effect and interaction plots windows(title="Main effects Full") #quartz(title="Main effects Full") MEPlot(timesFULL.aov,main="Main effects") windows(title="Interactions Full") #quartz(title="Interactions Full") IAPlot(timesFULL.aov,select=c(1,2,3,4),show.alias=T,main="Interactions") #Half fraction 1 times.FF1<-times[X.FULL[,16]==1,] #AoV and normal plot timesFF1.aov <- aov(Y~W*B*M*C,data=times.FF1) coeff.FF1<-coef(timesFF1.aov) X.FF1<-model.matrix(timesFF1.aov) XtX.FF1<-t(X.FF1)%*%X.FF1 #par(mfrow=c(1,1)) windows(title="QQplot FF1") #quartz(title="QQplot FF1") qqnorm.aov(timesFF1.aov,label=T) A.FF1<-Amatrix(X.FF1[,1:8],X.FF1[,9:16]) coeff.FF1A<-coeff.FULL[1:8]+A.FF1%*%coeff.FULL[9:16] #Half fraction 2 times.FF2<-times[X.FULL[,16]==-1,] #AoV and normal plot timesFF2.aov <- aov(Y~W*B*M*C,data=times.FF2) coeff.FF2<-coef(timesFF2.aov) X.FF2<-model.matrix(timesFF2.aov) XtX.FF2<-t(X.FF2)%*%X.FF2 #par(mfrow=c(1,1)) windows(title="QQplot FF2") #quartz(title="QQplot FF2") qqnorm.aov(timesFF2.aov,label=T) A.FF2<-Amatrix(X.FF2[,1:8],X.FF2[,9:16]) coeff.FF2A<-coeff.FULL[1:8]+A.FF2%*%coeff.FULL[9:16] #PB12 design PB12<-read.table("PB12b.txt",sep=" ") #PB122<-matrix(0,nrow=12,ncol=11) #for(i in 1:12) PB122[i,]<-2*as.numeric(PB12[i,])-3 #PB12<-PB122 times.PB<-times[match(data.frame(t(PB12)),data.frame(t(X.FULL[,2:5]))),] #AoV and normal plot timesPB.aov <- aov(Y~W+B+M+C,data=times.PB) coeff.PB<-coef(timesPB.aov) timesPB.aov2 <- aov(Y~W*B*M*C,data=times.PB) X.PB<-model.matrix(timesPB.aov2) XtX.PB<-t(X.PB)%*%X.PB #par(mfrow=c(1,1)) windows(title="QQplot PB") #quartz(title="QQplot PB") qqnorm.aov(timesPB.aov2,label=T) A.PB<-Amatrix(X.PB[,1:5],X.PB[,6:16]) coeff.PBA<-coeff.FULL[1:5]+A.PB%*%coeff.FULL[6:16] round(A.PB,dig=2)