#' Graphical illustration of selection uncertainty
#' @export vsc.ci
#' @param vsc.fit fitted results from VSC::vsc()
#' @param level significance level
#' @param type todo
#' @examples
#' X = matrix(rnorm(1000), nrow = 100)
#' Y = rowSums(X[,1:3])+rnorm(100)
#' mod.0 = vsc(X, Y, intercept = FALSE, q = 0, method.names = NULL)
#' print(vsc.ci(mod.0, level = 0.1, type = "original"))
vsc.ci<-function(vsc.fit, level, type = "original"){
mods = NULL
B = vsc.fit$param$B
mods = NULL
if(type == "PR"){
mods = vsc.fit$PR.mod.collection[,-1] # first is intercept
} else if(type == "original"){
mods = vsc.fit$mod.collection[,-1] # first is intercept
} else if(type == "conditional"){
mods = vsc.fit$mod.collection[,-1]
mods[which(mods==0)] = NA
}
p = ncol(mods)
li = ui = NULL
if(type == "conditional"){
m = nrow(mods)
li = sapply(1:p, function(i){
num.non.zero = sum(!is.na(mods[,i]))
ci.margin = num.non.zero-(1-level)*m
if(ci.margin<=1){
return (0)
} else{
j = floor(ci.margin/2)
return (sort(mods[,i],partial=j)[j])
}
})
ui = sapply(1:p, function(i){
num.non.zero = sum(!is.na(mods[,i]))
ci.margin = num.non.zero-(1-level)*m
if(ci.margin<=1){
return (0)
} else{
j = floor(ci.margin/2)
return (sort(mods[,i],decreasing = TRUE)[j])
}
})
} else{
li = sapply(1:p, function(i)
stats::quantile(mods[,i], level/2, na.rm = TRUE))
ui = sapply(1:p, function(i)
stats::quantile(mods[,i], 1-level/2, na.rm = TRUE))
}
names(li) = colnames(mods)
li[which(is.na(li))] = 0
ui[which(is.na(ui))] = 0
return (rbind(li = li,
ui = ui))
}
# ci.helper<-function(mods, alpha, type, b = NULL){
# p = ncol(mods)
#
# li = sapply(1:p, function(i) stats::quantile(mods[,i], alpha/2))
# ui = sapply(1:p, function(i) stats::quantile(mods[,i], 1-alpha/2))
# names(li) = colnames(mods)
#
# if(type == "quantile"){
# return (rbind(li = li,
# ui = ui))
# } else{
# if(is.null(b)){
# b = sapply(1:p, function(i) median(mods[,i]))
# }
# basic.li = 2*b-ui
# basic.ui = 2*b-li
# if(type == "basic"){
# return (rbind(li = basic.li,
# ui = basic.ui))
# } else if(type == "basic2"){
# i = which(b!=0)
# li[i] = basic.li[i]
# ui[i] = basic.ui[i]
# return (rbind(li = li,
# ui = ui))
# } else if(type == "basic3"){
# i = which(b!=0 & (sign(ui)==sign(li) | sign(ui)==0 | sign(li)==0))
# li[i] = basic.li[i]
# ui[i] = basic.ui[i]
# return (rbind(li = li,
# ui = ui))
# }
# }
# }
# # here the models does not include intercept
# get.ci.from.mod.collection<-function(mods, level = 0.1, vsc.sel = NULL){
# ci.interval = HDCI::ci(Beta = NULL, Beta_bootstrap = mods, type = "quantile", a = level)
# rownames(ci.interval) = c("li", "ui")
# if(!is.null(vsc.sel)){
# li = ci.interval["li",]
# ui = ci.interval["ui",]
# is.sel = sign(ui)==sign(li) & sign(li)!=0 & vsc.sel!=0
# ci.interval["ui", which(!is.sel & ui<0)] = 0
# ci.interval["li", which(!is.sel & li>0)] = 0
# }
# return (ci.interval)
# }
#
# ## ==== LPR helper ====
lm.PR<-function (x, y, intercept, est.b, ...)
{
PR.mod = HDCI::PartRidge(x = x, y = y, intercept = intercept,
varset = est.b[-1] != 0, lambda2 = 1/length(y), ...)
return (c(PR.mod$beta0, PR.mod$beta))
}
get.PR.mod.collection<-function(x, y, intercept, flds, vsc.mod.collection,
variable.freq, num.core = 1){
# # vsc.mod.collection includes intercept!!!
# LPR.mod = do.call(rbind, lapply(1:length(flds),
# function(fld.i){
# fld = flds[[fld.i]]
# return (lm.PR(x = x[-fld,], y = y[-fld], intercept = intercept,
# est.b = vsc.mod.collection[fld.i,]))
# }))
# return (LPR.mod)
result = NULL
if(num.core>1){
cl = makeCluster(min(detectCores()-1, num.core), outfile = "parallel_log.txt")
clusterExport(cl=cl, list("lm.PR"), envir=environment())
result = parLapply(cl, 1:length(flds), function(fld.i, param){
fld = param$flds[[fld.i]]
return (lm.PR(x = param$x[-fld,], y = param$y[-fld], intercept = param$intercept,
est.b = param$vsc.mod.collection[fld.i,]))
}, param = list(x = x, y = y, intercept = intercept, flds = flds, vsc.mod.collection = vsc.mod.collection))
stopCluster(cl)
rm(cl)
} else{
result = lapply(1:length(flds),
function(fld.i){
fld = flds[[fld.i]]
return (lm.PR(x = x[-fld,], y = y[-fld], intercept = intercept,
est.b = vsc.mod.collection[fld.i,]))
})
}
return (do.call(rbind, result))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.