permutation test

######### Permutation Test for WRR Model #########
permtest.wrr <- function(y, X, nPerm = 100, lambda = 1, 
                         center = TRUE, scale = TRUE, var.subset = NULL, ...){
  s <- stdize(X, center = center, scale = scale)
  X <- s$X
  meanX <- s$meanX
  stdX <- s$stdX
  if (is.null(var.subset)){
    var.subset <- 1:ncol(X)
  }
  
  res.all <- wrr(y, X, center = FALSE, scale = FALSE, lambda = lambda, ...)
  rmse.perm   <- matrix(NA, length(var.subset), nPerm + 1)
  effect.size <- rep(NA, length(var.subset))
  p.value     <- rep(NA, length(var.subset))
  
  
  ### Excute permutation test ###
  for (j in var.subset){
    
    for (p in 1:nPerm){
      
      ind <- sort(runif(nrow(X)), index.return = TRUE)$ix  
      X.perm   <- X
      X.perm[, j] <- X[ind, j]
      res.perm <- wrr(y, X.perm, center = FALSE, scale = FALSE, lambda = lambda, ...) ## Use all varibles (with permuted values in j-th variable) for weighted kernel ridge regression function
      rmse.perm[j, p] <- (sum(y - res.perm$yhat)^2)/length(y)^.5
    }
    
    res.minj <- wrr(y, X[, -j], center = FALSE, scale = FALSE, lambda = lambda, ...) ## Use varibles except j-th variable for weighted kernel ridge regression function
    effect.size[j] <- sum(y - res.all$yhat)^2 / sum(y - res.minj$yhat)^2
    rmse.perm[j, nPerm + 1] <- (sum(y - res.all$yhat)^2)/length(y)^.5
    p.value[j] <- 1 - sum(as.numeric(rmse.perm[j,] < rmse.perm[j, nPerm + 1]))/(nPerm + 1) ## Calculte the p-value
  }  
  
  mat <- cbind(effect.size, p.value)
  colnames(mat) <- c("effect.size","p-value")
  rownames(mat) <- colnames(X)[var.subset] 
  return(mat)
}