Nothing
#-----------------------------------------------------------------------------------------
#' Control parameters for \code{shapeselect()}
#'
#' Creates a list of parameters used by \code{shapeselect()} to control
#' basis construction, information-criterion calibration, simulation settings,
#' parallel computing etc.
#'
#' @param knots_lst Optional list specifying knots for spline basis functions,
#' one element per predictor. Use `NULL` to let \code{shapeselect()} (and its
#' helper routines) choose defaults.
#' @param boundary.knots Optional list specifiying knots for numeric predictors.
#' @param space Character string indicating the knots being equally spaced for spline
#' bases in underlying fitting routines (e.g., `"E"`).
#' @param nsim Integer; number of simulations used when estimating
#' mean degrees of freedom terms used in variable and shape selection criteria CBIC.
#' @param nfolds Integer; number of folds used for cross-validation when
#' \code{type = "classo"} (shrinkage selection). It is not in use for now.
#' @param eps Numeric tolerance used for detecting near-zero signal.
#' @param pnt Logical; whether to use penalized spline estimator.
#' @param multicore Logical/integer; parallel setting passed to \code{shapeselect()} via the
#' `parallel` control entry. Defaults to \code{getOption("cgam.parallel")}.
#' @param a Optional numeric; calibration constant used in the modified CBIC penalty term.
#' If `NULL`, \code{shapeselect()} chooses a value based on `n`, `p`, and the implied
#' total number of basis functions.
#' @param maxedf Logical; if `TRUE`, enforce an effective degrees of freedom cap in internal
#' forward-step routines (as implemented).
#' @param add_edfu Numeric; optional additive adjustment to the unconstrained EDF used in the
#' information criterion calculations.
#' @param add_kts Integer; optional adjustment to the number of knots used in spline bases.
#' If `NULL`, \code{shapeselect()} chooses a default depending on `pnt`.
#' @param shrinketa Logical; if `TRUE` and \code{type = "classo"}, shrink fitted component
#' contributions by the selected shrinkage parameters. It is not in use for now.
#' @param catnrep Logical; if `TRUE`, print the iteration counter during the selection path.
#' @param standardize Logical; if `TRUE`, standardize numeric predictors by their standard deviations.
#' @param center Logical; if `TRUE`, center numeric predictors when scaling them.
#' @param verbose Logical; if `TRUE`, print the model-path tracking matrix as the algorithm runs.
#' @param set4 Logical; optional flag passed through to internal routines used to restrict or
#' modify the candidate shape set to be increasing, decreasing, convex and concave; if false, then
#' combinations of monotonicity and convexity will be added
#'
#' @return A named list suitable for the \code{control} argument of \code{shapeselect()}.
#'
#' @examples
#' ctrl <- shapeselect_control(nsim = 100, nfolds = 10, verbose = TRUE)
#' ctrl
#'
#' @export
#'
shapeselect_control <- function(knots_lst = NULL, boundary.knots = NULL, space = "E", nsim = 200, nfolds = 5,
eps = 1e-8, pnt = TRUE, multicore = getOption("cgam.parallel"),
a = NULL, maxedf = FALSE, add_edfu = 0, add_kts = NULL,
shrinketa = FALSE, catnrep = FALSE, standardize = TRUE,
center = TRUE, verbose = FALSE, set4 = TRUE)
{
rslt <- list(knots_lst = knots_lst, boundary.knots = boundary.knots,
space = space, nsim = nsim, nfolds = nfolds,
eps = eps, pnt = pnt, parallel = multicore,
a = a, maxedf = maxedf, add_edfu = add_edfu, add_kts = add_kts,
shrinketa = shrinketa, catnrep = catnrep, standardize = standardize,
center = center, verbose = verbose, set4 = set4)
return (rslt)
}
#' Shape selection for additive components with shape or order constraints
#'
#' Fits an additive model by selecting variables with shape or order constraints
#' for each predictor using Cone Information Bayesian Criterion (CBIC).
#' The function can perform forward selection or stepwise (forward/backward) updates
#' over a set of candidate shapes appropriate to each predictor type.
#'
#' @param xmat A data frame of predictors. Columns may be numeric or factors.
#' Column names are used in output when available. If x is numeric, then the shape options are
#' increasing, decreasing, convex, concave or combination of monotonicity and convexity; if x is a
#' factor, then it can have monotonic or convex ordering options; if x is a factor, the option can also be "in" or "out-of" the model, or "in", "out-of" the model, or
#' following a "tree" ordering, depending on the `types` parameter.
#' @param y A response vector of length `nrow(xmat)`. It can be either numeric or binary.
#' @param types Optional character vector of length `ncol(xmat)` specifying the predictor type
#' for each column of `xmat`. This is used to determine which shape set is eligible for each
#' predictor. Values used for this parameter include `"num"`/`"numeric"` (continuous),
#' `"ord"`/`"ordinal"` (ordinal/categorical), `"tree"` (categorical),
#' and `"nom"`/`"nominal"` (categorical). Default is `NULL`, but the user must provide a value for each column of xmat;
#' otherwise the function will throw an error message.
#' @param direction Direction of the search. `"forward"` performs forward selection only.
#' `"both"` allows backward revisions (stepwise) after additions.
#' @param type Search strategy. `"stepwise"` adds predictors while allowing backward revisions.
#' `"classo"` then applies shrinkage/selection (e.g., via cross-validated lasso). The available option
#' is `"stepwise"` for now.
#' @param infocrit Information criterion used to decide adding or removing a predictor, and shape updates. Currently supports
#' `"CBIC"` and `"CAIC"`.
#' @param family A GLM family. May be a family function or a character string naming one. Available families are `gaussian` and `binomial`.
#' @param control A list of control parameters. Any entries not supplied fall back to
#' `shapeselect_control()`. Additional named arguments in `...` override both defaults and
#' `control`.
#' @param ... Additional control parameters (passed as overrides). These are merged into `control`
#' after defaults, so they take highest priority.
#'
#' @return An object of class `c("shapeselect","cgam")`, a list containing (among others):
#' \describe{
#' \item{mod_tried}{A matrix recording the selection path (chosen shapes, criterion values, adding or removing a predictor).}
#' \item{shps}{Numeric shape codes for each predictor (0 indicates excluded/flat).}
#' \item{shapes}{Character labels for the selected shapes.}
#' \item{etahat}{Fitted linear predictor values.}
#' \item{muhat}{Fitted mean response on the data scale (`family$linkinv(etahat)`).}
#' \item{etacomps}{Matrix of fitted component contributions by predictor (rows).}
#' \item{cic_vec}{Criterion values along the path.}
#' \item{cvec}{Shrinkage/selection vector (used when `type = "classo"`).}
#' \item{ddkeep_x_lst, ddkeep_z_lst}{Stored design/basis components for numeric and categorical predictors used for prediction and plotting.}
#' }
#'
#' @details
#' The algorithm begins with an intercept-only fit, then evaluates each inactive predictor under
#' a set of allowable shapes (depending on `types`). At each step it selects the predictor/shape
#' combination that improves the chosen information criterion. When `direction = "both"`, it
#' checks non-chosen shapes for already-selected predictors via backward steps.
#'
#' Shape and order choices in this implementation include
#' 1) smooth effects: monotone increasing/decreasing, convex/concave when `set4 = TRUE`, and combinations `set4 = FALSE`, if the corresponding value
#' in `types` is `"num"`/`"numeric"` (x is continuous).
#' 2) non-smooth effects: monotone increasing/decreasing and convex/concave (when `set4 = TRUE`), or their combinations (when `set4 = FALSE`),
#' if the corresponding value in `types` is `"ord"`/`"ordinal"` (i.e., x is ordinal, categorical, or numeric but treated as
#' discrete with an inherent ordering).
#' in `types` is `"ord"`/`"ordinal"` (x is ordinal/categorical/having numeric values but treated as discrete values with ordering).
#' 3) in or out-of the model, if the corresponding value in `types` is `"nom"`/`"nominal"` (x is categorical)
#' 4) in, out-of the model or following a tree ordering, if the corresponding value in `types` is `"tree"` (i.e., x is categorical or numeric but treated as
#' discrete with a tree ordering).
#' @examples
#' library(MASS)
#' data("Rubber", package = "MASS")
#' y <- Rubber[, 1]
#' xmat <- Rubber[, -1]
#' ans <- shapeselect(xmat, y, types = c("num", "num"), pnt = FALSE)
#'
#' print (ans$shapes)
#'
#' plot(ans)
#' @export
#-----------------------------------------------------------------------------------------
#coneA is used; one is included in bigmat
#main routine: shapeselect
#xmat must be a data frame
#tree: the first unique value of x is placebo
#add tree ordering when the type of x is 'tree': shps: 0 (out), 1 (tree), 2 (no order)
#-----------------------------------------------------------------------------------------
shapeselect <- function(xmat, y, types = NULL, direction = c("both", "forward"),
type = "stepwise", infocrit = c("CBIC", "CAIC"),
family = gaussian, control = list(),...)
{
cl <- match.call()
type <- match.arg(type)
infocrit <- match.arg(infocrit)
direction <- match.arg(direction)
if (is.null(types)) stop("'types' must be provided and have length ncol(xmat).")
defaults <- shapeselect_control()
extras <- list(...)
control <- modifyList(defaults, control)
control <- modifyList(control, extras)
#unpack the control list
knots_lst <- control$knots_lst
boundary.knots <- control$boundary.knots
space <- control$space
nsim <- control$nsim
nfolds <- control$nfolds
eps <- control$eps
pnt <- control$pnt
parallel <- control$parallel
a <- control$a
maxedf <- control$maxedf
add_edfu <- control$add_edfu
add_kts <- control$add_kts
shrinketa <- control$shrinketa
catnrep <- control$catnrep
standardize <- control$standardize
center <- control$center
verbose <- control$verbose
set4 <- control$set4
if (is.character(family))
family <- get(family, mode = "function", envir = parent.frame())
if (is.function(family))
family <- family()
if (is.null(family$family))
stop("'family' not recognized!")
#new
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
fml <- family$family
if(fml == "gaussian"){
wt.iter <- FALSE
} else {wt.iter <- TRUE}
xnms <- colnames(xmat)
if(is.null(xnms)){xnms <- paste('x', 1:NCOL(xmat), sep='')}
#new: use types to check if x is for spline or ordinal basis or just nominal
#use 18 for linear shape
#use 17 for flat
#shps <- 9:18
shps <- 9:17
#shps_char <- c('s.incr','s.decr','s.conv','s.conc','s.incr.conv','s.incr.conc','s.decr.conv','s.decr.conc', 'flat', 'linear')
shps_char <- c('s.incr','s.decr','s.conv','s.conc','s.incr.conv','s.incr.conc','s.decr.conv','s.decr.conc', 'flat')
mod_0 <- data.frame(x = character(), 'shape/order' = numeric(), cic = numeric(), check.names = FALSE)
n <- length(y)
cic_vec <- NULL
active <- NULL
x_rm <- NULL #keep track of redundant z predictors
im <- inactive <- 1:NCOL(xmat)
#--------------------------------
#step 0: project y onto 1
#--------------------------------
one <- rep(1, n)
xmat0 <- xmat #used xmat0 later in plot
#new
#xmat[] <- lapply(xmat, function(col) if(is.numeric(col)) scale(col, center = FALSE) else col)
p <- NCOL(xmat)
n <- NROW(xmat)
sdx <- rep(1, p)
mx <- rep(0, p)
#new: new rule for a
mult <- ifelse(pnt, 2, 1)
m_total <- 0
for(ix in 1:p){
if(is.numeric(xmat[,ix]) & types[ix] %in% c("num", "numeric")){
#if(is.numeric(xmat[,ix])){
#xmat[,i] <- (xmat[,i] - min(xmat[,i])) / (max(xmat[,i]) - min(xmat[,i]))
kts0 <- NULL
if(is.list(boundary.knots)){
kts0 <- boundary.knots[[ix]]
if (!is.numeric(kts0) | any(!is.finite(kts0)) | length(kts0) < 2){
warning(
"Invalid 'boundary.knots': must be a numeric vector of at least two finite values. ",
"Using range(x) instead.",
call. = FALSE
)
bounndary.knots[ix] <- list(NULL)
kts0 <- NULL
}
}
#new: check scale first
if(standardize){
xmat[,ix] <- scale(xmat[,ix], center=center)
#sdx <- c(sdx, attr(xmat[,ix], "scaled:scale"))
sdx[ix] <- attr(xmat[,ix], "scaled:scale")
if(center){
mx[ix] <- attr(xmat[, ix], 'scaled:center')
}
#new: scale user-defined boundary knots
if(length(kts0) >= 2){
if(center){
kts0 <- kts0 - mx[ix]
}
kts0 <- kts0 / sdx[ix]
boundary.knots[[ix]] <- kts0
}
}
#new:
n1 <- length(unique(xmat[,ix]))
m_total <- m_total + mult * (trunc(n1^(1/7)) + 6)
}
if(inherits(xmat[,ix], "factor") & types[ix] %in% "tree"){
xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
}
if(inherits(xmat[,ix], "factor") & types[ix] %in% "ord"){
xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
}
attr(xmat[, ix], "type") <- types[ix]
}
# for(ix in 1:p){
# if(is.numeric(xmat[,ix]) & types[ix] %in% c("num", "numeric")){
# #if(is.numeric(xmat[,ix])){
# #xmat[,i] <- (xmat[,i] - min(xmat[,i])) / (max(xmat[,i]) - min(xmat[,i]))
# xmat[,ix] <- scale(xmat[,ix], center=center)
# #sdx <- c(sdx, attr(xmat[,ix], "scaled:scale"))
# sdx[ix] <- attr(xmat[,ix], "scaled:scale")
# if(center){
# mx[ix] <- attr(xmat[, ix], 'scaled:center')
# }
# #new:
# n1 <- length(unique(xmat[,ix]))
# m_total <- m_total + mult * (trunc(n1^(1/7)) + 6)
# }
# if(inherits(xmat[,ix], "factor") & types[ix] %in% "tree"){
# xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
# }
# if(inherits(xmat[,ix], "factor") & types[ix] %in% "ord"){
# xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
# }
# attr(xmat[, ix], "type") <- types[ix]
# }
#new rule for a
if(is.null(a)){
a <- 1
# if(p >= n){
# a <- 2
# }
# #new: decide value of a according to the number of spline bases functions
# if(m_total > n & n > p){
# a <- 1
# }
# if(n > m_total){
# a <- 0.5
# }
}
#new rule for add_kts
if(is.null(add_kts)){
if(pnt){
add_kts <- 0
}else{
add_kts <- -2
}
}
#np will change later when including convex/concave x
np <- 1
zmat <- NULL
etahat <- muhat <- mean(y)
rss <- sum((y - etahat)^2)
#need a larger value to start cic/bic?
if(infocrit %in% c("CAIC", "CIC")) {
if(!wt.iter){
#test!
cic <- log(rss) + 2*1/n
} else {
#cat('fml:', fml, '\n')
llh <- llh.fun(y=y, muhat=rep(1/2, n), etahat=rep(0, n), n=n, fml=fml)
cic <- llh + 2*1/n
}
}
if(infocrit == "CBIC") {
if(!wt.iter){
#cic <- log(rss) + log(n) * 1 / n
#new
#cic <- log(rss) + a * log(n) * 1 / n^(6/7)
degree <- 6/7#1
cic <- log(rss) + a * log(n) * 1 / n^(degree)
} else {
m0 <- stats::glm(y ~ 1, family=family)
llh <- llh.fun(y=y, muhat=rep(1/2, n), etahat=rep(0, n), n=n, fml=fml)
#cic <- llh + log(n) * 1 / n
#new
#cic <- llh + a * log(n) * 1 / n^(6/7)
degree <- 6/7#1
cic <- llh + a * log(n) * 1 / n^(degree)
}
}
cic0 <- oldcic <- cic
cic_vec <- c(cic_vec, cic)
mod_0[1, 1] <- "1 vector"
mod_0[1, 2] <- "constant"
mod_0[1, 3] <- cic
y0 <- y
if(!wt.iter){
mu <- mean(y)
}
#----------------------------------------------------------------------
#step 1: regress y vs x1, y vs x2, ...., y vs xp
#----------------------------------------------------------------------
x_chosen <- shps_chosen <- NULL
#keep track of how many spline bases have been chosen
capm <- 0
ddkeep_z <- ddkeep_z_all <- ddkeep_x_all <- ddkeep_x <- NULL
#new: keep track of penalty matrix
dmat_keep <- NULL
ps_keep <- 0
capl <- NCOL(xmat) #capl includes z
#8 shapes for continuous x; 2 shapes for z
#include linear option
#1 is not in the projection space for residual
cic_st <- oldcic
rslt_step1 <- forward_cic(xmat, y, cic_st, infocrit, type, knots_lst, boundary.knots,
nsim, ddkeep_z_all, wt.iter, family, parallel, space, pnt, a,
maxedf, add_edfu, add_kts, set4)
#rslt_step1 <- forward_cic(xmat, y, cic_st, infocrit, type, knots_lst, nsim,
# ddkeep_z_all, wt.iter, family, parallel, space, pnt, a,
# maxedf, add_edfu, add_kts, set4)
cic_vals <- rslt_step1$cic_vals
#new: avoid wrong choice when e.g. conv and incr.conv has the diff of 1e-16
cic_vals <- round(cic_vals, 8)
#need to add a warning: is the smallest cic_val is for the flat shape for all predictors, then stop the algorithm
ddkeep_x <- rslt_step1$ddkeep_x
ddkeep_z <- rslt_step1$ddkeep_z
#etahatkeep is cumulative estimate; it's not smooth component chosen at each step
etahatkeep <- rslt_step1$etahatkeep
coefskeep <- rslt_step1$coefskeep
mskeep <- rslt_step1$mskeep
oldetahat <- etahatkeep
oldcoefs <- coefskeep
ktskeep <- rslt_step1$ktskeep
cic_st_next_step <- rslt_step1$cic_st_next_step
dmat_keep <- rslt_step1$dmat_keep
ps_keep <- rslt_step1$ps_keep
#NA's will be for categorical predictors
cic <- min(cic_vals, na.rm = TRUE)
#cic_vals = cic_vals[-9, ,drop = FALSE] #no need to remove -9; include flat
if(cic_st_next_step < oldcic){
x_index <- which(cic_vals == cic, arr.ind = TRUE)
} else {
#flat for all predictors
x_index <- which(cic_vals == round(oldcic, 8), arr.ind = TRUE)
}
x_index <- rbind(x_index)
#if(NROW(x_index) > 1 & any(x_index[, 1] == 9)){
if(NROW(x_index) > 1 & any(x_index[, 1] %in% c(5, 9))){
#if(NROW(x_index) > 1 & any(x_index[, 1] == nrow(cic_vals))){
#stop('check!')
#more than 1 x's best shape is flat; should we add them all at once?
x_index <- x_index[1, ]
} else if(NROW(x_index) > 1 & length(unique(x_index[, 2])) == 1){
#more than 1 shape of 1 x has the same cic_val; choose the one with smallest null dfmean
#dfmean_cbic_vals <- rslt_step1$dfmean_cbic_vals
col_ps <- unique(x_index[, 2])
dfmean_x_chosen <- dfmean_cbic_vals[, col_ps]
row_ps <- x_index[, 1]
x_index <- x_index[which.min(dfmean_x_chosen[row_ps]), ]
}
x_chosen <- c(x_chosen, inactive[x_index[2]])
#is_tree <- inherits(xmat_inactive[, x_index[2]], "factor") & attr(xmat_inactive[, x_index[2]], "type") %in% c("tree")
if(!inherits(xmat[, x_index[2]], "factor")){
if(attr(xmat[, x_index[2]], "type") %in% c("ord", "ordinal")){
shps <- c(1:8, 0)
if(wt.iter & infocrit == "CBIC"){
shps <- c(1:4, 0)
}
} else if (attr(xmat[, x_index[2]], "type") %in% c("tree")) {
shps <- c(1:2, 0)
} else {
shps <- 9:17
if(wt.iter & infocrit == "CBIC"){
shps <- c(9:12, 0)
}
}
shp_ch_i <- shps[x_index[1]]
shps_chosen <- c(shps_chosen, shps[x_index[1]])
} else {
if (cic == oldcic) {shp_ch_i <- 0}
if (cic < oldcic) {shp_ch_i <- 1}
shps_chosen <- c(shps_chosen, shp_ch_i)
}
#record chosen model's cic, x, shape
cic_vec <- c(cic_vec, cic)
colnames(cic_vals) <- colnames(xmat)
rownames(cic_vals) <- shps_char
#----------------------------------------------------------------------------------------------------
#step 2 to step k: regress y ~ x1 vs x2, ...., y ~ x1 vs xp
#stopping criteria: p = 0; oldcic < cic
#if classo = T, then stop when p = 0? go forward and backward to find the smallest CIC
#----------------------------------------------------------------------------------------------------
#mod_record_all <- list()
nrep <- 0
varlst <- varlst_z <- NULL #to keep track of delta's
#varlst <- NULL
#varlst_z <- 1 #start with one vector
#new keep track of predictors tested
#1 if x is in the model; 0 otherwise
totx <- NCOL(xmat)
mod_tried <- matrix(0, nrow = 1000, ncol = (totx + 3))
#colnames(mod_tried) = c(xnms, "# of x's", "cic", "new add")
xnms_with_id <- paste0(xnms, '|', 1:totx)
colnames(mod_tried) <- c(xnms_with_id, "# of x's", infocrit, "new add")#, "logrss", "edf0", "edf_max")
#used to create bigmat in the order cgam uses
is_z <- sapply(xmat, function(e) is.factor(e))
obs <- 1:ncol(xmat)
zid <- obs[is_z]
xid <- obs[!is_z]
#new:
mskeep_lst <- ddkeep_x_lst <- ddkeep_z_lst <- vector("list", length = totx)
dmat_keep_lst <- vector("list", length = totx)
ps_keep_vec <- rep(0, length = totx)
if(is.null(knots_lst)){
knots_lst <- vector("list", length = totx)
}
check1 <- cic_st_next_step < oldcic & length(inactive) > 0
rm_id <- NULL
x_rm <- NULL
ansi_back <- NULL
#new:
coefskeep_x_lst <- vector("list", length = totx)
coefskeep_z_lst <- vector("list", length = totx)
varlst_all <- NULL
elapsed_times <- NULL
ncols <- NULL
sttime <- proc.time()[3]
while(check1){
new_add <- inactive[x_index[2]]
x_rm_nrep <- inactive[rm_id]
x_rm <- c(x_rm, x_rm_nrep)
inactive <- inactive[-c(x_index[2], rm_id)]
active <- im[!im %in% inactive]
nrep <- nrep + 1
if(catnrep){
cat('nrep: ', nrep, '\n')
}
end_time <- proc.time()[3]
elapsed_time <- end_time - sttime
#cat('elapsed time: ', elapsed_time, '\n')
sttime <- end_time
elapsed_times <- c(elapsed_times, elapsed_time)
#new:
if(nrep > 1) {
mod_tried[nrep, ] <- mod_tried[nrep-1, ]
}
new_add_num <- id_ch <- new_add
#use she chosen shape for id_ch instead?
mod_tried[nrep, id_ch] <- rev(shps_chosen)[1]#1
mod_tried[nrep, totx + 1] <- length(active)
#------------
#test more!
#------------
#mod_tried[nrep, totx + 2] <- cic
mod_tried[nrep, totx + 2] = cic_st_next_step
mod_tried[nrep, totx + 3] <- id_ch
nac <- mod_tried[nrep, totx + 1]
#test more
if(!is.null(rm_id)){
nrep <- nrep + 1
mod_tried[nrep, ] <- mod_tried[nrep-1, ]
mod_tried[nrep, totx + 3] <- x_rm_nrep
}
if(verbose){
formatted_df <- mod_tried[1:nrep, , drop = FALSE]
cat(capture.output(print(formatted_df)), sep = "\n")
}
#etahat_all = etahat_all + etahatkeep
ddkeep_z_all <- cbind(ddkeep_z_all, ddkeep_z)
ddkeep_x_all <- cbind(ddkeep_x_all, ddkeep_x)#may not follow the order of x, may be 1,3,2.
np <- NCOL(ddkeep_z_all) #no intercept in it
capm <- NCOL(ddkeep_x_all)
if(!is.null(ddkeep_x)) {
ddkeep_x_lst[[new_add_num]] <- ddkeep_x
mskeep_lst[[new_add_num]] <- mskeep
if(!is.null(ktskeep)){
knots_lst[[new_add_num]] <- ktskeep
}else{
knots_lst[new_add_num] <- list(NULL)
}
#new:
dmat_keep_lst[[new_add_num]] <- dmat_keep
ps_keep_vec[new_add_num] <- ps_keep
#varlst not record z
varlst <- c(varlst, rep(new_add_num, NCOL(ddkeep_x)))
#new
#varlst_all <- c(varlst_all, rep(new_add_num, NCOL(ddkeep_x)))
#coefskeep_x_lst[[new_add_num]] <- (coefskeep[-1])[which(varlst_all %in% new_add_num)]
}
if(!is.null(ddkeep_z)) {
ddkeep_z_lst[[new_add_num]] <- ddkeep_z
mskeep_lst[[new_add_num]] <- mskeep
varlst_z <- c(varlst_z, rep(new_add_num, NCOL(ddkeep_z)))
}
#new: go backwards to see if other shapes of previously added predictors will decrease cic
if(direction == "both"){
if(nac >= 2 & (!shp_ch_i %in% c(0, 17))){
chs <- which(mod_tried[nrep-1, 1:totx] != 0)
#in backward step, use smaller oldcic?
oldcic <- cic_st_next_step
#new: used to determin degree of n in CBIC
shps_all <- mod_tried[nrep, 1:totx]
#cat('start backward step: ', '\n')
ansi_back <- backward_cic(xmat, y, varlst, varlst_z, ddkeep_x_all, ddkeep_x_lst, knots_lst,
boundary.knots,
ddkeep_z_all, ddkeep_z_lst, mskeep_lst, mod_tried, nrep, totx,
oldcic, infocrit, dfmean_old, dfmeanmax_old,
penalty_old, oldcoefs, oldetahat,
active, type, nsim, wt.iter, family,
parallel, space, pnt, dmat_keep_lst, ps_keep_vec,
a, add_edfu, add_kts, shps_all)
# ansi_back <- backward_cic(xmat, y, varlst, varlst_z, ddkeep_x_all, ddkeep_x_lst, knots_lst,
# ddkeep_z_all, ddkeep_z_lst, mskeep_lst, mod_tried, nrep, totx,
# oldcic, infocrit, dfmean_old, dfmeanmax_old,
# penalty_old, oldcoefs, oldetahat,
# active, type, nsim, wt.iter, family,
# parallel, space, pnt, dmat_keep_lst, ps_keep_vec,
# a, add_edfu, add_kts, shps_all)
#cat('finished backward step: ', '\n')
if(ansi_back$changed){
back_inactive <- ansi_back$back_inactive
if(length(back_inactive) > 0) {
inactive <- sort(c(inactive, back_inactive))
active <- im[!im %in% inactive]
}
#create new ddkeep_x_all, ddkeep_z_all
ddkeep_x_lst <- ansi_back$ddkeep_x_lst
#new:
coefskeep_x_lst <- ansi_back$coefskeep_x_lst
coefskeep_z_lst <- ansi_back$coefskeep_z_lst
#print (coefskeep_z_lst)
#new:
dmat_keep_lst <- ansi_back$dmat_keep_lst
ps_keep_vec <- ansi_back$ps_keep_vec
ddkeep_z_lst <- ansi_back$ddkeep_z_lst
#print (ddkeep_z_lst)
knots_lst <- ansi_back$knots_lst
mod_tried <- ansi_back$mod_tried
mskeep_lst <- ansi_back$mskeep_lst
cic <- ansi_back$cic_st
cic_st_next_step <- ansi_back$cic_st_next_step
#xtxkeep <- ansi_back$xtx
#cat('backward step xtx:', dim(xtxkeep), '\n')
#face <- ansi_back$face
etahatkeep <- ansi_back$etahatkeep
coefskeep <- ansi_back$coefskeep
oldcoefs <- coefskeep
oldetahat <- etahatkeep
seq_added <- mod_tried[1:nrep, totx + 3]
# ddkeep_x_lst_ord will only store dd's for chosen x
ddkeep_x_lst_ord <- ddkeep_x_lst[seq_added]
ddkeep_z_lst_ord <- ddkeep_z_lst[seq_added]
#keep the adding order of x, may not be 1,2,3; may be 1,3,2
ddkeep_x_all <- do.call(base::cbind, ddkeep_x_lst_ord)
ddkeep_z_all <- do.call(base::cbind, ddkeep_z_lst_ord)
#new:
#dmat_keep_all <- as.matrix(bdiag(Filter(Negate(is.null), dmat_keep_lst)))
#need to create varlst and varlst_z again
seq_added <- mod_tried[1:nrep, totx + 3]
obs <- 1:totx
obs_ord <- obs[seq_added]
if(!is.null(ddkeep_x_all)){
varlst_lst <- lapply(1:length(obs_ord), function(i, xlst) {if(NCOL(xlst[[i]]) >= 1) {rep(obs_ord[i], NCOL(xlst[[i]]))}}, xlst = ddkeep_x_lst_ord)
#print (varlst_lst)
varlst <- do.call(base::c, varlst_lst)
} else {varlst <- NULL}
if(!is.null(ddkeep_z_all)){
varlst_z_lst <- lapply(1:length(obs_ord), function(i, zlst) {if(NCOL(zlst[[i]]) >= 1) {rep(obs_ord[i], NCOL(zlst[[i]]))}}, zlst = ddkeep_z_lst_ord)
varlst_z <- do.call(base::c, varlst_z_lst)
} else {varlst_z <- NULL}
}
}
}
xmat_inactive <- xmat[, -active, drop=FALSE]
if(length(inactive) == 0) {break}
#oldcic <- cic
#cic_st <- cic
oldcic <- cic_st_next_step
cic_st <- cic_st_next_step
knots_lst_inactive <- NULL
boundary.knots_inactive <- NULL
if(!is.null(knots_lst)){
knots_lst_inactive <- knots_lst[-active]
}
if(!is.null(boundary.knots)){
boundary.knots_inactive <- boundary.knots[-active]
}
#cic_st should be the same as oldcic?
#cic_st is a little larger than cic_st_fixed_shps; it is cic val for 17; used in classo to make every x to be chosen??
#new: use shapes chosen so far to determine the denominator of CBIC
shps_all <- mod_tried[nrep, 1:totx]
shps_non_zero <- shps_all[!shps_all %in% c(0, 17)]
#new! dmat_keep_all should be of the same order as ddkeep_x_all
#cannot use this: may add some x twice if the x was out of the model in the previous step: 2 5 3 1 6 2
seq_added <- mod_tried[1:nrep, totx + 3]
#test!
seq_added <- seq_added[!duplicated(seq_added, fromLast = TRUE)]
ansi <- forward_cic_step2(xmat_inactive, y, cic_st, ddkeep_z_all, ddkeep_x_all, knots_lst_inactive,
boundary.knots_inactive,
oldcic, oldetahat, oldcoefs, infocrit, dfmean_old, dfmeanmax_old, penalty_old,
active, inactive, type, nsim, wt.iter, family, parallel, space, pnt,
dmat_keep_lst, ps_keep_vec, a, maxedf, add_edfu, add_kts, shps_non_zero, seq_added, set4)
#new! dmat_keep_all should be of the same order as ddkeep_x_all
#seq_added <- mod_tried[1:nrep, totx + 3]
#ansi <- forward_cic_step2(xmat_inactive, y, cic_st, ddkeep_z_all, ddkeep_x_all, knots_lst_inactive,
# oldcic, oldetahat, oldcoefs, infocrit, dfmean_old, dfmeanmax_old, penalty_old,
# active, inactive, type, nsim, wt.iter, family, parallel, space, pnt,
# dmat_keep_lst, ps_keep_vec, a, maxedf, add_edfu, add_kts, shps_non_zero, seq_added, set4)
#cat('finished 2nd step', ncol(ddkeep_x_all), '\n')
ncols <- c(ncols, ncol(ddkeep_x_all))
cic_vals <- ansi$cic_vals
#new: avoid tie of e.g. conv and s.incr.conv
cic_vals <- round(cic_vals, 8)
cic <- min(cic_vals, na.rm = TRUE)
ddkeep_x <- ansi$ddkeep_x
ddkeep_z <- ansi$ddkeep_z
etahatkeep <- ansi$etahatkeep
coefskeep <- ansi$coefskeep
oldetahat <- etahatkeep
oldcoefs <- coefskeep
rm_id <- ansi$rm_id
mskeep <- ansi$mskeep
#penalty_old <- ansi$penalty
ktskeep <- ansi$ktskeep
cic_st_next_step <- ansi$cic_st_next_step #used edf0 not edf_m
dmat_keep <- ansi$dmat_keep
ps_keep <- ansi$ps_keep
#loop through predictors in the inactive set
if(cic_st_next_step < oldcic) {
x_index <- which(cic_vals == min(cic_vals, na.rm = TRUE), arr.ind = TRUE)
} else {
x_index <- which(cic_vals == round(oldcic, 8), arr.ind = TRUE)
}
x_index <- rbind(x_index)
#test more: tie
#flat should be the 5th if y is binomial and CBIC
#if(NROW(x_index) > 1 & any(x_index[, 1] == 9)){
if(NROW(x_index) > 1 & any(x_index[, 1] %in% c(5, 9))){
#if(NROW(x_index) > 1 & any(x_index[, 1] == nrow(cic_vals))){
#stop('check!')
#more than 1 x's best shape is flat; should we add them all at once?
x_index <- x_index[NROW(x_index), ]
} else if(NROW(x_index) > 1 & length(unique(x_index[, 2])) == 1){
#more than 1 shape of 1 x has the same cic_val; choose the one with smallest null dfmean
dfmean_cbic_vals <- ansi$dfmean_cbic_vals
if(!anyNA(dfmean_cbic_vals)){
col_ps <- unique(x_index[, 2])
dfmean_x_chosen <- dfmean_cbic_vals[, col_ps]
row_ps <- x_index[, 1]
x_index <- x_index[which.min(dfmean_x_chosen[row_ps]), ]
} else {
#temp; no simulation for dfmean_cbic; choose any of the tie
x_index <- x_index[NROW(x_index), ]
}
}
x_index <- rbind(x_index)
x_chosen <- c(x_chosen, inactive[x_index[2]])
# print (x_chosen)
# print (head(xmat_inactive))
# print (x_index)
#is_tree <- inherits(xmat_inactive[, x_index[2]], "factor") & attr(xmat_inactive[, x_index[2]], "type") %in% c("tree")
if(!inherits(xmat_inactive[, x_index[2]], "factor")){
if(attr(xmat_inactive[, x_index[2]], "type") %in% c("ord", "ordinal")){
#cat('check shapes', '\n')
shps <- c(1:8, 0)
if(wt.iter & infocrit == "CBIC"){
shps <- c(1:4, 0)
}
} else if(attr(xmat_inactive[, x_index[2]], "type") %in% c("tree")) {
shps <- c(1:2, 0)
}else {
shps <- 9:17
if(wt.iter & infocrit == "CBIC"){
shps <- c(9:12, 0)
}
}
shp_ch_i <- shps[x_index[1]]
shps_chosen <- c(shps_chosen, shps[x_index[1]])
} else {
shp_ch_i <- ifelse(is.null(ddkeep_z), 0, 1)
shps_chosen <- c(shps_chosen, shp_ch_i)
}
#record chosen model's cic, x, shape
cic_vec <- c(cic_vec, cic)
colnames(cic_vals) <- colnames(xmat_inactive)
rownames(cic_vals) <- shps_char
#new
if(type == "classo") {
check1 <- length(inactive) > 0
}
if(type == "stepwise") {
check1 <- cic_st_next_step < oldcic & length(inactive) > 0 #works for diabetes together with cic_st = cic
#check1 = cic < oldcic & length(inactive) > 0
}
}
#cat('final step xtx:', dim(xtxkeep), '\n')
#new
ch1 <- is.null(ansi_back)
ch2 <- FALSE
if(!is.null(ansi_back)){
if(!ansi_back$changed) {
ch2 <- TRUE
}
}
if(ch1 | ch2){
if(length(varlst) > 0){
ux_added <- unique(varlst)
for(x_added in ux_added){
if(np == 0){
coefskeep_x_lst[[x_added]] <- (coefskeep[-1])[which(varlst %in% x_added)]
}
if(np > 0){
coefskeep_x_lst[[x_added]] <- (coefskeep[-c(1:(1+np))])[which(varlst %in% x_added)]
}
}
}
#new
if(length(varlst_z) > 0){
uz_added <- unique(varlst_z)
for(z_added in uz_added){
coefskeep_z_lst[[z_added]] <- (coefskeep[2:(1+np)])[which(varlst_z %in% z_added)]
}
}
}
#the order of gmat in ansi_back is not following the order of seq_added: need to get coefskeep here
#for etacomps
#new:
dv <- 0
if(pnt & nrep >= 1){
dmat <- Matrix::bdiag(Filter(Negate(is.null), dmat_keep_lst))
dmat <- as.matrix(dmat)
if(prod(dim(dmat)) != 0){
np_ix <- NCOL(ddkeep_z_all) + 1
if(np_ix > 0){
dm_zero <- matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat <- cbind(dm_zero, dmat)
}
dv <- crossprod(dmat)
}
}
#1 df for the intercept
#dfmeankeep <- dfmean_old #+ 1
#penalty <- penalty_old
if(nrep == 0) {
#mod_tried <- mod_0
mod_tried <- mod_tried[1, ,drop = FALSE]
rownames(mod_tried) <- paste("step", 0, sep = "")
mod_tried[, totx+2] <- mod_0$cic
yhat <- mu
} else if (nrep >= 1) {
mod_tried <- mod_tried[seq(nrep), ,drop = FALSE] #
#test!
#nrep <- min(totx, nrep)
mod_tried <- mod_tried[!duplicated(mod_tried[, totx + 3]), ,drop=FALSE]
nrep <- nrow(mod_tried)
rownames(mod_tried) <- paste("step", seq(nrep), sep = "")
seq_added <- mod_tried[, totx + 3]
#in case some z is redundant; rm_id is not NULL
yhat <- etahatkeep
}
#yhat = etahatkeep + mu
if(!wt.iter){
r2 <- 1 - sum((y0 - yhat)^2) / sum((y0 - mu)^2)
} else {
r2 <- NULL #change to be pseudo r2 later
}
#----------------------------------------------------------------------------------------------------
#etacomps store fitted components NOT orthogonal to 1 vector or [1,x] if convex/concave
#----------------------------------------------------------------------------------------------------
zcoefs <- xcoefs <- NULL
coefskeep_withone <- coefskeep #not ordered
coefskeep <- coefskeep[-1] #remove intercept coef
#new xcoefs and zcoefs should appear in the order of seq_add; otherwise etacomps is wrong!
if(!is.null(varlst_z)){
nz <- NCOL(ddkeep_z_all)
zcoefs <- coefskeep[1:nz]
xcoefs <- coefskeep[-c(1:nz)]
#test!
#zcoefs <- zcoefs[order(varlst_z)]
#xcoefs <- xcoefs[order(varlst)]
} else {
xcoefs <- coefskeep
#test!
#xcoefs <- xcoefs[order(varlst)]
}
rm_rep_x <- NULL
if(nrep >= 1) {
seq_added <- mod_tried[, totx + 3]
shps <- (mod_tried[nrep, 1:totx]) #not ordered
shps_added <- shps[seq_added]
#print (shps_added)
xnms_added <- xnms[seq_added]
xnms_cont <- xnms[!shps %in% 1]
xmat_cont <- xmat[, !shps %in% 1, drop = FALSE]
shps_cont <- shps[!shps %in% 1]
} else if (nrep == 0) {
xnms_cont <- xnms
xmat_cont <- xmat
#new:
shps <- rep(0, totx)
shps_cont <- shps
seq_added <- NULL
#print (shps_added)
}
#print (length(mskeep_lst))
#if a chosen shp is 17/flat, then let etacomp = 0; it will be removed at fmat
etacomps <- matrix(0, nrow = totx, ncol = n)
coefskeep_lst <- vector("list", length = totx)
#coefskeep_x_lst <- vector("list", length = totx)
#coefskeep_z_lst <- vector("list", length = totx)
if(nrep >= 1) {
totx_ch <- length(seq_added)
for(i in 1:totx_ch){
#print (i)
ps <- seq_added[i]
xnmi <- xnms_added[i]
#print (ps)
xi <- xmat[, ps]
shpi <- shps_added[i]
#if(shpi == 17 | ps %in% x_rm){#flat
if(shpi == 17 | shpi == 0){#flat
next
} else {
#one <- cbind(rep(1, n))
#xi = as.numeric(xi)
#if(shpi == 18) {
# etacomps[ps, ] <- ddkeep_z_all[, varlst_z %in% ps, drop = FALSE] %*% zcoefs[which(varlst_z %in% ps)]
# coefskeep_lst[[ps]] <- zcoefs[which(varlst_z %in% ps)]
#} else
if ((shpi == 1 & inherits(xi, "factor")) | (shpi == 2 & attr(xi, "type") == "tree")) {#need to use 1 for incr later
#etacomps[ps, ] <- (ddkeep_z_all[, which(varlst_z %in% ps), drop = FALSE]) %*% zcoefs[which(varlst_z %in% ps)]
etacomps[ps, ] <- ddkeep_z_lst[[ps]] %*% coefskeep_z_lst[[ps]]
#coefskeep_lst[[ps]] <- zcoefs[which(varlst_z %in% ps)]
#coefskeep_x_lst
#coefskeep_z_lst[[ps]] <- zcoefs[which(varlst_z %in% ps)]
#capk <- capk + length(zcoefs[which(varlst_z %in% ps)])
#zmat <- cbind(zmat, ddkeep_z_all[, which(varlst_z %in% ps), drop = FALSE])
} else if (shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5) {
#etacomps[ps, ] <- (ddkeep_x_all[, which(varlst %in% ps), drop = FALSE]) %*% xcoefs[which(varlst %in% ps)]
#etacomps[ps, ] <- etacomps[ps, ] + (xi-mean(xi)) * zcoefs[which(varlst_z %in% ps)]
etacomps[ps, ] <- ddkeep_x_lst[[ps]] %*% coefskeep_x_lst[[ps]]
#etacomps[ps, ] <- etacomps[ps, ] + ddkeep_z_lst[[ps]] %*% coefskeep_z_lst[[ps]]
etacomps[ps, ] <- etacomps[ps, ] + ddkeep_z_lst[[ps]] * coefskeep_z_lst[[ps]]
#coefskeep_lst[[ps]] <- c(xcoefs[which(varlst %in% ps)], zcoefs[which(varlst_z %in% ps)])
#coefskeep_x_lst[[ps]] <- xcoefs[which(varlst %in% ps)]
#coefskeep_z_lst[[ps]] <- zcoefs[which(varlst_z %in% ps)]
#print (ps)
#print (coefskeep_z_lst)
#zmat_x <- cbind(zmat_x, (xi-mean(xi)))
} else {
#etacomps[ps, ] <- (ddkeep_x_all[, which(varlst %in% ps), drop = FALSE]) %*% xcoefs[which(varlst %in% ps)]
etacomps[ps, ] <- ddkeep_x_lst[[ps]] %*% coefskeep_x_lst[[ps]]
#coefskeep_lst[[ps]] <- xcoefs[which(varlst %in% ps)]
#coefskeep_x_lst[[ps]] <- xcoefs[which(varlst %in% ps)]
}
}
}
}
#k-fold CV to find shrinkage for each component
if(nrep >= 1 & type == "classo"){
fmat <- t(etacomps)
colnames(fmat) <- xnms
if(anyNA(fmat)) {
print (length(coefskeep))
print (NCOL(ddkeep_x_all))
print (NCOL(ddkeep_z_all))
print (head(fmat))
stop('check fmat!')
}
one <- rep(1, n)
normx <- sqrt(drop(one %*% (fmat^2)))
nosignal <- normx / sqrt(n) < eps
}
nng.mod <- NULL
cv_lasso <- NULL
mod_final <- mod_0
if(nrep == 0) {
cvec <- rep(0, totx)
#mod_final = mod_0
cv_rslt <- list(cvec = cvec)
mod_tried <- mod_0
shps_added <- NULL
#colnames(mod_tried)[totx + 4] = "shrink parameter"
} else if(nrep >= 1) {
if(type == "classo") {
cvec0 <- rep(0, NCOL(fmat))
if(all(nosignal)) {
cvec <- cvec0
cv_rslt <- list(cvec = cvec)
} else {
#cv_rslt <- cv.classo(fmat, y0, nfolds = nfolds, type = "glmnet", family = family)
#cvec <- cv_rslt$cvec
if(any(nosignal)){
nosignal_id <- which(nosignal)
fmat <- fmat[, -nosignal_id, drop = FALSE]
}
cv_rslt <- cv.classo(fmat, y0, nfolds = nfolds, type = "glmnet", family = family)
cvec <- cv_rslt$cvec
if(any(nosignal)){
cvec0[-nosignal_id] <- cvec
cvec <- cvec0
}
}
cvec00 <- rep(0, NROW(mod_tried))
mod_tried <- cbind(mod_tried, cvec)
colnames(mod_tried)[NCOL(mod_tried)] <- "shrink parameter"
#new: change etacomps and etahatkeep here! each etacomps should be multiplied by its cvec
#need more test!
#use for predict.classo
totx_ch <- length(seq_added)
for(i in 1:totx_ch){
ps <- seq_added[i]
shpi <- shps_added[i]
if(shpi == 17 | shpi == 0){#flat
next
} else {
if(cvec[ps] == 0){
if((shpi == 1 & inherits(xmat[,ps], "factor")) | (shpi == 2 & attr(xmat[,ps], "type") == "tree") | shpi %in% c(11, 12) | shpi %in% c(3, 4)){
coefskeep_z_lst[ps] <- list(NULL)
}
if(shpi >= 9 & shpi <= 16 | shpi >= 1 & shpi <= 8){
coefskeep_x_lst[ps] <- list(NULL)
}
etacomps[ps, ] <- 0
etahatkeep <- colSums(etacomps) + coefskeep_withone[1]
yhat <- etahatkeep
}
if(cvec[ps] > 0 & shrinketa){
etacomps[ps, ] <- cvec[ps] * etacomps[ps, ]
etahatkeep <- colSums(etacomps) + coefskeep_withone[1]
yhat <- etahatkeep
}
}
}
} else {
seq_added <- mod_tried[1:nrep, totx + 3]
x_kept <- seq_added
cvec <- rep(0, totx)
if(any(shps_added %in% c(0, 17))) {
rm_ps <- which(shps_added %in% c(0, 17))
x_kept <- x_kept[-rm_ps]
}
cvec[x_kept] <- 1
cv_rslt <- NULL
}
if(type == "classo") {
shps[shps == 17] <- 0
shps_added_char <- sapply(shps_added, ShapeToChar_shapeselect)
cvec_added <- mod_tried[, NCOL(mod_tried)]
if(any(cvec_added < eps)) {
ps_out <- which(cvec_added < eps)
shps_added_char[ps_out] <- 'out-of-mod'
#new: test more
shps[ps_out] <- 0
}
} else {
shps[shps == 17] <- 0
shps_added_char <- sapply(shps_added, ShapeToChar_shapeselect)
}
}
#mod_final <- mod_tried[nrow(mod_tried), ,drop = FALSE]
#new: create some terms matching terms in cgam, to use vcov
#if(!wt.iter){
if(nrep >= 1){
ddx <- do.call(base::cbind, ddkeep_x_lst)
zkeep_id <- intersect(obs[sapply(ddkeep_z_lst, function(e) !is.null(e))], zid)
zmat <- do.call(base::cbind, ddkeep_z_lst[zkeep_id]) # no one, exclude conv/conc, used for plot.classo
capk <- NCOL(zmat)
xckeep_id <- setdiff(obs[sapply(ddkeep_z_lst, function(e) !is.null(e))], zid)
xcmat <- do.call(base::cbind, ddkeep_z_lst[xckeep_id]) #conv/conc/linear in the linear space# test linear
zmat <- cbind(zmat, xcmat)
} else {
ddx <- NULL
zmat <- NULL
capk <- 0
}
#match cgam
etahat <- yhat #added ybar if gaussian
muhat <- family$linkinv(etahat)
wt <- rep(1, n) #temp, replace it with what cgam has later
if(nrep >= 1){
bigmat <- t(cbind(one, zmat, ddx))
d0 <- NCOL(zmat) + 1 #add constant back
edf <- sum(abs(coefskeep) > 1e-12) #include 1
} else {
bigmat <- NULL
d0 <- 1
edf <- 0
}
#new: change shps to be words
tags <- rep('x', p)
tags[which(types %in% c('nom', 'nominal'))] <- 'z'
tags[which(types %in% 'tree')] <- 'tree'
#shapes in character
shapes <- mapply(ShapeToChar_shapeselect, shps, tags)
rslt <- list(mod_tried = mod_tried, etahatkeep = etahatkeep, ddkeep_x_all = ddkeep_x_all,
knots_lst = knots_lst, ddkeep_z_all = ddkeep_z_all,
mskeep_lst = mskeep_lst, varlst = varlst, varlst_z = varlst_z, coefskeep = coefskeep,
r2 = r2, cic_vec = cic_vec, d0 = d0, wt = wt, y = y0, edf = edf,
etahat = etahat, muhat = muhat, etacomps = etacomps, coefskeep_lst = coefskeep_lst,
cvec = cvec, cv_rslt = cv_rslt, family = family, dv = dv,
shps_added = shps_added, bigmat = bigmat, capk = capk, ddkeep_x_lst = ddkeep_x_lst,
ddkeep_z_lst = ddkeep_z_lst,
zcoefs = zcoefs, xcoefs = xcoefs, nfolds = nfolds, xmat = xmat0, xmat_cont = xmat_cont,
shps = shps, shapes = shapes, shps_cont = shps_cont, xnms = xnms, xnms_cont = xnms_cont, rm_id=rm_id,
ps_keep_vec = ps_keep_vec, type = type, space = space, pnt = pnt,
sdx = sdx, mx = mx, coefskeep_withone = coefskeep_withone, seq_added = seq_added, dmat_keep = dmat_keep,
coefskeep_x_lst = coefskeep_x_lst, coefskeep_z_lst = coefskeep_z_lst,
types = types, center = center, standardize=standardize, boundary.knots=boundary.knots,
elapsed_times = elapsed_times, a = a) #, ncols=ncols, a = a, m_total=m_total)
class(rslt) <- c("shapeselect", "cgam")
return (rslt)
}
#------------------------------------------
#create penalty matrix
#difference from quadratic or cubic splines
#for i-spl and c-spl
#2nd deri of c-spl is 1st deri of i-spl
#------------------------------------------
makedmat_1D_ic <- function(m, knots, shpi = 9) {
q <- 1
dmat <- matrix(0, nrow = (m - q), ncol = m)
nkts <- length(knots)
if(shpi %in% c(9, 10)){
vals0 <- iSpline(knots, knots = knots[-c(1, nkts)], degree = 1, derivs = 1L, scale=FALSE)
vals0 <- diag(vals0)
if(shpi == 9){
vals <- vals0
}
if(shpi == 10){
vals <- -vals0
}
}
if(shpi %in% c(11, 12)){
vals0 <- cSpline(knots, knots = knots[-c(1, nkts)], degree = 1, derivs = 2L, scale=FALSE)
vals0 <- diag(vals0)
if(shpi == 11){
vals <- vals0
}
if(shpi == 12){
vals <- -vals0
}
}
if(shpi == 13) {
vals0 <- cSpline(knots, knots = knots[-c(1, nkts)], degree = 1, derivs = 2L, scale=FALSE)
vals0 <- diag(vals0)
vals <- c(vals0, 0) #the last column in delta is x
}
if(shpi == 14) {
vals0 <- cSpline(knots, knots = knots[-c(1, nkts)], degree = 1, derivs = 2L, scale=FALSE)
vals0 <- -diag(vals0)
vals <- c(vals0, 0)
}
if(shpi == 15) {
vals0 <- cSpline(knots, knots = knots[-c(1, nkts)], degree = 1, derivs = 2L, scale=FALSE)
vals0 <- diag(vals0)
vals <- c(vals0, 0)
}
if(shpi == 16) {
vals0 <- cSpline(knots, knots = knots[-c(1, nkts)], degree = 1, derivs = 2L, scale=FALSE)
vals0 <- -diag(vals0)
vals <- c(vals0, 0)
}
#vals0 <- diag(vals0)
# if (shpi == 9 | shpi == 11) {
# vals <- vals0
# } else if (shpi == 10 | shpi == 12) {
# vals <- -vals0
# } else if (shpi == 13 | shpi == 14) {
# vals <- c(vals0, 0)
# } else if (shpi == 15 | shpi == 16) {
# vals <- c(-vals0, 0)
# }
for (i in 2:m) {
dmat[i-1, i-1] <- -1 * vals[i-1]; dmat[i-1, i] <- 1 * vals[i]
}
return (dmat)
}
#--------------------------------------------------------------------
#forward step used in classo or stepwise
#change something like cic_vals to be with 10 rows: 1 for tree ordering
#--------------------------------------------------------------------
forward_cic <- function(xmat, y, cic_st, infocrit = "CAIC", type = "stepwise",
knots_lst = NULL, boundary.knots = NULL,
nsim=50, ddkeep_z_all=NULL, wt.iter=FALSE,
family = gaussian(), parallel=FALSE, space="E", pnt=FALSE,
a=1, maxedf = FALSE, add_edfu=0, add_kts=0, set4=FALSE)
{
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
muhat.fun <- cicfamily$muhat.fun
ysim.fun <- cicfamily$ysim.fun
deriv.fun <- cicfamily$deriv.fun
dev.fun <- cicfamily$dev.fun
capl <- NCOL(xmat)
n <- NROW(xmat)
#new: add tree like z: 0 out, 1 tree, 2 unordered
#cic_vals <- matrix(NA, nrow = 10, ncol = capl)
#dfmean_cbic_vals <- matrix(NA, nrow = 10, ncol = capl)
#no 18
cic_vals <- matrix(NA, nrow = 9, ncol = capl)
dfmean_cbic_vals <- matrix(NA, nrow = 9, ncol = capl)
one <- cbind(rep(1, n))
ddkeep_x <- NULL
ddkeep_z <- NULL
dmat_keep <- NULL
ps_keep <- 0
dfmean <- dfmeankeep <- dfmeanmaxkeep <- 0 #shpi = 17
penalty <- penaltykeep <- 0
ktskeep <- NULL
mskeep <- NULL
etahatkeep <- NULL
coefskeep <- NULL
oldcic <- cic_st
cic_st_next_step <- oldcic
#dfmean_keep_cbic <- 1e+3 #?
#new:
#dfmean_caic_vec <- NULL
dm <- ps <- NULL
xtxkeep <- face <- facekeep <- NULL
for(ix in seq_len(capl))
{
#print (ix)
x <- xmat[, ix]
#if(inherits(x, "factor") & attr(x, 'type') != 'tree'){
if(inherits(x, "factor") & attr(x, 'type') == 'nom'){
dd <- model.matrix(~ x)[, -1, drop=FALSE]
ms <- apply(dd, 2, mean)
#ms <- one %*% rbind(ms)
dd <- apply(dd, 2, function(e) e - mean(e))
#compute cic
#xm <- dd
xm <- cbind(one, dd) #add intercept
#new: update np
np_ix <- NCOL(xm)
#ybar is removed from y
#add ybar back
zvec <- y
if(!wt.iter){
bhat <- solve(crossprod(xm), t(xm)) %*% zvec
etahat <- xm %*% bhat
} else {
#m0 <- glm(y ~ factor(x), family=family$family)
m0 <- stats::glm(y ~ xm - 1, family = family) #xm includes intercept, needs to remove 1 from glm
#bhat <- coef(m0)[-1]
bhat <- coef(m0)
etahat <- m0$linear.predictors
muhat <- family$linkinv(etahat)
}
rss <- sum((y - etahat)^2)
edf0 <- 0
dfmean <- length(bhat) + edf0
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml = family$family)
#test!
} else {llh <- log(rss)}
if (infocrit == "CBIC") {
#cic_val <- llh + log(n) * (dfmean) / n
#new
#cic_val <- llh + a * log(n) * (dfmean) / n^(6/7)
degree <- 6/7#1
#degree <- 1
cic_val <- llh + a * log(n) * (dfmean) / n^degree
}
if (infocrit == "CAIC") {
cic_val <- llh + 2 * (dfmean) / n
}
#if this x reduces cic, then let its etahat to be kept; update cic_st to be cic
#the new dd will be added to null space; np_add the the # of new columns added
cic_vals[1, ix] <- oldcic #start with the cic from the yhat = ybat
if(cic_val < cic_st){
etahatkeep <- etahat
coefskeep <- bhat
ddkeep_z <- dd
ddkeep_x <- NULL
#dmat_keep <- NULL
mskeep <- ms
cic_st <- cic_val
cic_st_next_step_0 <- llh + 2 * (dfmean) / n
if(cic_st_next_step_0 < oldcic) {
cic_st_next_step <- cic_st_next_step_0
}
np_add <- NCOL(dd)
capm_add <- 0
cic_vals[1, ix] <- cic_val
}
#cic_vals[1, ix] <- cic_val
} else {
if(attr(x, 'type') %in% c('num', 'numeric')){
if(set4){ #new
shps <- 9:12
}else{
shps <- 9:16
}
}
if(attr(x, 'type') %in% c('ord', 'ordinal')){
if(set4){#new
shps <- 1:4
}else{
shps <- 1:8
}
}
if(attr(x, 'type') %in% c('tree')){ #1: tree order #2: unordered
shps <- 1:2
}
#new: if we use cbic and wt.iter then we only have 4 options
if(wt.iter & infocrit == "CBIC"){
if(attr(x, 'type') %in% c('num', 'numeric')){
shps <- 9:12
}
if(attr(x, 'type') %in% c('ord', 'ordinal')){
shps <- 1:4
}
if(attr(x, 'type') %in% c('tree')){
shps <- 1:2
}
}
#shps <- c(9:16, 18) #include 17 later
#shps <- 9:18 #17 is flat
#shps <- c(9:16)
nobs <- length(shps)
#loop through shapes and compute cic
gmat_lst <- vector("list", length = nobs)
np_ix_lst <- vector("list", length = nobs)
#mu0_lst <- vector("list", length = nobs)
ms_lst <- vector("list", length = nobs)
dd_lst <- vector("list", length = nobs)
#new: penalty matrix
dmat_lst <- vector("list", length = nobs)
dm_lst <- vector("list", length = nobs)
#new:
knots_lst <- vector("list", length = nobs)
ps_vec <- rep(0, nobs)
cic_vals_max <- cic_vals_edf0 <- NULL
etahat_lst <- vector("list", length = nobs)
bhat_lst <- vector("list", length = nobs)
xtx_lst <- vector("list", length = nobs)
face_lst <- vector("list", length = nobs)
#record gmat, np_ix, etc
#no 17
for(i in 1:nobs)
{
#shpi <- c(9:16)[i]
#shpi <- c(9:16, 18)[i]
shpi <- shps[i]
#print (c(ix, shpi))
if(shpi == 18) {
#dd <- cbind(x)
ms <- mean(x)
#check more!
dd <- cbind(x - ms)
gmat <- cbind(one, dd)
np_ix <- 2
#test more
dm <- NULL
ps <- 0
ps_vec[nobs] <- ps
#dd <- x
} else {
if(shpi <= 8){
if(attr(x, 'type') %in% c('ord', 'ordinal')){
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = 0, space = space)
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
ps <- 0
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
if(shpi > 2 & shpi < 5){
gmat <- cbind(one, x - mean(x), dd)
np_ix <- 2
} else {
gmat <- cbind(one, dd)
np_ix <- 1
}
}
#new
if(attr(x, 'type') %in% c('tree')){
#ms <- NULL
dd <- t(tree.fun(x))
ms <- apply(dd, 2, mean)
dd <- apply(dd, 2, function(e) e - mean(e))
gmat <- cbind(one, dd)
if(shpi == 1){
#dd <- t(tree.fun(x))
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
ps <- 0
#ms <- apply(dd, 2, mean)
#dd <- apply(dd, 2, function(e) e - mean(e))
#gmat <- cbind(one, dd)
np_ix <- 1
}
if(shpi == 2){
#dd <- model.matrix(~ x)[, -1, drop=FALSE]
#dd <- t(tree.fun(x))
#ms <- apply(dd, 2, mean)
#dd <- apply(dd, 2, function(e) e - mean(e))
#gmat <- cbind(one, dd)
dm <- NULL
ps <- 0
np_ix <- 1 + ncol(dd)
#ms <- apply(dd, 2, mean)
#dd <- apply(dd, 2, function(e) e - mean(e))
}
}
#ps <- 0
#dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
#dm_lst[[ix]] <- dm #no penalty
}
if(shpi > 8){
mult <- ifelse(pnt, 2, 1)
spl_degree_i <- ifelse(shpi %in% c(9, 10), 2, 3) #use 2 for i-spl, 3 for c-spl
n1 <- length(unique(x))
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
#new: see if we can fix choosing 13 - 16 too much
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/7)) + 6 + add_kts, trunc(n1^(1/9)) + 6 + add_kts)
edfu <- ifelse(spl_degree_i == 2, nkts/mult, nkts/mult + 1) + add_edfu
kts <- seq.int(min(x), max(x), length.out = nkts)
#new: boundary knots
kts0 <- boundary.knots[[ix]]
if(length(kts0) >= 2){
x_rng <- range(x, na.rm = TRUE)
kts0_rng <- range(kts0, na.rm = TRUE)
tol <- 1e-10
reset_lwr <- kts0_rng[1] > x_rng[1] + tol
reset_upp <- kts0_rng[2] < x_rng[2] - tol
if(reset_lwr){
kts_lwr <- x_rng[1]
} else {
kts_lwr <- kts0_rng[1]
}
if(reset_upp){
kts_upp <- x_rng[2]
} else {
kts_upp <- kts0_rng[2]
}
kts[1] <- kts_lwr
kts[nkts] <- kts_upp
}
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = kts, space = space)
#dd is orth to 1
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
#kts <- dd_ans$knots
if(shpi > 10 & shpi < 13){
gmat <- cbind(one, x - mean(x), dd)
np_ix <- 2
} else {
gmat <- cbind(one, dd)
np_ix <- 1
}
#print (kts)
#new:
if(pnt){
nc <- NCOL(dd)
dm <- makedmat_1D_ic(m = nc, knots = kts, shpi = shpi) #different from makedmat_1D for testpar
qvi <- crossprod(dd)
dvi <- crossprod(dm)
#ps <- 1e-3
#9, 10, 11, 12, 13, 14, 15, 16
if(shpi %in% c(9, 11, 13)){
log_root_res <- uniroot(f = function(log_pen,...) .search_ps_shapeselect(exp(log_pen),...),
interval = c(-8, 8),
qv0 = qvi, dv = dvi, dd = dd, edfu = edfu,
extendInt = "yes",
tol = .Machine$double.eps^0.25)
# convert back to the original scale
#temp
if(abs(log_root_res$root) > 8) {
ps <- floor(exp(8))
} else {
ps <- exp(log_root_res$root)
}
#print (ps)
#print (log_root_res$root)
#temp!
#if(abs(ps) < 1e-5) {ps <- 1e-5}
ps_vec[which(shps %in% shpi)] <- ps
} else if (shpi == 10){
ps <- ps_vec[1]
ps_vec[which(shps %in% shpi)] <- ps
} else if (shpi == 12) {
ps <- ps_vec[3]
ps_vec[which(shps %in% shpi)] <- ps
} else {
ps <- ps_vec[5]
ps_vec[which(shps %in% shpi)] <- ps
}
dm <- sqrt(ps) * dm
}
}
}
dd_lst[[i]] <- dd
if(!is.null(ms)){
ms_lst[[i]] <- ms
}
gmat_lst[[i]] <- gmat
np_ix_lst[[i]] <- np_ix
if(pnt & !is.null(dm)){
dmat_lst[[i]] <- dm
dm_lst[[i]] <- dm
if(shpi > 8){
knots_lst[[i]] <- kts
}
}
}
#compute dfmean_max, which is the average of all dfmean_max over nsim simulations
dfmean_rslt <- compute_edf(family=family, gmat_lst, dmat_lst=dmat_lst,
mu0_lst=NULL, np_ix_lst, weights=NULL,
nsim=nsim, wt.iter=wt.iter, parallel=parallel,
pnt=pnt)
dfmean_caic <- dfmean_rslt$dfmean_caic
dfmean_cbic <- dfmean_rslt$dfmean_cbic
#new
#if(nobs > 4){
#dfmean_cbic_vals[-9, ix] <- dfmean_cbic
#}
#if(nobs == 4){
dfmean_cbic_vals[1:nobs, ix] <- dfmean_cbic
#}
#new:
#dfmean_caic_vec <- c(dfmean_caic_vec, dfmean_caic)
#go through all shapes to compute rss and cic/bic
#cic_vals[9, ix] <- oldcic #start with yhat = ybar
for(shpi in shps){
#ps_shpi <- which(c(9:16, 18) %in% shpi)
#ps_shpi <- which(c(9:16) %in% shpi)
ps_shpi <- which(shps %in% shpi)
gmat <- gmat_lst[[ps_shpi]]
ms <- ms_lst[[ps_shpi]]
dd <- dd_lst[[ps_shpi]]
np_ix <- np_ix_lst[[ps_shpi]]
#new
# if(shpi <= 8){
# pnt <- FALSE
# }
if(pnt){
dm <- dmat_lst[[ps_shpi]]
ps <- ps_vec[ps_shpi]
}
nc = ncol(gmat)
amat = diag(nc)
ps_zero = 1:np_ix
amat[ps_zero, ps_zero] = 0
bvec = cbind(rep(0, nc))
xtx = crossprod(gmat)
if(pnt & !is.null(dm)){
dm_zero <- matrix(0, nrow = nrow(dm), ncol = np_ix)
dmat <- cbind(dm_zero, dm)
dv <- crossprod(dmat) #absort penalty already
xtx <- xtx + dv
}
if (shpi == 18 | shpi == 2 & attr(x, "type") == 'tree') {
zvec <- y
if(!wt.iter){
bhat <- solve(crossprod(gmat), t(gmat)) %*% zvec
etahat <- gmat %*% bhat
} else {
#gmat has one
m0 <- stats::glm(y ~ gmat - 1, family=family)
bhat <- coef(m0)
etahat <- m0$linear.predictors
}
} else {
zvec <- y
if(!wt.iter){
if (!pnt | pnt & all(dm == 0)) {
dsend = gmat[, -c(1:np_ix), drop = FALSE]
zsend = gmat[, 1:np_ix, drop = FALSE]
ans = coneB(zvec, dsend, zsend)
etahat = ans$yhat
bhat = ans$coefs
face = ans$face
} else {
#ans <- coneB(zvec, dsend, zsend)
#etahat <- ans$yhat
#bhat <- ans$coefs
xty = crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
#ans = solve.QP(xtx, xty, t(amat), bvec)
ans = solve.QP(xtx, xty, amat, bvec)
bhat = ans$solution
bhat = round(bhat, 10)
etahat = gmat %*% bhat
face = ans$iact
}
} else {
gmat0 <- gmat
weights <- rep(1, n)
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
muhat <- family$linkinv(etahat)
diff <- 1
niter <- 0
sm <- 1e-5
#sm <- 1e-7
#mdiff <- TRUE
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
while (diff > sm & mdiff & niter < 10)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
#test
if(!pnt | pnt & all(dm == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ansi <- coneB(zvec, dsend, zsend)
bhat <- ansi$coefs
}
if(pnt & !is.null(dm) & !all(dm == 0)){
xtx = crossprod(gmat)
xtx = xtx + dv
xty = crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ansi = solve.QP(xtx, xty, amat, bvec)
#ansi = solve.QP(xtx, xty, t(amat), bvec)
bhat = ansi$solution
face = ansi$iact
}
bhat = round(bhat, 10)
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
#mdiff <- abs(max(muhat) - 1) > sm
}
}
}
rss <- sum((y - etahat)^2)
#new
etahat_lst[[ps_shpi]] <- etahat
bhat_lst[[ps_shpi]] <- bhat
xtx_lst[[ps_shpi]] <- xtx
face_lst[[ps_shpi]] <- face
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml=family$family)
} else {llh <- log(rss)}
if (infocrit == "CBIC") {
#cic_val <- llh + log(n) * (dfmean_caic) / n
#cic_val_edf0 <- llh + log(n) * (dfmean_cbic[ps_shpi]) / n
#new
#cic_val <- llh + a * log(n) * (dfmean_caic) / n^(6/7)
#cic_val_edf0 <- llh + a * log(n) * (dfmean_cbic[ps_shpi]) / n^(6/7)
if(shpi %in% c(9, 10)){
degree <- 6/7
} else if (shpi > 10 & shpi < 17){
degree <- 6/7#8/9 # 6/7#8/9
} else if (shpi > 0 & shpi < 9) {
degree <- 6/7#1 # 6/7#1
} else {
stop('check!')
}
cic_val <- llh + a * log(n) * (dfmean_caic) / n^(degree)
cic_val_edf0 <- llh + a * log(n) * (dfmean_cbic[ps_shpi]) / n^(degree)
}
if (infocrit == "CAIC") {
cic_val <- llh + 2 * (dfmean_caic) / n
cic_val_edf0 <- llh + 2 * (dfmean_cbic[ps_shpi]) / n
}
#new:
if(maxedf){
cic_vals_max <- c(cic_vals_max, cic_val)
} else {
cic_vals_max <- c(cic_vals_max, cic_val_edf0)
}
cic_vals_edf0 <- c(cic_vals_edf0, cic_val_edf0)
#test!
if(shpi > 8){
cic_vals[(shpi - 8), ix] <- cic_val
}
if(shpi <= 8){
cic_vals[shpi, ix] <- cic_val
#cic_vals[which(shps %in% shpi), ix] <- cic_val
}
#cic_vals[(shpi - 8), ix] <- cic_val_edf0
}
cic_vals[9, ix] <- oldcic #flat
#cic_vals[(1+nobs), ix] <- oldcic #flat
#new:
if(any(cic_vals_max < cic_st)){
#test!
#cic_vals[-9, ix] <- cic_vals_edf0 #no flat
if (wt.iter & infocrit == "CBIC") {
cic_vals[1:length(cic_vals_edf0), ix] <- cic_vals_edf0
} else {
cic_vals[1:length(cic_vals_edf0), ix] <- cic_vals_edf0
#cic_vals[-9, ix] <- cic_vals_edf0 #no flat, will not work when it is tree
#cic_vals[-(1+nobs), ix] <- cic_vals_edf0 #no flat
}
ps_shpi <- which(cic_vals_edf0 == min(cic_vals_edf0))
if(length(ps_shpi) > 1) {
ps_shpi <- rev(ps_shpi)[1]
}
etahatkeep <- etahat_lst[[ps_shpi]] #etahat #will include 1 for binomial; but etacomps will not use this
coefskeep <- bhat_lst[[ps_shpi]] # bhat
mskeep <- ms_lst[[ps_shpi]] # ms
ktskeep <- knots_lst[[ps_shpi]] # kts; will be NULL for ordinal
dd <- dd_lst[[ps_shpi]]
dm <- dmat_lst[[ps_shpi]]
ps <- ps_vec[ps_shpi]
shpi <- shps[ps_shpi]
if(shpi != 18) {
if(shpi == 2 & attr(x, 'type') == 'tree'){
ddkeep_x <- NULL
ddkeep_z <- dd
} else {
ddkeep_x <- dd
}
if(pnt){
dmat_keep <- dm
ps_keep <- ps
}
}
if(shpi != 2 | attr(x, 'type') != 'tree'){
ddkeep_z <- NULL
np_add <- 0
if (shpi > 10 & shpi < 13 | shpi == 18 | shpi > 2 & shpi < 5){
ddkeep_z <- cbind(x - mean(x))
np_add <- 1
}
}
cic_st_next_step_0 <- cic_vals_edf0[ps_shpi]
cic_st_next_step <- cic_st_next_step_0
#test
cic_st <- cic_st_next_step_0
#xtxkeep <- xtx_lst[[ps_shpi]]
#facekeep <- face_lst[[ps_shpi]]
}
}
}
rslt <- list(cic_vals = cic_vals, ddkeep_x = ddkeep_x, ddkeep_z = ddkeep_z,
dmat_keep = dmat_keep, ps_keep = ps_keep, mskeep = mskeep, ktskeep = ktskeep,
etahatkeep = etahatkeep, coefskeep = coefskeep, dfmean_cbic_vals = dfmean_cbic_vals,
cic_st = cic_st, cic_st_next_step = cic_st_next_step)
return (rslt)
}
#--------------------------------------------------------------------------------------------------------------------------------
forward_cic_step2 <- function(xmat_inactive, y, cic_st, ddkeep_z_all = NULL, ddkeep_x_all = NULL, knots_lst_inactive = NULL,
boundary.knots_inactive = NULL,
oldcic = NULL, oldetahat = NULL, oldcoefs = NULL, infocrit = "CAIC", dfmean_old = NULL,
dfmeanmax_old = NULL, penalty_old = NULL,
active, inactive, type = "stepwise", nsim=50, wt.iter=FALSE,
family = gaussian(), parallel=FALSE, space="E", pnt=FALSE, dmat_keep_lst=NULL,
ps_keep_vec=NULL, a=1, maxedf = FALSE,
add_edfu = 0, add_kts = 0, shps_non_zero, seq_added, set4 = FALSE) {
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
muhat.fun <- cicfamily$muhat.fun
ysim.fun <- cicfamily$ysim.fun
deriv.fun <- cicfamily$deriv.fun
dev.fun <- cicfamily$dev.fun
capl <- NCOL(xmat_inactive)
n <- NROW(xmat_inactive)
cic_st_0 <- cic_st #largetest cic_st before putting in any new x; cic val for 17 in classo
cic_st_next_step <- oldcic #for stepwise, cic_st and oldcic have the same starting val; cic_st will be replaced when decreased
#new: add tree
#cic_vals <- matrix(NA, nrow = 10, ncol = capl)
#dfmean_cbic_vals <- matrix(NA, nrow = 10, ncol = capl)
#no 18
cic_vals <- matrix(NA, nrow = 9, ncol = capl)
dfmean_cbic_vals <- matrix(NA, nrow = 9, ncol = capl)
etahatkeep <- oldetahat
coefskeep <- oldcoefs
ddkeep_z <- NULL
ddkeep_x <- NULL
ktskeep <- NULL
mskeep <- NULL
np_add <- 0
rm_id <- NULL
one <- cbind(rep(1, n))
#new:
dmat_keep <- dm <- dmat <- dmat_keep_all <- NULL
ps_keep <- ps <- 0
if(pnt){
dmat_keep_all <- Matrix::bdiag(Filter(Negate(is.null), dmat_keep_lst[seq_added]))
#dmat_keep_all <- Matrix::bdiag(dmat_keep_lst[seq_added])
dmat_keep_all <- as.matrix(dmat_keep_all)
if(prod(dim(dmat_keep_all)) == 0) {
dmat_keep_all <- NULL
}
}
for(ix in seq_len(capl)){
#print (ix)
x <- xmat_inactive[, ix]
#for categorical variable: oldcic is for 0; cic_val is for 1. they need to be compared with each other
if(inherits(x, "factor") & attr(x, 'type') == 'nom'){
#dd is orth to 1
dd <- model.matrix(~ x)[, -1, drop=FALSE]
#new:
#if(NCOL(ddkeep_z_all) == 0){
ms <- apply(dd, 2, mean)
dd <- apply(dd, 2, function(e) e - mean(e))
#}
#not sure how to do this in predict.classo....
# if(NCOL(ddkeep_z_all) > 0){
# dz_temp <- cbind(1, ddkeep_z_all)
# pmz <- dz_temp %*% MASS::ginv(crossprod(dz_temp)) %*% t(dz_temp)
# ms <- pmz %*% dd
# #coefms <- solve(crossprod(dz_temp)) %*% t(dz_temp) %*% dd
# dd <- dd - ms
# nosignal <- apply(dd, 2, function(e) all(round(e, 8) < 1e-8))
# if(all(nosignal)){
# rm_id <- c(rm_id, ix)
# cic_vals[1, ix] <- oldcic
# next
# } else if(any(!nosignal)){
# kept_id <- which(!nosignal)
# dd <- dd[, kept_id, drop = FALSE]
# ms <- ms[, kept_id, drop = FALSE]
# }
# }
#new: some column in dd may be the same as some column in ddkeep_z_all
#skip z if all columns in z are already in the space of ddkeep_z_all
# if(NCOL(ddkeep_z_all) > 0){
# cat('ms changed!', '\n')
# #pmz = ddkeep_z_all %*% solve(crossprod(ddkeep_z_all), t(ddkeep_z_all))
# pmz <- ddkeep_z_all %*% MASS::ginv(crossprod(ddkeep_z_all)) %*% t(ddkeep_z_all)
# ms_add <- pmz %*% dd
# dd <- dd - ms_add
# #test!
# ms_addt <- t(ms_add)
# ms <- ms_addt + ms
# ms <- t(ms)
# sgn_dd <- apply(dd, 2, function(e) all(round(e, 8) < 1e-8))
# if(all(sgn_dd)) {
# rm_id <- c(rm_id, ix)
# #new: make sure cic_vals is not all NA
# cic_vals[1, ix] <- oldcic
# next
# } else if(any(!sgn_dd)){
# #cat('ms changed!', '\n')
# kept_id <- which(!sgn_dd)
# dd <- dd[, kept_id, drop = FALSE]
# #ms <- one %*% rbind(ms) + ms_add
# #test!
# ms <- ms[, kept_id, drop = FALSE]
# }
# }
xm <- cbind(one, ddkeep_z_all, dd, ddkeep_x_all)
gmat <- xm
#zvec <- y - etahat_all #residual vector
zvec <- y
capm <- NCOL(ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + NCOL(dd) + 1
dmat <- dmat_keep_all
#vm <- gmat[, 1:np_ix]
if(capm == 0){
#no x has been chosen
edf0 = 0
dfmean = edf0 + ncol(xm)
dfmean_caic <- dfmean
dfmean_cbic <- dfmean
} else {
dfmean_rslt <- compute_edf(family = family, gmat_lst = list(gmat), dmat_lst=list(dmat),
mu0_lst = NULL,
np_ix_lst = list(np_ix), weights=NULL, nsim=nsim, wt.iter=wt.iter,
parallel=parallel, pnt=pnt) #oldcic is without z (edf is the average of the chosen shape without max)
dfmean_caic <- dfmean_rslt$dfmean_caic #only has 1 val; caic and cbic are the same
dfmean_cbic <- dfmean_rslt$dfmean_cbic #only has 1 val
}
#new:
#dfmean_caic_vec <- c(dfmean_caic_vec, dfmean_caic)
if(capm == 0){
if(!wt.iter){
bhat = solve(crossprod(xm), t(xm)) %*% zvec
etahat <- xm %*% bhat
} else {
m0 = stats::glm(y ~ xm - 1, family=family)
bhat = coef(m0)
etahat <- m0$linear.predictors
muhat <- family$linkinv(etahat)
}
}
if(capm > 0){ #no x
#gmat <- xm
nc = ncol(gmat)
amat = diag(nc)
ps_zero = 1:np_ix
amat[ps_zero, ps_zero] = 0
#amat[1:np_ix, 1:np_ix] = 0
bvec = cbind(rep(0, nc))
xtx = crossprod(gmat)
#test! avoid linear, shp = 18
if(pnt & !is.null(dmat)){
dm_lst = vector("list", length = 2)
dm0 = matrix(0, nrow = nrow(dmat), ncol = np_ix)
dm_lst[[1]] = dm0
dm_lst[[2]] = dmat
dmat = as.matrix(bdiag(dm_lst))
dv = crossprod(dmat)
xtx = xtx + dv
}
if(!wt.iter){
#test
if(!pnt | pnt & all(dmat == 0)){
dsend = gmat[, -c(1:np_ix), drop = FALSE]
zsend = gmat[, 1:np_ix, drop = FALSE]
ans = coneB(zvec, dsend, zsend)
bhat = ans$coefs
bhat = round(bhat, 10)
etahat = gmat %*% bhat
face = ans$face
} else {
xty = crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ans = solve.QP(xtx, xty, amat, bvec)
#ans = solve.QP(xtx, xty, t(amat), bvec)
bhat = ans$solution
bhat = round(bhat, 10)
#cat('finished bhat in check z:', '\n')
etahat = gmat %*% bhat
face = ans$iact
}
} else {
gmat0 <- gmat
weights <- rep(1, n)
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
muhat <- family$linkinv(etahat)
sm <- 1e-5
#sm <- 1e-7 #will cause solve error
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 10)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
#test
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ansi <- coneB(zvec, dsend, zsend)
bhat <- ansi$coefs
}
if(pnt & !is.null(dmat) & !all(dmat == 0)){
xtx = crossprod(gmat)
xtx = xtx + dv
xty = crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ansi = solve.QP(xtx, xty, amat, bvec)
#ansi = solve.QP(xtx, xty, t(amat), bvec)
bhat = ansi$solution
face = ansi$iact
}
bhat = round(bhat, 10)
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
}
#check
#vm <- cbind(one, ddkeep_z_all, dd)
rss <- sum((y - etahat)^2)
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml=family$family)
#test!
#llh <- 2 * (llh * n / 2) #2 times neg-loglike
} else {llh <- log(rss)}
#include 1 when computing cic
if (infocrit == "CBIC") {
#cic_val <- llh + log(n) * (dfmean_cbic) / n
#new:
#cic_val <- llh + a * log(n) * (dfmean_cbic) / n^(6/7)
all_shps <- c(shps_non_zero, 1)
if(any(all_shps %in% c(9, 10))){ #smooth incr/decr
degree <- 6/7
} else if (any(all_shps > 10 & all_shps < 17)){ #smooth conv/conc
degree <- 6/7#8/9 #
} else if (any(all_shps > 0 & all_shps < 9)) { #ord or parametric
degree <- 6/7#1 #
} else {
stop('check!')
}
cic_val <- llh + a * log(n) * (dfmean_cbic) / n^degree
}
if (infocrit == "CAIC") {
cic_val <- llh + 2 * (dfmean_caic) / n
}
#start with shape 0; oldcic
#cic_vals[1, ix] <- cic_st_0
cic_vals[1, ix] <- oldcic
if(cic_val < cic_st){
etahatkeep <- etahat
coefskeep <- bhat
ddkeep_z <- dd
mskeep <- ms
ddkeep_x <- NULL
dmat_keep <- NULL
ps_keep <- 0
ktskeep <- NULL
cic_st <- cic_val
cic_st_next_step_0 <- cic_val
cic_st_next_step <- cic_st_next_step_0
np_add <- NCOL(ddkeep_z)
cic_vals[1, ix] <- cic_val
#xtxkeep <- xtx
#facekeep <- face
}
#cat(ix, 'z', cic_val, '\n')
} else {
#test more
#cic_st <- oldcic
#shps <- 9:18 #17 is flat; 17 is included
if(attr(x, 'type') %in% c('num', 'numeric')){
#new
if(set4){
shps <- c(9:12, 17)
}else{
shps <- 9:17
}
#shps <- 9:17
}
if(attr(x, 'type') %in% c('ord', 'ordinal')){
#new
if(set4){
shps <- c(1:4, 0)
}else{
shps <- c(1:8, 0)
}
#shps <- c(1:8, 0)
}
if(attr(x, 'type') %in% c('tree')){
shps <- c(1, 2, 0)
}
#new: if we use cbic and wt.iter then we only have 4 options
if(wt.iter & infocrit == "CBIC"){
if(attr(x, 'type') %in% c('num', 'numeric')){
shps <- c(9:12, 17)
}
if(attr(x, 'type') %in% c('ord', 'ordinal')){
shps <- c(1:4, 0)
}
if(attr(x, 'type') %in% c('tree')){
shps <- c(1, 2, 0)
}
}
nobs <- length(shps)
gmat_lst <- vector("list", length = nobs)
np_ix_lst <- vector("list", length = nobs)
#mu0_lst <- vector("list", length = nobs)
ms_lst <- vector("list", length = nobs)
dd_lst <- vector("list", length = nobs)
kts_lst <- vector("list", length = nobs)
dmat_lst <- vector("list", length = nobs)
dm_lst <- vector("list", length = nobs)
ps_vec <- rep(0, length = nobs)
#dmat <- dmat_keep_all
#record gmat, np_ix, etc
for(i in 1:nobs)
{
shpi <- shps[i]
#print (c(ix, shpi))
if (shpi == 17 | shpi == 0) {
dd <- NULL
ms <- NULL
kts <- NULL
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + 1
dm <- NULL
ps <- 0
dmat <- dmat_keep_all
cic_val_edf0 <- oldcic
cic_val <- oldcic
} else if (shpi == 18) {
#dd <- cbind(x)
dd <- cbind(x - mean(x))
kts <- NULL
ms <- mean(x)
gmat <- cbind(one, ddkeep_z_all, dd, ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + 2
#dmat_lst[[i]] <- 0
#not_skip_cic <- TRUE
dm <- NULL
ps <- 0
dmat <- dmat_keep_all
} else {
if(shpi <= 8) {
if(attr(x, 'type') %in% c('ord')){
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = 0, space = space)
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
if(shpi > 2 & shpi < 5){
gmat <- cbind(one, ddkeep_z_all, x-mean(x), ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 2
} else {
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 1
}
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
}
if(attr(x, 'type') %in% c('tree')){
#ms <- NULL
dd <- t(tree.fun(x))
ms <- apply(dd, 2, mean)
dd <- apply(dd, 2, function(e) e - mean(e))
if(shpi == 1){
#dd <- t(tree.fun(x))
#ms <- apply(dd, 2, mean)
#dd <- apply(dd, 2, function(e) e - mean(e))
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 1
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
}
if(shpi == 2){
#dd <- model.matrix(~ x)[, -1, drop=FALSE]
#dd <- t(tree.fun(x))
#ms <- apply(dd, 2, mean)
#dd <- apply(dd, 2, function(e) e - mean(e))
gmat <- cbind(one, ddkeep_z_all, dd, ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + 1 + NCOL(dd)
dm <- NULL
}
}
ps <- 0
dmat <- dmat_keep_all
if(!is.null(dmat)){
dm_lst_tmp <- vector('list', length = 2)
dm_lst_tmp[[1]] <- dmat
dm_lst_tmp[[2]] <- dm
dmat <- as.matrix(Matrix::bdiag(dm_lst_tmp))
}else{
dmat <- dm
}
}
if(shpi > 8) {
mult <- ifelse(pnt, 2, 1)
spl_degree_i <- ifelse(shpi %in% c(9, 10), 2, 3) #use 2 for i-spl, 3 for c-spl
n1 <- length(unique(x))
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
#new
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/7)) + 6 + add_kts, trunc(n1^(1/9)) + 6 + add_kts)
edfu <- ifelse(spl_degree_i == 2, nkts/mult, nkts/mult + 1) + add_edfu
kts <- seq.int(min(x), max(x), length.out = nkts)
#new: boundary knots
kts0 <- boundary.knots_inactive[[ix]]
if(length(kts0) >= 2){
x_rng <- range(x, na.rm = TRUE)
kts0_rng <- range(kts0, na.rm = TRUE)
tol <- 1e-10
reset_lwr <- kts0_rng[1] > x_rng[1] + tol
reset_upp <- kts0_rng[2] < x_rng[2] - tol
if(reset_lwr){
kts_lwr <- x_rng[1]
} else {
kts_lwr <- kts0_rng[1]
}
if(reset_upp){
kts_upp <- x_rng[2]
} else {
kts_upp <- kts0_rng[2]
}
kts[1] <- kts_lwr
kts[nkts] <- kts_upp
}
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = kts, space = space)
#dd is orth to 1
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
#gmat <- cbind(zmat, ddkeep_z_all, ddkeep_x_all, dd)
if(shpi > 10 & shpi < 13){
gmat <- cbind(one, ddkeep_z_all, x-mean(x), ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 2
} else {
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 1
}
if(pnt){
nc <- NCOL(dd)
dm <- makedmat_1D_ic(m = nc, knots = kts, shpi = shpi) #different from makedmat_1D for testpar
qvi <- crossprod(dd)
dvi <- crossprod(dm)
#ps <- 1e-6
#9, 10, 11, 12, 13, 14, 15, 16
if(shpi %in% c(9, 11, 13)){
log_root_res <- uniroot(f = function(log_pen,...) .search_ps_shapeselect(exp(log_pen),...),
interval = c(-8, 8),
qv0 = qvi, dv = dvi, dd = dd, edfu = edfu,
extendInt = "yes",
tol = .Machine$double.eps^0.25)
# convert back to the original scale
#temp
if(abs(log_root_res$root) > 8) {
ps <- floor(exp(8))
} else {
ps <- exp(log_root_res$root)
}
#print (ps)
#ps <- exp(log_root_res$root)
ps_vec[which(shps %in% shpi)] <- ps
} else if (shpi == 10){
ps <- ps_vec[1]
ps_vec[which(shps %in% shpi)] <- ps
} else if (shpi == 12) {
ps <- ps_vec[3]
ps_vec[which(shps %in% shpi)] <- ps
} else {
ps <- ps_vec[5]
ps_vec[which(shps %in% shpi)] <- ps
}
dm <- sqrt(ps) * dm
dmat <- dmat_keep_all
if(!is.null(dmat)){
dm_lst_tmp <- vector('list', length = 2)
dm_lst_tmp[[1]] <- dmat
dm_lst_tmp[[2]] <- dm
dmat <- as.matrix(bdiag(dm_lst_tmp))
}else{
dmat <- dm
}
}
}
}
if(shpi == 17 | shpi == 0) {
dd_lst[i] <- list(NULL)
ms_lst[i] <- list(NULL)
kts_lst[i] <- list(NULL)
#dmat_lst[[i]] <- dmat
} else {
dd_lst[[i]] <- dd
if(!is.null(ms)){
ms_lst[[i]] <- ms
}
#dmat_lst[[i]] <- dmat
if(shpi == 18){
kts_lst[i] <- list(NULL)
} else {
#will be NULL for ordinal
if(shpi > 8){
kts_lst[[i]] <- kts
}
}
}
if(!is.null(dmat)){
dmat_lst[[i]] <- dmat #diff from forward_cic_step
}
if(pnt & !is.null(dm)){
dm_lst[[i]] <- dm
}
gmat_lst[[i]] <- gmat
np_ix_lst[[i]] <- np_ix
}
#new: compute llh first
#compute dfmean_max, which is the average of all dfmean_max over nsim simulations
#need to remove 17 when compute dfmean
#ps_17 <- which(shps %in% c(17, 0))
#dfmean_rslt <- compute_edf(family=family, gmat_lst=gmat_lst[-ps_17], dmat_lst=dmat_lst[-ps_17],
# mu0_lst=NULL, np_ix_lst=np_ix_lst[-ps_17],
# weights=NULL, nsim=nsim, wt.iter=wt.iter, parallel=parallel, pnt=pnt)
#dfmean_caic <- dfmean_rslt$dfmean_caic
#dfmean_cbic <- dfmean_rslt$dfmean_cbic
#dfmean_cbic_vals[-9, ix] <- dfmean_cbic
#new
# if(nobs > 5){
# dfmean_cbic_vals[-9, ix] <- dfmean_cbic
# }
#if(nobs == 5){
#no 17
#dfmean_cbic_vals[1:length(dfmean_cbic), ix] <- dfmean_cbic
#}
#test!
#dfmean_cbic_vals[1:length(dfmean_cbic), ix] <- dfmean_cbic
#new: compute the cic lower bound for each x first
# which uses the total number of predictors as dfmean
#dfmean_caic_vec <- c(dfmean_caic_vec, dfmean_caic)
#go through all shapes to compute rss and cic/bic
cic_vals_max <- cic_vals_edf0 <- NULL # rep(cic_st_0, nobs)
llh_vec <- NULL
etahat_lst <- vector("list", length = nobs)
bhat_lst <- vector("list", length = nobs)
face_lst <- vector("list", length = nobs)
xtx_lst <- vector("list", length = nobs)
for(shpi in shps){
#print (shpi)
#ps_shpi <- which(c(9:18) %in% shpi)
#ps_shpi <- which(c(9:17) %in% shpi)
ps_shpi <- which(shps %in% shpi)
gmat <- gmat_lst[[ps_shpi]]
ms <- ms_lst[[ps_shpi]]
dd <- dd_lst[[ps_shpi]]
dm <- dm_lst[[ps_shpi]]
np_ix <- np_ix_lst[[ps_shpi]]
kts <- kts_lst[[ps_shpi]]
if(pnt){
dmat <- dmat_lst[[ps_shpi]]
ps <- ps_vec[ps_shpi]
}
nc <- ncol(gmat)
amat <- diag(nc)
ps_zero <- 1:np_ix
amat[ps_zero, ps_zero] <- 0
bvec <- cbind(rep(0, nc))
xtx <- crossprod(gmat)
if(pnt & !is.null(dmat)){
dm_zero <- matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat <- cbind(dm_zero, dmat)
dv <- crossprod(dmat) #absort penalty already
xtx <- xtx + dv
}
if (shpi == 17 | shpi == 0) {
#etahat <- rep(mean(y), n)
etahat <- oldetahat
#rss <- sum((y - etahat)^2)
#cic_val <- oldcic
cic_val <- cic_st_0
cic_val_edf0 <- cic_st_0
bhat <- oldcoefs
#print (length(bhat))
dd <- NULL
ms <- NULL
#dm <- NULL
} else {
if(np_ix > 0 & NCOL(gmat) > np_ix){
zvec <- y
if(!wt.iter){
#test
if(!pnt | pnt & all(dmat == 0)){
dsend = gmat[, -c(1:np_ix), drop = FALSE]
zsend = gmat[, 1:np_ix, drop = FALSE]
ans = coneB(zvec, dsend, zsend)
bhat = ans$coefs
bhat = round(bhat, 10)
etahat = gmat %*% bhat
face = ans$face
} else {
xty = crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
#ans = solve.QP(xtx, xty, t(amat), bvec)
ans = solve.QP(xtx, xty, amat, bvec)
bhat = ans$solution
bhat = round(bhat, 10)
etahat = gmat %*% bhat
face = ans$iact
}
} else {
gmat0 <- gmat
weights <- rep(1, n)
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
muhat <- family$linkinv(etahat)
sm <- 1e-5
#sm <- 1e-7
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 10)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
#test
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ansi <- coneB(zvec, dsend, zsend)
bhat <- ansi$coefs
}
if(pnt & !is.null(dmat) & !all(dmat == 0)){
xtx = crossprod(gmat)
xtx = xtx + dv
xty = crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ansi = solve.QP(xtx, xty, amat, bvec)
#ansi = solve.QP(xtx, xty, t(amat), bvec)
bhat = ansi$solution
face = ansi$iact
}
bhat = round(bhat, 10)
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
} else if (np_ix > 0 & NCOL(gmat) == np_ix) { #shpi = 18 and no constrained x has been chosen
#zvec <- y
#bhat <- solve(crossprod(gmat), t(gmat)) %*% zvec
#etahat <- gmat %*% bhat
zvec <- y
if(!wt.iter){
#bhat <- solve(crossprod(gmat), t(gmat)) %*% zvec
bhat <- ginv(crossprod(gmat)) %*% t(gmat) %*% zvec
etahat <- gmat %*% bhat
} else {
m0 <- stats::glm(y ~ gmat - 1, family = family)
etahat <- m0$linear.predictors
bhat <- coef(m0)
muhat <- family$linkinv(etahat)
}
#etahat <- gmat %*% solve(crossprod(gmat), t(gmat)) %*% zvec
}
rss <- sum((y - etahat)^2)
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml = family$family)
#test!
#llh <- 2 * (llh * n / 2) #2 times neg-loglike
} else {llh <- log(rss)}
#new
llh_vec <- c(llh_vec, llh)
if(shpi != 17 & shpi != 0){ #remove this if?
#ps_shpi_bic = which(c(9:16, 18) %in% shpi)
if(shpi > 8){
ps_shpi_bic = which(c(9:16) %in% shpi)
}
if(shpi <= 8){
ps_shpi_bic = which(c(1:8) %in% shpi)
}
if (infocrit == "CBIC") {
#cic_val <- llh + log(n) * (dfmean_caic) / n
#cic_val_edf0 <- llh + log(n) * (dfmean_cbic[ps_shpi_bic]) / n
#new:
#cic_val <- llh + a * log(n) * (dfmean_caic) / n^(6/7)
#cic_val_edf0 <- llh + a * log(n) * (dfmean_cbic[ps_shpi_bic]) / n^(6/7)
all_shps <- c(shps_non_zero, shpi)
if(any(all_shps %in% c(9, 10))){ #smooth incr/decr
degree <- 6/7#
} else if (any(all_shps > 10 & all_shps < 17)){ #smooth conv/conc
degree <- 6/7#8/9
} else if (any(all_shps > 0 & all_shps < 9)) { #ord or parametric
degree <- 6/7#1
} else {
stop('check!')
}
#new
dfmean_lwr <- length(active) + 1
cic_val <- llh + a * log(n) * (dfmean_lwr) / n^(degree)
cic_val_edf0 <- llh + a * log(n) * (dfmean_lwr) / n^(degree)
}
if (infocrit == "CAIC") {
#new
dfmean_lwr <- length(active) + 1
cic_val <- llh + 2 * (dfmean_lwr) / n
cic_val_edf0 <- llh + 2 * (dfmean_lwr) / n
}
}
}
#cat('include x', '\n')
if(shpi > 8){
#cic_vals[(shpi - 8), ix] <- cic_val
#new: use the starting cic, bc we use lower bound for edf0 now, cic_vals could be much smaller than they should
cic_vals[(shpi - 8), ix] <- cic_st_0
}
if(shpi <= 8) {
if(shpi == 0) {
#cic_vals[9, ix] <- cic_val
#new
cic_vals[9, ix] <- cic_st_0
#cic_vals[nobs, ix] <- cic_val
} else {
#cic_vals[ps_shpi, ix] <- cic_val
#cic_vals[shpi, ix] <- cic_val
#new
cic_vals[shpi, ix] <- cic_st_0
}
}
#print (ps_shpi)
etahat_lst[[ps_shpi]] <- etahat
bhat_lst[[ps_shpi]] <- bhat
#face_lst[[ps_shpi]] <- face
#xtx_lst[[ps_shpi]] <- xtx
#cic_vals_max <- c(cic_vals_max, cic_val)
if(maxedf){ #not use anymore
cic_vals_max <- c(cic_vals_max, cic_val)
} else {
cic_vals_max <- c(cic_vals_max, cic_val_edf0)
}
cic_vals_edf0 <- c(cic_vals_edf0, cic_val_edf0)
#cic_vals[(shpi - 8), ix] <- cic_val_edf0
}
#new:
ps_to_sims <- NULL
dfmean_cbic_vals[1:(length(shps)-1), ix] <- length(active) + 1
shps_to_sims <- shps[-length(shps)]
if(any(cic_vals_max < cic_st)){
ps_shps_tosim <- which(cic_vals_max < cic_st)
shps_to_sims <- shps_to_sims[ps_shps_tosim]
ps_17 <- which(shps %in% c(17, 0))
#no need to worry about ps_17, because it will be the last one?
dfmean_rslt <- compute_edf(family=family, gmat_lst=gmat_lst[ps_shps_tosim], dmat_lst=dmat_lst[ps_shps_tosim],
mu0_lst=NULL, np_ix_lst=np_ix_lst[ps_shps_tosim],
weights=NULL, nsim=nsim, wt.iter=wt.iter,
parallel=parallel, pnt=pnt, shps_to_sims=shps_to_sims)
dfmean_caic <- dfmean_rslt$dfmean_caic
dfmean_cbic <- dfmean_rslt$dfmean_cbic
dfmean_cbic_vals[ps_shps_tosim, ix] <- dfmean_cbic
if (infocrit == "CBIC") {
cic_vals_max[ps_shps_tosim] <- llh_vec[ps_shps_tosim] + a * log(n) * (dfmean_cbic) / n^(degree)
cic_vals_edf0[ps_shps_tosim] <- llh_vec[ps_shps_tosim] + a * log(n) * (dfmean_cbic) / n^(degree)
}
if (infocrit == "CAIC") {
#new
cic_vals_max[ps_shps_tosim] <- llh_vec[ps_shps_tosim] + 2 * (dfmean_cbic) / n
cic_vals_edf0[ps_shps_tosim] <- llh_vec[ps_shps_tosim] + 2 * (dfmean_cbic) / n
}
}
#new
if(any(cic_vals_max < cic_st)){
if (wt.iter & infocrit == "CBIC") {
if(length(cic_vals_edf0) == 5){
cic_vals[1:length(cic_vals_edf0), ix] <- cic_vals_edf0
} else {
#tree
cic_vals[1:2, ix] <- cic_vals_edf0[1:2]
cic_vals[9, ix] <- cic_vals_edf0[3]
}
} else {
if(length(cic_vals_edf0) == 9){
cic_vals[, ix] <- cic_vals_edf0
} else if(length(cic_vals_edf0) == 5){ #new: set4
cic_vals[1:length(cic_vals_edf0), ix] <- cic_vals_edf0
} else {
#tree:
cic_vals[1:2, ix] <- cic_vals_edf0[1:2]
cic_vals[9, ix] <- cic_vals_edf0[3]
}
}
ps_shpi <- which(cic_vals_edf0 == min(cic_vals_edf0))
if(length(ps_shpi) > 1) {
ps_shpi <- rev(ps_shpi)[1]
}
shpi <- shps[ps_shpi]
etahatkeep <- etahat_lst[[ps_shpi]]
coefskeep <- bhat_lst[[ps_shpi]]
mskeep <- ms_lst[[ps_shpi]]
kts <- kts_lst[[ps_shpi]] #will be NULL for ordinal
dd <- dd_lst[[ps_shpi]]
dm <- dm_lst[[ps_shpi]]
ps <- ps_vec[ps_shpi]
if(shpi != 18) {
#new:
if(shpi == 2 & attr(x, 'type') == 'tree'){
ddkeep_x <- NULL
ddkeep_z <- dd
np_add <- ncol(dd)
} else {
ddkeep_x <- dd
ddkeep_z <- NULL
np_add <- 0
if (shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5) {
ddkeep_z <- cbind(x - mean(x))
np_add <- 1
}
}
if(pnt){
dmat_keep <- dm #new dmat, use dmat for all dmat's combined as block diag
ps_keep <- ps
}
ktskeep <- kts
}
if(shpi == 18) {
ddkeep_x <- NULL
if(pnt){
dmat_keep <- NULL
ps_keep <- 0
}
ktskeep <- NULL
#cat('include linear', '\n')
ddkeep_z <- cbind(x - mean(x))
np_add <- 1
}
cic_st_next_step_0 <- cic_vals_edf0[ps_shpi]
cic_st_next_step <- cic_st_next_step_0
cic_st <- cic_st_next_step_0
}
}
}
rslt <- list(cic_vals = cic_vals, ddkeep_x = ddkeep_x, ktskeep = ktskeep,
ddkeep_z = ddkeep_z, cic_st = cic_st, cic_st_next_step = cic_st_next_step,
etahatkeep = etahatkeep, coefskeep = coefskeep, mskeep = mskeep,
rm_id = rm_id, dmat_keep = dmat_keep, ps_keep = ps_keep,
dfmean_cbic_vals = dfmean_cbic_vals)
return (rslt)
}
#----------------------------------------------------------------------------------------------------------
#backward step used in classo or stepwise
#----------------------------------------------------------------------------------------------------------
backward_cic <- function(xmat, y, varlst, varlst_z, ddkeep_x_all,
ddkeep_x_lst, knots_lst=NULL, boundary.knots = NULL,
ddkeep_z_all, ddkeep_z_lst, mskeep_lst, mod_tried, nrep, totx,
oldcic = NULL, infocrit = "CAIC", dfmean_old, dfmeanmax_old, penalty_old, oldcoefs,
oldetahat, active, type = "stepwise", nsim = 100, wt.iter = FALSE,
family = gaussian(), parallel = FALSE, space = "E", pnt = FALSE,
dmat_keep_lst = list(NULL),
ps_keep_vec = NULL, a=0.5,
add_edfu=0, add_kts=0, shps_all)
{
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
muhat.fun <- cicfamily$muhat.fun
ysim.fun <- cicfamily$ysim.fun
deriv.fun <- cicfamily$deriv.fun
dev.fun <- cicfamily$dev.fun
chs <- which(mod_tried[nrep-1, 1:totx] != 0) #already removed z (shp = 0)
changed <- FALSE
rm_chi <- mod_tried[nrep-1, 1:totx]
shps <- 9:17
#new: if we use cbic and wt.iter then we only have 4 options
if(wt.iter & infocrit == "CBIC"){
shps <- c(9:12, 17)
}
ddkeep_x_all_0 <- ddkeep_x_all
#cat ('dim of x0', dim(ddkeep_x_all), '\n')
ddkeep_z_all_0 <- ddkeep_z_all
dmat_keep_lst_0 <- dmat_keep_lst
dmat_keep_all <- dmat_keep_all_0 <- dmat <- dm <- NULL
#new: avoid computing coefskeep again
varlst_0 <- varlst
varlst_z_0 <- varlst_z
#new! order dmat_keep_lst to get dmat_keep_all
#follow the same order of ddkeep_x_all
#cannot use this: may add some x twice if the x was out of the model in the previous step: 2 5 3 1 6 2
seq_added <- mod_tried[1:nrep, totx + 3]
#test!
seq_added <- seq_added[!duplicated(seq_added, fromLast = TRUE)]
if(pnt){
dmat_keep_all <- Matrix::bdiag(Filter(Negate(is.null), dmat_keep_lst[seq_added]))
dmat_keep_all <- as.matrix(dmat_keep_all)
if(prod(dim(dmat_keep_all)) == 0){
dmat_keep_all <- NULL
}
dmat_keep_all_0 <- dmat_keep_all
}
cic_st_next_step <- oldcic
cic_st <- oldcic
n <- NROW(xmat)
one <- cbind(rep(1, n))
#dfmeankeep <- dfmean_old
etahatkeep <- oldetahat
coefskeep <- oldcoefs
back_inactive <- NULL
back_inactive_check <- rep(FALSE, length(chs))
totx <- NCOL(xmat)
coefskeep_x_lst <- vector("list", length = totx)
coefskeep_z_lst <- vector("list", length = totx)
#xtxkeep <- oldxtx
#facekeep <- oldface
for(chi in chs){
#print (varlst_0)
#print (chi)
#-----------------------------------------------------------------------------------------
#change the following to include shapes loop
#go through all un tried shapes for the previously chosen x's to see if cic is smaller
#-----------------------------------------------------------------------------------------
x <- xmat[, chi, drop=TRUE]
#----------------------------------------------------
#won't check again for categorical predictor ?
#rm_chi[chi] records the chosen shape
if (rm_chi[chi] == 1 & inherits(x, "factor")) {
#should update varlst_0 and varlst_z_0 here
varlst_0 <- varlst
varlst_z_0 <- varlst_z
shps_not_tried <- 0
# ddkeep_z_all track changing z columns
# ddkeep_z_all_0 track final z columns whenever cic_st is decreased
ddkeep_z_all <- ddkeep_z_all_0[, -which(varlst_z == chi), drop=FALSE]
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all_0)
np_ix <- NCOL(ddkeep_z_all) + 1
#new
if(pnt){
dmat <- dmat_keep_all_0
}
dfmean_rslt <- compute_edf(family=family, gmat_lst=list(gmat),
dmat_lst = list(dmat),
mu0_lst=NULL, np_ix_lst=list(np_ix),
weights=NULL, nsim=nsim, wt.iter=wt.iter,
parallel=parallel, pnt=pnt)
dfmean_caic <- dfmean_rslt$dfmean_caic #val without z; val with z is oldcic
dfmean_cbic <- dfmean_rslt$dfmean_cbic
nc <- ncol(gmat)
amat <- diag(nc)
ps_zero <- 1:np_ix
amat[ps_zero, ps_zero] <- 0
bvec <- cbind(rep(0, nc))
xtx <- crossprod(gmat)
if(pnt & !is.null(dmat)){
dm_zero <- matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat <- cbind(dm_zero, dmat)
dv <- crossprod(dmat)
xtx <- xtx + dv
}
if(np_ix > 0 & NCOL(gmat) > np_ix) {# z and x
zvec <- y
if(!wt.iter){
#test
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ans <- coneB(zvec, dsend, zsend)
bhat <- ans$coefs
bhat = round(bhat, 10)
etahat <- gmat %*% bhat
} else {
xty <- crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ans <- solve.QP(xtx, xty, amat, bvec)
#ans = solve.QP(xtx, xty, t(amat), bvec)
bhat <- ans$solution
bhat = round(bhat, 10)
etahat <- gmat %*% bhat
}
} else {
gmat0 <- gmat
weights <- rep(1, n)
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
muhat <- family$linkinv(etahat)
#sm <- 1e-5
sm <- 1e-7
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 20)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
#test
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ansi <- coneB(zvec, dsend, zsend)
bhat <- ansi$coefs
}
if(pnt & !is.null(dmat) & !all(dmat == 0)){
xtx <- crossprod(gmat)
xtx <- xtx + dv
xty <- crossprod(gmat, zvec)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ansi <- solve.QP(xtx, xty, amat, bvec)
#ansi = solve.QP(xtx, xty, t(amat), bvec)
bhat <- ansi$solution
}
bhat = round(bhat, 10)
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
} else if (np_ix > 0 & NCOL(gmat) == np_ix) { #all z #shpi = 18 and no constrained x has been chosen
zvec <- y
if(!wt.iter){
#bhat <- solve(crossprod(gmat), t(gmat)) %*% zvec
bhat <- ginv(crossprod(gmat)) %*% t(gmat) %*% zvec
bhat = round(bhat, 10)
etahat <- gmat %*% bhat
} else {
m0 <- stats::glm(y ~ gmat - 1, family = family)
etahat <- m0$linear.predictors
bhat <- coef(m0)
muhat <- family$linkinv(etahat)
}
}
# #mu0 is used for CIC simulation later
rss <- sum((y - etahat)^2)
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml = family$family)
} else {llh <- log(rss)}
if (infocrit == "CBIC") {
#cic_val <- llh + log(n) * (dfmean_cbic) / n
#new
#cic_val <- llh + a * log(n) * (dfmean_cbic) / n^(6/7)
all_shps <- shps_all
all_shps[chi] <- 0
if(any(all_shps %in% c(9, 10))){ #smooth incr/decr
degree <- 6/7
} else if (any(all_shps > 10 & all_shps < 17)){ #smooth conv/conc
degree <- 6/7#8/9 #6/7#
} else if (any(all_shps > 0 & all_shps < 9)) { #ord or parametric
degree <- 6/7#1 #6/7#
} else if (all(all_shps %in% c(0, 17))) { #won't happen
stop('check!')
}
cic_val <- llh + a * log(n) * (dfmean_cbic) / n^(degree)
}
if (infocrit == "CAIC") {
cic_val <- llh + 2 * (dfmean_caic) / n
}
# #new: compare without z vs with z; cic_val[1] is with z
bool <- (cic_val - cic_st > 0) & (cic_val - cic_st) < 1e-8
if(cic_val < cic_st | bool){
#cat('z changed in backward!', 'nrep:', nrep, '\n')
changed <- TRUE
new_shp <- shps_not_tried #will always = 0
mod_tried[nrep, chi] <- new_shp
#new: check more
mod_tried[nrep, totx + 2] <- cic_val
etahatkeep <- etahat
coefskeep <- bhat
ddkeep_z <- NULL
ddkeep_z_lst[chi] <- list(NULL)
#new:
varlst_z_0 <- varlst_z_0[-which(varlst_z_0 == chi)]
#new:
#coefskeep_z_lst[chi] <- list(NULL)
mskeep_lst[chi] <- list(NULL)
np_add <- 0
knots_lst[chi] <- list(NULL)
ddkeep_x <- NULL
ddkeep_x_lst[chi] <- list(NULL) #is already null bc it is factor?
#new: need to update ddkeep_z_all_0 and ddkeep_x_all_0 here in case >= 2 predictors' shapes are changed
seq_added <- mod_tried[1:nrep, totx + 3]
ddkeep_z_lst_ord <- ddkeep_z_lst[seq_added]
ddkeep_z_all_0 <- do.call(base::cbind, ddkeep_z_lst_ord)
obs <- 1:totx
obs_ord <- obs[seq_added]
#new
coefs_nointercept <- coefskeep[-1]
#print (varlst_z_0)
if(length(varlst_z_0) > 0){
zid <- 1:length(varlst_z_0)
zcoefs <- coefs_nointercept[zid]
xcoefs <- coefs_nointercept[-zid]
}
if(length(varlst_z_0) == 0){
zcoefs <- NULL
xcoefs <- coefs_nointercept
}
#new: we should update ddkeep_x_all and ddkeep_z_all here
ddkeep_x_all_temp <- do.call(base::cbind, ddkeep_x_lst)
ddkeep_z_all_temp <- do.call(base::cbind, ddkeep_z_lst)
accx <- NCOL(ddkeep_x_all_temp)
accz <- NCOL(ddkeep_z_all_temp)
coefskeep_x_lst <- vector("list", length = totx)
coefskeep_z_lst <- vector("list", length = totx)
if(accx > 0){
#new
#print (varlst_0)
ux_added <- unique(varlst_0)
for(x_added in ux_added){
coefskeep_x_lst[[x_added]] <- xcoefs[which(varlst_0 %in% x_added)]
}
#order to be put back in coefskeep
ord_x <- order(varlst_0)
xcoefs <- xcoefs[ord_x]
}
if(accz > 0){
varlst_z_lst <- lapply(1:length(obs_ord), function(i, zlst) {if(NCOL(zlst[[i]]) >= 1) {rep(obs_ord[i], NCOL(zlst[[i]]))}}, zlst = ddkeep_z_lst_ord)
varlst_z <- do.call(base::c, varlst_z_lst)
#new
uz_added <- unique(varlst_z_0)
for(z_added in uz_added){
coefskeep_z_lst[[z_added]] <- zcoefs[which(varlst_z_0 %in% z_added)]
}
ord_z <- order(varlst_z_0)
zcoefs <- zcoefs[ord_z]
} else {varlst_z <- NULL}
#new: put b0 back
coefskeep <- c(coefskeep[1], zcoefs, xcoefs)
cic_st <- cic_val
if (infocrit == "CBIC"){
#cic_st_next_step_0 <- llh + log(n) * (dfmean_cbic) / n
#new
#cic_st_next_step_0 <- llh + a * log(n) * (dfmean_cbic) / n^(6/7)
cic_st_next_step_0 <- llh + a * log(n) * (dfmean_cbic) / n^(degree)
}
if (infocrit == "CAIC"){
cic_st_next_step_0 <- llh + 2 * (dfmean_caic) / n
}
if(cic_st_next_step_0 < oldcic){
cic_st_next_step <- cic_st_next_step_0
}
#new: put it back to active set
back_inactive_check[which(chs %in% chi)] <- TRUE
}
#}
} else {
#shps <- 9:18
#should update varlst_0 and varlst_z_0 here
varlst_0 <- varlst
varlst_z_0 <- varlst_z
if(attr(x, 'type') %in% c('num', 'numeric')){
shps <- 9:17
}
if(attr(x, 'type') %in% c('ord', 'ordinal')){
shps <- c(1:8, 0)
}
if(attr(x, 'type') %in% c('tree')){
shps <- c(1:2, 0)
}
#new: if we use cbic and wt.iter then we only have 4 options
if(wt.iter & infocrit == "CBIC"){
if(attr(x, 'type') %in% c('num', 'numeric')){
shps <- c(9:12, 17)
}
if(attr(x, 'type') %in% c('ord', 'ordinal')){
shps <- c(1:4, 0)
}
if(attr(x, 'type') %in% c('tree')){
shps <- c(1:2, 0)
}
}
shps_not_tried <- shps[!shps %in% rm_chi[chi]]
#new: speed up the code
if(attr(x, 'type') %in% c('num', 'numeric')){
shps_not_tried <- 17
} else {
shps_not_tried <- 0
}
#cic_st = mod_tried[nrep, totx + 2]
#test more!
nobs <- length(shps_not_tried) + 1 #include the chosen shape too
obs <- 1:nobs
gmat_lst <- vector("list", length = nobs) #accumulative
np_ix_lst <- vector("list", length = nobs)
#mu0_lst <- vector("list", length = nobs)
ms_lst <- vector("list", length = nobs)
dd_lst <- vector("list", length = nobs) #for new
kts_lst <- vector("list", length = nobs)
dmat_lst <- vector("list", length = nobs) #accumulative
dm_lst <- vector("list", length = nobs) #for new
ps_vec <- rep(0, nobs)
#new
#st_new_vec <- rep(0, nobs)
#ed_new_vec <- rep(0, nobs)
#add the chosen shp first
if(rm_chi[chi] %in% c(17, 0)){
dd_lst[1] <- list(NULL)
#dmat_lst[1] <- list(NULL)
} else if (rm_chi[chi] %in% c(18) | rm_chi[chi] == 2 & attr(x, 'type') == 'tree'){
if(rm_chi[chi] %in% c(18)){
dd_lst[[1]] <- cbind(x - mean(x))
} else {
dd_lst[[1]] <- ddkeep_z_all_0[, which(varlst_z == chi), drop=FALSE]
}
#dmat_lst[1] <- list(NULL)
} else {
dd_lst[[1]] <- ddkeep_x_all_0[, which(varlst == chi), drop=FALSE]
if(pnt){
dm_lst[[1]] <- dmat_keep_lst_0[[chi]]
}
#dmat_lst[1] <- dmat_keep_all_0[, which(varlst == chi), drop=FALSE]
}
gmat_lst[[1]] <- cbind(one, ddkeep_z_all_0, ddkeep_x_all_0)
np_ix_lst[[1]] <- NCOL(ddkeep_z_all_0) + 1
#new:
if(pnt & !is.null(dmat_keep_all_0)){
dmat_lst[[1]] <- dmat_keep_all_0
}
#record gmat, np_ix, etc
for(i in 2:nobs)
{
shpi <- shps_not_tried[i-1]
if(rm_chi[chi] %in% c(17, 18, 0) | rm_chi[chi] == 2 & attr(x, 'type') == 'tree'){
ddkeep_x_all <- ddkeep_x_all_0
if(pnt){
dmat_keep_all <- dmat_keep_all_0
}
} else {
ddkeep_x_all <- ddkeep_x_all_0[, -which(varlst == chi), drop=FALSE]
if(pnt){
dmat_keep_all <- Matrix::bdiag(Filter(Negate(is.null), dmat_keep_lst_0[seq_added[-which(seq_added %in% chi)]]))
dmat_keep_all <- as.matrix(dmat_keep_all)
if(prod(dim(dmat_keep_all)) == 0){
dmat_keep_all <- NULL
}
}
}
#if the chosen shape is conv or conc, then need to remove from ddkeep_z_all
if(rm_chi[chi] > 10 & rm_chi[chi] < 13 | rm_chi[chi] == 18 | rm_chi[chi] > 2 & rm_chi[chi] < 5 | rm_chi[chi] == 2 & attr(x, 'type') == 'tree'){
ddkeep_z_all <- ddkeep_z_all_0[, -which(varlst_z == chi), drop=FALSE]
} else {
ddkeep_z_all <- ddkeep_z_all_0
}
if (shpi == 17 | shpi == 0) {
ms <- NULL
dd <- NULL
kts <- NULL
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + 1
if(pnt){
dmat <- dmat_keep_all
#wrong! shp not tried is 17, then cic_val is not oldcic
#cic_val <- oldcic
#not_skip_cic <- FALSE
dm <- NULL
ps <- 0
}
} else if (shpi == 18) {
#dd <- cbind(x)
ms <- mean(x)
dd <- cbind(x - mean(x))
kts <- NULL
gmat <- cbind(one, ddkeep_z_all, dd, ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + 2
if(pnt){
dmat <- dmat_keep_all
dm <- NULL
ps <- 0
}
#not_skip_cic <- TRUE
} else {
if(shpi <= 8){
if(attr(x, 'type') %in% c('ord')){
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = 0, space = space)
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
if(shpi > 2 & shpi < 5){
gmat <- cbind(one, ddkeep_z_all, x-mean(x), ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 2
} else {
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 1
}
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
}
if(attr(x, 'type') %in% c('tree')){
dd <- t(tree.fun(x))
ms <- apply(dd, 2, mean)
dd <- apply(dd, 2, function(e) e - mean(e))
if(shpi == 1){
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 1
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
}
if(shpi == 2){
gmat <- cbind(one, ddkeep_z_all, dd, ddkeep_x_all)
np_ix <- NCOL(ddkeep_z_all) + 1 + NCOL(dd)
dm <- NULL
}
}
ps <- 0
dmat <- dmat_keep_all
if(!is.null(dmat)){
dm_lst_tmp <- vector('list', length = 2)
dm_lst_tmp[[1]] <- dmat
dm_lst_tmp[[2]] <- dm
dmat <- as.matrix(bdiag(dm_lst_tmp))
} else {
dmat <- dm
}
}
if(shpi > 8){
mult <- ifelse(pnt, 2, 1)
spl_degree_i <- ifelse(shpi %in% c(9, 10), 2, 3) #use 2 for i-spl, 3 for c-spl
n1 <- length(unique(x))
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
#new
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/7)) + 6 + add_kts, trunc(n1^(1/9)) + 6 + add_kts)
edfu <- ifelse(spl_degree_i == 2, nkts/mult, nkts/mult + 1) + add_edfu
kts <- seq.int(min(x), max(x), length.out = nkts)
#new: boundary knots
kts0 <- boundary.knots[[chi]]
if(length(kts0) >= 2){
x_rng <- range(x, na.rm = TRUE)
kts0_rng <- range(kts0, na.rm = TRUE)
tol <- 1e-10
reset_lwr <- kts0_rng[1] > x_rng[1] + tol
reset_upp <- kts0_rng[2] < x_rng[2] - tol
if(reset_lwr){
kts_lwr <- x_rng[1]
} else {
kts_lwr <- kts0_rng[1]
}
if(reset_upp){
kts_upp <- x_rng[2]
} else {
kts_upp <- kts0_rng[2]
}
kts[1] <- kts_lwr
kts[nkts] <- kts_upp
}
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = kts, space = space)
#dd is orth to 1
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
if(shpi > 10 & shpi < 13){
gmat <- cbind(one, ddkeep_z_all, x-mean(x), ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 2
} else {
gmat <- cbind(one, ddkeep_z_all, ddkeep_x_all, dd)
np_ix <- NCOL(ddkeep_z_all) + 1
}
if(pnt){
nc <- NCOL(dd)
dm <- makedmat_1D_ic(m = nc, knots = kts, shpi = shpi) #different from makedmat_1D for testpar
qvi <- crossprod(dd)
dvi <- crossprod(dm)
#ps <- 1e-6
#9, 10, 11, 12, 13, 14, 15, 16
#if(shpi %in% c(9, 11, 13)){
#ps <- uniroot(f = .search_ps, qv0 = qvi, dv = dvi, dd = dd,
# edfu = edfu, interval = c(1e-10, 1e+10), tol = .Machine$double.eps^0.25)$root
#max_log_pen <- log(NROW(dd))
log_root_res <- uniroot(f = function(log_pen,...) .search_ps_shapeselect(exp(log_pen),...),
interval = c(-8, 8),
qv0 = qvi, dv = dvi, dd = dd, edfu = edfu,
extendInt = "yes",
tol = .Machine$double.eps^0.25)
# convert back to the original scale
#temp
if(abs(log_root_res$root) > 8) {
ps <- floor(exp(8))
} else {
ps <- exp(log_root_res$root)
}
#ps <- exp(log_root_res$root)
ps_shpi <- which(shps_not_tried %in% shpi) + 1 #1st is for shp chosen
#print (ps_shpi)
ps_vec[ps_shpi] <- ps
#} else if (shpi == 10){
# ps <- ps_vec[2]
# ps_shpi <- which(shps_not_tried %in% shpi) + 1 #1st is for shp chosen
# ps_vec[ps_shpi] <- ps
#} else if (shpi == 12) {
# ps <- ps_vec[4]
# ps_shpi <- which(shps_not_tried %in% shpi) + 1 #1st is for shp chosen
# ps_vec[ps_shpi] <- ps
#} else {
# ps <- ps_vec[6]
# ps_shpi <- which(shps_not_tried %in% shpi) + 1 #1st is for shp chosen
# ps_vec[ps_shpi] <- ps
#}
dm <- sqrt(ps) * dm
dmat <- dmat_keep_all
if(!is.null(dmat)){
dm_lst_tmp <- vector('list', length = 2)
dm_lst_tmp[[1]] <- dmat
dm_lst_tmp[[2]] <- dm
dmat <- as.matrix(bdiag(dm_lst_tmp))
} else {
dmat <- dm
}
}
#not_skip_cic <- TRUE
}
}
if(shpi == 17 | shpi == 0) {
dd_lst[i] <- list(NULL)
ms_lst[i] <- list(NULL)
kts_lst[i] <- list(NULL)
} else {
if(shpi == 18) {
#dd_lst[i] <- list(NULL)
dd_lst[[i]] <- dd
ms_lst[[i]] <- ms
kts_lst[i] <- list(NULL)
} else {
dd_lst[[i]] <- dd
ms_lst[[i]] <- ms
#NULL for ordinal
if(shpi > 8){
kts_lst[[i]] <- kts
}
}
}
if(pnt & !is.null(dm)){
dm_lst[[i]] <- dm
}
gmat_lst[[i]] <- gmat
np_ix_lst[[i]] <- np_ix
#mu0_lst[[i]] <- mu0
if(pnt & !is.null(dmat)){
dmat_lst[[i]] <- dmat
}
}
shp17_id <- NULL
if(any(shps_not_tried %in% c(17, 0))){
#shp17_id = which(shps_not_tried %in% 17)
shp17_id <- which(shps_not_tried %in% c(17, 0)) + 1 #1st one is for chosen shp
}
dfmean_rslt <- compute_edf(family=family, gmat_lst=gmat_lst, dmat_lst=dmat_lst,
mu0_lst=NULL, np_ix_lst=np_ix_lst,
weights=NULL, shp17_id=shp17_id, nsim=nsim, wt.iter=wt.iter,
parallel=parallel, pnt=pnt)
dfmean_caic <- dfmean_rslt$dfmean_caic #if 17 is in shps not tried, the caic has two vals, the 2nd one is for 17
dfmean_cbic <- dfmean_rslt$dfmean_cbic
#print (dfmean_cbic)
#go through all shapes to compute rss and cic/bic
for(shpi in shps_not_tried){
#print (shpi)
ps_shpi <- which(shps_not_tried %in% shpi) + 1 #1st is for shp chosen
#ps_shpi = which(9:18 %in% shpi)
#new
#st_new <- st_new_vec[ps_shpi]
#ed_new <- ed_new_vec[ps_shpi]
np_ix <- np_ix_lst[[ps_shpi]]
gmat <- gmat_lst[[ps_shpi]]
ms <- ms_lst[[ps_shpi]]
if(shpi %in% c(17, 0)) {
dd <- NULL
dm <- NULL
} else {
dd <- dd_lst[[ps_shpi]]
if(pnt){
dm <- dm_lst[[ps_shpi]]
}
}
kts <- kts_lst[[ps_shpi]]
if(pnt){
dmat <- dmat_lst[[ps_shpi]]
ps <- ps_vec[ps_shpi]
}
nc <- ncol(gmat)
amat <- diag(nc)
#amat[1:np_ix, 1:np_ix] = 0
ps_zero <- 1:np_ix
amat[ps_zero, ps_zero] <- 0
bvec <- cbind(rep(0, nc))
xtx <- crossprod(gmat)
if(pnt & !is.null(dmat)){
dm_zero <- matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat <- cbind(dm_zero, dmat)
dv <- crossprod(dmat)
xtx <- xtx + dv
}
if(np_ix > 0 & NCOL(gmat) > np_ix) {# z and x
zvec <- y
if(!wt.iter){
#test
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ans <- coneB(zvec, dsend, zsend)
bhat <- ans$coefs
bhat = round(bhat, 10)
etahat <- gmat %*% bhat
face <- ans$face
} else {
xty <- crossprod(gmat, zvec)
#test
# if(qr(xtx)$rank < NCOL(xtx)){
# xtx <- xtx + diag(1e-8, NCOL(xtx))
# }
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ans <- solve.QP(xtx, xty, amat, bvec)
#ans = solve.QP(xtx, xty, t(amat), bvec)
bhat <- ans$solution
bhat = round(bhat, 10)
etahat <- gmat %*% bhat
face <- ans$iact
}
} else {
gmat0 <- gmat
weights <- rep(1, n)
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
muhat <- family$linkinv(etahat)
#sm <- 1e-5
sm <- 1e-7
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 20)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
#test
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ansi <- coneB(zvec, dsend, zsend)
bhat <- ansi$coefs
}
if(pnt & !is.null(dmat) & !all(dmat == 0)){
xtx <- crossprod(gmat)
xtx <- xtx + dv
xty <- crossprod(gmat, zvec)
#test
# if(qr(xtx)$rank < NCOL(xtx)){
# xtx <- xtx + diag(1e-8, NCOL(xtx))
# }
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
ansi <- solve.QP(xtx, xty, amat, bvec)
#ansi = solve.QP(xtx, xty, t(amat), bvec)
bhat <- ansi$solution
face <- ansi$iact
}
bhat = round(bhat, 10)
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
#new
# coef_xc_new <- NULL
# coef_x_new <- NULL
# if(shpi > 10 & shpi < 13){
# coef_xc_new <- bhat[np_ix]
# coef_x_new <- bhat[st_new:ed_new]
# } else {
# #coef_xc_new <- NULL
# coef_x_new <- bhat[st_new:ed_new]
# }
} else if (np_ix > 0 & NCOL(gmat) == np_ix) { #all z #shpi = 18 and no constrained x has been chosen
zvec <- y
if(!wt.iter){
#bhat <- solve(crossprod(gmat), t(gmat)) %*% zvec
bhat <- ginv(crossprod(gmat)) %*% t(gmat) %*% zvec
bhat = round(bhat, 10)
etahat <- gmat %*% bhat
} else {
m0 <- stats::glm(y ~ gmat - 1, family = family)
etahat <- m0$linear.predictors
bhat <- coef(m0)
muhat <- family$linkinv(etahat)
}
}
rss <- sum((y - etahat)^2)
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml = family$family)
#test!
} else {llh <- log(rss)}
ps_shpi_bic <- which(shps_not_tried %in% shpi) + 1 #1st is for chosen shp
if (infocrit == "CBIC" & shpi %in% c(17, 0)) {
#cic_val <- llh + log(n) * (dfmean_caic[2]) / n
#new
#cic_val <- llh + a * log(n) * (dfmean_caic[2]) / n^(6/7)
all_shps <- shps_all
all_shps[chi] <- shpi
if(any(all_shps %in% c(9, 10))){ #smooth incr/decr
degree <- 6/7
} else if (any(all_shps > 10 & all_shps < 17)){ #smooth conv/conc
degree <- 6/7 #8/9#6/7#
} else if (any(all_shps > 0 & all_shps < 9)) { #ord or parametric; e.g., all_shps = c(8,17)
degree <- 6/7 #1# 6/7#
} else if (all(all_shps %in% c(0, 17))) {
stop('check!')
}
cic_val <- llh + a * log(n) * (dfmean_caic[2]) / n^(degree)
}
if (infocrit == "CBIC" & shpi != 17 & shpi != 0) {
#cic_val <- llh + log(n) * (dfmean_cbic[ps_shpi_bic]) / n
#new:
#cic_val <- llh + a * log(n) * (dfmean_cbic[ps_shpi_bic]) / n^(6/7)
all_shps <- shps_all
all_shps[chi] <- shpi
if(any(all_shps %in% c(9, 10))){ #smooth incr/decr
degree <- 6/7
} else if (any(all_shps > 10 & all_shps < 17)){ #smooth conv/conc
degree <- 6/7# 8/9#6/7#
} else if (any(all_shps > 0 & all_shps < 9)) { #ord or parametric
degree <- 6/7# 1#6/7#
} else {
stop('check!')
}
cic_val <- llh + a * log(n) * (dfmean_cbic[ps_shpi_bic]) / n^(degree)
}
if (infocrit == "CAIC" & shpi %in% c(17, 0)) {
cic_val <- llh + 2 * (dfmean_caic[2]) / n
}
if (infocrit == "CAIC" & shpi != 17 & shpi != 0) {
#test
cic_val <- llh + 2 * (dfmean_cbic[ps_shpi_bic]) / n
}
#test more!
bool <- (cic_val - cic_st > 0) & (cic_val - cic_st) < 1e-8
#new: if etahat = oldetahat, then we should remove x even though cic_val > cic_st?
bool2 <- isTRUE(all.equal(etahat, etahatkeep, tolerance = 1e-8))
bool <- bool | bool2
if(cic_val < cic_st | bool){
#cat(paste0('x', chi), 'changed in backward!', 'nrep:', nrep, '\n')
changed <- TRUE
new_shp <- shpi
#new:
if(new_shp == 17){
new_shp <- 0
}
#new: put chi back to active set
if(shpi == 17 | shpi == 0) {
back_inactive_check[which(chs %in% chi)] <- TRUE
#cat('x', chi, 'put back in inactive set!', 'nrep:', nrep, '\n')
}
mod_tried[nrep, chi] <- new_shp
#new: check more
mod_tried[nrep, totx+2] <- cic_val
etahatkeep <- etahat
#print (head(etahatkeep))
coefskeep <- bhat
if(shpi != 18) {
#new:
if(shpi != 2 | attr(x, 'type') != 'tree'){
ddkeep_x <- dd
} else {
ddkeep_z <- dd
ddkeep_x <- NULL
}
if(pnt){
dmat_keep <- dm
ps_keep <- ps
}
ktskeep <- kts
if(shpi == 17 | shpi == 0) {
ddkeep_x_lst[chi] <- list(NULL)
#new
if(any(varlst_0 == chi)){
varlst_0 <- varlst_0[-which(varlst_0 == chi)]
}
#test more
ddkeep_z_lst[chi] <- list(NULL)
#new
if(any(varlst_z_0 == chi)){
varlst_z_0 <- varlst_z_0[-which(varlst_z_0 == chi)]
}
if(pnt){
dmat_keep_lst[chi] <- list(NULL)
ps_keep_vec[chi] <- 0
}
mskeep_lst[chi] <- list(NULL)
kts_lst[chi] <- list(NULL)
} else {
if(shpi == 2 & attr(x, 'type') == 'tree'){
ddkeep_x_lst[chi] <- list(NULL)
if(any(varlst_0 == chi)){
varlst_0 <- varlst_0[-which(varlst_0 == chi)]
}
if(pnt){
dmat_keep_lst[chi] <- list(NULL)
ps_keep_vec[chi] <- 0
}
mskeep_lst[[chi]] <- ms
} else {
ddkeep_x_lst[[chi]] <- dd
#new
if(any(varlst_0 == chi)){
varlst_0 <- varlst_0[-which(varlst_0 == chi)]
}
varlst_0 <- c(varlst_0, rep(chi, ncol(dd)))
#new
#coefskeep_x_lst[[chi]] <- coef_x_new
if(pnt){
dmat_keep_lst[[chi]] <- dm
ps_keep_vec[chi] <- ps
}
mskeep_lst[[chi]] <- ms
#NULL for ordinal
if(shpi > 8){
knots_lst[[chi]] <- kts
}
}
}
if(shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5 | shpi == 2 & attr(x, 'type') == 'tree'){
if(shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5){
ddkeep_z <- cbind(x - mean(x))
ddkeep_z_lst[[chi]] <- cbind(x - mean(x))
#new
if(any(varlst_z_0 == chi)){
varlst_z_0 <- varlst_z_0[-which(varlst_z_0 == chi)]
}
varlst_z_0 <- c(varlst_z_0, rep(chi, 1))
np_add <- 1
} else {
ddkeep_z <- dd
ddkeep_z_lst[[chi]] <- dd
#new
if(any(varlst_z_0 == chi)){
varlst_z_0 <- varlst_z_0[-which(varlst_z_0 == chi)]
}
varlst_z_0 <- c(varlst_z_0, rep(chi, NCOL(dd)))
np_add <- NCOL(dd)
}
} else {
ddkeep_z <- NULL
ddkeep_z_lst[chi] <- list(NULL)
#new
if(any(varlst_z_0 == chi)){
varlst_z_0 <- varlst_z_0[-which(varlst_z_0 == chi)]
}
#new
#coefskeep_z_lst[chi] <- list(NULL)
np_add <- 0
}
#varlst_lst[[chi]] <- rep(chi, NCOL(dd))
}
if (shpi == 18) {
ddkeep_z <- cbind(x - mean(x))
ddkeep_z_lst[[chi]] <- cbind(x - mean(x))
#new
if(any(varlst_z_0 == chi)){
varlst_z_0 <- varlst_z_0[-which(varlst_z_0 == chi)]
}
varlst_z_0 <- c(varlst_z_0, rep(chi, 1))
mskeep_lst[[chi]] <- mean(x)
np_add <- 1
knots_lst[chi] <- list(NULL)
ddkeep_x <- NULL
ddkeep_x_lst[chi] <- list(NULL)
#new
if(any(varlst_0 == chi)){
varlst_0 <- varlst_0[-which(varlst_0 == chi)]
}
if(pnt){
dmat_keep <- NULL
ps_keep <- 0
dmat_keep_lst[chi] <- list(NULL)
ps_keep_vec[chi] <- 0
}
}
#new: need to update ddkeep_z_all_0 and ddkeep_x_all_0 here in case >= 2 predictors' shapes are changed
seq_added <- mod_tried[1:nrep, totx + 3]
# ddkeep_x_lst_ord will only store dd's for chosen x
ddkeep_x_lst_ord <- ddkeep_x_lst[seq_added]
ddkeep_z_lst_ord <- ddkeep_z_lst[seq_added]
if(pnt){
dmat_keep_lst_0 <- dmat_keep_lst
#dmat_keep_all_0 <- as.matrix(bdiag(Filter(Negate(is.null), dmat_keep_lst_0)))
#new!
#dmat_keep_all_0 <- as.matrix(bdiag(dmat_keep_lst_0[seq_added]))
dmat_keep_all_0 <- as.matrix(bdiag(Filter(Negate(is.null), dmat_keep_lst_0[seq_added])))
}
ddkeep_x_all_0 <- do.call(base::cbind, ddkeep_x_lst_ord)
ddkeep_z_all_0 <- do.call(base::cbind, ddkeep_z_lst_ord)
obs <- 1:totx
obs_ord <- obs[seq_added]
#new: not in the order of x1,x2...; the chi one is at the end.
coefs_nointercept <- coefskeep[-1]
#print (varlst_z_0)
if(length(varlst_z_0) > 0){
zid <- 1:length(varlst_z_0)
zcoefs <- coefs_nointercept[zid]
xcoefs <- coefs_nointercept[-zid]
}
if(length(varlst_z_0) == 0){
zcoefs <- NULL
xcoefs <- coefs_nointercept
}
#cat('chi:', chi, '\n')
#cat('shpi:', shpi, '\n')
#print (varlst_0)
#print (zcoefs)
#new: we should update ddkeep_x_all and ddkeep_z_all here
ddkeep_x_all_temp <- do.call(base::cbind, ddkeep_x_lst)
ddkeep_z_all_temp <- do.call(base::cbind, ddkeep_z_lst)
accx <- accz <- 0
accx <- NCOL(ddkeep_x_all_temp)
accz <- NCOL(ddkeep_z_all_temp)
totx <- NCOL(xmat)
#new
coefskeep_x_lst <- vector("list", length = totx)
coefskeep_z_lst <- vector("list", length = totx)
if(accx > 0){
varlst_lst <- lapply(1:length(obs_ord), function(i, xlst) {if(NCOL(xlst[[i]]) >= 1) {rep(obs_ord[i], NCOL(xlst[[i]]))}}, xlst = ddkeep_x_lst_ord)
#print (varlst_lst)
varlst <- do.call(base::c, varlst_lst)
#new
#print (varlst_0)
ux_added <- unique(varlst_0)
for(x_added in ux_added){
coefskeep_x_lst[[x_added]] <- xcoefs[which(varlst_0 %in% x_added)]
}
#test!?
ord_x <- order(varlst_0)
xcoefs <- xcoefs[ord_x]
} else {varlst <- NULL}
if(accz > 0){
varlst_z_lst <- lapply(1:length(obs_ord), function(i, zlst) {if(NCOL(zlst[[i]]) >= 1) {rep(obs_ord[i], NCOL(zlst[[i]]))}}, zlst = ddkeep_z_lst_ord)
varlst_z <- do.call(base::c, varlst_z_lst)
#new
uz_added <- unique(varlst_z_0)
for(z_added in uz_added){
coefskeep_z_lst[[z_added]] <- zcoefs[which(varlst_z_0 %in% z_added)]
}
#test!? correct; match order of columns in xmat
ord_z <- order(varlst_z_0)
zcoefs <- zcoefs[ord_z]
} else {varlst_z <- NULL}
#new
# totx <- NCOL(xmat)
# coefskeep_x_lst <- vector("list", length = totx)
# coefskeep_z_lst <- vector("list", length = totx)
# if(length(varlst_0) > 0){
# ux_added <- unique(varlst_0)
# for(x_added in ux_added){
# coefskeep_x_lst[[x_added]] <- xcoefs[which(varlst_0 %in% x_added)]
# }
# }
# if(length(varlst_z_0) > 0){
# uz_added <- unique(varlst_z_0)
# for(z_added in uz_added){
# coefskeep_z_lst[[z_added]] <- zcoefs[which(varlst_z_0 %in% z_added)]
# }
# }
#new: put b0 back
coefskeep <- c(coefskeep[1], zcoefs, xcoefs)
cic_st <- cic_val
if (infocrit == "CBIC") {
#degree is determined earlier
cic_st_next_step_0 <- llh + a * log(n) * (dfmean_cbic[ps_shpi_bic]) / n^(degree)
}
if (infocrit == "CAIC") {
#cic_val <- llh + 2 * (dfmean_caic) / n
cic_st_next_step_0 <- llh + 2 * (dfmean_cbic[ps_shpi_bic]) / n
}
if(cic_st_next_step_0 < oldcic){
cic_st_next_step <- cic_st_next_step_0
}
}
}
}
}
#new:
if(any(back_inactive_check)){
back_inactive <- chs[back_inactive_check]
}
rslt = list(ddkeep_x_lst = ddkeep_x_lst, ddkeep_z_lst = ddkeep_z_lst,
dmat_keep_lst = dmat_keep_lst, ps_keep_vec = ps_keep_vec,
knots_lst = knots_lst,
mskeep_lst = mskeep_lst,
mod_tried = mod_tried, changed = changed,
cic_st = cic_st, cic_st_next_step = cic_st_next_step,
back_inactive = back_inactive,
etahatkeep = etahatkeep, coefskeep = coefskeep,
coefskeep_x_lst = coefskeep_x_lst, coefskeep_z_lst = coefskeep_z_lst)
return (rslt)
}
#--------------------------------------------------------------------
.compute_edf0 = function(isim, gmat_lst, mu0_lst=NULL, np_ix_lst,
dmat_lst=NULL, weights=NULL,
wt.iter=FALSE, family=gaussian(), pnt=FALSE,
xtx_lst=NULL, uinv_lst=NULL,
xtxinv_lst=NULL, qv0_lst=NULL, dv_lst=NULL,
amat_lst=NULL, bvec_lst=NULL,
atil_lst=NULL, edfi_base_vec=NULL, shps_to_sims=NULL)
{
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
muhat.fun <- cicfamily$muhat.fun
ysim.fun <- cicfamily$ysim.fun
deriv.fun <- cicfamily$deriv.fun
dev.fun <- cicfamily$dev.fun
nobs <- length(gmat_lst)
#9:16, 18 #17 is in the list in the backward step
obs <- 1:nobs #9 or 10 shapes, including linear, will be fit in each isim?
eps <- 1e-8
edfs <- rep(0, nobs)
# ysim.fun <- function(n, mu0 = NULL, fml = object$family, shp0 = NULL) {
# if (fml == "binomial") {
# ysim <- 1:n*0
# ysim[runif(n) < .5] <- 1
# }
# if (fml == "poisson") {
# if (!is.null(mu0)) {
# ysim <- rpois(n, mu0)
# }
# }
# if (fml == "gaussian") {
# ysim <- rnorm(n)
# }
# if (fml == "Gamma") {
# ysim <- rgamma(n, shape=1)
# }
# ysim
# }
#new:
if(nobs > 1) {
if(is.null(shps_to_sims)){
if(nobs == 9) {
obs <- c(1, 3, 5, 9)
}
if(nobs == 8) {
obs <- c(1, 3, 5)
}
#new for bin and cbic
if(nobs == 4) {
obs <- c(1, 3)
}
if(nobs == 5) { #backward step for bin and cbic
obs <- c(1, 3, 5)
}
#if(nobs == 10) {
# obs <- c(1, 3, 5, 9, 10)
#}
} else {
if(nobs == 2 | nobs == 3){
obs <- 1:nobs
} else {
mono_ps <- mono_ps_other <- mono_ps_1 <- NULL
conv_ps <- conv_ps_other <- conv_ps_1 <- NULL
if(any(shps_to_sims %in% c(1, 2, 9, 10))){
mono_ps <- sort(which(shps_to_sims %in% c(1, 2, 9, 10)))
mono_ps_1 <- mono_ps[1]
if(length(mono_ps) > 1) {
mono_ps_other <- mono_ps[-1]
}
}
if(any(shps_to_sims > 2 & shps_to_sims < 9) | any(shps_to_sims > 10 & shps_to_sims < 17)){
conv_ps <- sort(which(shps_to_sims > 2 & shps_to_sims < 9 | shps_to_sims > 10 & shps_to_sims < 17))
conv_ps_1 <- conv_ps[1]
if(length(conv_ps) > 1){
conv_ps_other <- conv_ps[-1]
}
}
obs <- c(mono_ps_1, conv_ps_1)
}
}
}
for(i in obs){
#print(i)
gmat <- gmat_lst[[i]]
#mu0 <- mu0_lst[[i]]
#mu0 <- rep(1, NROW(gmat))
np_ix <- np_ix_lst[[i]]
#new
amat <- amat_lst[[i]]
bvec <- bvec_lst[[i]]
xtx <- xtx_lst[[i]]
xtxinv <- xtxinv_lst[[i]]
uinv <- uinv_lst[[i]]
qv0 <- qv0_lst[[i]]
atil <- atil_lst[[i]]
edfi_base <- edfi_base_vec[i]
#new
if(pnt){
dmat <- dmat_lst[[i]] #ignored np_ix
#kts <- knots_lst[[i]] #used for bin
dv <- dv_lst[[i]] #considered np_ix
} else {
dmat <- NULL
}
if(NCOL(gmat) == np_ix) {
#print ('no sim')
#if(!wt.iter){
edfi <- np_ix
#} else {
# edfi <- np_ix + 1 #include intercept and remove later
#}
}
if (NCOL(gmat) > np_ix) {
n <- NROW(gmat)
#test
#mu0 <- rep(0, n)
#if(!wt.iter){
# ysim <- mu0 + rnorm(n)
#ysim <- rnorm(n)
#} else {
ysim <- ysim.fun(n, mu0 = NULL, fml = family$family, shp0 = NULL)
#}
#-------------------------------------------
#use qprog
xtx <- qv0
if(pnt & !is.null(dmat)){
#test
dmat_zero <- matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat2 <- cbind(dmat_zero, dmat)
dv <- crossprod(dmat2)
xtx <- qv0 + dv
}
#-------------------------------------------
#if(np_ix > 0){
if(!wt.iter){
#ans <- suppressMessages(try(coneB(ysim, dsend, zsend), silent = TRUE))
#-------------------------------------------
#use qprog
#gmat0 = gmat
if (!pnt | pnt & all(dmat == 0)) {
zvec <- ysim
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
#test!
ans <- suppressMessages(try(coneB(zvec, dsend, zsend), silent = TRUE))
if (inherits(ans, "try-error")) {
break
} else {
face <- ans$face
bhat <- ans$coefs
bhat = round(bhat, 10)
}
} else {
xty <- crossprod(gmat, ysim)
#test
#amat is symmetric, no need to transpose it
ans <- try(solve.QP(xtx, xty, amat, bvec), silent = TRUE)
if (!inherits(ans, "try-error")) {
bhat <- ans$solution
bhat = round(bhat, 10)
face <- ans$iact
}
}
} else {
gmat0 <- gmat
#weights <- rep(1, n)
etahat <- etahat.fun(n, ysim, fml = family$family)
gr <- gr.fun(ysim, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(ysim, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, ysim, fml = family$family)
#muhat <- binomial(link = "logit")$linkinv(etahat)
muhat <- family$linkinv(etahat)
sm <- 1e-5
#sm <- 1e-7
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 10)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(ysim, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(ysim, etahat, n, weights=NULL, fml = family$family)
#ch = muhat > 1e-5 & muhat < 1-(1e-5)
#wt[ch] = muhat[ch]*(1-muhat[ch])
#wt[!ch] = 1e-5
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, ysim, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
#dsend <- gmat[, -c(1:np_ix), drop = FALSE]
#zsend <- gmat[, 1:np_ix, drop = FALSE]
#test!
#ans <- suppressMessages(try(coneB(zvec, dsend, cbind(onewt, zsend)), silent = TRUE))
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
#test!
ans <- suppressMessages(try(coneB(zvec, dsend, zsend), silent = TRUE))
if (inherits(ans, "try-error")) {
break
} else {
face <- ans$face
bhat <- ans$coefs
}
}
if(pnt & !is.null(dmat) & !all(dmat == 0)){
xtx <- crossprod(gmat)
xtx <- xtx + dv
#test!
xty <- crossprod(gmat, zvec)
ans <- try(solve.QP(xtx, xty, amat, bvec), silent = TRUE)
if (inherits(ans, "try-error")) {
break
} else {
bhat <- ans$solution
face <- ans$iact
}
}
bhat = round(bhat, 10)
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
#ans <- tryCatch({coneB(ysim, dsend, zsend)}, error = function(e) {e})
#}
if(inherits(ans, "try-error")) {
edfi <- NULL
} else {
#bhat <- ans$solution
#face <- ans$iact
if(!pnt | pnt & all(dmat == 0)){
edfi <- sum(abs(na.omit(bhat)) > eps)
#edfi <- sum(abs(na.omit(ans$solution)) > eps)
} else {
#new: if !wt.iter, use pre-computed edfi_base; otherwise, use pre-computed xtxinv
if(wt.iter | !wt.iter){
umat <- chol(xtx)
uinv <- backsolve(umat, diag(ncol(xtx)))
atil <- amat %*% uinv
xtxinv <- tcrossprod(uinv)
bigp_base <- xtxinv %*% crossprod(gmat)
edfi_base <- sum(diag(bigp_base))
}
if (all(face == 0) | length(face) == 0) {
edfi <- edfi_base
} else {
smat <- -t(atil[face, ,drop = FALSE])
#new: avoid solve error
qr_smat <- qr(smat)
smat <- qr.Q(qr_smat)
pmat_polar <- tcrossprod(smat)
if(!wt.iter){
bigp_polar <- uinv %*% tcrossprod(pmat_polar, uinv) %*% qv0
}else{
bigp_polar <- uinv %*% tcrossprod(pmat_polar, uinv) %*% crossprod(gmat)
}
edfi <- edfi_base - sum(diag(bigp_polar))
}
}
}
}
if(!is.null(edfi)){
edfs[i] <- edfi
}
}
#new:
if(nobs > 1){
if(nobs >= 4){
if(is.null(shps_to_sims)){
edfs[2] <- edfs[1]
edfs[4] <- edfs[3]
if(nobs > 5){
edfs[6:8] <- edfs[5]
}
} else {
if(!is.null(mono_ps_other)){
edfs[mono_ps_other] <- edfs[mono_ps_1]
}
if(!is.null(conv_ps_other)){
edfs[conv_ps_other] <- edfs[conv_ps_1]
}
}
}
}
return (edfs)
}
#--------------------------------------------------------------------
#simulate edf by lapply or mclapply
#--------------------------------------------------------------------
compute_edf <- function(family = gaussian(link="identity"), gmat_lst, dmat_lst=NULL,
mu0_lst=NULL, np_ix_lst, weights=NULL,
shp17_id=NULL, nsim=50, wt.iter=FALSE,
parallel=FALSE, pnt=FALSE, shps_to_sims=NULL){
#not faster
#totcores = detectCores()
#new:
#library(RhpcBLASctl)
# This is the "Magic Fix" for matrix math in parallel?
#blas_set_num_threads(1)
#omp_set_num_threads(1)
ng <- length(gmat_lst)
xtx_lst <- vector('list', length = ng)
xtxinv_lst <- vector('list', length = ng)
uinv_lst <- vector('list', length = ng)
amat_lst <- vector('list', length = ng)
bvec_lst <- vector('list', length = ng)
qv0_lst <- vector('list', length = ng)
dv_lst <- vector('list', length = ng)
atil_lst <- vector('list', length = ng)
edfi_base_vec <- numeric(ng)
for(i in 1:ng){
gmat <- gmat_lst[[i]]
np_ix <- np_ix_lst[[i]]
dmat <- dmat_lst[[i]]
nc <- ncol(gmat)
amat <- diag(nc)
ps_zero <- 1:np_ix
amat[ps_zero, ps_zero] <- 0
bvec <- cbind(rep(0, nc))
qv0 <- crossprod(gmat)
xtx <- qv0
if(pnt & !is.null(dmat)){
#test
dmat_zero <- matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat2 <- cbind(dmat_zero, dmat)
dv <- crossprod(dmat2)
dv_lst[[i]] <- dv #need this for binomial
#test, compute this term later
#xtx <- qv0 + dv
#xtx <- qv0
}
#test
#qv0 <- Matrix(qv0, sparse = TRUE)
qv0_lst[[i]] <- qv0
#test
#xtx <- Matrix(xtx, sparse = TRUE)
# calculate eigenvalues
ev <- eigen(xtx, only.values = TRUE)$values
# check if the smallest eigenvalue is too close to zero or negative
if (min(ev) < 1e-8) {
diag(xtx) <- diag(xtx) + 1e-8
}
xtx_lst[[i]] <- xtx
amat_lst[[i]] <- amat
bvec_lst[[i]] <- bvec
umat <- chol(xtx)
uinv <- backsolve(umat, diag(ncol(xtx)))
atil <- amat %*% uinv
uinv_lst[[i]] <- uinv
atil_lst[[i]] <- atil
xtxinv <- tcrossprod(uinv)
xtxinv_lst[[i]] <- xtxinv
if(!wt.iter){
bigp_base <- xtxinv %*% qv0
edfi_base <- sum(diag(bigp_base))
edfi_base_vec[i] <- edfi_base
}
}
if(parallel){
cores <- min(8, parallel::detectCores(logical = FALSE) - 1)
#test!
#cl <- makeCluster(cores)
#edfs <- parLapply(cl, 1:nsim, function(isim){
# .compute_edf0(isim, gmat_lst, mu0_lst=NULL, np_ix_lst, dmat_lst, weights=NULL, wt.iter, family, pnt)
#})
#stopCluster(cl)
edfs <- suppressWarnings(parallel::mclapply(1:nsim,
.compute_edf0, gmat_lst, mu0_lst=NULL,
np_ix_lst, dmat_lst,
weights=NULL, wt.iter,
family, pnt, xtx_lst, uinv_lst,
xtxinv_lst, qv0_lst, dv_lst,
amat_lst, bvec_lst, atil_lst,
edfi_base_vec, shps_to_sims, mc.cores = cores))
if(all(sapply(edfs, is.null))){
edfs <- parallel::mclapply(1:nsim, .compute_edf0, gmat_lst, mu0_lst=NULL,
np_ix_lst, dmat_lst,
weights=NULL, wt.iter,
family, pnt, xtx_lst, uinv_lst,
xtxinv_lst, qv0_lst, dv_lst,
amat_lst, bvec_lst, atil_lst,
edfi_base_vec, shps_to_sims, mc.cores = 1)
}
}else{
edfs <- lapply(1:nsim, .compute_edf0, gmat_lst, mu0_lst=NULL,
np_ix_lst, dmat_lst,
weights=NULL, wt.iter, family, pnt, xtx_lst, uinv_lst,
xtxinv_lst, qv0_lst, dv_lst, amat_lst, bvec_lst,
atil_lst, edfi_base_vec, shps_to_sims)
}
edfs <- Filter(Negate(is.null), edfs)
edfs <- simplify2array(edfs)
#print (edfs)
#edfs is a # of shapes by nsim matrix
#row averages will be used for CBIC
#col maxes' average will be used for CAIC
nobs <- length(gmat_lst)
if(nobs > 1){
if(is.null(shp17_id)){
dfmean_cbic <- rowMeans(edfs)
dfmean_caic <- max(dfmean_cbic)
} else {
dfmean_cbic <- rowMeans(edfs)
dfmean_caic_no17 <- max(dfmean_cbic[-shp17_id])
dfmean_caic_17 <- dfmean_cbic[shp17_id]
dfmean_caic <- c(dfmean_caic_no17, dfmean_caic_17)
}
#save(edfs, file='edfs.Rda')
rslt <- list(dfmean_caic = dfmean_caic, dfmean_cbic = dfmean_cbic)
} else {
dfmean <- mean(edfs)
rslt <- list(dfmean_caic = dfmean, dfmean_cbic = dfmean)
}
return (rslt)
}
#--------------------------------
#visualization
#--------------------------------
#plot.shapeselect <- function(object, pages = 1, ci = FALSE, mn = TRUE, mu = FALSE,...)
plot.shapeselect <- function(x,...)
{
object <- x
extras <- list(...)
if(exists("pages", where = extras)){
pages <- extras$pages
} else {
pages <- 1
}
if(exists("ci", where = extras)){
ci <- extras$ci
} else {
ci <- FALSE
}
if(exists("mn", where = extras)){
mn <- extras$mn
} else {
mn <- TRUE
}
if(exists("mu", where = extras)){
mu <- extras$mu
} else {
mu <- FALSE
}
layout_for_cgam <- function(n_plots, pages = 1) {
ppp <- ceiling(n_plots / max(pages, 1L))
nr <- floor(sqrt(ppp))
nc <- ceiling(ppp / nr)
return (c(nr, nc))
}
xmat <- object$xmat
p <- ncol(xmat)
types <- object$types
for(ix in 1:p){
if(inherits(xmat[,ix], "factor") & types[ix] %in% "tree"){
xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
}
if(inherits(xmat[,ix], "factor") & types[ix] %in% "ord"){
xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
}
}
family <- object$family
etacomps <- object$etacomps
#theme_use <- ggthemes::theme_solarized()
#theme_use <- ggthemes::theme_economist()
shps <- object$shps
xnms <- object$xnms
not_selected <- xnms[which(shps %in% c(0, 17))]
zmat <- object$zmat
np <- object$d0
capk <- object$capk
varlist <- object$varlst #used varlst in classo
varlist <- sort(varlist) #follow xmat order, not the x's added order
varlist_z <- object$varlst_z
varlist_z <- sort(varlist_z)
#dd <- object$dd #include 1
dd <- t(object$bigmat) #only contains basis functions of selected predictors
#if(onlySelect){
#kept <- !shps %in% c(0, 1, 17) #remove z
kept <- !shps %in% c(0, 17)
xmat <- xmat[, kept, drop = FALSE]
types <- types[kept]
etacomps <- etacomps[kept, ,drop = FALSE]
shps <- shps[kept]
xnms <- xnms[kept]
#}
p <- ncol(xmat)
if(pages == 1){
oldpar <- par(mfrow = layout_for_cgam(n_plots = p, pages = pages))
}
if(ci){
cov_beta <- vcov.cgam(object)
} else {
cov_beta <- NULL
}
extras <- list(...)
uvar_x <- unique(varlist)
uvar_z <- unique(varlist_z)
ncon <- 1 #for conv/conc
ziter <- 1 #for z
#types <- object$types
shapes_char <- object$shapes[kept]
#p <- ncol(xmat) #need to define it again! only use chosen xs
for (k in 1:p) {
x_k <- xmat[, k]
xnm_k <- xnms[k]
#shp_k <- ShapeToChar_shapeselect(shps[k], types[k])
shp_k <- shapes_char[k]
eta_k <- etacomps[k, ]
if(mu){
eta_k <- family$linkinv(eta_k)
}
#ps_k <- which(varlist == k)
if(is.factor(x_k)){
#x_k <- as.numeric(levels(x_k))[x_k]
ps_k <- which(varlist_z == uvar_z[ziter]) + 1 #consider constant
ziter <- ziter + 1
} else {
ps_k <- which(varlist == uvar_x[k])
ps_k <- ps_k + np
}
#range_eta <- range(eta_k, na.rm = TRUE)
#if(!is.null(cov_beta) & !xnm_k %in% not_selected){
if(!is.null(cov_beta)){
if(shps[k] > 2 & shps[k] < 5 | shps[k] > 10 & shps[k] < 13){
ps_k <- c(1 + capk + ncon, ps_k) #constant, nominal, conv/conc
ncon <- ncon + 1
}
cov_k <- cov_beta[ps_k, ps_k]
dd_k <- dd[, ps_k, drop = FALSE]
sd_eta_k <- diag(dd_k %*% cov_k %*% t(dd_k))^.5
lwr_eta_k <- eta_k - 2 * sd_eta_k
upp_eta_k <- eta_k + 2 * sd_eta_k
if(mu){
lwr_eta_k <- family$linkinv(lwr_eta_k)
upp_eta_k <- family$linkinv(upp_eta_k)
}
range_eta <- range(lwr_eta_k, upp_eta_k, na.rm = TRUE)
} else {
range_eta <- range(eta_k, na.rm = TRUE)
}
padding <- 0.05 * diff(range_eta)
ylim_use <- c(range_eta[1] - padding, range_eta[2] + padding)
#--------------------------
#vis each component
ord_k <- order(x_k)
tp_k <- ifelse(is.factor(x_k), 'p', 'l')
xaxt_k <- ifelse(is.factor(x_k), 'n', 's')
x_k_tmp <- x_k
if(is.factor(x_k_tmp)){
#x_k_tmp <- as.numeric(levels(x_k))[x_k]
x_k_tmp <- as.integer(x_k) #- 1
}
if(mn){
plot(x_k_tmp[ord_k], eta_k[ord_k], xlab = bquote(~ .(xnm_k)), ylab = bquote(hat(eta)[.(k)]),
ylim = ylim_use,
main = bquote(~ hat(eta)[.(k)] ~ "vs" ~ .(xnm_k) ~ "(" * .(shp_k) * ")"),
#main = bquote(~ .(shp_k)),
xaxt = xaxt_k,
type = tp_k, lwd = 1)
}else{
plot(x_k_tmp[ord_k], eta_k[ord_k], xlab = bquote(~ .(xnm_k)), ylab = "",
ylim = ylim_use,
#main = bquote(~ hat(eta)[.(k)] ~ "vs" ~ .(xnm_k) ~ "(" * .(shp_k) * ")"),
main = NULL,
xaxt = xaxt_k,
type = tp_k, lwd = 1)
}
if(!exists("rug", where = extras)){
rug <- TRUE
} else {
rug <- extras$rug
}
if(rug){
rug(x_k, ticksize = 0.03, side = 1)
}
#if(!is.null(cov_beta) & !xnm_k %in% not_selected){
if(!is.null(cov_beta) & !is.factor(x_k)){
lines(x_k[ord_k], lwr_eta_k[ord_k], col = 1, lty = 2)
lines(x_k[ord_k], upp_eta_k[ord_k], col = 1, lty = 2)
}
#print (attr(x_k, 'type') == 'tree')
if(is.factor(x_k) | types[k] == 'tree'){
points(x_k_tmp, eta_k, pch=20)
if(is.factor(x_k)){
lines(x_k_tmp[ord_k], eta_k[ord_k], lty = 2)
axis(1, at = unique(x_k_tmp), labels = unique(levels(x_k))) # custom
}
}
if(pages > 1){
readline("Press <Enter> to continue...")
}
}
if(pages == 1){
on.exit(par(oldpar)) # restores original settings
}
}
#---------------------------------------------------------------------------------------------------------------------------
predict.shapeselect <- function(object, newdata, interval = c("none", "confidence", "prediction"),
type = c("response", "link"), level = 0.95, n.mix = 100,...)
{
interval <- match.arg(interval)
type <- match.arg(type)
family <- object$family
space <- object$space
pnt <- object$pnt
center <- object$center
standardize <- object$standardize
boundary.knots <- object$boundary.knots
#cicfamily <- CicFamily(family)
#muhat.fun <- cicfamily$muhat.fun
if (!inherits(object, "classo")) {
warning("calling predict.classo(<fake-classo-object>) ...")
}
if (missing(newdata) | is.null(newdata)) {
etahat <- object$etahat
muhat <- object$muhat
#muhat <- muhat.fun(etahat, fml = family$family)
#muhat <- family$linkinv(etahat)
ans <- list(fit = muhat, etahat = etahat)
return (ans)
}
if (!is.data.frame(newdata)) {
stop ("newdata must be a data frame!")
}
if(interval != "none"){
cmat <- vcov.cgam(object)
} else {
cmat <- NULL
}
#xmatpr = object$xmat
#match newdata by the xnms
knots_lst <- object$knots_lst
mskeep_lst <- object$mskeep_lst
cvec <- object$cvec
sdx <- object$sdx
mx <- object$mx
coefskeep_withone <- object$coefskeep_withone #will lead to wrong prediction when we have conc/conv; bigmat is not ordered the same as cbind(one, ddkeepz_all, ddkeepx_all)
coefs_x_lst <- object$coefskeep_x_lst
coefs_z_lst <- object$coefskeep_z_lst
#temp
xnms <- object$xnms
shps <- object$shps
seq_added <- object$seq_added
xnms_kept <- xnms[which(!shps %in% c(0))]
#if(!all(colnames(newdata) %in% xnms_kept) | !all(xnms_kept %in% colnames(newdata))){
if(!all(names(newdata) %in% xnms)){
stop ("newdata must include at least all selected predictors!")
} else {
if(NCOL(newdata) == length(xnms_kept)){
newdata <- newdata[, xnms_kept, drop = F]
} else if(NCOL(newdata) == length(xnms)){
newdata <- newdata[, xnms, drop = F]
}
}
#######################
#local helper function#
#######################
my_line <- function(xp = NULL, y, x, end, start) {
slope <- NULL
intercept <- NULL
yp <- NULL
slope <- (y[end] - y[start]) / (x[end] - x[start])
intercept <- y[end] - slope * x[end]
yp <- intercept + slope * xp
ans <- new.env()
ans$slope <- slope
ans$intercept <- intercept
ans$yp <- yp
ans
}
if(any(shps > 0)){
xmat <- object$xmat #not scaled
types <- object$types
#xmat0 <- xmat
n <- nrow(xmat)
p <- length(xnms)
for(i in 1:p){
if(is.numeric(xmat[,i]) & types[i] %in% c('num', 'numeric')){
#if(center){
xmat[,i] <- (xmat[,i] - mx[i]) / sdx[i]
newdata[,i] <- (newdata[,i] - mx[i]) / sdx[i]
#} else {
# xmat[,i] <- xmat[,i] / sdx[i]
# newdata[,i] <- newdata[,i] / sdx[i]
#}
#new: scale user-defined boundary knots
# if(is.list(boundary.knots)){
# kts0 <- boundary.knots[[i]]
# }
#
# if(length(kts0) >= 2){
# #if(center){
# kts0 <- (kts0 - mx[i]) / sdx[i]
# #} else {
# # kts0 <- kts0 / sdx[i]
# #}
# }
#xmat[,i] <- xmat[,i] / sdx[i]
#xmat[,ix] <- scale(xmat[,ix], center=center)
}
#test
if(inherits(xmat[,i], "factor") & types[i] %in% "tree"){
xmat[,i] <- as.numeric(as.character(xmat[,i]))
newdata[,i] <- as.numeric(as.character(newdata[,i]))
}
if(inherits(xmat[,i], "factor") & types[i] %in% "ord"){
xmat[,i] <- as.numeric(as.character(xmat[,i]))
newdata[,i] <- as.numeric(as.character(newdata[,i]))
}
}
#print (newdata)
#create xmatpr for new data
zmat <- NULL
xcmat <- NULL #x's when shp is conv or conc
delta <- NULL
ddkeep_x_lst <- vector("list", length = p)
ddkeep_z_lst <- vector("list", length = p)
one <- rep(1, nrow(newdata))
varlst <- varlst_z <- NULL
#shps_added <- shps[seq_added]
#xnms_added <- xnms[seq_added]
for(ix in 1:p){
#print (ix)
type_ix <- types[ix]
shpi <- shps[ix]
if(shpi %in% c(0, 17)){
next
}
x <- newdata[, ix]
ms <- mskeep_lst[[ix]]
if(inherits(x, "factor")){
#dd <- model.matrix(~ x)[, -1, drop=FALSE]
#dd <- apply(dd, 2, function(e) e - mean(e))
lvli <- levels(x)
xui <- unique(x)
ztbi <- levels(xmat[, ix])
dd <- NULL
zcoefid <- NULL
if (!any(lvli %in% ztbi)) {
stop ("new factor level must be among factor levels in the fit!")
} else {
#baseline: ztbi[1]
#if(any(lvli != ztbi[1])){
klvls <- length(ztbi)
rn <- nrow(newdata)
#remove baseline?
dd <- matrix(0, nrow = rn, ncol = (klvls - 1))
#new: test!
#if(nrow(ms) == 1) {
for (i1 in 1:rn) {
if (x[i1] != ztbi[1]) {
id_col <- which(ztbi %in% x[i1]) - 1
#print (id_col)
dd[i1, id_col] <- 1
dd[i1, id_col] <- dd[i1, id_col] - ms[id_col] #ms[i1, id_col]
}
}
#temp
for(ic in 1:ncol(dd)){
ps_base <- which(dd[, ic] == 0)
dd[ps_base, ic] <- dd[ps_base, ic] - ms[ic] #ms[ps_base, ic]
}
#}
# if(nrow(ms) > 1) {
# oldx <- xmat[, ix]
# ddold <- model.matrix(~ oldx)[, -1, drop=FALSE]
# #xs <- sort(oldx)
# #ord <- order(oldx)
# nr <- nrow(dd)
# #nc <- length(x)
# nc <- ncol(dd)
# ms2 <- matrix(0, nrow = nr, ncol = nc)
# for (i1 in 1:nc) {
# xs <- sort(ddold[, i1])
# ord <- order(ddold[, i1])
# for (i2 in 1:nr) {
# ms2[i2, i1] <- my_line(xp = x[i1], y = ms[i2, ][ord], x = xs, end = n, start = 1)$yp
# }
# }
#
# }
zmat <- cbind(zmat, dd)
#test
ddkeep_z_lst[[ix]] <- dd
varlst_z <- c(varlst_z, rep(ix, NCOL(dd)))
}
#}
} else {
oldx <- xmat[, ix]
#if(any(x > max(oldx)) | any(x < min(oldx))){
# stop ("No extrapolation is allowed in classo prediction!")
#}
shpi <- shps[ix]
ms <- mskeep_lst[[ix]]
if(type_ix %in% c("num", "numeric")){
kts <- knots_lst[[ix]]
# if(type_ix %in% c("ord", "ordinal")){
# kts <- 0
# }
#print (kts)
dd_ans <- suppressWarnings(
makedelta(x, sh = shpi, numknots = 0, knots = kts, space = space, suppre = TRUE, interp = TRUE))
dd <- dd_ans$amat
#cat('finished delta')
if(shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5){
ms <- t(ms)
}
#need to make basis orthogonal to one or [1, x]
if (shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5) {
xs <- sort(oldx)
ord <- order(oldx)
nr <- nrow(dd)
nc <- length(x)
ms2 <- matrix(0, nrow = nr, ncol = nc)
for (i1 in 1:nc) {
for (i2 in 1:nr) {
ms2[i2, i1] <- my_line(xp = x[i1], y = ms[i2, ][ord], x = xs, end = n, start = 1)$yp
}
}
dd <- dd - ms2
} else {
dd <- dd - ms
}
delta <- cbind(delta, t(dd))
}
if(type_ix %in% c("ord", "ordinal")){
dd <- pred_del(oldx, shpi, x, ms)
}
if(type_ix %in% c("tree")){
#easier than predict.cgam
uoldx <- unique(oldx)
#pl <- uoldx[1]
if (is.numeric(oldx)) {
if (0 %in% oldx) {
pl <- 0
} else {pl <- min(oldx)}
} else {
pl <- uoldx[1]
}
nx <- uoldx[uoldx != pl]
is_pl <- x %in% pl
#use the order of unique oldx to determine the treatment level of new: x
dd <- matrix(0, nrow = length(nx), ncol = length(x))
#if new x only has placebo, then just center each row of dd by ms
for(i in 1:nrow(dd)){
xi <- nx[i]
dd[i, ] <- dd[i, ] - ms[which(nx %in% xi)]
}
if(any(!is_pl)){
for(i in 1:nrow(dd)){
xi <- nx[i]
if(any(x %in% xi)){
dd[i, which(x %in% xi)] <- 1 - ms[which(nx %in% xi)]
}
}
}
}
if(shpi == 2 & type_ix %in% c("tree")) {
ddkeep_z_lst[[ix]] <- t(dd)
varlst_z <- c(varlst_z, rep(ix, NROW(dd)))
} else {
ddkeep_x_lst[[ix]] <- t(dd)
varlst <- c(varlst, rep(ix, NROW(dd)))
}
if(shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5){
#zmat <- cbind(zmat, x - mean(oldx))
xcmat <- cbind(xcmat, x - mean(oldx))
ddkeep_z_lst[[ix]] <- x - mean(oldx)
varlst_z <- c(varlst_z, rep(ix, 1))
}
}
}
#temp
# if(shrinketa){
# shps_added <- shps[seq_added]
# xnms_added <- xnms[seq_added]
# for(i in 1:length(seq_added)){
# ps <- seq_added[i]
# xnmi <- xnms_added[i]
# shpi <- shps_added[i]
# if(shpi == 0){
# next
# }
# if(shpi == 1 | shpi == 2 & type_ix %in% c("tree")) {
# coefs_z_lst[[ps]] <- coefs_z_lst[[ps]] * cvec[ps]
# } else {
# coefs_x_lst[[ps]] <- coefs_x_lst[[ps]] * cvec[ps]
# if(shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5){
# coefs_z_lst[[ps]] <- coefs_z_lst[[ps]] * cvec[ps]
# }
# }
# }
# }
#xmatpr <- cbind(one, zmat, xcmat, delta)
#print (dim(xmatpr))
#get new fit
#fit <- xmatpr %*% coefskeep_withone[1:ncol(xmatpr)]
#coefsuse <- c(onecoef, zcoefs_nox, xccoefs, xcoefs)
#ddkeep_x_lst_ord <- ddkeep_x_lst[seq_added]
#ddkeep_z_lst_ord <- ddkeep_z_lst[seq_added]
delta <- do.call(base::cbind, ddkeep_x_lst)
zmat <- do.call(base::cbind, ddkeep_z_lst) #include conv/conc
xmatpr <- cbind(one, zmat, delta)
#coefs_x_lst_ord <- coefs_x_lst[seq_added]
#coefs_z_lst_ord <- coefs_z_lst[seq_added]
coefsuse_z <- do.call(base::c, coefs_z_lst)
coefsuse_x <- do.call(base::c, coefs_x_lst)
coefsuse <- c(coefskeep_withone[1], coefsuse_z, coefsuse_x)
fit <- xmatpr %*% coefsuse
#fit <- switch(type, link = fit, response = object$family$linkinv(fit))
etacomps <- matrix(0, nrow = p, ncol = length(fit))
etacomps_lwr <- etacomps_upp <- etacomps
#shps_added <- shps[seq_added]
#xnms_added <- xnms[seq_added]
if(!is.null(cmat)){
cmat_noone <- cmat[2:nrow(cmat),2:nrow(cmat)]
if(NCOL(zmat) > 0){
ps_z <- 1:NCOL(zmat)
cmat_z <- cmat_noone[ps_z, ps_z] #include conv/conc
cmat_x <- cmat_noone[-ps_z, -ps_z]
} else {
cmat_x <- cmat_noone
}
} else {
cmat_z <- cmat_x <- NULL
}
for(k in 1:p){
#ps <- seq_added[i]
# xnmi <- xnms_added[i]
# shpi <- shps_added[i]
shpi <- shps[k]
if(shpi == 0){#flat
next
} else {
if (shpi == 1 | shpi == 2) {#need to use 1 for incr later
#etacomps[ps, ] <- ddkeep_z_lst[[ps]] %*% coefs_z_lst[[ps]]
if(types[k] %in% c('nom', 'nominal') | shpi == 2 & types[k] %in% 'tree'){
etacomps[k, ] <- ddkeep_z_lst[[k]] %*% coefs_z_lst[[k]]
if(!is.null(cmat_z)){
ps_k <- which(varlst_z == k)
dd_k <- zmat[, ps_k, drop = FALSE]
cov_k <- cmat_z[ps_k, ps_k]
}
}
if(types[k] %in% c('ord', 'ordinal') | shpi == 1 & types[k] %in% 'tree'){
etacomps[k, ] <- ddkeep_x_lst[[k]] %*% coefs_x_lst[[k]]
if(!is.null(cmat_x)){
ps_k <- which(varlst == k)
dd_k <- delta[, ps_k, drop = FALSE]
cov_k <- cmat_x[ps_k, ps_k]
}
}
} else if (shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5) {
etacomps[k, ] <- ddkeep_x_lst[[k]] %*% coefs_x_lst[[k]]
etacomps[k, ] <- etacomps[k, ] + ddkeep_z_lst[[k]] * coefs_z_lst[[k]]
if(!is.null(cmat)){
ps_k <- c(which(varlst_z == k), which(varlst == k))
dd_k <- cbind(zmat[, which(varlst_z == k), drop = FALSE], delta[, which(varlst == k), drop = FALSE])
cov_k <- cmat_noone[ps_k, ps_k]
}
} else {
etacomps[k, ] <- ddkeep_x_lst[[k]] %*% coefs_x_lst[[k]]
if(!is.null(cmat_x)){
ps_k <- which(varlst == k)
dd_k <- delta[, ps_k, drop = FALSE]
cov_k <- cmat_x[ps_k, ps_k]
}
}
if(!is.null(cmat)){
eta_k <- etacomps[k, ]
#ps_k <- which(varlst == k)
#dd_k <- xmat[obs, ps_k, drop = FALSE]
#cov_k <- cov_beta[ps_k, ps_k]
sd_eta_k <- diag(dd_k %*% cov_k %*% t(dd_k))^.5
mult <- qnorm((1 - level)/2, lower.tail=FALSE)
lwr_eta_k <- eta_k - mult * sd_eta_k
upp_eta_k <- eta_k + mult * sd_eta_k
etacomps_lwr[k, ] <- lwr_eta_k
etacomps_upp[k, ] <- upp_eta_k
}
}
}
if(interval != "none"){
coveta <- diag(xmatpr %*% cmat %*% t(xmatpr))
se.fit <- sqrt(coveta)
# new: C.I. level
mult <- qnorm((1 - level)/2, lower.tail=FALSE)
upp <- fit + mult*se.fit
lwr <- fit - mult*se.fit
fit <- switch(type, link = fit, response = object$family$linkinv(fit))
lwr <- switch(type, link = lwr, response = object$family$linkinv(lwr))
upp <- switch(type, link = upp, response = object$family$linkinv(upp))
#se.fit <- switch(type, link = se.fit, response = object$family$linkinv(se.fit))
ans <- list(fit = fit, lwr = lwr, upp = upp, se.fit = se.fit, xmatpr = xmatpr, newddx = delta, newddz = zmat,
etacomps = etacomps, etacomps_lwr = etacomps_lwr, etacomps_upp = etacomps_upp)
} else {
#print (type)
fit <- switch(type, link = fit, response = object$family$linkinv(fit))
ans <- list(fit = fit, xmatpr = xmatpr, newddx = delta, newddz = zmat, etacomps = etacomps)
}
}
if(all(shps == 0)){
y <- object$y
m0 <- glm(y ~ 1, family = family$family)
fit <- predict(m0, newdata)
xmatpr <- predict(m0, newdata, type = 'terms')
ans <- list(fit = fit, xmatpr = xmatpr)
}
return (ans)
}
#############################
#predict delta for ordinal x#
#############################
# pred_del = function(x, sh, xp, ms) {
# n = length(xp)
# #x = (x - min(x)) / (max(x) - min(x))
# xu = sort(unique(x))
# n1 = length(xu)
# sigma = NULL
# #######################
# #local helper function#
# #######################
# my_line = function(xp = NULL, y, x, end, start) {
# slope = NULL
# intercept = NULL
# yp = NULL
# slope = (y[end] - y[start]) / (x[end] - x[start])
# intercept = y[end] - slope * x[end]
# yp = intercept + slope * xp
# ans = new.env()
# ans$slope = slope
# ans$intercept = intercept
# ans$yp = yp
# ans
# }
# # increasing or decreasing
# if (sh < 3) {
# sigma = matrix(0, nrow = n1 - 1, ncol = n)
# for (i in 1: (n1 - 1)) {
# sigma[i, xp > xu[i]] = 1
# }
# if (sh == 2) {sigma = -sigma}
# for (i in 1:(n1 - 1)) {sigma[i, ] = sigma[i, ] - ms[i]}
# }
# if (sh == 3 | sh == 4) {
# # convex or concave
# sigma = matrix(0, nrow = n1 - 2, ncol = n)
# #for (i in 1: (n1 - 2)) {
# # sigma[i, x > xu[i]] = x[x > xu[i]] - xu[i]
# #}
# for (i in 1: (n1 - 2)) {
# sigma[i, xp > xu[i+1]] = xp[xp > xu[i+1]] - xu[i+1]
# }
# if (sh == 4) {sigma = -sigma}
# #xm = cbind(1:n*0+1, xp)
# #xpx = solve(t(xm) %*% xm)
# #pm = xm %*% xpx %*% t(xm)
# #sigma = sigma - sigma %*% t(pm)
# xs = sort(x)
# ord = order(x)
# nx = length(x)
# obs = 1:nx
# m = nrow(ms)
# ms0 = matrix(0, nrow = m, ncol = n)
# for (i1 in 1:n) {
# for (i2 in 1:m) {
# ms0[i2, i1] = my_line(xp = xp[i1], y = ms[i2, ][ord], x = xs, end = nx, start = 1)$yp
# }
# }
# sigma = sigma - ms0
# }
# if (sh > 4 & sh < 9) {
# sigma = matrix(0, nrow = n1 - 1, ncol = n)
# if (sh == 5) { ### increasing convex
# for (i in 1:(n1 - 1)) {
# sigma[i, xp > xu[i]] = (xp[xp > xu[i]] - xu[i]) / (max(x) - xu[i])
# }
# for (i in 1:(n1 - 1)) {sigma[i,] = sigma[i,] - ms[i]}
# } else if (sh == 6) { ## decreasing convex
# for (i in 1:(n1 - 1)) {
# sigma[i, xp < xu[i + 1]] = (xp[xp < xu[i + 1]] - xu[i + 1]) / (min(x) - xu[i + 1])
# }
# for (i in 1:(n1 - 1)) {sigma[i,] = sigma[i,] - ms[i]}
# #print (ms)
# } else if (sh == 7) { ## increasing concave
# for (i in 1:(n1 - 1)) {
# sigma[i, xp < xu[i + 1]] = (xp[xp < xu[i + 1]] - xu[i + 1]) / (min(x) - xu[i + 1])
# }
# for (i in 1:(n1 - 1)) {sigma[i,] = -sigma[i,] + ms[i]}
# } else if (sh == 8) {## decreasing concave
# for (i in 1:(n1 - 1)) {
# sigma[i, xp > xu[i]] = (xp[xp > xu[i]] - xu[i]) / (max(x) - xu[i])
# }
# for (i in 1:(n1 - 1)) {sigma[i,] = -sigma[i,] + ms[i]}
# }
# }
# return (sigma)
# }
#--------------------------------------------------------------------------------
# zvec.fun <- function(cvec = NULL, wt = NULL, y, sm = 1e-7, fml = object$family) {
# n <- length(y)
# if (fml == "gaussian") {
# #zvec <- y
# zvec <- wt^(1/2) * y
# }
# if (fml == "poisson") {
# #zvec <- cvec / wt
# zvec <- cvec / sqrt(wt)
# }
# if (fml == "binomial") {
# zvec = 1:n*0
# zvec[wt == 0] <- 1 / sm
# zvec[wt > 0] <- cvec[wt > 0] / sqrt(wt[wt > 0])
# }
# if (fml == "Gamma") {
# zvec <- cvec / sqrt(wt)
# }
# zvec
# }
# wt.fun <- function(y, etahat = NULL, n = NULL, weights = NULL, fml = object$family){
# if (is.null(weights)) {
# weights <- 1:n*0 + 1
# }
# w <- weights
# if (fml == "poisson") {
# wt <- w * exp(etahat)
# }
# if (fml == "binomial") {
# if (all(etahat == 0)){
# #wt <- 1:n*0 + 1/4
# wt <- w * (1:n*0 + 1/4)
# } else {
# wt <- 1:n*0
# wt <- (1 - 1 / (1 + exp(etahat))) * (1 / (1 + exp(etahat))) #pi_ij(1 - pi_ij)
# wt <- c(wt)
# # for (i in 1:n) {
# # if (etahat[i] > 100) {
# # wt[i] <- 0
# # } else {
# # wt[i] <- w[i] * exp(etahat[i]) / ((1 + exp(etahat[i]))^2)
# # }
# # }
# }
# }
# if (fml == "gaussian") {
# wt <- w # (1:n*0 + 1) / w
# }
# if (fml == "Gamma") {
# wt <- w * y * exp(-etahat)
# }
# wt <- as.vector(wt)
# wt
# }
#
# etahat.fun <- function(n, y, fml = object$family){
# if (fml == "poisson") {
# etahat <- 1:n*0 + log(mean(y))
# }
# if (fml == "binomial") {
# etahat <- 1:n*0
# }
# if (fml == "Gamma") {
# etahat <- 1:n*0 + log(mean(y))
# }
# etahat
# }
# gr.fun <- function(y, etahat = NULL, weights = NULL, fml = object$family){
# n <- length(y)
# if (is.null(weights)) {
# weights <- 1:n*0 + 1
# }
# w <- weights
# if (fml == "poisson") {
# gr <- w * (exp(etahat) - y)
# }
# if (fml == "binomial") {
# if (all(etahat == 0)) {
# gr <- w * (1/2 - y)
# } else {
# gr <- 1:n*0
# #ignore weights
# gr <- 1 - 1 / (1 + exp(etahat)) - y
#
# # for (i in 1:n) {
# # if (etahat[i] > 100) {
# # gr[i] <- w[i] * (1 - y[i])
# # } else {gr[i] <- w[i] * (exp(etahat[i]) / (1 + exp(etahat[i])) - y[i])}
# # }
# }
# }
# if (fml == "Gamma") {
# gr <- w * (1 - y * exp(-etahat))
# }
# gr
# }
# ysim.fun <- function(n, mu0 = NULL, fml = object$family, shp0 = NULL) {
# if (fml == "binomial") {
# ysim <- 1:n*0
# ysim[runif(n) < .5] <- 1
# }
# if (fml == "poisson") {
# if (!is.null(mu0)) {
# ysim <- rpois(n, mu0)
# }
# }
# if (fml == "gaussian") {
# ysim <- rnorm(n)
# }
# if (fml == "Gamma") {
# ysim <- rgamma(n, shape=1)
# }
# ysim
# }
# llh.fun <- function(y, muhat = NULL, etahat = NULL, phihat = NULL, n = NULL, weights = NULL, fml = object$family){
# #sm <- 1e-7
# sm <- 1e-5
# if (is.null(weights)) {
# weights <- 1:n*0 + 1
# }
# w <- weights
# #new: avoid Inf
# if (fml == "poisson") {
# llh <- 2 * sum(w[w!=0] * (muhat[w!=0] - y[w!=0] * etahat[w!=0])) / n
# }
# if (fml == "binomial") {
# llh <- 0
# if (all(0 <= y) & all(y <= 1)) {
# for (i in 1:n) {
# if (muhat[i] > 0 & muhat[i] < 1) {
# llh <- llh + w[i] * (y[i] * log(muhat[i]) + (1 - y[i]) * log(1 - muhat[i]))
# }
# }
# #ch <- muhat > 0 & muhat < 1
# #llh <- sum(w[ch] * (log(1 - muhat[ch]) - y[ch] * etahat[ch]))
# llh <- (-2/n) * llh
# } else {
# stop ("y values must be 0 <= y <= 1!")
# }
# }
# if (fml == "gaussian") {
# if (all(w == 1)) {
# llh <- log(sum((y - etahat)^2))
# } else {
# llh <- log(sum(w[w!=0] * (y[w!=0] - etahat[w!=0])^2)) - sum(log(w[w!=0])) / n
# }
# }
# if (fml == "Gamma") {
# vuhat <- 1 / phihat
# #print (vuhat)
# #llh <- 2 * sum(w[w!=0] * (etahat[w!=0] + y[w!=0] * exp(-etahat[w!=0]))) / n
# llh <- 2 / n * (vuhat * sum(w[w!=0] * (etahat[w!=0] + y[w!=0] * exp(-etahat[w!=0]))) + n * (log(gamma(vuhat)) - vuhat * log(vuhat)) - (vuhat-1) * sum(log(y)))
# #print (llh)
# }
# llh
# }
.search_ps_shapeselect = function(pen, qv0, dv, dd, edfu, covmat_inv=NULL) {
# 1. Calculate qv
qv = qv0 + pen * dv
# 2. Stronger Regularization
# If 1e-10 failed, your matrix is extremely sensitive.
# Let's use a slightly larger nudge or a relative nudge.
diag(qv) <- diag(qv) + (max(diag(qv)) * 1e-12) + 1e-8
#diag(qv) <- diag(qv) + 1e-10
# 3. Use tryCatch for the whole block to catch "Singular" errors
val <- tryCatch({
umat = chol(qv)
# Efficiently calculate the trace:
# instead of solve(qv, t(dd)), use backsolve/forsolve with the Cholesky factor
# This is numerically much more stable
# Let A = dd %*% qv_inv %*% t(dd). We want trace(A).
tmp <- backsolve(umat, t(dd), transpose = TRUE)
trace_val <- sum(tmp^2)
return(trace_val - .8 * edfu)
}, error = function(e) {
# If the math breaks, return a value that tells uniroot to keep looking
# Usually a large positive value if we are looking for a root from above
return(1e10)
})
return(val)
}
#--------------------------------------------------------------------
#cross validation to find the c vector
#remember to acknowledge
#--------------------------------------------------------------------
create_folds <- function(n_data, nfolds = 5){
folds <- rep(list(numeric(0)), nfolds)
available_samples <- 1:n_data
fold_size <- floor(n_data/nfolds)
if(fold_size == n_data/nfolds){
for(fold_id in 1:nfolds){
folds[[fold_id]] <- sample(available_samples, fold_size)
available_samples <- available_samples[!(available_samples %in% folds[[fold_id]])]
}
} else{
for(fold_id in 1:(nfolds-1)){
folds[[fold_id]] <- sample(available_samples, min(fold_size + rbinom(1, 1, 0.5), length(available_samples)))
available_samples <- available_samples[!(available_samples %in% folds[[fold_id]])]
}
folds[[nfolds]] <- available_samples
}
return(folds)
}
#------------------------------------------------------------------------------------------
cv.classo = function(fmat, y, nfolds = 5, type = c("glmnet", "solve.qp"), family=gaussian()) {
#sdxs = apply(fmat, 2, sd) #used later to go back to the original scale
#mxs = apply(fmat, 2, mean)
#fmat = scale(fmat, center = FALSE, scale = TRUE) #scale new x
type = match.arg(type)
#temp:
if(NCOL(fmat) == 1) {
#type = "solve.qp"
cvec = 1
rslt = list(cvec = cvec)
#return (rslt)
} else {
if(type == "solve.qp") {
n = length(y)
#vec = 1:n
#shuffle the observations
#vec_sample = sample(vec, size = n, replace = FALSE)
#split this data into k folds with roughly equal sample sizes
#folds = split(vec_sample, 1:nfolds)
folds = create_folds(n, nfolds)
p = NCOL(fmat)
amat = diag(p)
#new: each c value <= 1
amat = rbind(amat, -diag(p))
amat = rbind(amat, c(rep(-1, p)))
#amat = cbind(0, amat)
#what's the steps for s?
s.max = 1*p
#s.max = 1.5*p
s.min = 0.01
s.step = (s.max - s.min)/49
s.grid = seq(from=s.min, to=s.max, by=s.step)
ns = length(s.grid)
#record all 5-fold cv error for s on the s.grid
cv_errors_all_s = vector("numeric", length = ns)
cv_errors_mat = matrix(NA, nrow = ns, ncol = nfolds)
ck_mat = matrix(NA, nrow = ns, ncol = (p))
#y is already centered
df = cbind(y, fmat)
for(is in 1:ns){
#print (is)
si = s.grid[is]
test_errors_s = vector("numeric", nfolds) #initiate an error vector
for(i in 1:nfolds){
test_set_id = folds[[i]] #take out the ith fold as the test set ID
test_set = df[test_set_id, ,drop = FALSE] #specify the test set
train_set = df[-test_set_id, ,drop = FALSE] #specify the training set
xs_train = train_set[, -1, drop = FALSE] #put y at the 1st column
xs_test = test_set[, -1, drop = FALSE] #put y at the 1st column
y_train = train_set[, 1, drop = FALSE]
y_test = test_set[, 1, drop = FALSE]
#apply the trained model to the test set and get yhat_test
#trained model: for si, get a set of ck's, fk's are already estimated in step 1
qv = crossprod(xs_train)
cv = crossprod(xs_train, y_train)
#new: each c value >= 0 and <=1
bvec = c(rep(0, p), rep(-1, p), -si)
ansi = solve.QP(qv, cv, t(amat), bvec)
ck_vec = ansi$solution
ck_mat[is, ] = ck_vec #no intercept
yhat_test = xs_test %*% ck_vec
test_errors_s[i] = mean((y_test - yhat_test)^2) #compute the test error for the ith iteration
}
cv_errors_all_s[is] = mean(test_errors_s) # sum(k-fold cv)/k
cv_errors_mat[is, ] = test_errors_s
}
#print (cv_errors_all_s)
#print (which.min(cv_errors_all_s))
#if there is a tie, then choose the smallest s
best.s_ps = which(cv_errors_all_s == min(cv_errors_all_s))[1]
best.s = s.grid[best.s_ps]
cvec_kfold = ck_mat[best.s_ps, ,drop = TRUE]
#cvec_kfold = (cvec_kfold / sdxs) |> round(5)
#final model: use all data points
qv = crossprod(fmat)
cv = crossprod(fmat, y)
#best.s is not cv.glmnet$lambda.min
#new: each c value >= 0 and <=1
bvec = c(rep(0, p), rep(-1, p), -best.s)
ans.qp = solve.QP(qv, cv, t(amat), bvec)
cvec = ans.qp$solution
cvec = cvec |> round(7)
#back to the original scale
#cvec = (cvec / sdxs) |> round(5)
#for(ic in seq(NCOL(ck_mat))) {ck_mat[,ic] = ck_mat[,ic]/sdxs[ic]}
rslt = list(ck_mat = ck_mat, cvec = cvec, best.s_ps = best.s_ps, best.s = best.s)
#return (rslt)
}
#y is not centered.
if(type == "glmnet") {
cv_lasso = glmnet::cv.glmnet(fmat, y, alpha=1, intercept=TRUE, lower.limits=0,
standardize=FALSE, nfolds=nfolds, family=family, type.measure = "deviance")
#cv_lasso = glmnet::cv.glmnet(fmat, y, alpha=1, intercept=TRUE,
# standardize=FALSE, nfolds=nfolds, family=family, type.measure = "deviance")
cvec = coef(cv_lasso, s = "lambda.min")[-1]
rslt = list(cvec = cvec)#, nng.mod = nng.mod, cv_lasso = cv_lasso)
#return (rslt)
}
}
return (rslt)
}
##############################
#tranform shapes back to char
##############################
ShapeToChar_shapeselect = function(shp, tag = "x") {
#if (max(shp) > 16 | min(shp) < 0) {
# stop ('No such a shape! A shape value must be between 0 and 16.')
#}
if (tag == "x") {
if (shp == 0) {
#shp = 18
shp = 19 #this should include z too...need to add tag later
}
switch(shp,
cs1 = {ch = 'incr'},
#need to change later:
#cs1 = {ch = 'in the model'},
cs2 = {ch = 'decreasing'},
cs3 = {ch = 'convex'},
cs4 = {ch = 'concave'},
cs5 = {ch = 'increasing convex'},
cs6 = {ch = 'decreasing convex'},
cs7 = {ch = 'increasing concave'},
cs8 = {ch = 'decreasing concave'},
cs9 = {ch = 'smooth increasing'},
cs10 = {ch = 'smooth deceasing'},
cs11 = {ch = 'smooth convex'},
cs12 = {ch = 'smooth concave'},
cs13 = {ch = 'smooth increasing convex'},
cs14 = {ch = 'smooth increasing concave'},
cs15 = {ch = 'smooth decreasing convex'},
cs16 = {ch = 'smooth decreasing concave'},
#cs17 = {ch = 's'},
cs17 = {ch = 'flat'},
cs18 = {ch = 'linear'},
cs19 = {ch = 'out of the model'}
#{print ('No such a shape')}
)
} else if (tag == "z") {
if (shp == 0) {
shp = 2
}
switch(shp,
cs1 = {ch = 'in the model'},
cs2 = {ch = 'out of the model'}
)
} else if (tag == "tree") {
if (shp == 0) {
shp = 3
}
switch(shp,
cs1 = {ch = 'tree'},
cs2 = {ch = 'unordered'},
cs3 = {ch = 'out of the model'}
)
}
return (ch)
}
#--------------------------------------------------------------------
#print.classo: show the final choices
#--------------------------------------------------------------------
print.shapeselect = function(x,...)
{
print(x$family)
cat("Call:\n")
#print(object$mod_final)
if(inherits(x, "shapeselectAll")){
cat("Top models:\n")
print(head(x$mod_tried))
} else {
print(x$mod_tried)
}
cat("\n")
invisible(x)
}
#--------------------------------
#all subsets
#--------------------------------
shapeselectAll <- function(xmat, y, family = gaussian, nsim = 100,
eps = 1e-8, space = 'E', infocrit = c("CBIC", "CAIC"), pnt = TRUE, parallel = TRUE,
message = FALSE, a = .5, add_edfu = 0, decrease_edfu = 0, types = NULL, center = FALSE)
{
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
#linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
#muhat.fun <- cicfamily$muhat.fun
#ysim.fun <- cicfamily$ysim.fun
#deriv.fun <- cicfamily$deriv.fun
#dev.fun <- cicfamily$dev.fun
infocrit <- match.arg(infocrit)
if (is.character(family))
family <- get(family, mode = "function", envir = parent.frame())
if (is.function(family))
family <- family()
if (is.null(family$family))
stop("'family' not recognized!")
fml <- family$family
if(fml == "gaussian"){
wt.iter <- FALSE
} else {wt.iter <- TRUE}
xnms <- colnames(xmat)
if(is.null(xnms)){xnms <- paste('x', 1:NCOL(xmat), sep='')}
#new: use types to check if x is for spline or ordinal basis or just nominal
p <- NCOL(xmat)
for(ix in 1:p){
attr(xmat[, ix], "type") <- types[ix]
}
n <- length(y)
#--------------------------------
#step 0: project y onto 1
#--------------------------------
one <- rep(1, n)
etahat <- muhat <- mean(y)
rss <- sum((y - etahat)^2)
#scale numeric x first:
p <- NCOL(xmat)
sdx <- mx <- rep(0, p)
for(ix in 1:p){
if(is.numeric(xmat[,ix]) & types[ix] %in% c("num", "numeric")){
#xmat[,ix] <- scale(xmat[,ix], center=F)
#sdx[ix] <- attr(xmat[,ix], "scaled:scale")
xmat[,ix] <- scale(xmat[,ix], center=center)
#sdx <- c(sdx, attr(xmat[,ix], "scaled:scale"))
sdx[ix] <- attr(xmat[,ix], "scaled:scale")
if(center){
mx[ix] <- attr(xmat[, ix], 'scaled:center')
}
}
if(inherits(xmat[,ix], "factor") & types[ix] %in% "tree"){
xmat[,ix] <- as.numeric(as.character(xmat[,ix]))
}
attr(xmat[, ix], "type") <- types[ix]
}
#need a larger value to start cic/bic?
if(infocrit %in% c("CAIC", "CIC")) {
if(!wt.iter){
cic <- log(rss) + 2*1/n
} else {
llh <- llh.fun(y=y, muhat=rep(1/2, n), etahat=rep(0, n), n=n, fml=fml)
#test!
#llh <- 2 * (llh * n / 2) #2 times neg-loglike
cic <- llh + 2*1/n
}
}
if(infocrit == "CBIC") {
if(!wt.iter){
#cic <- n*log(rss) + log(n) * 1
cic <- log(rss) + a * log(n) * 1 / n^(6/7)
#cic <- log(rss) + log(n) * 1 / n^(6/7)
} else {
#m0 <- glm(y ~ 1, family=family)
#cic <- n*log(m0$deviance) + log(n) * 1
#test
#cic <- cic / n
llh <- llh.fun(y=y, muhat=rep(1/2, n), etahat=rep(0, n), n=n, fml=fml)
#cic <- llh + log(n) * 1 / n
#cic <- llh + a * log(n) * 1 / n^(6/7)
cic <- llh + a * log(n) * (1+add_edfu) / n^(6/7)
#print (cic)
}
}
#----------------------------------------------------------------------
#step 1: list all shape combs
#----------------------------------------------------------------------
totx <- NCOL(xmat)
tot_combs <- 1
shps_lst <- vector("list", length = totx)
for(ic in 1:totx){
if(inherits(xmat[[ic]], 'factor')){
acc <- 2 #include flat
#shps_lst[[ic]] <- c(0, 1)
#new
#if(attr(xmat[[ic]], 'type') %in% c('nom', 'nominal')){
shps_lst[[ic]] <- c(0, 1)
#}
} else {
#acc <- 9 #include flat: 17
if(attr(xmat[[ic]], 'type') %in% c('ord', 'ordinal')){
shps_lst[[ic]] <- c(1:8, 0)
}
if(attr(xmat[[ic]], 'type') %in% c('tree')){
shps_lst[[ic]] <- c(1:2, 0)
}
if(attr(xmat[[ic]], 'type') %in% c('num', 'numeric')){
shps_lst[[ic]] <- 9:17
}
acc <- length(shps_lst[[ic]])
#new: if we use cbic and wt.iter then we only have 4 options
if(wt.iter & infocrit == "CBIC"){
if(attr(xmat[[ic]], 'type') %in% c('ord', 'ordinal')){
shps_lst[[ic]] <- c(1:4, 0)
}
if(attr(xmat[[ic]], 'type') %in% c('num', 'numeric')){
shps_lst[[ic]] <- c(9:12, 17)
}
if(attr(xmat[[ic]], 'type') %in% c('tree')){
shps_lst[[ic]] <- c(1:2, 0)
}
acc <- length(shps_lst[[ic]])
}
}
tot_combs <- tot_combs * acc
}
msg <- paste("We will have a total of", tot_combs, "models.", "Continue?")
if(message){
if(tot_combs > 160){
#library(svDialogs)
res <- dlgMessage(msg, "yesno")$res
if (res == "yes") {
cat("Computing CBIC for each model....", "\n")
} else {
stop("You can call shapeselect instead.", "\n")
}
}
}
mod_tried <- matrix(0, nrow = (tot_combs+1), ncol = (totx + 5))
colnames(mod_tried) <- c(xnms, infocrit, "edf0", "-2/n*loglike", "penalty", "penalty for spline")
#the 1st row is for etahat = ybar
mod_tried[1, totx+1] <- cic
mod_tried[1, totx+2] <- 1
#----------------------------------------------------------------------
#step 2: make delta, compute edf and cic
#----------------------------------------------------------------------
shp_combs <- expand.grid(shps_lst)
colnames(shp_combs) <- xnms
zmat_all <- varlst_all <- varlst_z_all <- delta_all <- coefs_lst <- vector("list", length = tot_combs)
for(icomb in 1:tot_combs){
#print (icomb)
shp_comb_i <- as.numeric(shp_combs[icomb, ])
if(all(shp_comb_i %in% c(0, 17))){
mod_tried[icomb+1, ] <- mod_tried[1, ,drop = TRUE]
} else {
non_flat_ps <- which(!shp_comb_i %in% c(0, 17))
#new: for ordinal
#if(all(shp_comb_i[non_flat_ps] < 9)){
# pnt = FALSE
#}
gmat <- NULL
np_ix <- 1
zmat <- NULL
delta <- NULL
dmat <- NULL
varlst <- NULL
varlst_z <- NULL
dm_lst <- vector("list", length = length(non_flat_ps))
for(ix in non_flat_ps){
x <- xmat[, ix]
if(inherits(x, "factor") & attr(x, 'type') != 'tree'){
#new:
#if(attr(x, 'type') %in% c('nom', 'nominal')){
dd <- model.matrix(~ x)[, -1, drop=FALSE]
ms <- apply(dd, 2, mean)
dd <- apply(dd, 2, function(e) e - mean(e))
zmat <- cbind(zmat, dd) #no intercept
varlst_z <- c(varlst_z, rep(ix, ncol(dd)))
#}
} else {
shpi <- shp_comb_i[ix]
#new: include ordinal
if(shpi <= 8){
if(attr(x, 'type') %in% c('ord', 'ordinal')){
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = 0, space = space)
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
ps <- 0
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
}
#new
if(attr(x, 'type') %in% c('tree')){
#ms <- NULL
dd <- t(tree.fun(x))
ms <- apply(dd, 2, mean)
dd <- apply(dd, 2, function(e) e - mean(e))
#gmat <- cbind(one, dd)
if(shpi == 1){
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
ps <- 0
}
if(shpi == 2){
dm <- NULL
ps <- 0
}
}
if(shpi != 2 | attr(x, 'type') != 'tree'){
delta <- cbind(delta, dd)
}
if(shpi != 2 | attr(x, 'type') != 'tree'){
varlst <- c(varlst, rep(ix, ncol(dd)))
}
if(shpi > 2 & shpi < 5 | shpi == 2 & attr(x, 'type') == 'tree'){
if(shpi > 2 & shpi < 5){
zmat <- cbind(zmat, x - mean(x)) #remove mean(x)?
varlst_z <- c(varlst_z, ix)
}
if(shpi == 2 & attr(x, 'type') == 'tree'){
zmat <- cbind(zmat, dd)
varlst_z <- c(varlst_z, rep(ix, ncol(dd)))
}
}
ps <- 0
if(shpi != 2 | attr(x, 'type') != 'tree'){
dm <- matrix(0, nrow=ncol(dd), ncol=ncol(dd))
dm_lst[[ix]] <- dm #no penalty
}
}
if(shpi > 8){
mult <- ifelse(pnt, 2, 1)
spl_degree_i <- ifelse(shpi %in% c(9, 10), 2, 3) #use 2 for i-spl, 3 for c-spl
n1 <- length(unique(x))
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5, trunc(n1^(1/7)) + 5)
nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/7)) + 6, trunc(n1^(1/9)) + 6)
#new
#nkts <- mult * ifelse(spl_degree_i == 2, trunc(n1^(1/5)) + 5 - decrease_edfu,
# trunc(n1^(1/7)) + 5 + add_edfu)
edfu <- ifelse(spl_degree_i == 2, nkts/mult, nkts/mult + 1)
#print (c(nkts, edfu))
kts <- seq.int(min(x), max(x), length.out = nkts)
dd_ans <- makedelta(x, sh = shpi, numknots = 0, knots = kts, space = space)
#dd is orth to 1
ms <- dd_ans$ms
dd <- t(dd_ans$amat)
#kts <- dd_ans$knots
delta <- cbind(delta, dd) #splines excluding conv/conc identity term
varlst <- c(varlst, rep(ix, ncol(dd)))
if(shpi > 10 & shpi < 13){
zmat <- cbind(zmat, x - mean(x))
varlst_z <- c(varlst_z, ix)
}
if(pnt){
nc <- NCOL(dd)
dm <- makedmat_1D_ic(m = nc, knots = kts, shpi = shpi)
qvi <- crossprod(dd)
dvi <- crossprod(dm)
ps <- uniroot(f = .search_ps_shapeselect, qv0 = qvi, dv = dvi, dd = dd,
edfu = edfu, interval = c(1e-10, 1e+10), tol = .Machine$double.eps^0.25)$root
#print (ps)
#ps <- .5
#test!
#if(shpi > 10) {ps = 0}
#if(shpi < 11) {ps = 1}
dm <- sqrt(ps) * dm
dm_lst[[ix]] <- dm
}
}
}
}
#---------------
#compute etahat
#---------------
gmat <- cbind(one, zmat, delta)
if(!is.null(zmat)){
zmat_all[[icomb]] <- zmat
varlst_z_all[[icomb]] <- varlst_z
}
if(!is.null(delta)){
delta_all[[icomb]] <- delta
varlst_all[[icomb]] <- varlst
}
gmat0 <- gmat
np_ix <- 1 + NCOL(zmat)
dmat <- as.matrix(bdiag(Filter(Negate(is.null), dm_lst)))
#new: change dmat to be 0 if it is 1 by 1 0
# if(dim(dmat)[1] == 1 & dim(dmat)[2] == 1){
# dmat <- dmat[1, 1]
# }
qv0 <- crossprod(gmat)
xtx <- qv0
nc <- ncol(gmat)
amat <- diag(nc)
ps_zero <- 1:np_ix
amat[ps_zero, ps_zero] = 0
bvec <- cbind(rep(0, nc))
if(pnt & !is.null(dmat)){
dmat_zero = matrix(0, nrow = nrow(dmat), ncol = np_ix)
dmat2 = cbind(dmat_zero, dmat)
dv = crossprod(dmat2)
print (dim(qv0))
print (dim(dv))
xtx = qv0 + dv
}
if(!wt.iter){
xty = crossprod(gmat, y)
# if(!pnt){
# dsend <- gmat[, -c(1:np_ix), drop = FALSE]
# zsend <- gmat[, 1:np_ix, drop = FALSE]
# ans <- suppressMessages(try(coneB(y, dsend, zsend), silent = TRUE))
# }
#if(pnt & !is.null(dmat)){
ans = suppressMessages(try(solve.QP(xtx, xty, t(amat), bvec), silent = TRUE))
#}
if (!inherits(ans, "try-error")) {
#if(pnt){
bhat = ans$solution
face = ans$iact
#}
# if(!pnt){
# face = ans$face
# bhat = ans$coefs
# }
etahat = gmat %*% bhat
muhat = etahat
}
} else {
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
muhat <- family$linkinv(etahat)
#sm <- 1e-8
sm <- 1e-5
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 10)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(y, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(y, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
ans <- suppressMessages(try(coneB(zvec, dsend, zsend), silent = TRUE))
if (inherits(ans, "try-error")) {
break
} else {
face <- ans$face
bhat <- ans$coefs
}
}
if(pnt & !is.null(dmat)){
xtx = crossprod(gmat)
xtx = xtx + dv
xty = crossprod(gmat, zvec)
ans = suppressMessages(try(solve.QP(xtx, xty, t(amat), bvec), silent = TRUE))
if (inherits(ans, "try-error")) {
break
} else {
bhat = ans$solution
face = ans$iact
}
}
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
coefs_lst[[icomb]] <- bhat
#---------------
#compute cic
#---------------
cores = min(8, parallel::detectCores(logical = FALSE) - 1)
if(parallel){
edfs = parallel::mclapply(1:nsim,
.compute_edf0_fix_shps, gmat=gmat0, np_ix, dmat,
pnt, qv0, amat, bvec, dv,
family, ps, mc.cores = (cores))
}else{
edfs = lapply(1:nsim,
.compute_edf0_fix_shps, gmat=gmat0, np_ix, dmat,
pnt, qv0, amat, bvec, dv, family, ps)
}
edfs = Filter(Negate(is.null), edfs)
edfs = simplify2array(edfs)
#cat('shpi comb', icomb, '\n')
#print (edfs)
if(length(edfs) > 0){
dfmean <- mean(edfs, na.rm = TRUE)
}
#temp
if(length(edfs) == 0) {
dfmean <- np_ix
}
if(wt.iter){
llh <- llh.fun(y, muhat, etahat, phihat=NULL, n, weights=NULL, fml=family$family)
#test!
#llh <- llh * n / 2
#llh <- 2 * (llh * n / 2) #2 times neg-loglike
} else {
rss <- sum((y - etahat)^2)
llh <- log(rss)
}
if (infocrit == "CBIC") {
#cic_val <- n * llh + log(n) * (dfmean)
#test
#cic_val <- llh + log(n) * (dfmean) / n
#new
#cic_val <- llh + a * log(n) * (dfmean) / n^(6/7)
if(any(non_flat_ps %in% c(9, 10))){ #smooth incr/decr
degree <- 6/7
} else if (any(non_flat_ps > 10 & non_flat_ps < 17)){ #smooth conv/conc
degree <- 6/7#8/9
} else if (all(non_flat_ps > 0 & non_flat_ps < 9)) { #ord or parametric
degree <- 6/7#1
} else {
stop('check!')
}
cic_val <- llh + a * log(n) * (dfmean) / n^degree
#new
# if(shpi %in% c(9, 10)){
# bic_p <- 3
# }
# if(shpi %in% c(11, 12, 13, 14, 15, 16)){
# bic_p <- 3
# }
# if(shpi %in% 18){
# bic_p <- 2
# }
# cic_val <- llh + a * log(n) * (dfmean) * n^(-1) #* n^(-2bic_p/(2*bic_p+1))
}
if (infocrit == "CAIC") {
cic_val <- llh + 2 * (dfmean) / n
}
mod_tried[icomb + 1, 1:totx] <- shp_comb_i
mod_tried[icomb + 1, totx + 1] <- cic_val
mod_tried[icomb + 1, totx + 2] <- dfmean
mod_tried[icomb + 1, totx + 3] <- llh
mod_tried[icomb + 1, totx + 4] <- cic_val - llh
mod_tried[icomb + 1, totx + 5] <- ps
}
}
#----------------------------------------------------------------------
#step 3: choose the best one
#----------------------------------------------------------------------
#mod_tried <- unique(mod_tried)
cics <- as.vector(mod_tried[, totx+1])
optim_id <- which.min(cics)[1]
optim_shp_comb <- mod_tried[optim_id, 1:totx]
optim_shp_comb <- as.vector(optim_shp_comb)
ordcic <- order(cics)
mod_tried_unordered <- mod_tried
mod_tried <- mod_tried[ordcic, ]
#----------------------------------------------------------------------
#step 4: compute etacomps for the best
#----------------------------------------------------------------------
mod_best <- mod_tried[1, ,drop = FALSE]
n <- length(y)
etacomps <- matrix(0, nrow = totx, ncol = n)
etahatkeep <- ifelse(family$family == 'gaussian', rep(mean(y), n), rep(0, n))
muhatkeep <- family$linkinv(etahatkeep)
if(optim_id > 1) {
optim_id_noflat <- optim_id-1
coefskeep <- coefs_lst[[optim_id_noflat]]
ddkeep_z_all <- zmat_all[[optim_id_noflat]]
ddkeep_x_all <- delta_all[[optim_id_noflat]]
varlst <- varlst_all[[optim_id_noflat]]
varlst_z <- varlst_z_all[[optim_id_noflat]]
coefskeep0 <- coefskeep
coefskeep <- coefskeep[-1] #remove intercept
capk <- length(varlst_z)
if(capk > 0){
zcoefs <- coefskeep[1:capk]
xcoefs <- coefskeep[-c(1:capk)]
} else {
zcoefs <- NULL
xcoefs <- coefskeep
}
#new
etahatkeep <- cbind(one, ddkeep_z_all, ddkeep_x_all) %*% coefskeep0
muhatkeep <- family$linkinv(etahatkeep)
for(i in 1:totx){
ps <- i
xi <- xmat[, i]
shpi <- optim_shp_comb[i]
#if(shpi == 17 | ps %in% x_rm){#flat
if(shpi == 17 | shpi == 0){#flat
next
} else {
#one <- cbind(rep(1, n))
#xi = as.numeric(xi)
if(shpi == 18) {
etacomps[ps, ] <- ddkeep_z_all[, varlst_z %in% ps, drop = FALSE] %*% zcoefs[which(varlst_z %in% ps)]
} else if (shpi == 1 | shpi == 2) {#need to use 1 for incr later
if(attr(xi, 'type') %in% c('nom', 'nominal') | attr(xi, 'type') %in% 'tree' & shpi == 2){
etacomps[ps, ] <- (ddkeep_z_all[, which(varlst_z %in% ps), drop = FALSE]) %*% zcoefs[which(varlst_z %in% ps)]
}
if(attr(xi, 'type') %in% c('ord', 'ordinal') | attr(xi, 'type') %in% 'tree' & shpi == 1){
etacomps[ps, ] <- (ddkeep_x_all[, which(varlst %in% ps), drop = FALSE]) %*% xcoefs[which(varlst %in% ps)]
}
} else if (shpi > 10 & shpi < 13 | shpi > 2 & shpi < 5) {
etacomps[ps, ] <- (ddkeep_x_all[, which(varlst %in% ps), drop = FALSE]) %*% xcoefs[which(varlst %in% ps)]
etacomps[ps, ] <- etacomps[ps, ] + (xi-mean(xi)) * zcoefs[which(varlst_z %in% ps)]
} else {
etacomps[ps, ] <- (ddkeep_x_all[, which(varlst %in% ps), drop = FALSE]) %*% xcoefs[which(varlst %in% ps)]
}
}
}
}
#----------------------------------------------------------------------
#step 5: Lasso shrink? no need
#----------------------------------------------------------------------
#----------------------------------------------------------------------
#return result
#----------------------------------------------------------------------
if(any(optim_shp_comb == 17)) {
ps_17 = which(optim_shp_comb == 17)
optim_shp_comb[ps_17] = 0
}
tags <- rep('x', p)
tags[which(types %in% c('nom', 'nominal'))] <- 'z'
tags[which(types %in% 'tree')] <- 'tree'
#shapes in character
shapes <- mapply(ShapeToChar_Classo, optim_shp_comb, tags)
rslt <- list(mod_tried = mod_tried, etacomps = etacomps, mod_best = mod_best, shps = optim_shp_comb, shapes=shapes,
family = family, etahat = etahatkeep, muhat = muhatkeep)
class(rslt) <- c("shapeselectAll", "shapeselect")
return (rslt)
}
#####################################
#simulate edf for a fixed shape comb#
#####################################
.compute_edf0_fix_shps <- function(isim, gmat, np_ix, dmat = NULL, pnt = FALSE,
qv0 = NULL, amat = NULL, bvec = NULL, dv = NULL,
family = gaussian(link = 'identity'), ps=0)
{
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
muhat.fun <- cicfamily$muhat.fun
ysim.fun <- cicfamily$ysim.fun
deriv.fun <- cicfamily$deriv.fun
dev.fun <- cicfamily$dev.fun
eps = 1e-8
wt.iter = ifelse(family$family == 'gaussian', FALSE, TRUE)
if(NCOL(gmat) == np_ix) {
edfi <- np_ix
}
if (NCOL(gmat) > np_ix) {
n <- NROW(gmat)
#ysim <- ysim.fun(n, mu0 = NULL, fml = family$family, shp0 = NULL)
if(family$family == 'binomial'){
ysim <- rbinom(n, size = 1, prob = .5)
}else{
ysim <- ysim.fun(n, mu0 = NULL, fml = family$family, shp0 = NULL)
}
#-------------------------------------------
#use qprog
#nc = ncol(gmat)
#amat = diag(nc)
#ps_zero = 1:np_ix
#amat[ps_zero, ps_zero] = 0
#bvec = cbind(rep(0, nc))
gmat0 = gmat
#qv0 = crossprod(gmat0)
xtx = qv0
#use this to test edfi
#if(round(ps, 8) == 0){
# pnt = FALSE
#}
if(pnt & !is.null(dmat)){
#dmat_zero = matrix(0, nrow = nrow(dmat), ncol = np_ix)
#dmat2 = cbind(dmat_zero, dmat)
#dv = crossprod(dmat2)
xtx = qv0 + dv
}
#-------------------------------------------
if(!wt.iter){
xty = crossprod(gmat, ysim)
#if(pnt){
ans = suppressMessages(try(solve.QP(xtx, xty, t(amat), bvec), silent = TRUE))
#}
# if(!pnt){
# dsend <- gmat[, -c(1:np_ix), drop = FALSE]
# zsend <- gmat[, 1:np_ix, drop = FALSE]
# ans <- suppressMessages(try(coneB(ysim, dsend, zsend), silent = TRUE))
# }
if (!inherits(ans, "try-error")) {
#if(pnt){
bhat = ans$solution
face = ans$iact
#}
# if(!pnt){
# face <- ans$face
# bhat <- ans$coefs
# }
}
} else {
etahat <- etahat.fun(n, ysim, fml = family$family)
gr <- gr.fun(ysim, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(ysim, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, ysim, fml = family$family)
muhat <- family$linkinv(etahat)
#sm <- 1e-8
sm <- 1e-5
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- 1
niter <- 0
while (diff > sm & mdiff & niter < 10)
{
oldmu <- muhat
niter <- niter + 1
gr <- gr.fun(ysim, etahat, weights=NULL, fml = family$family)
wt <- wt.fun(ysim, etahat, n, weights=NULL, fml = family$family)
cvec <- wt * etahat - gr
zvec <- zvec.fun(cvec, wt, ysim, fml = family$family)
gmat <- sweep(gmat0, MARGIN = 1, sqrt(wt), FUN = '*')
if(!pnt | pnt & all(dmat == 0)){
dsend <- gmat[, -c(1:np_ix), drop = FALSE]
zsend <- gmat[, 1:np_ix, drop = FALSE]
#test!
ans <- suppressMessages(try(coneB(zvec, dsend, zsend), silent = TRUE))
if (inherits(ans, "try-error")) {
break
} else {
face <- ans$face
bhat <- ans$coefs
}
}
if(pnt & !is.null(dmat)){
xtx = crossprod(gmat)
xtx = xtx + dv
#test!
xty = crossprod(gmat, zvec)
ans = suppressMessages(try(solve.QP(xtx, xty, amat, bvec), silent = TRUE))
#ans = suppressMessages(try(solve.QP(xtx, xty, t(amat), bvec), silent = TRUE))
if (inherits(ans, "try-error")) {
break
} else {
bhat = ans$solution
face = ans$iact
}
}
etahat <- gmat0 %*% bhat
muhat <- family$linkinv(etahat)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
diff <- mean((muhat - oldmu)^2)
}
}
if(inherits(ans, "try-error")) {
edfi <- NULL
} else {
if(!pnt | pnt & all(dmat == 0)){
#test
bhat <- round(bhat, 10)
edfi <- sum(abs(na.omit(bhat)) > eps)
} else {
#edfi <- sum(abs(na.omit(bhat)) > eps)
umat = chol(xtx)
uinv = backsolve(umat, diag(ncol(xtx)))
atil = amat %*% uinv
xtxinv = tcrossprod(uinv)
if(!wt.iter){
bigp_base = xtxinv %*% qv0
}else{
bigp_base = xtxinv %*% crossprod(gmat)
}
edfi_base = sum(diag(bigp_base))
if (all(face == 0) | length(face) == 0) {
edfi = edfi_base
} else {
smat = -t(atil[face, ,drop = FALSE])
#new: avoid solve error
qr_smat = qr(smat)
smat = qr.Q(qr_smat)
#smat = qr.Q(qr_smat)[, 1:(qr_smat$rank), drop = FALSE]
pmat_polar = tcrossprod(smat)
if(!wt.iter){
bigp_polar = uinv %*% tcrossprod(pmat_polar, uinv) %*% qv0
}else{
bigp_polar = uinv %*% tcrossprod(pmat_polar, uinv) %*% crossprod(gmat)
}
edfi = edfi_base - sum(diag(bigp_polar))
}
}
}
}
return (edfi)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.