R/page.test.unbalanced.R

In MultNonParam: Multivariate Nonparametric Methods

page.test.unbalanced<-function(x,trt,blk,sides=2){
#' Perform Page test for unbalanced two-way design
#' @param x A vector of responses
#' @param trt A vector of consecutive integers starting at 1 indicating treatment
#' @param blk A vector of consecutive integers starting at 1 indicating block
#' @param sides A single integer indicating sides.  Defaults to 2.
#' @return P-value for Page test.
#' @examples
#' page.test.unbalanced(rnorm(15),rep(1:3,5),rep(1:5,rep(3,5)))
#' @export
#' @importFrom stats var
#' @importFrom stats pnorm
#' @importFrom stats rnorm
   trt<-as.numeric(as.factor(trt))
   blk<-as.numeric(as.factor(blk))
   makecmat<-function(trt,blk){
      nmat<-table(trt,blk)
      cmat<-array(NA,rep(length(unique(trt)),2))
      nt<-length(unique(trt))
      mvec<-rep(NA,nt)
      nb<-length(unique(blk))
      nsl<-apply(nmat,2,sum)
      nsk<-apply(nmat,1,sum)
      mvec<-rep(0,nt)
      for(k in seq(nt)){
         for(l in seq(nb)) mvec[k]<-mvec[k]+nmat[k,l]*(nsl[l]+1)
      }
      mvec<-mvec/(2*nsk)
      for(k in seq(nt-1)) for(m in ((k+1):nt)) {
         cmat[k,m]<-0
         for(l in seq(nb))
            cmat[k,m]<-cmat[k,m]-nmat[k,l]*nmat[m,l]*(nsl[l]+1)
         cmat[k,m]<-cmat[k,m]/(12*nsk[k]*nsk[m])
         cmat[m,k]<-cmat[k,m]
      }
      for(k in seq(nt)){
         cmat[k,k]<-0
         for(l in seq(length(unique(blk))))
            cmat[k,k]<-cmat[k,k]+nmat[k,l]*sum(nmat[-k,l])*(nsl[l]+1)
         cmat[k,k]<-cmat[k,k]/(12*nsk[k]^2)
      }
      return(list(mvec=mvec,cmat=cmat))
   }
   pagestatu<-function(x,trt,blk){
      ranksums<-rep(0,length(unique(trt)))
      names(ranksums)<-unique(trt)
      for(oneblk in unique(blk)){
         rr<-rank(x[blk==oneblk])
         thesetrt<-trt[blk==oneblk]
         for(onet in unique(trt)){
            ranksums[onet]<-ranksums[onet]+sum(rr[thesetrt==onet])
         }
      }
      return(ranksums/table(trt))
   }
   fun.checkvar<-function(trt,blk,nsamp=10000){
      ranksums<-rep(0,length(unique(trt)))
      names(ranksums)<-unique(trt)
      rankmeans<-array(NA,c(nsamp,length(unique(trt))))
      for(m in 1:nsamp){
         x<-rnorm(length(trt))
         rankmeans[m,]<-pagestatu(x,trt,blk)
      }
      return(list(variance=var(rankmeans),mean=apply(rankmeans,2,mean)))
   }
# First executable
   rm<-pagestatu(x,trt,blk)
   nt<-length(rm)
   moms<-makecmat(trt,blk)
   z<-seq(nt)%*%(rm-moms$mvec)/sqrt(seq(nt)%*%moms$cmat%*%seq(nt))
   pval<-pnorm(z,lower.tail=F)
   if(sides==2) {
      pval<-min(pval,1-pval)*2
      alts<-"two-sided"
   }else{alts<-"greater"}
   nullval<-0
   names(nullval)<-"median"
   parmout<-c(seq(nt)%*%moms$mvec,seq(nt)%*%moms$cmat%*%seq(nt))
   names(parmout)<-c("Expectation","Variance")
   out<-list(null.value=nullval,alternative=alts,method="Page",
      statistic=seq(nt)%*%rm,z=z,p.value=pval,parameters=parmout,
      rankmeans=rm,rankmeanexpect=moms$mvec,rankmeanvar=moms$cmat)
   class(out)<-"htest"
   return(out)
}
#fun.checkvar(trt,blk)