permzh <-function(k,P1,P2,perms) {#permutation distribution
r=2  #number of doses
c=4 #number of categories (same for all variables)
P=rbind(P1,P2)
Pperms=cbind(P[1:k,],P[(k+1):(k*r),]) #this matrix has doses across rows
W0=matrix(-99,perms,k)
A=kronecker(diag(k),matrix(c(1,2,3,4),1,c))
for (perm in 1:perms)
    {
    #P=rbind(Pperms[,v[allperms[[sim]]]],Pperms[,v[-allperms[[sim]]]]) #permuted matrix (if sim=1==original matrix)
    y=c(compy(t(P[(1:k),])),compy(t(P[(k+1):(2*k),])))
    y=matrix(y,c*k,r,byrow=FALSE)    # holds the counts for each category for each AE, colums are different dose levels
    n1=dim(P1)[2]; n2=dim(P2)[2]
    p=y/c(n1,n2)
    s=A%*%(p[,2]-p[,1])

    if (perm==1)  #computes covariance matrix, the same for all perms since pooled
        {
        ppooled=rowSums(y)/(n1+n2)
        CovAE=vector("list",k)
        for (j in 1:k)
            {
            CovAE[[j]]=(1/n1+1/n2) * ( diag(ppooled[(c*j-(c-1)):(c*j)]) - ppooled[(c*j-(c-1)):(c*j)]%*%t(ppooled[(c*j-(c-1)):(c*j)]) )
            }
        Covd=bdiag(CovAE) # ignore correlation among endpoints
        Sigma = (A  %*% Covd %*%t(A))          # statistic W_0'
        Sigmainv = diag(Sigma)^(-0.5)
        index=(Sigmainv==Inf)
        Sigmainv[index==1]=1
        }
    W0[perm,]=t(Sigmainv*s)  # differences in weigthed mean scores
    P=Pperms[,sample(1:(n1+n2),replace=FALSE)]
    P=rbind(P[,1:n1],P[,(n1+1):(n1+n2)])
    }  #end permutations
meanW0=apply(W0,1,mean)
sig.meanT=mean(meanW0>meanW0[1])+0.5*mean(meanW0==meanW0[1]) 
maxs=apply(W0,1,max)
sig.maxT=mean(maxs>max(W0[1,]))
#permpvalue=matrix(NA,k,1)
#for (i in 1:k)
#    {
#    permpvalue[i]=mean(W0[,i]>W0[1,i])+0.5*mean(W0[,i]==W0[1,i])
#    }
#permpvalue
adjpvalue=matrix(99,k,1)
ind=order(W0[1,],decreasing=TRUE)
j=1
for (i in ind) {
    adjpvalue[i]=mean(maxs>W0[1,i])+0.5*mean(maxs==W0[1,i])  
    if (j<(k-1)) maxs=apply(W0[, -ind[1:j]],1,max) else maxs=W0[, -ind[1:j]]
    if (j>1 && adjpvalue[ind[j]]<adjpvalue[ind[j-1]]) adjpvalue[ind[j]]=adjpvalue[ind[j-1]]
    j=j+1
    }

return(list(global.sig.meanT=sig.meanT,global.sig.maxT=sig.maxT,adjp=adjpvalue))
}