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