Genbasis <- function(sseq, df, degree, type = c("bs", "fourier", "OBasis")) {
type <- match.arg(type)
interval <- range(sseq)
nknots <- df - (degree + 1)
if(nknots > 0) {
knots <- ((1:nknots) / (nknots + 1)) * diff(interval) + interval[1]
} else {
knots <- NULL
}
basis <- switch(type,
"bs" = bs(x = sseq, df = df, degree = degree,
Boundary.knots = interval, intercept = TRUE),
"fourier" = eval.basis(sseq,
basisobj = create.fourier.basis(rangeval = interval, nbasis = df)
),
"OBasis" = evaluate(OBasis(expand.knots(c(interval[1], knots, interval[2])),
order = degree + 1),
sseq)
)
return(basis)
}
cv.test <- function(outlist, y, X, foldid, lam, trim = 0, keep = FALSE) {
nlam <- length(lam)
n <- length(y)
predmat <- matrix(NA, n, nlam)
nfolds <- max(foldid)
for (i in seq(nfolds)) {
which <- foldid == i
#pred <- X_test %*% outlist[[i]]$path
#pred <- apply(outlist[[i]], 1 ,function(beta,X_test) X_test %*% beta, X_test=X[which, ])
pred <- predict(outlist[[i]], X[which, , drop = FALSE])
nlami <- length(outlist[[i]]$lam)
predmat[which, seq(nlami)] <- pred }
cvraw <- (y - predmat)^2
N <- length(y) - apply(is.na(predmat), 2, sum)
cvm <- apply(cvraw, 2, mean, na.rm = TRUE)
cvsd <- sqrt(apply(scale(cvraw, cvm, FALSE)^2, 2, mean, na.rm = TRUE)/(N - 1))
if(trim > 0) {
cv.trim <- apply(cvraw, 2, function(x) {
x <- x[!is.na(x)]
x.boundary <- quantile(x, probs = c(trim / 2, 1 - trim / 2))
x <- x[x < x.boundary[2]]
x <- x[x >= x.boundary[1]]
x.mean <- mean(x)
x.sd <- sd(x)
return(c(MEAN = x.mean, SD = x.sd))
})
} else {cv.trim <- NULL}
output <- list()
output$cvm <- cvm
output$cvsd <- cvsd
output$cvmtrim <- cv.trim[1, ]
output$cvsdtrim <- cv.trim[2, ]
if(keep) output$fit.preval <- predmat
return(output)
}
getmin <- function(lam, cvm, cvsd, digits = 5) {
idx <- !is.na(cvm)
cvm1 <- cvm[idx]
cvsd1 <- cvsd[idx]
cvm1 <- round(cvm1, digits = digits)
cvmin <- min(cvm1, na.rm = TRUE)
idmin <- cvm1 <= cvmin
lam.min <- max(lam[idmin])
idmin <- match(lam.min, lam)
semin <- round((cvm1 + cvsd1), digits = digits)[idmin]
idmin <- cvm1 <= semin
lam.1se <- max(lam[idmin])
output <- list(lam.min = lam.min, lam.1se = lam.1se)
return(output)
}
ggetmin <- function(lam, cvm, cvsd, digits = 5, k_list) {
#cvmin <- apply(cvm, 1, min, na.rm = TRUE)
cvm1 <- round(cvm, digits = digits)
cvmin <- min(cvm1, na.rm = TRUE)
lam.min <- apply(cvm1 <= cvmin, 1, function(x, lam)
ifelse(length(lam[x]) == 0, 0, max(lam[x], na.rm = TRUE)), lam=lam)
# sapply(apply(cvm1 <= cvmin, 1, function(x, lam) lam[x], lam=lam),
# function(x) ifelse(length(x) > 0, max(x, na.rm = TRUE), 0))
lam.min_k <- which.max(lam.min)
lam.min <- max(lam.min)
idmin <- match(lam.min, lam)
semin <- round((cvm + cvsd), digits = digits)[lam.min_k, idmin]
lam.1se <- apply(cvm1 <= semin, 1, function(x, lam)
ifelse(length(lam[x]) == 0, 0, max(lam[x], na.rm = TRUE)), lam=lam)
# sapply(apply(cvm1 <= semin, 1, function(x, lam) lam[x], lam=lam),
# function(x) ifelse(length(x) > 0, max(x, na.rm = TRUE), 0))
lam.1se_k <- which.max(lam.1se)
lam.1se <- max(lam.1se)
output <- list(lam.min = c("lam" = lam.min, "df" = k_list[lam.min_k]),
lam.1se = c("lam" = lam.1se, "df" = k_list[lam.1se_k]))
return(output)
}
vet <- function(beta, p, k) {
p1 <- p * k
coef <- matrix(beta[1:p1], byrow = TRUE, nrow = p)
result<- list(C = coef)
coef <- beta[(p1+1):length(beta)]
result$b <- coef
return(result)
}
Nzero <- function(beta, p, k, tol) {
coef <- vet(beta, p = p, k = k)$C
group <- apply(coef, 1, function(x) ifelse(max(abs(x)) > tol , TRUE, FALSE))
group <- (1:p)[group]
return(group)
}
ERROR_fun <- function(beta_fit, beta_true,
basis_fit, basis_true, sseq,
m, p, Nzero_group, tol = 0) {
error.list <- list()
pos_group <- 1:Nzero_group
neg_group <- (1:p)[-pos_group]
df_fit <- (length(beta_fit) - 1 - m) / p
df_true <- (length(beta_true) -1 -m) / p
Non.zero_select <- Nzero(beta = beta_fit, p = p, k = df_fit, tol = tol)
error.list$class_error <- Class_error(pos_select = Non.zero_select,
pos_group = pos_group, neg_group = neg_group)
error.list$coef_error <- Coef_error(coef_true = beta_true, coef_esti = beta_fit,
p = p, k_fit = df_fit, k_true = df_true)
error.list$curve_error <- Curve_error(coef_true = beta_true, coef_esti = beta_fit, p = p,
k_fit = df_fit, k_true = df_beta,
basis_true = basis_true, basis_fit = basis_fit, sseq = sseq)
error.list$Non.zero <- Non.zero_select
return(error.list)
}
##classification error
Class_error <- function(pos_select, pos_group, neg_group) {
pos <- length(pos_group)
neg <- length(neg_group)
FP_select <- pos_select[which( pos_select %in% neg_group )]
FN_select <- pos_group[which( !(pos_group %in% pos_select) )]
TP_select <- pos_select[which( pos_select %in% pos_group )]
TN_select <- neg_group[which( !(neg_group %in% pos_select) )]
FP <- length(FP_select)
FN <- length(FN_select)
TP <- length(TP_select)
TN <- length(TN_select)
FPR <- FP / neg # false positive rate
FNR <- FN / pos # false negative rate
Sensitivity <- TP / pos # true positive rate
Specificity <- TN / neg # true negative rate
PPV <- TP / (TP + FP) # precision/positive predictive value
NPV <- TN/ (TN + FN) # negative predictive value
FDR <- 1 - PPV # false discovery rate
ACC <- (TP + TN) / p
F1_score <- (2 * TP) / (2 * TP + FP + FN)
MCC <- (TP * TN - FP * FN) / sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN)) # Matthews correlation coefficient
Informedness <- Sensitivity + Specificity - 1 # Youden's J statistic
Markedness <- PPV + NPV - 1
error <- c(FP = FP, FN = FN, TP = TP, TN = TN,
FPR = FPR, FNR = FNR, Sensi = Sensitivity, Speci = Specificity,
PPV = PPV, NPV = NPV, FDR = FDR, ACC = ACC, F1_score = F1_score,
MCC = MCC, Informedness = Informedness, Markedness = Markedness)
return(error)
}
##Coefficients estimation Norm error
Coef_error <- function(coef_true, coef_esti, p, k_fit, k_true) {
coef_esti.comp <- vet(beta = coef_esti, p = p, k = k_fit)$C
coef_esti.cont <- vet(beta = coef_esti, p = p, k = k_fit)$b
coef_true.comp <- vet(beta = coef_true, p = p, k = k_true)$C
coef_true.cont <- vet(beta = coef_true, p = p, k = k_true)$b
#cat("1")
error <- vector()
error <- c(L2.cont = sqrt(sum( (coef_esti.cont - coef_true.cont)^2 )) )
if(k_fit == k_true) {
L2_diff.beta <- sqrt(sum( (coef_true - coef_esti)^2 ))
L2_diff.comp <- apply(coef_esti.comp - coef_true.comp, 1,
function(x) sqrt(sum(x^2)) )
L2_diff_inf.comp <- max(L2_diff.comp)
L2_diff_L1.comp <- sum(L2_diff.comp)
L2_diff.comp <- sqrt(sum( (as.vector(coef_esti.comp - coef_true.comp))^2 ))
error <- c(error, L2.beta = L2_diff.beta, L2.comp = L2_diff.comp,
L2_inf.comp = L2_diff_inf.comp, L2_L1.comp = L2_diff_L1.comp)
}
#cat('2')
return(error)
}
Curve_error <- function(coef_true, coef_esti, p,
k_fit, k_true, basis_true, basis_fit, sseq) {
coef_esti.comp <- vet(beta = coef_esti, p = p, k = k_fit)$C
coef_true.comp <- vet(beta = coef_true, p = p, k = k_true)$C
curve_true <- coef_true.comp %*% t(basis_true)
curve_esti <- coef_esti.comp %*% t(basis_fit)
curve_diff <- abs(curve_esti - curve_true) ## p by length(sseq)
ns <- length(sseq)
time_diff <- sseq[2] - sseq[1]
extra_sum <- rowSums(curve_diff[, c(1, ns)]^2) * time_diff / 2## crossprod(curve_diff[, c(1, ns)]) * time_diff/2
add_sum <- apply(curve_diff, 1, function(x) sum(x^2)) * time_diff
ITG <- add_sum - extra_sum
L2 <- sqrt(ITG)
L2_L1 <- sum(L2)
L2_inf <- max(L2)
extra_sum <- rowSums(curve_diff[, c(1, ns)]) * time_diff / 2## crossprod(curve_diff[, c(1, ns)]) * time_diff/2
add_sum <- rowSums(curve_diff) * time_diff
ITG <- add_sum - extra_sum
L1 <- ITG
L1_L1 <- sum(L1)
L1_inf <- max(L1)
L1_each.max <- max(curve_diff)
error <- c(L1_each.max = L1_each.max,
L1_inf = L1_inf, L1_L1 = L1_L1,
L2_inf = L2_inf, L2_L1 = L2_L1 )
return(error)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.