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]]