multcomp<-function(a,cont.matrix="allpairs",type="scheffe",conf.level=.95) {

fact=attr(a$terms,"term.labels")[1]
levels=a$xlevels[[1]]
nlevels=length(levels)

#get sample sizes
n=numeric()
for(i in 1:nlevels) n[i]=length(which(a$model[fact]==levels[i]))

#get level means
#levelmeans=numeric()
#for(i in 1:nlevels) levelmeans[i]=mean(a$model[[1]],[fact])
levelmeans=tapply(a$model[[1]],a$model[fact],mean)


if (class(cont.matrix)=="character" && cont.matrix=="allpairs") {
   ncompare=choose(nlevels,2)
   cont.matrix=matrix(0,nlevels,ncompare)
   currcol=1
   for(i in 1:(nlevels-1)) { 
      for(j in (i+1):nlevels) {
         cont.matrix[i,currcol]=1
         cont.matrix[j,currcol]=-1
         currcol=currcol+1
      }
   }
}

#set up strings to describe contrasts for easier reading
cont.desc=character()
for (i in 1:ncol(cont.matrix)) {
   cont.desc[i]=""
   for (j in 1:nlevels) {
      if (cont.matrix[j,i] != 0) {
         if (cont.matrix[j,i]>0) st="+"
         else st="-"
         if (abs(cont.matrix[j,i])==1)
            cont.desc[i]=paste(cont.desc[i],st,sep="")
         else cont.desc[i]=paste(cont.desc[i],st,"(",round(abs(cont.matrix[j,i]),3),")",sep="")
         cont.desc[i]=paste(cont.desc[i],levels[j],sep="")
      }
   }
}

mse=sum(a$residuals^2/a$df.residual)
if (type=="tukey") coeff=qtukey(conf.level,nlevels,a$df.residual)/sqrt(2)
else if (type=="bonf") coeff=qt(1-(1-conf.level)/(2*ncol(cm)),a$df.residual)
else coeff=sqrt((nlevels-1)*qf(conf.level,nlevels-1,a$df.residual))

lhs= t(cont.matrix)%*%levelmeans
rhs=sqrt(coeff*mse*colSums(cont.matrix^2/n))

lb=lhs-rhs
ub=lhs+rhs
data.frame(contrast=cont.desc,estimate=lhs,lb=lb,ub=ub)

}