bootzh <-function(k,P1,P2,boots) {#gives bootstrap results
r=2  #number of doses
c=4 #number of categories (same for all variables)
n1=dim(P1)[2]; n2=dim(P2)[2]
A=kronecker(diag(k),matrix(c(1,2,3,4),1,c)) #matrix of scores

y=c(compy(t(P1)),compy(t(P2)))
y=matrix(y,c*k,r,byrow=FALSE)    # holds the counts for each category for each AE, colums are different dose levels
p=y/c(n1,n2)
s.obs=A%*%(p[,2]-p[,1])
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))          # for statistic W_0'
Sigmainv = diag(Sigma)^(-0.5)
index=(Sigmainv==Inf)
Sigmainv[index==1]=1
T2.obs=t(Sigmainv*s.obs)

T2=index=matrix(-99,boots,k)
#set.seed(123)
for (boot in 1:boots) {
    P1.boot=P1[,sample(1:n1,replace=TRUE)]
    P2.boot=P2[,sample(1:n2,replace=TRUE)]
    y=c(compy(t(P1.boot)),compy(t(P2.boot)))
    y=matrix(y,c*k,r,byrow=FALSE)    # holds the counts for each category for each AE, colums are different dose levels
    p=y/c(n1,n2)
    s=A%*%(p[,2]-p[,1])-s.obs
    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))          # for statistic W_0'
    Sigmainv = diag(Sigma)^(-0.5)
    index[boot,]=(Sigmainv==Inf)
    Sigmainv[index[boot,]==1]=1
    T2[boot,]=t(Sigmainv*s)
    }
meanT2=apply(T2,1,mean,na.rm=TRUE)
sig.meanT=mean(meanT2>mean(T2.obs))+0.5*mean(meanT2==mean(T2.obs))
maxT2=apply(T2,1,max,na.rm=TRUE)
sig.maxT=mean(maxT2>max(T2.obs))
#bootpvalue=matrix(NA,k,1)
#for (i in 1:k)
#    {
#    bootpvalue[i]=mean(T[,i]>T2.obs[i], na.rm=TRUE)+0.5*mean(T[,i]==T2.obs[i], na.rm=TRUE)
#    }
#cbind(bootpvalue,colMeans(index))
adjpvalue=matrix(99,k,1)
ind=order(T2.obs,decreasing=TRUE)
j=1
for (i in ind) {
    adjpvalue[i]=mean(maxT2>T2.obs[i])+0.5*mean(maxT2==T2.obs[i])  
    if (j<(k-1)) maxT2=apply(T2[, -ind[1:j]],1,max,na.rm=TRUE) else maxT2=T2[, -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))
}