Nothing
#' @importFrom stats plogis
J.BipLog.BIN <- function(xt, par, k, lambda, ...) {
#---- Para iniciar.
x = as.matrix(xt)
n = nrow(x)
p = ncol(x)
aik = n * k
bjk = p * (k + 1)
dTheta = aik + bjk
A = matrix(par[1:aik], n, k)
B = matrix(par[(aik + 1):dTheta], p, k + 1)
lin = cbind(rep(1, n), A) %*% t(B)
hX = plogis(lin)
J = -sum(x * log(hX) + (1 - x) * log(1 - hX), na.rm = TRUE) + lambda *
sum(A^2, na.rm = TRUE)/2 + lambda * sum(B[, -1]^2, na.rm = TRUE)/2
return(J)
}
Grad.BipLog.BIN <- function(xt, par, k, lambda, ...) {
x = as.matrix(xt)
n = nrow(x)
p = ncol(x)
aik = n * k
bjk = p * (k + 1)
dTheta = aik + bjk
A = matrix(par[1:aik], n, k)
B = matrix(par[(aik + 1):dTheta], p, k + 1)
lin = cbind(rep(1, n), A) %*% t(B)
hX = plogis(lin)
dA = (hX - x) %*% B[, -1] + lambda * A
dB = t(hX - x) %*% cbind(rep(1, n), A) + lambda * cbind(rep(0, p), B[,
-1])
gradient = c(c(dA), c(dB))
return(gradient)
}
Indsup.BIN <- function(xs, B, par, k, lambda, ...) {
x = as.matrix(xs)
n = nrow(x)
p = ncol(x)
aik = n * k
dTheta = aik
A = matrix(par[1:aik], n, k)
Bs = as.matrix(B)
lin = cbind(rep(1, n), A) %*% t(Bs)
hX = plogis(lin)
L = -sum(x * log(hX) + (1 - x) * log(1 - hX), na.rm = TRUE) + lambda *
sum(A^2, na.rm = TRUE)/2 + lambda * sum(Bs[, -1]^2, na.rm = TRUE)/2
return(L)
}
Indsup.GradBIN <- function(xs, B, par, k, lambda, ...) {
x = as.matrix(xs)
Bs = as.matrix(B)
n = nrow(x)
p = ncol(x)
aik = n * k
dTheta = aik
A = matrix(par[1:aik], n, k)
lin = cbind(rep(1, n), A) %*% t(Bs)
hX = plogis(lin)
dA = (hX - x) %*% Bs[, -1] + lambda * A
gradient = c(dA)
return(gradient)
}
train_miss_pattern <- function(data, K = 7){
X <- as.matrix(data)
out <- sapply(1:K, function(x) seq(x, length(X), by = K), simplify = FALSE)
train <- list()
for(i in 1:K){
temp <- X
temp[out[[i]]] <- NA
train[[i]] <- temp
}
lout <- list(missWold = out, Xtrain = train)
return(lout)
}
thresholds <- function(x, P, ncuts = 100){
if (!requireNamespace("tidyr", quietly = TRUE)) {
stop("Package \"tidyr\" needed for this function to work. Please install it.",
call. = FALSE)
}
if (!requireNamespace("dplyr", quietly = TRUE)) {
stop("Package \"dplyr\" needed for this function to work. Please install it.",
call. = FALSE)
}
P <- as.matrix(P)
if(is.null(colnames(x))) colnames(x) <- paste0('V', 1:ncol(x))
N1 <- ceiling(sum(x, na.rm = T))
N0 <- length(as.matrix(x)) - N1 - sum(is.na(x))
TE <- lapply(seq(0, 1, length.out = ncuts), function(z){
Pr <- ifelse(P>=z, 1, 0)
c1 <- 1 - apply((Pr == 1) & (x == 1), 2, sum, na.rm=TRUE)/apply(x == 1, 2, sum, na.rm = T)
c2 <- 1 - apply((Pr == 0) & (x == 0), 2, sum, na.rm=TRUE)/apply(x == 0, 2, sum, na.rm = T)
TE <- 100/2 * (c1 + c2)
lista <- list(TE)
})
TEp <- data.frame(dplyr::bind_rows(TE), threshold = seq(0, 1, length.out = ncuts))
thresholds <- TEp |>
tidyr::pivot_longer(-threshold, names_to = "variable", values_to = "BACC") |>
dplyr::group_by(variable) |>
dplyr::mutate(merror = min(BACC)) |>
dplyr::filter(BACC == merror) |>
dplyr::mutate(row = dplyr::row_number(), nclass = max(row)) |>
dplyr::mutate(sel = ifelse(nclass > 1 & row ==2, 1,
ifelse(nclass == 1 & row ==1, 1, 0))) |>
dplyr::filter(sel == 1) |> dplyr::ungroup() |>
dplyr::select(variable, threshold, BACC)
thresholds <- thresholds[match(colnames(x), thresholds$variable),]
Pr <- matrix(NA, nrow(P), ncol(P))
for(p in 1:ncol(P)){
Pr[,p] <- ifelse( P[,p] >= thresholds$threshold[p], 1, 0)
}
c1 <- 1 - sum(apply((Pr == 1) & (x == 1), 2, sum, na.rm=TRUE))/N1
c2 <- 1 - sum(apply((Pr == 0) & (x == 0), 2, sum, na.rm=TRUE))/N0
BACC <- round(100/2 * (c1 + c2), 2)
out <- list(pred = Pr, thres = thresholds, BACC = BACC)
}
crossval <- function(x, k = 2, K = 7, thres = NULL, method = NULL, type = NULL){
folds <- train_miss_pattern(x, K = K)
missVal <- folds$missWold
train <- folds$Xtrain
err_T = list()
cv_errD = list()
for(i in 1:K){
if(method == "CG" & k > 0){
missBip <- BiplotML::LogBip(train[[i]], k = k,
method = method, type = type, plot = FALSE)
P <- fitted_LB(missBip, type = "response")
}else if(method == "BFGS" & k > 0){
missBip <- BiplotML::LogBip(train[[i]], k = k,
method = method, plot = FALSE)
P <- fitted_LB(missBip, type = "response")
}else if(method == "MM" & k > 0){
missBip <- BiplotML::LogBip(train[[i]], k = k,
method = "MM", random_start = FALSE, plot = FALSE)
P <- fitted_LB(missBip, type = "response")
}else{
P <- rep(1, nrow(x)) %*% t(as.matrix(colMeans(train[[i]], na.rm=TRUE)))
}
Xhat <- matrix(NA, nrow(P), ncol(P))
for(p in 1:ncol(P)){
Xhat[,p] <- ifelse(P[,p] >= thres$threshold[p], 1, 0)
}
Xhat[is.na(Xhat)] <- 0
Xhat_pred <- Xhat[missVal[[i]]]
xReal <- x[missVal[[i]]]
n1 <- sum(xReal); n1t <- sum(x)
n0 <- length(xReal) - n1; n0t <- length(x) - n1t
err0 <- sum(ifelse(xReal == 0 & Xhat_pred == 1, 1, 0))/n0
err1 <- sum(ifelse(xReal == 1 & Xhat_pred == 0, 1, 0))/n1
err0t <- sum(ifelse(x == 0 & Xhat == 1, 1, 0))/n0t
err1t <- sum(ifelse(x == 1 & Xhat == 0, 1, 0))/n1t
err_T[[i]] = 100/2 * (err0t + err1t)
cv_errD[[i]] <- 100/2 * (err0 + err1)
}
cvT <- round(mean(sapply(err_T, mean, na.rm = TRUE), na.rm = TRUE), 2)
cvD <- round(mean(sapply(cv_errD, mean, na.rm = TRUE), na.rm = TRUE), 2)
out <- list(cvT = cvT, cvD = cvD)
return(out)
}
log_like <- function(x, w = NULL, theta){
n <- nrow(x)
p <- ncol(x)
if (!is.null(w)) {
w = matrix(rep(1, n*p), n, p)
}
e_log_like <- matrix(NA, nrow = n, ncol = p)
P <- plogis(theta)
for (i in 1:nrow(x)) {
for (j in 1:ncol(x)) {
e_log_like[i,j] <- w[i,j] * (x[i,j] * log(P[i, j]) + (1 - x[i,j]) * log(1 - P[i, j]))
}
}
log_like <- -sum(e_log_like)
return(log_like)
}
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.