#' @export
#' @method coef milr
coef.milr <- function(object, ...) {
return(object$best_model$coeffiecents)
}
#' Fitted Response of milr Fits
#'
#' @param object A fitted obejct of class inheriting from \code{"milr"}.
#' @param type The type of fitted response required. Default is \code{"bag"}, the fitted labels of bags.
#' The \code{"instance"} option returns the fitted labels of instances.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @method fitted milr
#' @importFrom stats fitted
fitted.milr <- function(object, type = "bag", ...) {
stopifnot(length(type) == 1)
if (type == "bag") {
return(object$best_model$fitted$bag)
} else if (type == "instance") {
return(object$best_model$fitted$instance)
}
}
#' Predict Method for milr Fits
#'
#' @param object A fitted obejct of class inheriting from \code{"milr"}.
#' @param newdata Default is \code{NULL}. A matrix with variables to predict.
#' @param bag_newdata Default is \code{NULL}. A vector. The labels of instances to bags.
#' If \code{newdata} and \code{bag_newdata} both are \code{NULL}, return the fitted result.
#' @param type The type of prediction required. Default is \code{"bag"}, the predicted labels of bags.
#' The \code{"instance"} option returns the predicted labels of instances.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @method predict milr
#' @importFrom stats predict
predict.milr <- function(object, newdata = NULL, bag_newdata = NULL, type = "bag", ...) {
stopifnot(length(type) == 1, type %in% c("bag", "instance"))
if (is.null(newdata) && is.null(bag_newdata))
return(fitted(object, type = type))
if (is.null(newdata) && !is.null(bag_newdata))
stop("newdata cannot be NULL!")
if (!is.null(newdata) && is.null(bag_newdata))
stop("bag_newdata cannot be NULL!")
if (type == "bag") {
return(coef(object) %>>% getMilrProb(cbind(1, newdata), bag_newdata) %>>%
`>`(0.5) %>>% as.numeric)
} else if (type == "instance") {
return(logit(cbind(1, newdata), coef(object)) %>>% `>`(0.5) %>>% as.numeric)
}
}
#' @export
#' @method print milr
print.milr <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {
if (length(coef(x)) > 0) {
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2, quote = FALSE)
cat("\n")
}
cat("Residual Deviance: ", format(signif(x$best_model$deviance, digits)), "\n")
cat(" BIC: ", format(signif(x$best_model$BIC, digits)), "\n")
invisible(x)
}
#' @export
#' @method summary milr
#' @importFrom stats pnorm
summary.milr <- function(object, ...) {
if (object$best_model$lambda_chosen == 0) {
summary <- list(loglik = object$best_model$loglik,
beta = object$best_model$coeffiecents,
se = sqrt(object$best_model$var),
z = object$best_model$coeffiecents / sqrt(object$best_model$var),
lambda = object$best_model$lambda_chosen)
summary$pvalue <- pnorm(abs(summary$z), 0, 1, FALSE) * 2
} else {
summary <- list(loglik = object$best_model$loglik,
beta = object$best_model$coeffiecents,
lambda = object$best_model$lambda_chosen)
}
class(summary) <- "summary.milr"
return(summary)
}
#' @export
#' @method print summary.milr
#' @importFrom stats printCoefmat
print.summary.milr <- function(x, digits = max(3L, getOption("digits") - 2L), ...) {
message(sprintf("Log-Likelihood: %.3f.", x$loglik))
if (x$lambda == 0) {
outMat <- cbind(x$beta, x$se, x$z, x$pvalue) %>>%
`colnames<-`(c("Estimate", "Std.Err", "Z value", "Pr(>z)"))
message("Estimates:")
printCoefmat(outMat, digits = digits)
} else {
message(sprintf("Chosen Penalty: %.3f.", x$lambda))
outMat <- cbind(x$beta) %>>% `colnames<-`("Estimate")
message("Estimates:")
printCoefmat(outMat, digits = digits)
}
}
# function of generation of cv index
cvIndex_f <- function(n, fold) {
fold_n <- n %/% fold
rem <- n - fold_n * fold
index <- rep(1:fold, fold_n)
if (rem > 0)
index <- c(index, 1:rem)
return(sample(index, length(index)))
}
#' Maximum likelihood estimation of multiple-instance logistic regression with LASSO penalty
#'
#' Please refer to \link{milr-package}.
#'
#' @param y a vector. Bag-level binary labels.
#' @param x the design matrix. The number of rows of \code{x} must be equal to the length of \code{y}.
#' @param bag a vector, bag id.
#' @param lambda the tuning parameter for LASSO-penalty. If \code{lambda} is a real value number, then the \code{milr}
#' fits the model based on this lambda value. Second, if \code{lambda} is vector, then the optimal lambda value would be
#' be chosen based on the optimality criterion, \code{lambdaCriterion}.
#' Finally, if \code{lambda = -1}, then the optimal lambda value would be chosen automatically.
#' The default is 0.
#' @param numLambda An integer, the maximum length of LASSO-penalty. in atuo-tunning mode
#' (\code{lambda = -1}). The default is 20.
#' @param lambdaCriterion a string, the used optimality criterion for tuning the \code{lambda} value.
#' It can be specified with \code{lambdaCriterion = "BIC"} or \code{lambdaCriterion = "deviance"}.
#' @param nfold an integer, the number of fold for cross-validation to choose the optimal \code{lambda} when
#' \code{lambdaCriterion = "deviance"}.
#' @param maxit an integer, the maximum iteration for the EM algorithm. The default is 500.
#' @return An object with S3 class "milr".
#' \itemize{
#' \item{lambda}{a vector of candidate lambda values.}
#' \item{cv}{a vector of predictive deviance via \code{nfold}-fold cross validation
#' when \code{lambdaCriterion = "deviance"}.}
#' \item{deviance}{a vector of deviance of candidate model for each candidate lambda value.}
#' \item{BIC}{a vector of BIC of candidate model for each candidate lambda value.}
#' \item{best_index}{an integer, indicates the index of the best model among candidate lambda values.}
#' \item{best_model}{a list of the information for the best model including deviance (not cv deviance),
#' BIC, chosen lambda, coefficients, fitted values, log-likelihood and variances of coefficients.}
#' }
#' @examples
#' set.seed(100)
#' beta <- runif(5, -5, 5)
#' trainData <- DGP(40, 3, beta)
#' testData <- DGP(5, 3, beta)
#' # default (not use LASSO)
#' milr_result <- milr(trainData$Z, trainData$X, trainData$ID)
#' coef(milr_result) # coefficients
#' fitted(milr_result) # fitted bag labels
#' fitted(milr_result, type = "instance") # fitted instance labels
#' summary(milr_result) # summary milr
#' predict(milr_result, testData$X, testData$ID) # predicted bag labels
#' predict(milr_result, testData$X, testData$ID, type = "instance") # predicted instance labels
#'
#' # use BIC to choose penalty (not run)
#' #milr_result <- milr(trainData$Z, trainData$X, trainData$ID,
#' # exp(seq(log(0.01), log(50), length = 30)))
#' #coef(milr_result) # coefficients
#' #fitted(milr_result) # fitted bag labels
#' #fitted(milr_result, type = "instance") # fitted instance labels
#' #summary(milr_result) # summary milr
#' #predict(milr_result, testData$X, testData$ID) # predicted bag labels
#' #predict(milr_result, testData$X, testData$ID, type = "instance") # predicted instance labels
#'
#' # use auto-tuning (not run)
#' #milr_result <- milr(trainData$Z, trainData$X, trainData$ID, lambda = -1, numLambda = 20)
#' #coef(milr_result) # coefficients
#' #fitted(milr_result) # fitted bag labels
#' #fitted(milr_result, type = "instance") # fitted instance labels
#' #summary(milr_result) # summary milr
#' #predict(milr_result, testData$X, testData$ID) # predicted bag labels
#' #predict(milr_result, testData$X, testData$ID, type = "instance") # predicted instance labels
#'
#' # use cv in auto-tuning (not run)
#' #milr_result <- milr(trainData$Z, trainData$X, trainData$ID,
#' # lambda = -1, numLambda = 20, lambdaCriterion = "deviance")
#' #coef(milr_result) # coefficients
#' #fitted(milr_result) # fitted bag labels
#' #fitted(milr_result, type = "instance") # fitted instance labels
#' #summary(milr_result) # summary milr
#' #predict(milr_result, testData$X, testData$ID) # predicted bag labels
#' #predict(milr_result, testData$X, testData$ID, type = "instance") # predicted instance labels
#' @importFrom numDeriv hessian
#' @importFrom glmnet glmnet
#' @name milr
#' @rdname milr
#' @export
milr <- function(y, x, bag, lambda = 0, numLambda = 20L, lambdaCriterion = "BIC",
nfold = 10L, maxit = 500L) {
# if x is vector, transform it to matrix
if (is.vector(x))
x <- matrix(x, ncol = 1)
if (!is.matrix(x))
x <- as.matrix(x)
# if column names of x is missing, assign xi
if (is.null(colnames(x)))
colnames(x) <- paste0("x", 1L:ncol(x))
if (!all(y %in% c(0, 1)))
stop("y must be 0 and 1.")
bag <- as.integer(factor(bag))
# input check
alpha <- 1
stopifnot(length(unique(y)) == 2, length(y) == nrow(x),
all(is.finite(y)), is.numeric(y), all(is.finite(x)), is.numeric(x),
length(lambda) >= 1, all(is.finite(lambda)), is.numeric(lambda),
is.finite(alpha), is.numeric(alpha), is.finite(maxit), is.numeric(maxit),
abs(maxit - floor(maxit)) < 1e-4, nfold < nrow(x),
lambdaCriterion %in% c("deviance", "BIC"))
if (length(lambda) == 1 && lambda == -1) {
message("The penalty term is selected automatically with ", numLambda, " candidates.")
m <- table(bag)
zi <- tapply(y, bag, function(x) sum(x) > 0) %>>% as.numeric
lambdaMax <- sqrt(sum(m-1)) * sqrt(sum(m**(1-2*zi)))
lambda <- exp(seq(log(lambdaMax/1000), log(lambdaMax), length = numLambda))
} else if (length(lambda) == 1 && lambda == 0) {
message("Lasso-penalty is not used.")
} else {
message("Use the user-defined lambda vector.")
}
# initial value for coefficients
#init_beta <- coef(glm(y~x))
init_beta <- glmnet(x, y, standardize = T, alpha = 0, lambda = lambda[1]) %>>% coef %>>% as.vector
beta_history <- matrix(NA, ncol(x) + 1, length(lambda) + 1)
beta_history[ , 1] <- init_beta
unique_bag <- unique(bag)
n_bag <- length(unique_bag)
cv <- NULL
if (length(lambda) > 1) {
# save the history of betas along lambda
message(sprintf("Using the criterion %s to choose optimized penalty.", lambdaCriterion))
if (lambdaCriterion == "BIC") {
dev <- vector("numeric", length(lambda))
BIC <- vector("numeric", length(lambda))
for (i in seq_along(lambda)) {
beta_history[, i + 1] <- milr_cpp(y, cbind(1, x), bag, beta_history[, i], lambda[i], alpha, maxit)
dev[i] <- -2 * getLogLikMilr(beta_history[, i + 1], y, cbind(1, x), bag)
BIC[i] <- dev[i] + sum(beta_history[, i + 1] != 0) * log(n_bag)
}
locMin <- which.min(BIC)
} else if (lambdaCriterion == "deviance") {
# calculate the bag label
bagLvl <- tapply(y, bag, function(x) any(x > 0))
# avoid the sample size of one class is less than nfold to retrun error in cvIndex_f.
if (nfold > min(table(bagLvl))) {
message("nfold is bigger than the number of size of one class.")
message(sprintf("Therefore, nfold is changed to %i.", min(table(bagLvl))))
nfold <- min(table(bagLvl))
}
# stratified sampling for cv on bag labels
cvIndex_bag_p <- cvIndex_f(sum(bagLvl == 1), nfold)
cvIndex_bag_n <- cvIndex_f(sum(bagLvl == 0), nfold)
# convert the cv index on bag labels to cv index on samples
cvIndex_sample <- lapply(1:nfold, function(x){
tmp <- c(which(bagLvl == 1)[which(cvIndex_bag_p == x)], which(bagLvl == 0)[which(cvIndex_bag_n == x)])
which(bag %in% unique_bag[tmp])
})
dev_cv <- matrix(NA, length(lambda), nfold)
dev <- vector("numeric", length(lambda))
BIC <- vector("numeric", length(lambda))
cv_betas <- replicate(nfold, beta_history, simplify = FALSE)
for (i in seq_along(lambda)) {
# CV to calculate the deviances
for (j in 1:nfold) {
trainSetIndex <- do.call(c, cvIndex_sample[setdiff(1:nfold, j)])
cv_betas[[j]][ , i + 1] <- milr_cpp(y[trainSetIndex], cbind(1, x[trainSetIndex, ]),
bag[trainSetIndex], cv_betas[[j]][ , i], lambda[i],
alpha, maxit)
dev_cv[i, j] <- -2 * getLogLikMilr(cv_betas[[j]][ , i + 1], y[cvIndex_sample[[j]]],
cbind(1, x[cvIndex_sample[[j]], ]), bag[cvIndex_sample[[j]]])
}
# calculate the beta under lambda
beta_history[, i + 1] <- milr_cpp(y, cbind(1, x), bag, beta_history[, i], lambda[i], alpha, maxit)
dev[i] <- -2 * getLogLikMilr(beta_history[, i + 1], y, cbind(1, x), bag)
BIC[i] <- dev[i] + sum(beta_history[, i + 1] != 0) * log(n_bag)
}
cv <- rowMeans(dev_cv)
locMin <- which.min(cv)
}
lambda_out <- lambda[locMin]
message(sprintf("The chosen penalty throuth %s is %.4f.", lambdaCriterion, lambda[locMin]))
beta <- beta_history[ , locMin + 1]
} else {
beta_history[,2] <- milr_cpp(y, cbind(1, x), bag, init_beta, lambda, alpha, maxit)
beta <- beta_history[,2]
dev <- -2 * getLogLikMilr(beta, y, cbind(1, x), bag)
BIC <- dev + sum(beta != 0) * log(n_bag)
lambda_out <- lambda
locMin <- 1
}
# return beta with names
beta <- as.vector(beta) %>>% `names<-`(c("intercept", colnames(x)))
rownames(beta_history) <- c("intercept", colnames(x))
# fitted response
fit_yij <- (beta %>>% (logit(cbind(1, x), .)) > 0.5) %>>% as.numeric(.)
fit_y <- getMilrProb(beta, cbind(1, x), bag) %>>% `>`(0.5) %>>% as.numeric
# calculate the hessian matrix when without using lasso
if (lambda_out == 0)
bateVar <- -diag(solve(hessian(function(b) getLogLikMilr(b, y, cbind(1, x), bag), beta)))
else
bateVar <- NULL
out <- list(lambda = lambda, cv = cv, deviance = dev, BIC = BIC,
beta = beta_history[, 2L:ncol(beta_history)],
best_index = locMin,
best_model = list(lambda_chosen = lambda_out,
deviance = dev[locMin],
BIC = BIC[locMin],
coeffiecents = beta,
fitted = list(bag = fit_y, instance = fit_yij),
loglik = getLogLikMilr(beta, y, cbind(1, x), bag),
var = bateVar)
)
class(out) <- "milr"
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.