pychen.ping@gmail.com

zw12356@gmail.com

chunhaoyang@ufl.edu

smchang@mail.ncku.edu.tw

kuojunglee@mail.ncku.edu.tw

Abstract

The purpose of the proposed package **milr** is to analyze multiple-instance data. Ordinary multiple-instance data consists of many independent bags, and each bag is composed of several instances. The statuses of bags and instances are binary. Moreover, the statuses of instances are not observed, whereas the statuses of bags are observed. The functions in this package are applicable for analyzing multiple-instance data, simulating data via logistic regression, and selecting important covariates in the regression model. To this end, maximum likelihood estimation with an expectation-maximization algorithm is implemented for model estimation, and a lasso penalty added to the likelihood function is applied for variable selection. Additionally, an `milr`

object is applicable to generic functions `fitted`

, `predict`

and `summary`

. Simulated data and a real example are given to demonstrate the features of this package.

Multiple-instance learning (MIL) is used to model the class labels which are associated with bags of observations instead of the individual observations. This technique has been widely used in solving many different real-world problems. In the early stage of the MIL application, Dietterich, Lathrop, and Lozano-Pérez (1997) studied the drug-activity prediction problem. A molecule is classified as a good drug if it is able to bind strongly to a binding site on the target molecule. The problem is: one molecule can adopt multiple shapes called the conformations and only one or a few conformations can bind the target molecule well. They described a molecule by a bag of its many possible conformations whose binding strength remains unknown. An important application of MIL is the image and text categorization, such as in Maron and Ratan (1998), Andrews, Tsochantaridis, and Hofmann (2003), J. Zhang et al. (2007), Zhou, Sun, and Li (2009), W. Li et al. (2011), Kotzias et al. (2015), to name a few. An image (bag) possessing at least one particular pattern (instance) is categorized into one class; otherwise, it is categorized into another class. For example, Maron and Ratan (1998) treated the natural scene images as bags, and, each bag is categorized as the scene of waterfall if at least one of its subimages is the waterfall. Whereas, Zhou, Sun, and Li (2009) studied the categorization of collections (bags) of posts (instances) from different newsgroups corpus. A collection is a positive bag if it contains 3% posts from a target corpus category and the remaining 97% posts, as well as all posts in the negative bags, belong to the other corpus categories. MIL is also used in medical researches. The UCSB breast cancer study (Kandemir, Zhang, and A. (2014)) is such a case. Patients (bags) were diagnosed as having or not having cancer by doctors; however, the computer, initially, had no knowledge of which patterns (instances) were associated with the disease. Furthermore, in manufacturing processes (R.-B. Chen et al. (2016)), a product (bag) is defective as long as one or more of its components (instances) are defective. In practice, at the initial stage, we only know that a product is defective, and we have no idea which component is responsible for the defect.

Several approaches have been offered to analyze datasets with multiple instances, e.g., Maron (1998), S. Ray and Craven (2005), X. Xu and Frank (2004), Q. Zhang and Goldman (2002). From our point of view, the statuses of these components are missing variables, and thus, the Expectation-Maximization (EM) algorithm (Dempster, Laird, and Rubin (1977)) can play a role in multiple-instance learning. By now the toolboxes or libraries available for implementing MIL methods are developed by other computer softwares. For example, J. Yang (2008) and Tax and Cheplygina (2016) are implemented in MATLAB software, but neither of them carries the methods based on logistic regression model. Settles, Craven, and Ray (2008) provided the Java codes including the method introduced in S. Ray and Craven (2005). Thus, for R users, we are first to develop a MIL-related package based on logistic regression modelling which is called multiple-instance logistic regression (MILR). In this package, we first apply the logistic regression defined in S. Ray and Craven (2005) and X. Xu and Frank (2004), and then, we use the EM algorithm to obtain maximum likelihood estimates of the regression coefficients. In addition, the popular lasso penalty (R. Tibshirani (1996)) is applied to the likelihood function so that parameter estimation and variable selection can be performed simultaneously. This feature is especially desirable when the number of covariates is relatively large.

To fix ideas, we firstly define the notations and introduce the construction of the likelihood function. Suppose that the dataset consists of \(n\) bags and that there are \(m_i\) instances in the \(i\)th bag for \(i=1,\dots, n\). Let \(Z_i\) denote the status of the \(i\)th bag, and let \(Y_{ij}\) be the status of the \(j\)th instance in the \(i\)th bag along with \(x_{ij} \in \Re^p\) as the corresponding covariates. We assume that the \(Y_{ij}\) follow independent Bernoulli distributions with defect rates of \(p_{ij}\), where \(p_{ij}=g\left(\beta_0+x_{ij}^T\beta\right)\) and \(g(x) = 1/\left(1+e^{-x}\right)\). We also assume that the \(Z_i\) follow independent Bernoulli distributions with defect rates of \(\pi_i\). Therefore, the bag-level likelihood function is

\[\begin{equation}\label{eq:L} L\left(\beta_0,\beta\right)=\prod_{i=1}^n\pi_i^{z_i}\left(1-\pi_i\right)^{1-z_i}. \end{equation}\]To associate the bag-level defect rate \(\pi_i\) with the instance-level defect rates \(p_{ij}\), several methods have been proposed. The bag-level status is defined as \(Z_i=I\left(\sum_{j=1}^{m_i}Y_{ij}>0\right)\). If the independence assumption among the \(Y_{ij}\) holds, the bag-level defect rate is \(\pi_i=1-\prod_{j=1}^{m_i}(1-p_{ij})\). On the other hand, if the independence assumption might not be held, X. Xu and Frank (2004) and S. Ray and Craven (2005) proposed the softmax function to associate \(\pi_i\) to \(p_{ij}\), as follows:

\[\begin{equation}\label{eq:softmax} s_i\left(\alpha\right)=\sum_{j=1}^{m_i}p_{ij}\exp{\left\{\alpha p_{ij}\right\}} \Big/ \sum_{j=1}^{m_i}\exp{\left\{\alpha p_{ij}\right\}}, \end{equation}\]where \(\alpha\) is a pre-specified nonnegative value. X. Xu and Frank (2004) used \(\alpha=0\), therein modeling \(\pi_i\) by taking the average of \(p_{ij}\), \(j=1,\ldots,m_i\), whereas S. Ray and Craven (2005) suggested \(\alpha=3\). We observe that the likelihood () applying neither the \(\pi_i\) function nor the \(s_i(\alpha)\) function results in effective estimators.

Below, we begin by establishing the E-steps and M-steps required for the EM algorithm and then attach the lasso penalty for the estimation and feature selection. Several computation strategies applied are the same as those addressed in Friedman, Hastie, and Tibshirani (2010). Finally, we demonstrate the functions provided in the **milr** package via simulations and on a real dataset.

If the instance-level statuses, \(y_{ij}\), are observable, the complete data likelihood is \[\prod_{i=1}^n\prod_{j=1}^{m_i}p_{ij}^{y_{ij}}q_{ij}^{1-y_{ij}}~,\] where \(q_{ij}=1-p_{ij}\). An ordinary approach, such as the Newton method, can be used to solve this maximal likelihood estimate (MLE). However, considering multiple-instance data, we can only observe the statuses of the bags, \(Z_i=I\left(\sum_{j=1}^{m_j}Y_{ij}>0\right)\), and not the statuses of the instances \(Y_{ij}\). As a result, we apply the EM algorithm to obtain the MLEs of the parameters by treating the instance-level labels as the missing data.

In the E-step, two conditional distributions of the missing data given the bag-level statuses \(Z_i\) are \[Pr\left(Y_{i1}=0,\ldots,Y_{im_i}=0\mid Z_i=0\right)=1\] and \[ Pr\left(Y_{ij}=y_{ij}, \quad j=1,\dots, m_i \mid Z_i=1\right) = \frac{ \prod_{j=1}^{m_i}p_{ij}^{y_{ij}}q_{ij}^{1-y_{ij}}\times I\left(\sum_{j=1}^{m_i}y_{ij}>0\right) }{1-\prod_{l=1}^{m_i}q_{il}}. \] Thus, the conditional expectations are

\[\begin{equation*} E\left(Y_{ij}\mid Z_i=0\right)=0 \quad \mbox{ and } \quad E\left(Y_{ij}\mid Z_i=1\right)=\frac{p_{ij}}{1-\prod_{l=1}^{m_i}q_{il}}\equiv\gamma_{ij}. \end{equation*}\]The \(Q\) function at step \(t\) is \(Q\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right) = \sum_{i=1}^nQ_i\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right)\), where \(Q_i\) is the conditional expectation of the complete log-likelihood for the \(i\)th bag given \(Z_i\), which is defined as

\[\begin{align*} Q_i\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right) & = E\left(\sum_{j=1}^{m_i}y_{ij}\log{\left(p_{ij}\right)}+\left(1-y_{ij}\right)\log{\left(q_{ij}\right)} ~\Bigg|~ Z_i=z_i,\beta_0^t,\beta^t\right) \\ & = \sum_{j=1}^{m_i}z_i\gamma_{ij}^t\left(\beta_0+x_{ij}^T\beta\right)-\log{\left(1+e^{\beta_0+x_{ij}^T\beta}\right)}. \end{align*}\]Note that all the \(p_{ij}\), \(q_{ij}\), and \(\gamma_{ij}\) are functions of \(\beta_0\) and \(\beta\), and thus, we define these functions by substituting \(\beta_0\) and \(\beta\) by their current estimates \(\beta_0^t\) and \(\beta^t\) to obtain \(p_{ij}^t\), \(q_{ij}^t\), and \(\gamma_{ij}^t\), respectively.

In the M-step, we maximize this \(Q\) function with respect to \(\left(\beta_0, \beta\right)\). Since the maximization of the nonlinear \(Q\) function is computationally expensive, following Friedman, Hastie, and Tibshirani (2010), the quadratic approximation to \(Q\) is applied. Taking the second-order Taylor expansion about \(\beta_0^t\) and \(\beta^t\), we have \(Q\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right) =Q_Q\left(\beta_0,\beta\mid \beta_0^t,\beta^t\right) + C + R_2\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right)\), where \(C\) is a constant in terms of \(\beta_0\) and \(\beta\), \(R_2\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right)\) is the remainder term of the expansion and \[
Q_Q\left(\beta_0,\beta\mid \beta_0^t,\beta^t\right) =
-\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^{m_i}w_{ij}^t\left[u_{ij}^t-\beta_0-x_{ij}^T\beta\right]^2,
\] where \(u_{ij}^t=\beta_0+x_{ij}^T\beta^t+\left(z_i\gamma^t_{ij}-p_{ij}^t\right)\Big/\left(p_{ij}^tq_{ij}^t\right)\) and \(w_{ij}^t=p_{ij}^tq_{ij}^t\). In the **milr** package, instead of maximizing \(Q\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right)\), we maximize its quadratic approximation, \(Q_Q\left(\beta_0,\beta\mid\beta_0^t,\beta^t\right)\). Since the objective function is quadratic, the roots of \(\partial Q_Q / \partial \beta_0\) and \(\partial Q_Q / \partial \beta\) have closed-form representations.

We adopt the lasso method (R. Tibshirani (1996)) to identify active features in this MILR framework. The key is to add the \(L_1\) penalty into the objective function in the M-step so that the EM algorithm is capable of performing estimation and variable selection simultaneously. To this end, we rewrite the objective function as

\[\begin{equation}\label{eq:lasso} \underset{\beta_0,\beta}{\min}\left\{-Q_Q\left(\beta_0,\beta\mid \beta_0^t,\beta^t\right)+\lambda\sum_{k=1}^p\left|\beta_k\right|\right\}. \end{equation}\]Note that the intercept term \(\beta_0\) is always kept in the model; thus, we do not place a penalty on \(\beta_0\). In addition, \(\lambda\) is the tuning parameter, and we will introduce how to determine this parameter later. We applied the shooting algorithm (Fu (1998), milr_paper) to update \(\left(\beta^t_0,\beta^t\right)\).

The **milr** package contains a data generator, `DGP`

, which is used to generate the multiple-instance data for the simulation studies, and two estimation approaches, `milr`

and `softmax`

, which are the main tools for modeling the multiple-instance data. In this section, we introduce the usage and default setups of these

The function `DGP`

is the generator for the multiple-instance-type data under the MILR framework.

To use the `DGP`

function, the user needs to specify an integer `n`

as the number of bags, a vector `m`

of length \(n\) as the number of instances in each bag, and a vector `beta`

of length \(p\), with the desired number of covariates, and the regression coefficients, \(\beta\), as in `DGP(n, m, beta)`

. Note that one can set `m`

as an integer for generating the data with an equal instance size `m`

for each bag. Thus, the total number of observations is \(N=\sum_{i=1}^n m_i\). The `DGP`

simulates the labels of bags through the following steps:

In the **milr** package, we provide two approaches to model the multiple-instance data: the proposed `milr`

(R.-B. Chen et al. (2016)) and the `softmax`

approach (X. Xu and Frank (2004)). To implement these two approaches, we assume that the number of observations and covariates are \(N\) and \(p\), respectively. The input data for both `milr`

and `softmax`

are separated into three parts: the bag-level statuses, `y`

, as a vector of length \(N\); the \(N\times p\) design matrix, `x`

; and `bag`

, the vector of indices of length \(N\), representing the indices of the bag to which each instance belongs.

```
milr(y, x, bag, lambda, numLambda, lambdaCriterion, nfold, maxit)
softmax(y, x, bag, alpha, ...)
```

For the `milr`

function, specifying `lambda`

in different ways controls whether and how the lasso penalty participates in parameter estimation. The default value of `lambda`

is \(0\). With this value, the ordinary MLE is applied, i.e., no penalty term is considered. This is the suggested choice when the number of covariates \(p\) is small. When \(p\) is large or when variable selection is desired, users can specify a \(\lambda\) vector of length \(\kappa\); otherwise, by letting `lambda = -1`

, the program automatically provides a \(\lambda\) vector of length \(\kappa=\)`numLambda`

as the tuning set. Following Friedman, Hastie, and Tibshirani (2010), the theoretical maximal value of \(\lambda\) in () is

The automatically specified sequence of \(\lambda\) values ranges from \(\lambda_{min}=\lambda_{max}/1000\) to \(\lambda_{max}\) in ascending order.

The default setting for choosing the optimal \(\lambda\) among these \(\lambda\) values is the Bayesian information criterion (BIC), \(-2\log{(likelihood)} + p^*\times\log{(n)}\), where \(p^*\) is the number of nonzero regression coefficients. Alternatively, the user can use the options `lambdaCriterion = "deviance"`

and `nfold = K`

with an integer `K`

to obtain the best \(\lambda\) that minimizes the predictive deviance through ‘bag-wise’ K-fold cross validation. The last option, `maxit`

, indicates the maximal number of iterations of the EM algorithm; its default value is 500.

For the `softmax`

function, the option `alpha`

is a nonnegative real number for the \(\alpha\) value in (). The maximum likelihood estimators of the regression coefficients are obtained by the generic function `optim`

. Note that no variable selection approach is implemented for this method.

Two generic accessory functions, `coef`

and `fitted`

, can be used to extract the regression coefficients and the fitted bag-level labels returned by `milr`

and `softmax`

. We also provide the significance test based on Wald’s test for the `milr`

estimations without the lasso penalty through the `summary`

function. In addition, to predict the bag-level statuses for the new data set, the `predict`

function can be used by assigning three items: `object`

is the fitted model obtained by `milr`

or `softmax`

, `newdata`

is the covariate matrix, and `bag\_newdata`

is the bag indices of the new dataset. Finally, the MIL model can be used to predict the bag-level labels and the instances-level labels. The option `type`

in `fitted`

and `predicted`

functions controls the type of output labels. The default option is `type = "bag"`

which results the bag-level prediction. Otherwise, by setting `type = "instance"`

, the instances-level labels will be presented.

```
fitted(object, type)
predict(object, newdata, bag_newdata, type)
```

We illustrate the usage of the **milr** package via simulated and real examples.

We demonstrate how to apply the `milr`

function for model estimation and variable selection. We simulate data with \(n=50\) bags, each containing \(m=3\) instances and regression coefficients \(\beta = (-2, -1, 1, 2, 0.5, 0, 0, 0, 0, 0)\). Specifically, the first four covariates are important.

```
library(milr)
library(pipeR)
set.seed(99)
# set the size of dataset
numOfBag <- 50
numOfInstsInBag <- 3
# set true coefficients: beta_0, beta_1, beta_2, beta_3
trueCoefs <- c(-2, -2, -1, 1, 2, 0.5, 0, 0, 0, 0, 0)
trainData <- DGP(numOfBag, numOfInstsInBag, trueCoefs)
colnames(trainData$X) <- paste0("X", 1:ncol(trainData$X))
(instanceResponse <- as.numeric(with(trainData, tapply(Z, ID, any))))
```

```
## [1] 1 1 1 1 1 1 0 1 0 0 0 1 0 1 1 1 1 1 1 1 0 1 1 1 0 0 0 1 0 1 1 0 0 1 1
## [36] 1 1 1 1 1 1 1 0 1 0 0 1 0 0 1
```

Since the number of covariates is small, we then use the `milr`

function to estimate the model parameters with `lambda = 0`

. One can apply `summary`

to produce results including estimates of the regression coefficients and their corresponding standard error, testing statistics and the P-values under Wald’s test. The regression coefficients are returned by the function `coef`

.

```
# fit milr model
milrFit_EST <- milr(trainData$Z, trainData$X, trainData$ID, lambda = 0)
```

`## Lasso-penalty is not used.`

```
# call the Wald test result
summary(milrFit_EST)
```

`## Log-Likelihood: -3.969.`

`## Estimates:`

```
## Estimate Std.Err Z value Pr(>z)
## intercept -5.74156 2.55833 -2.2443 0.02482 *
## X1 -5.75495 2.61345 -2.2021 0.02766 *
## X2 -0.81066 0.98933 -0.8194 0.41255
## X3 5.15616 2.41108 2.1385 0.03247 *
## X4 5.14240 2.45991 2.0905 0.03657 *
## X5 2.26233 1.22981 1.8396 0.06583 .
## X6 -2.55779 1.66536 -1.5359 0.12457
## X7 -0.83306 1.25562 -0.6635 0.50703
## X8 0.38746 0.82876 0.4675 0.64013
## X9 1.18565 0.98515 1.2035 0.22878
## X10 -1.42732 1.14157 -1.2503 0.21118
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```
# call the regression coefficients
coef(milrFit_EST)
```

```
## intercept X1 X2 X3 X4 X5
## -5.7415608 -5.7549500 -0.8106618 5.1561629 5.1424035 2.2623284
## X6 X7 X8 X9 X10
## -2.5577931 -0.8330627 0.3874582 1.1856494 -1.4273210
```

The generic function `table`

builds a contingency table of the counts for comparing the true bag-level statuses and the fitted bag-level statuses (obtained by the option `type = "bag"`

) and the `predict`

function is used to predict the labels of each bag with corresponding covariate \(X\). On the other hand, The fitted and predicted instance-level statuses can also be found by setting `type = "instance"`

in the `fitted`

and `predict`

functions.

`fitted(milrFit_EST, type = "bag") `

```
## [1] 1 1 1 1 1 1 0 1 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1 0 0 1 1
## [36] 1 1 1 1 1 1 1 0 1 0 0 1 0 0 1
```

```
# fitted(milrFit_EST, type = "instance") # instance-level fitted labels
table(DATA = instanceResponse, FITTED = fitted(milrFit_EST, type = "bag"))
```

```
## FITTED
## DATA 0 1
## 0 16 1
## 1 1 32
```

```
# predict for testing data
testData <- DGP(numOfBag, numOfInstsInBag, trueCoefs)
colnames(testData$X) <- paste0("X", 1:ncol(testData$X))
(instanceResponseTest <- as.numeric(with(trainData, tapply(Z, ID, any))))
```

```
## [1] 1 1 1 1 1 1 0 1 0 0 0 1 0 1 1 1 1 1 1 1 0 1 1 1 0 0 0 1 0 1 1 0 0 1 1
## [36] 1 1 1 1 1 1 1 0 1 0 0 1 0 0 1
```

```
pred_EST <- with(testData, predict(milrFit_EST, X, ID, type = "bag"))
# predict(milrFit_EST, testData$X, testData$ID,
# type = "instance") # instance-level prediction
table(DATA = instanceResponseTest, PRED = pred_EST)
```

```
## PRED
## DATA 0 1
## 0 6 11
## 1 13 20
```

Next, the \(n < p\) cases are considered. We generate a data set with \(n=50\) bags, each with 3 instances and \(p=100\) covariates. Among these covariates, only the first five of them, \(X_1,\ldots,X_5\), are active and their nonzero coefficients are the same as the previous example. First, we manually specify 50 \(\lambda\) values manually from 0.01 to 50. The `milr`

function chooses the best tuning parameter which results in the smallest BIC. For this dataset, the chosen model is a constant model.

```
set.seed(99)
# Set the new coefficienct vector (large p)
trueCoefs_Lp <- c(-2, -2, -1, 1, 2, 0.5, rep(0, 95))
# Generate the new training data with large p
trainData_Lp <- DGP(numOfBag, numOfInstsInBag, trueCoefs_Lp)
colnames(trainData_Lp$X) <- paste0("X", 1:ncol(trainData_Lp$X))
# variable selection by user-defined tuning set
lambdaSet <- exp(seq(log(0.01), log(50), length = 50))
milrFit_VS <- with(trainData_Lp, milr(Z, X, ID, lambda = lambdaSet))
# grep the active factors and their corresponding coefficients
coef(milrFit_VS) %>>% `[`(abs(.) > 0)
```

```
## intercept X1 X4 X20 X63 X84
## -1.51632737 -1.05288885 1.18563408 -0.04499043 0.16079512 -0.25585103
```

Second, we try the auto-tuning feature implemented in `milr`

by assigning `lambda = -1`

. The total number of tuning \(\lambda\) values is indicated by setting `nlambda`

. The following example shows the result of the best model chosen among 50 \(\lambda\) values. The slice `$lambda`

shows the auto-tuned \(\lambda\) candidates and the slice `$BIC`

returns the corresponding value of BIC for every candidate \(\lambda\) value. Again, the chosen model is a constant model.

```
# variable selection using auto-tuning
milrFit_auto_VS <- milr(trainData_Lp$Z, trainData_Lp$X, trainData_Lp$ID,
lambda = -1, numLambda = 50)
# the auto-selected lambda values
milrFit_auto_VS$lambda
```

```
## [1] 0.08366600 0.09633265 0.11091697 0.12770929 0.14704389
## [6] 0.16930566 0.19493775 0.22445043 0.25843120 0.29755649
## [11] 0.34260517 0.39447402 0.45419557 0.52295869 0.60213223
## [16] 0.69329228 0.79825355 0.91910546 1.05825380 1.21846856
## [21] 1.40293909 1.61533761 1.85989230 2.14147143 2.46568036
## [26] 2.83897302 3.26878047 3.76365880 4.33345942 4.98952524
## [31] 5.74491641 6.61467032 7.61610098 8.76914363 10.09675163
## [36] 11.62535337 13.38537839 15.41186309 17.74514826 20.43168207
## [41] 23.52494473 27.08651314 31.18728661 35.90889832 41.34534032
## [46] 47.60483463 54.81198758 63.11027032 72.66487490 83.66600265
```

```
# the values of BIC under each lambda value
milrFit_auto_VS$BIC
```

```
## [1] 204.39082 169.23484 161.52169 142.07592 138.33626 142.48906 134.91614
## [8] 131.26204 127.67262 124.14601 124.59860 117.29284 110.06191 106.83584
## [15] 107.63389 108.56954 101.70868 83.26213 88.55978 86.33643 80.77398
## [22] 71.63801 70.10472 68.81315 71.82667 62.77118 61.82642 57.43703
## [29] 61.35707 56.40626 59.60851 63.82736 65.51461 72.05761 71.21319
## [36] 71.21319 71.21319 71.21319 71.21319 71.21319 71.21319 71.21319
## [43] 71.21319 71.21319 71.21319 71.21319 71.21319 71.21319 71.21319
## [50] 71.21319
```

```
# grep the active factors and their corresponding coefficients
coef(milrFit_auto_VS) %>>% `[`(abs(.) > 0)
```

```
## intercept X1 X4 X84
## -1.3055617 -0.8400178 0.8850002 -0.1157281
```

Instead of using BIC, a better way to choose the proper \(\lambda\) is using the cross validation by setting `lambdaCriterion = "deviance"`

. The following example shows the best model chosen by minimizing the predictive deviance via ‘bag-wise’ 10-fold cross validation. The results of the predictive deviance for every candidate \(\lambda\) can be found in the slice `$cv`

. Twenty-nine covariates were identified including the first four true active covariates, \(X_1,\ldots,X_4\).

```
# variable selection using auto-tuning with cross validation
milrFit_auto_CV <- milr(trainData_Lp$Z, trainData_Lp$X, trainData_Lp$ID,
lambda = -1, numLambda = 50,
lambdaCriterion = "deviance", nfold = 10)
# the values of predictive deviance under each lambda value
milrFit_auto_CV$cv
```

```
## [1] 6.943834 7.405697 7.534674 7.508246 7.697270 7.806679 7.667741
## [8] 7.489969 7.263860 7.040731 6.824969 6.598542 6.391749 6.121951
## [15] 5.956360 5.850343 5.804736 5.429982 5.333415 4.802979 4.426573
## [22] 4.252164 3.948187 3.799091 3.763219 3.844775 3.998868 4.224706
## [29] 4.441538 4.714765 5.069669 5.616141 6.471764 6.732428 6.730117
## [36] 6.730117 6.730117 6.730117 6.730117 6.730117 6.730117 6.730117
## [43] 6.730117 6.730117 6.730117 6.730117 6.730117 6.730117 6.730117
## [50] 6.730117
```

```
# grep the active factors and their corresponding coefficients
coef(milrFit_auto_CV) %>>% `[`(abs(.) > 0)
```

```
## intercept X1 X4 X20 X29
## -1.8280887083 -1.3420366754 1.5776332899 -0.0345905122 0.1278154706
## X63 X65 X74 X84 X95
## 0.4075393849 0.0477684451 0.0042718113 -0.4275612105 -0.0062782663
## X96 X97
## -0.0006204478 0.1415587818
```

According to another simulation study which is not shown in this paper, in contrast to cross-validation, BIC does not perform well for variable selection in terms of multiple-instance logistic regressions. However, it can be an alternative when performing cross-validation is too time consuming.

Hereafter, we denote the proposed method with the lasso penalty by MILR-LASSO for brevity. In the following, we demonstrate the usage of MILR-LASSO and the `softmax`

approach on a real dataset, called MUSK1. The MUSK1 data set consists of 92 molecules (bags) of which 47 are classified as having a musky smell and 45 are classified to be non-musks. The molecules are musky if at least one of their conformers (instances) were responsible for the musky smell. However, knowledge about which conformers are responsible for the musky smell is unknown. There are 166 features that describe the shape, or conformation, of the molecules. The goal is to predict whether a new molecules is musk or non-musk. This dataset is one of the popular benchmark datasets in the field of multiple-instance learning research and one can download the dataset from the following weblink.

```
dataName <- "MIL-Data-2002-Musk-Corel-Trec9.tgz"
dataUrl <- "http://www.cs.columbia.edu/~andrews/mil/data/"
```

We use the `untar`

function to decompress the downloaded file and extract the `MUSK1`

dataset. Then, with the following data preprocessing, we reassemble the `MUSK1`

dataset in a `"data.frame"`

format. The first 2 columns of the `MUSK1`

dataset are the bag indices and the bag-level labels of each observation. Starting with the third column, there are \(p=166\) covariates involved in the `MUSK1`

dataset.

```
filePath <- file.path(getwd(), dataName)
# Download MIL data sets from the url
if (!file.exists(filePath))
download.file(paste0(dataUrl, dataName), filePath)
# Extract MUSK1 data file
if (!dir.exists("MilData"))
untar(filePath, files = "MilData/Musk/musk1norm.svm")
# Read and Preprocess MUSK1
library(data.table)
MUSK1 <- fread("MilData/Musk/musk1norm.svm", header = FALSE) %>>%
`[`(j = lapply(.SD, function(x) gsub("\\d+:(.*)", "\\1", x))) %>>%
`[`(j = c("bag", "label") := tstrsplit(V1, ":")) %>>%
`[`(j = V1 := NULL) %>>% `[`(j = lapply(.SD, as.numeric)) %>>%
`[`(j = `:=`(bag = bag + 1, label = (label + 1)/2)) %>>%
setnames(paste0("V", 2:(ncol(.)-1)), paste0("V", 1:(ncol(.)-2))) %>>%
`[`(j = paste0("V", 1:(ncol(.)-2)) := lapply(.SD, scale),
.SDcols = paste0("V", 1:(ncol(.)-2)))
X <- paste0("V", 1:(ncol(MUSK1) - 2), collapse = "+") %>>%
(paste("~", .)) %>>% as.formula %>>% model.matrix(MUSK1) %>>% `[`( , -1L)
Y <- as.numeric(with(MUSK1, tapply(label, bag, function(x) sum(x) > 0)))
```

To fit an MIL model without variable selection, the **milr** package provides two functions. The first is the `milr`

function with `lambda = 0`

. The second approach is the `softmax`

function with a specific value of `alpha`

. Here, we apply the approaches that have been introduced in X. Xu and Frank (2004) and S. Ray and Craven (2005), called the \(s(0)\) (`alpha=0`

) and \(s(3)\) (`alpha=3`

) methods, respectively. The optimization method in `softmax`

is chosen as the default settings of the generic function `optim`

, that is, the method.

As suggested by one reviewer, it is relevant to compare the computational efficiencies and convergence rates of the `milr`

and `softmax`

functions implemented in this package. Note that, the `milr`

approach is written in C++ and so is the objective function in `softmax`

, and, we only consider their performance affected by their common tuning parameter, `maxit`

, the total number of iterations. For each approach, the total number of iterations are set from 5,000 to 25,000, and, the computation task was performed by a laptop with Intel Core M-5Y71 CPU 1.4 GHz and 8GB RAM. Moreover, the performance in model fitting is assessed based on the classification accuracy. We use the generic function `table`

to produce the contingency tables and calculate the classification accuracy values accordingly.

The left panel of Figure 1 shows the computational cost of each approach along with the increment of the total number of iterations. As expected, the computational cost increases with the number of iterations linearly for both functions. However, the slope for the `milr`

function is much flatter than the slope for the `softmax`

function. A further result of MILR not shown here suggests that, for this dataset, the coefficient estimate of the MILR approach converges between 15,000 and 16,000 iterations. The resulting accuracy of each model is shown in the right panel of Figure 1 which indicates that the MILR approach requires fewer iterations to achieve the best fit.

```
# set the iterations from 5000 to 25000
itSet <- seq(5000, 20000, 3000)
outDT <- data.table(iteration = rep(itSet, each = 3),
method = rep(c("s_0", "s_3", "milr"), length(itSet)),
time = NA_real_, acc = NA_real_)
zzz <- sapply(seq(1L, nrow(outDT), 3L), function(i){
j <- (i-1)/3 + 1
# record the computation time
tmp <- system.time(
softmaxFit_0 <- softmax(MUSK1$label, X, MUSK1$bag, alpha = 0,
control = list(maxit = itSet[j]))
)[3]
set(outDT, i, 3L, tmp)
tmp <- system.time(
softmaxFit_3 <- softmax(MUSK1$label, X, MUSK1$bag, alpha = 3,
control = list(maxit = itSet[j]))
)[3]
set(outDT, i + 1L, 3L, tmp)
tmp <- system.time(
# use a very small lambda so that milr do the estimation
# without evaluating the Hessian matrix
milrFit <- milr(MUSK1$label, X, MUSK1$bag, lambda = 1e-7, maxit = itSet[j])
)[3]
set(outDT, i + 2L, 3L, tmp)
# calculate the accuracy
tmp <- table(DATA = Y, FIT_s0 = fitted(softmaxFit_0, type = "bag"))
set(outDT, i, 4L, sum(diag(tmp))/sum(tmp))
tmp <- table(DATA = Y, FIT_s3 = fitted(softmaxFit_3, type = "bag"))
set(outDT, i + 1L, 4L, sum(diag(tmp))/sum(tmp))
tmp <- table(DATA = Y, FIT_MILR = fitted(milrFit, type = "bag"))
set(outDT, i + 2L, 4L, sum(diag(tmp))/sum(tmp))
invisible(NULL)
})
outDT[ , `:=`(iteration = iteration / 1000, acc = acc * 100)]
library(plyr)
outDT2 <- melt(outDT, 1:2, 3:4) %>>%
`[`(j = variable := mapvalues(variable, c("time", "acc"),
c("CPU Time (sec)", "Accuracy (%)")))
```

For variable selection, we apply the MILR-LASSO approach. First, the tuning parameter set is chosen automatically by setting \(\lambda = -1\), and the best \(\lambda\) value is obtained by minimizing the predictive deviance with 10-fold cross validation among `nlambda = 100`

candidates. In total it costs about 130 seconds to choose the optimal \(\lambda\) value and there are 19 active covariates detected by the MILR-LASSO approach. Using these active covariates, the reduced MILR model performs 89.13% classification accuracy.

```
# MILR-LASSO
milrSV <- milr(MUSK1$label, X, MUSK1$bag, lambda = -1, numLambda = 100,
lambdaCriterion = "deviance", maxit = 16000)
# show the detected active covariates
sv_ind <- names(which(coef(milrSV)[-1L] != 0)) %>>%
(~ print(.)) %>>% match(colnames(X))
```

```
## [1] "V31" "V36" "V37" "V76" "V83" "V105" "V106" "V108" "V109" "V116"
## [11] "V118" "V124" "V126" "V129" "V132" "V136" "V147" "V162" "V163"
```

```
# use a very small lambda so that milr do the estimation
# without evaluating the Hessian matrix
milrREFit <- milr(MUSK1$label, X[ , sv_ind], MUSK1$bag,
lambda = 1e-7, maxit = 16000)
table(DATA = Y, FIT_MILR = fitted(milrREFit, type = "bag"))
```

```
## FIT_MILR
## DATA 0 1
## 0 39 6
## 1 4 43
```

Following the discussion above, we use 10-fold cross validation and compare the prediction accuracy among four MIL models which are \(s(0)\), \(s(3)\), the MILR model with all covariates, and, the MILR model fitted by the selected covariates via MILR-LASSO. The resulting prediction accuracies are 83.70%, 77.17%, 75.00% and 81.52%, respectively.

```
set.seed(99)
predY <- matrix(0, length(Y), 4L) %>>%
`colnames<-`(c("s0","s3","milr","milr_sv"))
folds <- 5
foldBag <- rep(1:folds, floor(length(Y) / folds) + 1,
length = length(Y)) %>>% sample(length(.))
foldIns <- rep(foldBag, table(MUSK1$bag))
for (i in 1:folds) {
# prepare training and testing sets
ind <- which(foldIns == i)
# train models
fit_s0 <- softmax(MUSK1[-ind, ]$label, X[-ind, ], MUSK1[-ind, ]$bag,
alpha = 0, control = list(maxit = 25000))
fit_s3 <- softmax(MUSK1[-ind, ]$label, X[-ind, ], MUSK1[-ind, ]$bag,
alpha = 3, control = list(maxit = 25000))
# milr, use a very small lambda so that milr do the estimation
# without evaluating the Hessian matrix
fit_milr <- milr(MUSK1[-ind, ]$label, X[-ind, ], MUSK1[-ind, ]$bag,
lambda = 1e-7, maxit = 16000)
fit_milr_sv <- milr(MUSK1[-ind, ]$label, X[-ind, sv_ind], MUSK1[-ind, ]$bag,
lambda = 1e-7, maxit = 16000)
# store the predicted labels
ind2 <- which(foldBag == i)
# predict function returns bag response in default
predY[ind2, 1L] <- predict(fit_s0, X[ind, ], MUSK1[ind, ]$bag)
predY[ind2, 2L] <- predict(fit_s3, X[ind, ], MUSK1[ind, ]$bag)
predY[ind2, 3L] <- predict(fit_milr, X[ind, ], MUSK1[ind, ]$bag)
predY[ind2, 4L] <- predict(fit_milr_sv, X[ind, sv_ind], MUSK1[ind, ]$bag)
}
table(DATA = Y, PRED_s0 = predY[ , 1L])
```

```
## PRED_s0
## DATA 0 1
## 0 37 8
## 1 10 37
```

`table(DATA = Y, PRED_s3 = predY[ , 2L])`

```
## PRED_s3
## DATA 0 1
## 0 27 18
## 1 5 42
```

`table(DATA = Y, PRED_MILR = predY[ , 3L])`

```
## PRED_MILR
## DATA 0 1
## 0 29 16
## 1 6 41
```

`table(DATA = Y, PRED_MILR_SV = predY[ , 4L])`

```
## PRED_MILR_SV
## DATA 0 1
## 0 35 10
## 1 12 35
```

This article introduces the usage of the R package **milr** for analyzing multiple-instance data under the framework of logistic regression. In particular, the package contains two approaches: summarizing the mean responses within each bag using the softmax function (X. Xu and Frank (2004), S. Ray and Craven (2005)) and treating the instance-level statuses as hidden information as well as applying the EM algorithm for estimation (R.-B. Chen et al. (2016)). In addition, to estimate the MILR model, a lasso-type variable selection technique is incorporated into the latter approach. The limitations of the developed approaches are as follows. First, we ignore the potential dependency among instance statuses within one bag. Random effects can be incorporated into the proposed logistic regression to represent the dependency. Second, according to our preliminary simulation study, not shown in this paper, the maximum likelihood estimator might be biased when the number of instances in a bag is large, say, \(m_i=100\) or more. Bias reduction methods, such as Firth (1993) and Quenouille (1956), can be applied to alleviate this bias. These attempts are deferred to our future work.

Andrews, Stuart, Ioannis Tsochantaridis, and Thomas Hofmann. 2003. “Support Vector Machines for Multiple-Instance Learning.” In *Advances in Neural Information Processing Systems 15*, 561–68.

Chen, Ray-Bing, Kuang-Hung Cheng, Sheng-Mao Chang, Shuen-Lin Jeng, Ping-Yang Chen, Chun-Hao Yang, and Chi-Chun Hsia. 2016. “Multiple-Instance Logistic Regression with Lasso Penalty.” *ArXiv Preprint ArXiv:1607.03615*.

Dempster, Arthur P, Nan M Laird, and Donald B Rubin. 1977. “Maximum Likelihood from Incomplete Data via the EM Algorithm.” *Journal of the Royal Statistical Society B* 39 (1). JSTOR: 1–38.

Dietterich, Thomas G, Richard H Lathrop, and Tomás Lozano-Pérez. 1997. “Solving the Multiple Instance Problem with Axis-Parallel Rectangles.” *Artificial Intelligence* 89 (1). Elsevier: 31–71.

Firth, David. 1993. “Bias Reduction of Maximum Likelihood Estimates.” *Biometrika* 80 (1). Biometrika Trust: 27–38.

Friedman, Jerome, Trevor Hastie, and Rob Tibshirani. 2010. “Regularization Paths for Generalized Linear Models via Coordinate Descent.” *Journal of Statistical Software* 33 (1). NIH Public Access: 1–22.

Fu, Wenjiang J. 1998. “Penalized Regressions: The Bridge Versus the Lasso.” *Journal of Computational and Graphical Statistics* 7 (3). Taylor & Francis: 397–416.

Kandemir, M., C. Zhang, and Hamprecht F. A. 2014. “Empowering Multiple Instance Histopathology Cancer Diagnosis by Cell Graphs.” In *Medical Image Computing and Computer Assisted Intervention 17 (Pt 2)*, 228–35.

Kotzias, Dimitrios, Misha Denil, Nando de Freitas, and Padhraic Smyth. 2015. “From Group to Individual Labels Using Deep Features.” In *ACM Sigkdd International Conference on Knowledge Discovery and Data Mining 21*, 597–606.

Li, Wen, Lixin Duan, Dong Xu, and Ivor Wai-Hung Tsang. 2011. “Text-Based Image Retrieval Using Progressive Multi-Instance Learning.” In *International Conference on Computer Vision ’11*, 2049–55.

Maron, Oded. 1998. “Learning from Ambiguity.” PhD thesis, Massachusetts Institute of Technology.

Maron, Oded, and Aparna Lakshmi Ratan. 1998. “Multiple-Instance Learning for Natural Scene Classification.” In *International Conference on Machine Learning 15*, 341–49.

Quenouille, Maurice H. 1956. “Notes on Bias in Estimation.” *Biometrika* 43 (3/4). JSTOR: 353–60.

Ray, Soumya, and Mark Craven. 2005. “Supervised Versus Multiple Instance Learning: An Empirical Comparison.” In *International Conference on Machine Learning 22*, 697–704.

Settles, B., M. Craven, and S. Ray. 2008. “Multiple-Instance Active Learning.” In *Advances in Neural Information Processing Systems*, 1289–96.

Tax, D. M. J., and V. Cheplygina. 2016. “MIL, a Matlab Toolbox for Multiple Instance Learning.” http://prlab.tudelft.nl/david-tax/mil.html.

Tibshirani, Robert. 1996. “Regression Shrinkage and Selection via the Lasso.” *Journal of the Royal Statistical Society B* 58 (1). JSTOR: 267–88.

Xu, Xin, and Eibe Frank. 2004. “Logistic Regression and Boosting for Labeled Bags of Instances.” In *Pacific-Asia Conference on Knowledge Discovery and Data Mining ’04*, 272–81.

Yang, J. 2008. “MILL: A Multiple Instance Learning Library.” http://www.cs.cmu.edu/~juny/MILL/.

Zhang, Jianguo, Marcin Marszałek, Svetlana Lazebnik, and Cordelia Schmid. 2007. “Local Features and Kernels for Classification of Texture and Object Categories: A Comprehensive Study.” *International Journal of Computer Vision* 73 (2). Springer-Verlag: 213–38.

Zhang, Qi, and Sally A Goldman. 2002. “EM-DD: An Improved Multiple-Instance Learning Technique.” In *Advances in Neural Information Processing Systems 14*, 1073–80.

Zhou, Zhi-Hua, Yu-Yin Sun, and Yu-Feng Li. 2009. “Multi-Instance Learning by Treating Instances as Non-Iid Samples.” In *Annual International Conference on Machine Learning 26*, 1249–56.