yAR <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,20,21,22)  
yARout <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,20,210,22)  
levelARa <- c(rep(1,9),rep(2,9))
levelARb <- c(1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3)


twowayhyp2  <- function(p,q)
     {
          amat <- t(contr.sum(p)) %x% t(rep(1,q))
          bmat <- t(rep(1,p)) %x% t(contr.sum(q))
          abmat <- t(contr.sum(p)) %x% t(contr.sum(q)) 
          list(amat=amat,bmat=bmat,abmat=abmat)
     }

 # read in the ww.r file

source("http://www.stat.wmich.edu/mckean/HMC/Rcode/AppendixB/ww.r")

library(quantreg)


twowayfrt2 = function(y,a,b,inda,indb,type="Fisher",alpha=.05)
  {
# This returns the Type III Wilcoxon tests for a two-way
# analysis. It also returns the associated main effect MCPs
# (Fisher’s or bonferroni). These are easy to interpret if
# interaction is not present.
# twowayfrt2 uses twowayhyp2
#
cell = (inda - 1)*b + indb
wmat = cellmnxy(cell)
r1 = c(1,rep(0,(a*b - 1)))
c1 = rep(1,(a*b - 1))
c2 = diag(c1)
et = cbind(c1,c2)
emat = rbind(r1,et)
xmat = wmat%*%emat
x1 = xmat[,2:(a*b)]
tempw = wwest(x1,y,"WIL",print.tbl=F)
beta = tempw$tmp1$coefficients
tauhat = tempw$tmp2$tau
temph = twowayhyp2(a,b)
ahmat0 = temph$amat%*%emat
ahmat = ahmat0[,2:(a*b)]
bhmat0 = temph$bmat%*%emat
bhmat = bhmat0[,2:(a*b)]
abhmat0 = temph$abmat%*%emat
abhmat = abhmat0[,2:(a*b)]
print("Wilcoxon Test on Main Effects (Type 3) for Factor A")
dropa = droptest(x1,y,ahmat,delta=.80,param=2,print.tbl=T)
print("Wilcoxon Test on Main Effects (Type 3) for Factor B")
dropb = droptest(x1,y,bhmat,delta=.80,param=2,print.tbl=T)
print("Wilcoxon Test on Interaction Between Factors A and B")
dropab = droptest(x1,y,abhmat,delta=.80,param=2,print.tbl=T)
muhat = emat%*%beta
n = length(y)
cellss = t(t(rep(1,n))%*%wmat)
cellhat = matrix(rep(0,(a*b)),ncol=b)
secell = matrix(rep(0,(a*b)),ncol=b)
ic = 1
for(i in 1:a){
  for(j in 1:b){
    cellhat[i,j] = muhat[ic]
    secell[i,j] = 1/cellss[ic]
    ic = ic + 1
  }
  }
colmean = t(t(rep(1,a))%*%cellhat)/a
colse = t(t(rep(1,a))%*%secell)/a^2
rowmean = t(t(rep(1,b))%*%t(cellhat))/b
rowse = t(t(rep(1,b))%*%t(secell))/b^2
df = n - a*b
if(type=="Fisher"){crit = qt(1-(alpha/2),df)}
numba = choose(a,2)
numb = choose(b,2)
numall = numba+numb
if(type=="Bonferroni"){crit = qt(1-(alpha/(2*numall)),df)}
am1 = a - 1
ic = 1
firstlev = rep(0,am1)
secondlev = rep(0,am1)
estdiff = rep(0,numba)
sediff = rep(0,numba)
lcidiff = rep(0,numba)
ucidiff = rep(0,numba)
for(i in 1:am1){
  ip1 = i + 1
  for(j in ip1:a){
    firstlev[ic] = i
    secondlev[ic] = j
    estdiff[ic] = rowmean[i] - rowmean[j]
    sediff[ic] = crit*tauhat*sqrt(rowse[i] + rowse[j])
    lcidiff[ic] = estdiff[ic] - sediff[ic]
    ucidiff[ic] = estdiff[ic] + sediff[ic]
    ic = ic + 1
   }
 }
rowmcp = cbind(firstlev,secondlev,estdiff,sediff,lcidiff,ucidiff)
bm1 = b - 1
ic = 1
firstlev = rep(0,bm1)
secondlev = rep(0,bm1)
estdiff = rep(0,numb)
sediff = rep(0,numb)
lcidiff = rep(0,numb)
ucidiff = rep(0,numb)
for(i in 1:bm1){
  ip1 = i + 1
  for(j in ip1:b){
    firstlev[ic] = i
    secondlev[ic] = j
    estdiff[ic] = colmean[i] - colmean[j]
    sediff[ic] = crit*tauhat*sqrt(colse[i] + colse[j])
    lcidiff[ic] = estdiff[ic] - sediff[ic]
    ucidiff[ic] = estdiff[ic] + sediff[ic]
    ic = ic + 1
   }
 }
colmcp = cbind(firstlev,secondlev,estdiff,sediff,lcidiff,ucidiff)
print(rowmcp)
print(colmcp)
# print(dim(x1))
# print(dim(tempw$ans))
# print(tempw$ans)
# print(tempw$ans[1:12])
# print(xmat)
pred1 <- xmat %*% tempw$ans[1:nrow(tempw$ans)]
par(mfrow=c(2,2))
boxplot(y ~ as.factor(inda)*as.factor(indb),varwidth=TRUE,col=1:a)
plot(pred1,tempw$tmp1$residuals)
qqnorm(tempw$tmp1$residuals)
qqline(tempw$tmp1$residuals)
hist(tempw$tmp1$residuals,freq=FALSE)
lines(density(tempw$tmp1$residuals))
par(mfrow=c(1,1))
# res1  <- y - pred1
# windows()
# plot(res1,tempw$tmp1$residuals)
list(dropa=dropa,dropb=dropb,dropab=dropab,rowmcp=rowmcp,colmcp=colmcp )
}

 # ----------------------------------------------------------------------------------

 # view the data

boxplot(yAR ~ as.factor(levelARa)*as.factor(levelARb))
boxplot(yARout ~ as.factor(levelARa)*as.factor(levelARb))

 # Note that the lm - based ANOVA requires the factor levels as factors

anova(lm(yAR ~ as.factor(levelARa)*as.factor(levelARb)))
anova(lm(yARout ~ as.factor(levelARa)*as.factor(levelARb)))

 # the RGLM code, on the other hand, requires the factor levels as numeric objects

levelARa <- as.numeric(levelARa)
levelARb <- as.numeric(levelARb)

 # the twowayfrt2 function is adapted from a function by J. Terpstra
 # the arguments are:
 # (response vector, number of levels of factor a, number of levels of factor b,
 #     the level of factor a, the level of factor b, 
 #     type of multiple comparison method, alpha level)

yAR.RGLMfit <- twowayfrt2(yAR,2,3,levelARa,levelARb,type="Fisher",alpha=.05)

yARout.RGLMfit <- twowayfrt2(yARout,2,3,levelARa,levelARb,type="Fisher",alpha=.05)