# This file contains external functions from other packages, trimmed down or optimized for the simulation.
# All authors will be credited.
#' @title .
#' Copyright: https://github.com/cvborkulo/IsingFit
#' @export
IsingFitEssential <- function(x, family ='binomial', AND = TRUE, gamma = 0.25, lowerbound.lambda = NA) {
if (family != "binomial")
stop("This procedure is currently only supported for binary (family='binomial') data")
nvar <- ncol(x)
p <- nvar - 1
intercepts <- betas <- lambdas <- list(vector, nvar)
nlambdas <- rep(0, nvar)
for (i in 1:nvar) {
a <- glmnet::glmnet(x[, -i], x[, i], family = family)
intercepts[[i]] <- a$a0
betas[[i]] <- a$beta
lambdas[[i]] <- a$lambda
nlambdas[i] <- length(lambdas[[i]])
}
P <- logl <- sumlogl <- J <- matrix(0, max(nlambdas), nvar)
for (i in 1:nvar) {
J[1:ncol(betas[[i]]), i] <- colSums(betas[[i]] != 0)
}
logl_M <- P_M <- array(0, dim = c(nrow(x), max(nlambdas), nvar))
N <- nrow(x)
for (i in 1:nvar) {
betas.ii <- as.matrix(betas[[i]])
int.ii <- intercepts[[i]]
y <- matrix(0, nrow = N, ncol = ncol(betas.ii))
xi <- x[, -i]
NB <- nrow(betas.ii)
for (bb in 1:NB) {
y <- y + betas.ii[rep(bb, N), ] * xi[, bb]
}
y <- matrix(int.ii, nrow = N, ncol = ncol(y), byrow = TRUE) + y
n_NA <- max(nlambdas) - ncol(y)
if (n_NA > 0) {
for (vv in 1:n_NA) {
y <- cbind(y, NA)
}
}
P_M[, , i] <- exp(y * x[, i])/(1 + exp(y))
logl_M[, , i] <- log(P_M[, , i])
}
logl_Msum <- colSums(logl_M, 1, na.rm = FALSE)
sumlogl <- logl_Msum
sumlogl[sumlogl == 0] = NA
penalty <- J * log(nrow(x)) + 2 * gamma * J * log(p)
EBIC <- -2 * sumlogl + penalty
lambda.mat <- matrix(NA, nrow(EBIC), ncol(EBIC))
for (i in 1:nvar) {
lambda.mat[, i] <- c(lambdas[[i]], rep(NA, nrow(EBIC) - length(lambdas[[i]])))
}
if (!is.na(lowerbound.lambda)) {
EBIC <- EBIC/(lambda.mat >= lowerbound.lambda) * 1
}
lambda.opt <- apply(EBIC, 2, which.min)
lambda.val <- rep(NA, nvar)
thresholds <- 0
for (i in 1:length(lambda.opt)) {
lambda.val[i] <- lambda.mat[lambda.opt[i], i]
thresholds[i] <- intercepts[[i]][lambda.opt[i]]
}
weights.opt <- matrix(, nvar, nvar)
for (i in 1:nvar) {
weights.opt[i, -i] <- betas[[i]][, lambda.opt[i]]
}
asymm.weights <- weights.opt
diag(asymm.weights) = 0
if (AND == TRUE) {
adj <- weights.opt
adj <- (adj != 0) * 1
EN.weights <- adj * t(adj)
EN.weights <- EN.weights * weights.opt
meanweights.opt <- (EN.weights + t(EN.weights))/2
diag(meanweights.opt) <- 0
}
else {
meanweights.opt <- (weights.opt + t(weights.opt))/2
diag(meanweights.opt) <- 0
}
# Return the main results.
return(
list(
weiadj = meanweights.opt,
thresholds = thresholds
)
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.