R/mbmdr.PermTest.R

Defines functions mbmdr.PermTest

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)){
	 cat("Output readed from file: ",x$call$output,"\n")
 	 x$result <- read.table(x$call$output,header=TRUE,sep=";",stringsAsFactors=FALSE)
 }
  
 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$ADJ.P)
   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){
   	cat("Model \"",m,"\" not found in the mbmdr object\n")
    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"
 # colnames(out)[7+length(model)] <- "Perm_ADJ.P"
 # out$Perm_PH <- ph
 # out$Perm_PL <- pl
 # aux <- c(ph,pl)
 # out$Perm_ADJ.P <- 2*min(aux[!is.na(aux)])
 }
 
 class(object) <- "mbmdr.PermTest"
 return(object)
}

Try the mbmdr package in your browser

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

mbmdr documentation built on May 30, 2017, 7:12 a.m.