# R/get.perm.mat.R In MRPP: Multiresponse permutation procedure and its variable importance and variable selection methods.

```if (FALSE) {  ## old code
get.perm2.mat <-
function(trt, B=100L) ## trt needs have exactly 2 levels; other situations are not implemented yet
{
if(!exists('.Random.seed', envir=globalenv())) runif(1L)
save.seed=get('.Random.seed', envir=globalenv())
n=table(trt)
if(length(n) != 2L) stop('"get.perm2.mat" only works when "trt" has exactly two levels.')
N=sum(n); m=min(n)
#perms=replicate(B, sort(sample(N, m)))
perms=combn2R(N, m, R=B, sample.method="noReplace")
idx1=sort(which(as.factor(trt)==names(n[which.min(n)])))
colid=which(apply(perms == idx1 , 2L, all))
if (length(colid)==0L) {
perms[,1L]=idx1
}else if(colid!=1L){
perms[, colid] = perms[,1L]
perms[,1L] = idx1
}
attr(perms,'.Random.seed')=save.seed
perms
}

permuteTrt <-
function(trt, B=100L) ## permutation matrix for one way design
## returns a length(trt) by B matrix, if B is not larger than multinomial coefficient.
{
if(!exists('.Random.seed', envir=globalenv())) runif(1)
save.seed=get('.Random.seed', envir=globalenv())
n=table(trt)
N=sum(n);
#    multcoeff=exp(lfactorial(N)-sum(lfactorial(n)))
#    if(multcoeff>=B){   ## list all effective permutations
#        B=multcoeff
#        get.combn=function(vec, n){
#            if(length(vec)==n[1]) return(matrix(vec))
#            this.ans=combn(length(vec), n[1])
#            ans=apply(this.ans, 2, function(x){
#                    tmp=get.combn(vec[-x],n[-1])
#                    rbind(matrix(rep(vec[x],ncol(tmp)),n[1]), tmp)
#                })
#            matrix(ans,length(vec))
#        }
#        ans.mat=get.combn(N,n)
#    }else{  ## generate ramdom samples from all permutations
get.choose=function(N,n)exp(lfactorial(N)-sum(lfactorial(n)))
get.combn=function(vec, n, R){
if(length(vec)==n[1]) return(matrix(vec))
if(FALSE){  ## this is old implementation that is buggy, but still useful
k=choose(length(vec), n[1])
rest=get.choose(length(vec)-n[1],n[-1])
if(k<=R) {
this.ans=combn(length(vec), n[1])
if(rest*k>R){
sampled=table(sort(sample(rest*k, R)-1)%/%rest+1)
remainder=numeric(k)
remainder[as.numeric(names(sampled))]=sampled
}else   remainder=rep(rest,k)
}else {
this.ans=combn2R(length(vec), n[1],R=R) ## BUG: this does not allow duplicates
remainder=rep(1L, R)
}
}else{
remainder=table(bdSample(seq(choose(length(vec), n[1L])),
R,
get.choose(length(vec)-n[1L], n[-1L])
)
)
this.ans=combn(length(vec), n[1L])[,as.integer(names(remainder)),drop=FALSE]
}

ans.mat=matrix(NA_real_, length(vec), R)
cur.col=0
for(i in 1:ncol(this.ans)){
if(remainder[i]==0) next
tmp=Recall(vec[-this.ans[,i]], n[-1], remainder[i])
tmpans=rbind(matrix(rep(vec[this.ans[,i]],ncol(tmp)),n[1]), tmp)
ans.mat[,cur.col+1:ncol(tmpans)]=tmpans
cur.col=cur.col+ncol(tmpans)
}
ans.mat[,1:cur.col,drop=FALSE]
}
ans.mat=get.combn(1:N, n, B)
#    }
identity=which(colSums(abs(ans.mat-1:N))==0)
if(length(identity)==0){
ans.mat[,1]=1:N
}else{
ans.mat[,identity]=ans.mat[,1]
ans.mat[,1]=1:N
}
perms=ans.mat[order(trt),]
attr(perms,'.Random.seed')=save.seed
attr(perms,'trt')=trt
perms
}

bdSample=function(x, n, ubounds)
# sample n elements of vector x with replacement subject to the number of each element being sampled to be no larger than the corresponding upper bound ubounds
{
N=length(x)
x0=seq(N)
if(length(ubounds)!=N) ubounds=rep(ubounds, length=N)
stopifnot(sum(ubounds)>n)
counts=integer(N)
idx=rep(TRUE, N)
m=n
ans=c()
while(m>0L){
tmp=x0[idx][sample(sum(idx), m, replace=TRUE)]
tmpCt=table(tmp); unq=as.integer(names(tmpCt))
outOfBd= counts[unq]+tmpCt > ubounds[unq]
if(any(outOfBd)){
tmpCt[outOfBd]=ubounds[unq[outOfBd]] - counts[unq[outOfBd]]
idx[unq[outOfBd]]=FALSE
}
ans=c(ans, rep(unq, tmpCt))
counts[unq]=counts[unq]+tmpCt
m=m-sum(tmpCt)
}
if(length(ans)==1) x[ans] else x[sample(ans)]
}

mchooseZ=function(N, n)  ## multinomial coef for a vector of n
{
ans = chooseZ(N, n[1L])
if(length(n)==1L) ans else ans*Recall(N-n[1L], n[-1L])
}
}

nparts=function(n)
factorialZ(sum(n))/prod(c(factorialZ(n), factorialZ(table(n))))  ## total number of distinct trt assignments

permuteTrt <-
function(trt, B=100L, idxOnly = FALSE) ## permutation matrices for one way design
## returns a length(trt) by B matrix, if B is not larger than multinomial coefficient.
{
n=table(trt)
cn=cumsum(n)
ntrts=length(n)
N=cn[ntrts]
#mc=mchooseZ(N, n)
ordn=order(n, decreasing=TRUE)

SP=nparts(n)
part0=split(seq_len(N),trt);

if(B>=SP){ # list all partitions
sp=setparts(n)
B=ncol(sp)
#        for(i in seq(ntrts))  ans[[i]]=matrix( apply(sp==i,2L,function(xx)sort(which(xx)) ), ncol=B)
ans=split(row(sp),sp); ## this line and the next replace the previous line
for(i in seq(ntrts)) dim(ans[[i]])=c(length(ans[[i]])/B,B)
names(ans)=names(n)[ordn]

# swapping the original assignment to the first permutation
flag=TRUE
for(b in seq(B))
if(setequal(part0, lapply(ans, '[', , b))){ idx=b; flag=FALSE; break }
if(isTRUE(flag)){
warning("The first permutaiton may not be the original assignment.")
}
for(i in seq(ntrts)) {tmp=ans[[i]][,b]; ans[[i]][,b]=ans[[i]][,1L]; ans[[i]][,1L]=tmp}

if(isTRUE(idxOnly)){
warning("'idxOnly=TRUE' has not been implemented yet. Full results are returned.")
}#else{
attr(ans, 'idx') = NA_character_
class(ans)='permutedTrt'
#}
}else{   #sample from all permutations using factoradic number. Ideally, a sample from 1:SP should work, but how to do this without enumerating all SP possibilities using setparts?

decfr=HSEL.bigz(factorialZ(N), B)
idx=which(decfr==0L)
if(length(idx)>0L) decfr[idx]=decfr[1L]
decfr[1L]=as.bigz(0L)

decfrCC=drop(as.character(decfr))  ## for speed only
if(isTRUE(idxOnly)) {	## save memory
ans = lapply(part0, as.matrix)
attr(ans, 'idx') = decfrCC
class(ans) = 'permutedTrt'
return(ans)
}

ans=lapply(sapply(part0,length), matrix, data=NA_integer_, ncol=B)
buff = integer(N); buff[1L]
for(b in seq(B)){
# perm=dec2permvec(decfr[b],N)  ## This subsetting decfr[b] is the slowest part!
perm=dec2permvec(decfrCC[b],N)  ## This change speeds up for about 8~9X.
# for(i in seq(ntrts)) ans[[i]][,b]=sort.int(perm[part0[[i]]])
for(i in seq(ntrts)) ans[[i]][,b]=.Call(radixSort_prealloc, perm[part0[[i]]], buff)  ## radix sort with pre-allocated buffer space
}
names(ans)=names(part0)
attr(ans, 'idx') = NA_character_
class(ans)='permutedTrt'
}
ans
}

nperms.permutedTrt=function(permutedTrt)
{
if(is.na(attr(permutedTrt, 'idx')[1L])) {
ncol(permutedTrt[[1L]])
}else   length(attr(permutedTrt, 'idx'))
}

ntrt.permutedTrt=function(permutedTrt)
{
sapply(permutedTrt, NROW)
}

trt.permutedTrt=function(permutedTrt)
{ ## previous implementation was bugged; this is the corrected version.
ans=rep(NA_integer_, sum(sapply(permutedTrt,NROW)))#, levels=seq_along(permutedTrt))
for(i in seq_along(permutedTrt)) ans[permutedTrt[[i]][,1L]] = i
class(ans)='factor'
levels(ans)=names(permutedTrt)
ans
}
```

## Try the MRPP package in your browser

Any scripts or data that you put into this service are public.

MRPP documentation built on May 2, 2019, 4:46 p.m.