# R/mbmdr.PermTest.R In mbmdr: Model Based Multifactor Dimensionality Reduction

#### Documented in mbmdr.PermTest

```mbmdr.PermTest <-
function(x, n, model=NULL, sig.level=1){

if(!class(x)=="mbmdr"){
cat("Argument must be a mbmdr object generated by MBMDR function.\n")
return(NULL)
}
if(sig.level<0 || sig.level>1){
cat("Confidence level must be in [0,1] interval.\n")
return(NULL)
}

if(!is.null(x\$call\$output)){
}

if(is.null(model)){
ind <- order(pmax(x\$result\$WH,x\$result\$WL,na.rm=TRUE),decreasing=TRUE)
aux <- x\$result[ind,]
ind <- order(aux\$MIN.P)
model <- sapply(aux[ind[1],1:x\$call\$order],as.character)
}

model <- as.matrix(rbind(model))

object <- list()
object\$n <- n
object\$mbmdr <- x

for(i in 1:nrow(model)){
m <- model[i,]
m <- colnames(x\$data[1,m])
ind <- rep(TRUE, nrow(x\$result))
for(k in 1:length(m)) ind <- ind & x\$result[,k]==m[k]
if(sum(ind)==0){
next
}

result <- NULL
toCall <- x\$call
toCall\$data <- x\$data[,m]
toCall\$printStep1 <- NULL
toCall\$output <- NULL
toCall\$list.models <- NULL
toCall\$first.model <- NULL
for(k in 1:n){
toCall\$y <- sample(x\$y)
aux <- eval(toCall)\$result
result <- rbind(result,aux)
}
result\$Wmax <- apply(result[,c("WH","WL")],1,max,na.rm=TRUE)
wmax <- max(x\$result[ind,c("WH","WL")],na.rm=TRUE)
perm_p <- sum(result\$Wmax>wmax)/n
## out <- x\$result[ind,!colnames(x\$result) %in% c("PH","PL","MAX.P")]
out <- x\$result[ind,!colnames(x\$result) %in% c("PH","PL","MIN.P")]
out\$Wmax <- wmax
out\$Perm.P <- perm_p

if(sig.level<1){
if(perm_p < (5/n)){
out\$CI.lower <- NA
out\$CI.upper <- NA
}
else{
error <- qnorm(1-sig.level/2)*sqrt(perm_p*(1-perm_p)/n)
out\$CI.lower <- perm_p - error
out\$CI.upper <- perm_p + error
}
}

object\$PermTest <- rbind(object\$PermTest,out)

# ph <- ifelse(x\$result[ind,]\$NH>0,sum(result\$WH[!is.na(result\$WH)]>x\$result[ind,]\$WH)/n, NA)
# pl <- ifelse(x\$result[ind,]\$NL>0,sum(result\$WL[!is.na(result\$WL)]>x\$result[ind,]\$WL)/n, NA)
# out <- x\$result[ind,]
# colnames(out)[3+length(model)] <- "Perm_PH"
# colnames(out)[6+length(model)] <- "Perm_PL"
# out\$Perm_PH <- ph
# out\$Perm_PL <- pl
# aux <- c(ph,pl)