\section{Introduction}
\input{sections/intro.tex}

\section{Method}
\input{sections/methods.tex}

\section{Experiments and Results}
\input{sections/results.tex}

\section{Conclusions}
\input{sections/conclusion.tex}

\bibliographystyle{splncs18}

\input{references.bbl}

\end{document}

\subsection{Segmentation of CVC using modified U-Net}
\label{method:u_net}
To identify regions of interest, we adapted the commonly used U-Net~\cite{ronneberger2015u} for CVC segmentation. While the original loss function works well for large structures, because CVCs are small structures, we used the exponential logarithmic loss function proposed recently~\cite{wong20183d} to address the highly imbalanced label sizes. When the CXR has no CVC of interest, the segmentation output is a blank image, representing no interesting region.

\subsection{Feature Extraction}
\label{sec:feature_ip}
Given the segmentation of possible CVCs in the CXR from the U-Net, we design image processing features describing the segmentation, its relation to prior knowledge of CVC contours, and its relation to key anatomical features of the chest as below. These describe the overall properties of the potential CVCs, even if the initial detected region is imperfect.

\textbf{Spatial prior for each class}: To obtain shape priors, our clinicians annotated the contours of CVCs in hundreds of training images. The CVCs are traced from their anatomical origin of insertion to their tip to give contours that are reflective for different types of CVCs. We averaged the manual annotations per-pixel for each class~(Fig.~\ref{fig:example}(e)) and blurred them spatially to obtain signature spatial priors for each of the CVC classes: left/right PICC, IJ, subclavian and Swan-Ganz catheters as shown in Fig.~\ref{fig:example}(e).

We then performed a pixel-wise multiplication of the segmentation output with the prior for each class and characterized this overlap using an $n$-bin intensity histogram, and histogram of oriented gradients (HoG) features for each class. This informs us how well a CVC segmentation aligns with priors of particular CVC classes.

\textbf{Relation to anatomical features}:
We extracted the segmentation of chest anatomical structures including the clavicles, lungs, heart and mediastinum using an independent U-Net again trained on clinician annotations. We obtained the Euclidean distance distributions of the segmentation relative to the center of these different chest anatomies.
These provide contextual information and distinguish confusing classes such as PICC and subclavian lines.

\textbf{Size and shape properties}: We also characterized the shape and size properties of the positive regions of segmentation: the area, length and width characterize the overall presence/absence of CVCs, the histogram of oriented gradients (HoG) describe the shape contours which are crucial for type classification since different CVC types have different shape signatures, as shown in Fig.~\ref{fig:example}(e).

\subsection{Classification using Random Forests}
The extracted features are used for classification on the presence of CVCs, and identifying the type of CVCs.
We employ a random forest (RF) for each task. The first RF yields a binary presence/absence label, and the second provides a multi-label output, with 4 indicators, one for each type of CVC: PICC, IJ, subclavian and Swan-Ganz. \\

We chose the current architecture after experimenting with different end-to-end deep learning architectures. Feeding spatial priors and anatomical segmentation as additional channels into a VGG-like network failed to learn distinguishing relations. For our specific problem characterized by small area footprint, long and indistinguishable tubular structures with uneven sample sizes, random forests were found to have better generalization.
\subsection{Data}

We worked with the NIH dataset comprising of 112,000 CXRs~\cite{NIHcxr_wang2017chestx}. A subset of these were labelled and annotated for different tasks:

\textbf{CVC Segmentation - Pixel level annotation:}
A random sample of 1500 AP CXRs was selected from the NIH dataset. The CXRs with CVCs were annotated at the pixel level, to provide annotations for
359 IJ lines, 78 subclavian lines, 277 PICC lines, and 32 Swan-Ganz catheters, yielding a total of 608 annotated images of size 512 x 512. The remaining images have no lines. An overlay of the different annotations, based on CVC types, is shown in Fig.~\ref{fig:example}(e).

\textbf{CVC Presence - Whole image level annotation:}
A subset of 3000 images was chosen from the NIH dataset, which a radiologist labeled globally for presence of external medical devices. This resulted in 2381 CXRs with some device present, and 619 CXRs with an absence of any device. Since devices are usually associated with catheters, this label extends to CVC presence.

\textbf{CVC Type - Whole image level annotation:}
A subset of around 16,000 CXRs was sampled from the NIH dataset and a group of radiologists reported on these using a semi-structured template, where one section was dedicated for device reporting.
A recently proposed NLP sentence clustering algorithm~\cite{bechmarkISBI2019} was applied to the device section and a sentence from each resultant cluster was validated manually by two clinicians to derive the global labels for the presence of different CVCs.

This resulted in 10,746 CXRs with at least one type of externally inserted catheter, with 4249 PICC lines, 1651 IJ lines, 201 subclavian lines, 192 Swan-Ganz catheters, and 4453 CXRs with other catheters including airway and drainage tubes.
The related dataset has been released as part of the MICCAI 2019 Multimodal Learning for Clinical Decision Support (ML-CDS) Challenge.

\subsection{U-Net Segmentation}
\label{sec:unet_results}

\begin{figure}[t]
\centering
    \includegraphics[width=\textwidth]{images/Unet_2examples.png}
    \caption{Some example CXR images, overlay of clinician annotations of CVCs, and overlay of U-Net segmentations on the original CXR.}
\label{fig:unet}
\end{figure}

We trained a U-Net for identifying the lines as discussed in Sec.~\ref{method:u_net}, treating the all clinician annotated CVCs as belonging to a positive class without specific distinction on CVC type. We split the pixel-level annotated images into 80\% training and 20\% validation.
We trained the U-Net until convergence using the Adam optimizer with learning rate of $5e^{-5}$ and the exponential loss with the best-performing weights ($w_{\text{Dice}} = 0.8 $ and $w_{\text{Cross}} = 0.2 $)  from~\cite{wong20183d}.

Some segmentation results are shown in Fig.~\ref{fig:unet}.
We observe that sections of the CVCs are missed in some cases, and there are false positives in some confusing areas of the CXRs.

We evaluate the quality of our segmentation by computing the extent of overlap between the ground truth annotations and the U-Net segmentation output on the held-out set. Since CVCs are thin structures, for reliable overlap estimation between automatic and manual CVC segmentation, we enlarged the binary manual segmentations via a dilation operation. With the radius of 2-pixel dilation, 75\% cases have $>$50\% overlap and 84\% cases have $>$40\% overlap and a 5-pixel dilation radius resulted in 80\% and 90\% cases with greater than 50\% and 40\% overlap, respectively.

\setlength{\tabcolsep}{3.5pt}
\begin{table}[t]
\centering
\caption{Results for CVC detection (P: Precision, R: Recall, Acc: Acccuracy, AUC: Area under ROC, DN: DenseNet, VGG:VGG16, SP: spatial prior). Best values for each column are in bold.}
\resizebox{\textwidth}{!}{
\begin{tabular}{l|cccc|l|cccc}
\hline
Method & P & R & Acc & AUC & Method & P & R & Acc & AUC\\ \hline
1. DN & {20.0} & {20.0} & {20.0} & {50.0} & 8. Seg-VGG-RF & {80.2} & {95.0} & {78.2} & {53.6} \\
2. VGG & {20.0} & {20.0} & {20.0} & {50.0} & 9. + CXR & {79.2} & {97.2} & {79.8} & {49.8} \\
3. DN+VGG & {81.6} & {95.0} & {75.3} & {63.0} & 10. Mask-DN-RF & {84.0} & {92.4} & {79.4} & {62.0} \\
4. DN-RF & 79.4 & \textbf{98.8} & 78.8 & 50.2 & 11. Mask-VGG-RF & 84.2 & 91.6 & 79.6 & 62.8 \\
5. VGG-RF & 79.4 & 98.0 & 78.4 & 50.0 & 12. Seg-SP-RF (Ours) & 89.6 & {89.6} & {83.8} & 75.0 \\
6. Seg-DN-RF & 80.0 & 95.8 & 77.8 & 58.6 & 13. + HoG (Ours) & 91.4 & 89.4 & \textbf{85.2} & 78.8 \\
7. + CXR & 79.6 & 97.2 & 78.4 & 51.4 & 14. + Anatomy (Ours) & \textbf{91.6} & {89.6} & \textbf{85.2} & \textbf{79.4} \\ \hline
\end{tabular}
}
\label{table:presence}
\end{table}

\subsection{CVC Presence Identification}
We next focus on the task of identifying CVC presence in CXRs. We performed a 5-fold cross-validation using a 60-20-20 split for training-validation-testing using the 3000 CXRs labeled for device presence. The classifier outputs a binary label indicating the presence of at least one CVC in the CXR (label: 1), or the absence of any CVCs (label: 0). Results for this task are presented in Table~\ref{table:presence}. The parameters (number of trees, and depth of trees) for all random forests were chosen with hyper-parameter tuning for validation performance, for each fold.

To set the baseline, we trained state-of-the-art VGG16~\cite{simonyan2014very} and  DenseNet~\cite{huang2017densely} networks
for the classification task directly on the CXRs, fine-tuning their weights pre-trained on ImageNet. We observed that this yields poor classifiers, with less than 50\% accuracy (Items 1-2).
Concatenating the features from DenseNet and VGG16 and performing heavy hyper-parameter tuning results in a moderately improved performance (Item 3). Overall, the networks were unable to recognize the discriminative regions and performed poorly, due to the small area footprint, long tubular structures that blend into the background, and uneven sample sizes.
 Thus, in further experiments, we treated the ImageNet pre-trained neural networks as feature extractors, feeding the pre-final layer outputs to random forest (RF) classifiers, henceforth abbreviated NN-RF.
Using the original CXR as input to the NN-RF improved accuracy, while the area under ROC (AUC) still remained at 50\% (Items 4-5).

Next, we processed all the CXRs through the U-Net generating segmentation images. We fed combinations of the segmentation (Seg) and the original CXR image to NN-RFs: Seg alone (Items 6, 8), Seg with CXR as one of the image channels (Items 7, 9) and the original CXR masked to zero out the lowest values of the segmentation output to create a masked CXR image focused on regions of potential CVCs (Items 10-11). These showed considerable improvements in the AUC, while the other metrics (precision, recall, F-score and accuracy) remained primarily unchanged.

Finally, we extracted the image-processing features describing the size, shape, likelihood based on CVC spatial priors, and relation to chest anatomical elements, as presented in Section~\ref{sec:feature_ip}. The simplest set of features comprising the spatial prior information itself yielded a 12\% increase in AUC (Item 12). The further addition of size and HoG shape features (Item 13) and anatomical relation information (Item 14) improved the classifier performance further, giving 85.2\% accuracy at a precision of 91.6\% and a recall of 89.6\%.

\subsection{CVC Type Identification}

For the classification problem of identifying the type of CVC present, we again perform a 5-fold cross-validation using a 60-20-20 split for training-validation-testing. We used the 10,746 CXRs which have at least one catheter.
We report the support-weighted average performance metrics, and the precision and recall for each CVC type in Table~\ref{table:type} using the different methods presented for CVC presence identification (except for DenseNet+VGG concatenation).

As in the previous classification problem, pre-trained DenseNet and VGG16 are unable to perform well (Item 1-2), though treating these as feature extractors for a RF classification improves performance substantially (Item 3-4). Utilizing the features extracted from fine-tuned VGG performs similarly achieving an average 79.6\% precision, 49.2\% recall and 59.5\% accuracy, while DenseNet performs worse.
Segmentations and CXRs input to the networks for feature extraction marginally improve the average accuracy and AUC (Item 5-10), leaving the other metrics unchanged.
The image-processing based features perform the best (Item 11-13), reaching an average classification accuracy of 78.2\% at a precision of 95.2\%.
Our method also achieves lower metric variation as brought out by the standard deviation values of the average accuracy.

From our best-performing classifier (Item 13), it can be observed that recall for CVCs other than PICC is under 50\%.
This reveals that the RF classifier performs best on PICC lines, primarily due to the fact that PICC lines make up about 40\% of the CVC images while having the simplest contours. More complex contours like the Swan-Ganz are under-represented and comprise only 2\% of all CVCs, and make the task challenging. This could be overcome by obtaining more data for these classes.
Some example CXRs with the U-Net output, and the final class prediction are presented in Fig.~\ref{fig:final_examples}. It can be observed that our approach identifies the type of CVC correctly despite the incomplete segmentations from the U-Net. Our method can also handle the presence of multiple CVCs in the same CXR.

\begin{figure*}[t]
\centering
    \includegraphics[width=\textwidth]{images/final_examples.png}
  \caption{Some examples of CXRs with U-Net output overlayed. Despite the incomplete segmentation by the U-Net, our method correctly labels each of the catheters.}
\label{fig:final_examples}
\end{figure*}

\begin{table}[t]
\caption{Results for CVC type identification (Mean $\pm$ standard deviation, P: Precision, R: Recall, AUC: Area under ROC, DN: DenseNet, SP: spatial prior). Best values for each column are in bold. Our algorithm: Rows 11-13.}
\resizebox{\textwidth}{!}{
\begin{tabular}{l|cc|cc|cc|cc||cccc}
 & \multicolumn{2}{c|}{PICC} & \multicolumn{2}{c|}{IJ} & \multicolumn{2}{c|}{SC} & \multicolumn{2}{c|}{SG} & \multicolumn{4}{c}{Weighted Average} \\ \hline
Method & P & R & P & R & P & R & P & R & P & R & Accuracy & AUC \\ \hline
1. DN & 34.6 & 11.6 & 13.2 & 22.2 & 0.0 & 0.0 & 2.6 & 51.4 & 34.2 & 20.4 & 20.0 $\pm$ 0.71 & 51.0 $\pm$ 1.00 \\
2. VGG & 36.6 & 2.8 & 16.6 & 8.8 & 1.2 & 16.8 & 2.8 & 6.2 & 36.0 & 27.6 & 28.0 $\pm$ 0.00 & 48.6 $\pm$ 0.89 \\ \hline
3. DN-RF & 78.4 & 55.0 & 96.6 & 32.4 & \textbf{100.0} & \textbf{33.8} & \textbf{100.0} & \textbf{25.8} & 84.6 & 47.4 & 66.2 $\pm$ 1.10 & 67.0 $\pm$ 0.71 \\
4. VGG-RF & 77.2 & 54.6 & 96.4 & 32.2 & \textbf{100.0} & 33.2 & \textbf{100.0} & \textbf{25.8} & 83.8 & 47.2 & 65.6 $\pm$ 0.89 & 66.8 $\pm$ 0.84 \\ \hline
5. Seg-DN-RF & 77.4 & 64.8 & 84.2 & \textbf{38.2} & \textbf{100.0} & 33.2 & \textbf{100.0} & 24.6 & 80.8 & 55.4 & 68.6 $\pm$ 0.89 & 68.6 $\pm$ 0.89 \\
6. +CXR & 78.2 & 52.2 & 98.2 & 33.2 & \textbf{100.0} & \textbf{33.8} & \textbf{100.0} & 25.2 & 85.0 & 45.8 & 65.4 $\pm$ 1.14 & 66.8 $\pm$ 0.45 \\ \hline
7. Seg-VGG-RF & 76.4 & 63.8 & 82.4 & 37.6 & \textbf{100.0} & 32.6 & \textbf{100.0} & \textbf{25.8} & 79.2 & 54.8 & 67.4 $\pm$ 0.55 & 68.2 $\pm$ 0.84 \\
8. +CXR & 76.2 & 52.6 & 96.2 & 32.8 & \textbf{100.0} & 33.2 & \textbf{100.0} & 23.8 & 83.2 & 45.8 & 64.8 $\pm$ 1.30 & 66.6 $\pm$ 0.55 \\ \hline
9. Mask-DN-RF & 75.6 & 60.0 & 88.6 & 35.6 & \textbf{100.0} & \textbf{33.8} & \textbf{100.0} & 25.2 & 80.6 & 51.6 & 66.6 $\pm$ 0.89 & 67.8 $\pm$ 0.45 \\
10. Mask-VGG-RF & 74.6 & 58.0 & 87.2 & 34.8 & \textbf{100.0} & \textbf{33.8} & \textbf{100.0} & 25.2 & 79.2 & 50.2 & 65.4 $\pm$ 1.14 & 67.2 $\pm$ 0.84 \\ \hline
11. Seg-SP-RF  & 88.6 & 75.4 & 95.8 & 33.4 & \textbf{100.0} & 31.6 & \textbf{100.0} & \textbf{25.8} & 91.2 & 61.6 & 75.6  $\pm$ 0.55 & 70.0  $\pm$ 0.71 \\
12. + HoG  & 93.0 & 78.8 & 99.0 & 32.4 & \textbf{100.0} & 33.2 & \textbf{100.0} & \textbf{25.8} & 95.2 & \textbf{63.8} & 78.2 $\pm$ 0.84 & 70.8 $\pm$ 0.84 \\
13. + Anatomy  & \textbf{93.0} & \textbf{78.8} & \textbf{99.6} & 32.2 & \textbf{100.0} & \textbf{33.8} & \textbf{100.0} & \textbf{25.8} & \textbf{95.6} & 63.6 & \textbf{78.4} $\pm$ \textbf{0.55} & \textbf{71.0} $\pm$ \textbf{0.71} \\ \hline
\end{tabular}}
\label{table:type}
\end{table}


\section{Introduction}
\label{sec:intro}
		Multi-parameter regression (MPR) modelling refers to the approach whereby covariates are allowed to enter the model through multiple distributional parameters simultaneously. This is in contrast to the standard approaches where covariates enter through a single parameter (e.g., a location parameter) while holding the remaining parameter(s) (e.g., a scale or shape parameters) constant. Multi-parameter regression approaches have been considered in many areas of applied statistics. One of the earliest examples appears in econometrics literature where \citet{park1966estimation} proposed a log-linear model for the scale parameter in the normal linear model when the assumption of homoscedasticity is violated and describes a two stage process for the estimation of the parameters. \citet{harvey1976estimating} develops a maximum likelihood estimation procedure for these so called heteroscedastic regression models. \citet{stirling1985heteroscedastic} and \citet{aitkin1987modelling} illustrate the procedure using examples and GLIM (Generalized Linear Interactive Modelling) macros for the implementation of these models. Other examples of multi-parameter regression models include an extension of generalized linear models to allow for the joint modelling of the mean and dispersion \citep{smyth1989generalized, nelder1991generalized, lee1998generalized}, the generalized additive models for location, scale and shape (GAMLSS) which goes beyond the exponential family to general parameters (not necessarily location and scale) for a variety of distributions \citep{rigby2005generalized}. \citet{taylor2004joint} model the location and scale parameters of the t-distribution jointly to rectify the dependence of the scale parameter on the location.
		\par
		In this article, we consider MPR modelling in the setting of survival analysis. Examples of the MPR modelling framework include \citet{anderson1991nonproportional} who extended the Weibull accelerated failure time model such that both the location and dispersion parameters depend on covariates. A set of models known as the threshold regression models, model survival time through an underlying diffusion process (e.g., Weiner process) whose drift and initial state parameters both depend on covariates \citep{lee2006threshold, aalen2008survival}. More recently, \citet{burke2017multi} explored the general use of parameteric hazard MPR models, and demonstrated the favourable interpretation and flexibility afforded by jointly modelling the scale and shape parameters of a Weibull distribution.
		\par
		The variable selection problem is omnipresent in most, if not all, statistical applications. This problem arises due to the uncertainty associated with the selection of a subset of explanatory variables to model some variable of interest \citep{thompson1978selection, george2000variable}. With the growing availability of data, this problem has received a significant amount of attention in recent years. Traditionally, variable selection has been performed using methods such as backward elimination, forward selection, and stepwise regression; \citet{miller2002subset} provides a comprehensive discussion of these methods. Due to their inherent discreteness {(i.e., covariates are either ``in'' or ``out'')}, these methods can be unstable \citep{breiman1996heuristics}. Furthermore, they are known to be computationally intensive with a power explosion of the number of possible submodels ($2^{p}$) to be considered for $p$ predictors. Modern approaches have been focused on penalized regression whereby parameter estimation and (continuous) model selection are both carried simultaneously. Methods such as the least absolute shrinkage and selection operator (LASSO) \citep{tibshirani1996regression}, smoothly clipped absolute deviation (SCAD) \citep{fan2001variable}, the elastic net \citep{zou2005regularization}, and adaptive lasso (ALASSO) \citep{zou2006adaptive} are used to simultaneously select variables and estimate their regression coefficients.
		\par
		Standard regression approaches only have one regression component, and, therefore, variable selection literature mainly focuses on this, i.e., selection of covariates in a location parameter (cf.~\citet{fan2010selective}). Beyond this standard setting, MPR models require variable selection in multiple distributional parameters. However, little work has been done in this area. \citet{wu2012variable} considered penalized likelihood for variable selection in the case of an inverse Gaussian joint mean and dispersion model, while \citet{wu2013variable} and \citet{wu2014variable} applied this procedure to joint location-scale skew-normal and skew-t-normal models. None of the aforementioned papers consider censored data, and, furthermore, these papers only consider the case of a common tuning parameter for the two regression components. Therefore, the main objective of this article is to develop penalized variable selection in MPR models more generally. Using the Weibull MPR model as an example, we investigate the need for a separate tuning parameter for each regression component.  We select tuning parameters based on the BIC function, and, because this is multi-modal, we propose the use of a differential evolution ``global'' optimization procedure.
		\par
		The remainder of the paper is structured as follows. In Section 2, the Weibull multi-parameter regression model, the penalty functions and the penalized likelihood estimation procedure are introduced. Section 3 describes the model fitting process and the algorithm used in the selection of the tuning parameter(s). Simulation studies to investigate the characteristics of the algorithm and evaluate its performance in both variable selection and parameter estimation are given in Section 4. The proposed methodologies are then illustrated on a real data example, a lung cancer dataset, in Section 5. A discussion of the proposed methods and some concluding remarks are given in the last section, Section 6.
\section{Model Formulation}
\label{sec:Formulation}
 \subsection{Weibull MPR Model}
		Although the variable selection methods we consider in this article can be applied to any parametric MPR model, it is helpful to focus on a specific example. We therefore consider the Weibull MPR model as the Weibull distribution is one of the most popular parametric survival distributions. {The Weibull hazard function for survival time $\tilde{T}_i$ corresponding to the $i$th individual is given by }
		\begin{equation}\nonumber
		h(t|x_i,z_i)=\tau_i\gamma_i t^{\gamma_i-1} \ ,
		\end{equation}
		for $i = 1, 2, \dots, n$, where $\tau_i>0$ is the scale parameter and $\gamma_i>0$ is the shape parameter. The Weibull MPR model is obtained by letting both distributional parameters depend on covariates as follows:
		$$\log(\tau_i)=x^{T}_{i}\beta, \qquad \log(\gamma_i)=z^{T}_{i}\alpha,$$
		where $x_i=(1,x_{i1},\dots,x_{ip})^T$ and $z_i=(1,z_{i1},\dots,z_{iq})^T$ are scale and shape covariate vectors which may or may not have covariates in common, $\beta=(\beta_0,\beta_1,\dots,\beta_p)^T$ and $\alpha=(\alpha_0,\alpha_1,\dots,\alpha_q)^T$ are the corresponding regression coefficients, and the log link is used to ensure positivity of the parameters.
		\par
		Parameter estimation within the unpenalized MPR model can be carried out in a standard fashion using maximum likelihood. {First, let $T_i = \text{min}(\tilde{T}_i, C_i) $ be the observed survival time for the $i$th individual. Then the associated log-likelihood function given by  }
	 	\begin{equation}\label{eq:NonPen_lik}
		\ell(\theta)= \sum^{n}_{i=1}{\delta_i \{\log{\tau_i} + \log{\gamma_i} + (\gamma_i-1) \log{t_i}\}}- \tau_i t^{\gamma_i}_i,
		\end{equation}
		where $\theta=(\beta^T,\alpha^T)^T$ is the full parameter vector, {$t_i$ is the realization of $T_i$}, and $\delta_i$ is the censoring indicator which takes the value 0 for censored survival times and 1 for uncensored survival times. Beyond the Weibull case we consider here, the likelihood function is $\sum_{i=1}^n \delta_i \log h(t_i|x_i,z_i) - H(t_i|x_i,z_i)$ where $H(t|x_i,z_i) = \int_0^t h(u|x_i,z_i) du$ is the cumulative hazard function.

	\subsection{Penalized Likelihood\label{sec:penalizedlike}}
		Penalized MPR estimation can be developed on the basis of maximising a penalized log-likelihood given by
		\begin{equation}\label{eq:Model1}
		\ell_\lambda(\theta)= \ell(\theta)- n \sum^{p}_{j = 0} J_{\lambda_{\beta_j}}(\lvert\beta_j\rvert) - n \sum^{q}_{j = 0} J_{\lambda_{\alpha_j}}(\lvert\alpha_j\rvert),
		\end{equation}
		where $\ell(\theta)$ is the unpenalized likelihood, $\lambda = (\lambda_{\beta_0}, \lambda_{\beta_1}, \ldots,\lambda_{\beta_p}, \lambda_{\alpha_0}, \lambda_{\alpha_1}, \ldots,\lambda_{\alpha_q})$ is a vector of coefficient-specific tuning parameters, and $J_{\lambda_{\beta_j}}(\cdot)$ and $J_{\lambda_{\alpha_j}}(\cdot)$ are scale and shape penalty functions which we assume have the same functional form (but differ with respect to the tuning parameter). As is standard practice, we assume that the intercepts are not penalized, and, therefore, define $\lambda_{\beta_0} \equiv \lambda_{\alpha_0} \equiv 0$ (rather than, for example, assuming the intercepts are zero as other authors do); we also assume that the covariates are standardized. Although we have defined $\lambda$ quite generally, we will in fact impose constraints on this vector (beyond fixing $\lambda_{\beta_0} \equiv \lambda_{\alpha_0} \equiv 0$) by considering the following possibilities (for $j \ne 0$):
 \begin{enumerate}[label=\roman*)]
   \item single penalty,
   $$\lambda_{\beta_j} = \lambda_{\alpha_j} = \lambda,$$
   \item single adaptive penalty,
   $$\lambda_{\beta_j} = \lambda w_{\beta_j}, \qquad \lambda_{\alpha_j} = \lambda w_{\alpha_j},$$ where $w_{\beta_j}$ and $w_{\alpha_j}$ are predefined weights,
   \item separate non-adaptive penalties,
   $$\lambda_{\beta_j} = \lambda_{\beta}, \qquad \lambda_{\alpha_j} = \lambda_\alpha,$$
   \item separate adaptive penalties,
   $$\lambda_{\beta_j} = \lambda_{\beta} w_{\beta_j}, \qquad \lambda_{\alpha_j} = \lambda_{\alpha} w_{\alpha_j}.$$
		\end{enumerate}
		Note that the only ``adaptive" penalty considered for the purpose of this article is the ALASSO penalty. (i) and (ii) are standard approaches where a single penalty, $\lambda$, applies to the whole vector of parameters. This is reasonable in more standard setting where there is only a $\beta$ vector.  However, in this particular MPR setting, we have two separate distributional parameters, which exist on different scales. For this reason we investigate methods (iii) and (iv) which apply different penalties to the two regression vectors via $\lambda_\beta$ and $\lambda_\alpha$.
		\par
		For the purpose of this article, we consider the most commonly used penalties, namely: the least absolute shrinkage and selection operator (LASSO) \citep{tibshirani1996regression},
		\begin{equation*}
		J_{\lambda_{\theta_j}}(\lvert\theta_j\rvert) = \lambda_{\theta_j}\lvert\theta_j\rvert,
		\end{equation*}
		which, although popular, is known to select too many variables \citep{radchenko2008variable}; the non-convex smoothly clipped absolute deviation (SCAD) \citep{fan2001variable},
		\begin{equation*}
			J_{\lambda_{\theta_j}}(\lvert\theta_j\rvert) =  \left\{
				\begin{array}{ll}
				\lambda_{\theta_j}(\lvert\theta_j\rvert) & \qquad \text{if }\ \lvert\theta_j\rvert \le  \lambda_{\theta_j},\\
				\frac{2a \lambda_{\theta_j} \lvert\theta_j\rvert - \theta_{j}^2 -  \lambda_{\theta_j}^2}{2(a-1)} & \qquad \text{if }\  \lambda_{\theta_j} <  \lvert\theta_j\rvert < a \lambda_{\theta_j},\\
				\frac{\lambda_{\theta_j}^2(a+1)}{2} & \qquad \text{if }\ \lvert\theta_j\rvert \ge a \lambda_{\theta_j},\\
				\end{array}
				\right.
				\end{equation*}
				where $a = 3.7$, and the adaptive LASSO (ALASSO) \citep{zou2006adaptive},
				\begin{equation*}
					J_{\lambda_{\theta_j}}(\lvert\theta_j\rvert) = \lambda_{\theta_j} w \lvert\theta_j\rvert\ ,
				\end{equation*}
				where, typically, {$w=1/\lvert\hat\theta_{0,j}\rvert$ and $\hat\theta_{0,j}$} is a unpenalized estimate of $\theta_j$. These so called adaptive weights are used to apply different penalties to different regression coefficients such that a larger amount of shrinkage is applied to the unimportant variables. Note that, here, we use $\theta_j$ to denote a generic regression coefficient, and $\lambda_{\theta_j}$ is the corresponding tuning parameter. The latter two penalties (i.e., SCAD and ALASSO) are known to possess the oracle property, i.e., the procedure asymptotically identifies the right subset model and estimates the coefficients and covariance matrix as though the true model were known in advance \citep{fan2001variable}.
		\section{Penalized Estimation Procedure}
		\label{sec:Estimation}
 \subsection{Model Fitting}
 We define
		\begin{equation}\label{eq:Est_theta}
			\hat{\theta}_{\lambda} = \underset{\theta}{\text{argmax}}\ \ell_{\lambda}(\theta),
		\end{equation}
		where $\ell_\lambda(\theta)$ is given by (\ref{eq:Model1}). The corresponding score functions are given by
     \begin{align} \label{eq:Pen_Lik}
            	\begin{split}
            	\frac{\partial \ell_\lambda}{\partial \beta} &= \frac{\partial \ell}{\partial \beta} - nV_{\beta} = X^T U_{\beta} - nV_{\beta},\\
            	\frac{\partial \ell_\lambda}{\partial \alpha} &= \frac{\partial \ell}{\partial \alpha} - nV_{\alpha} = Z^T U_{\alpha} - nV_{\alpha},
            	\end{split}
         \end{align}
     where $X$ is an $n\times(p+1)$ matrix whose $i$th row is $x_i$, $Z$ is an $n\times(q+1)$ matrix whose $i$th row is $z_i$; $U_{\beta}$ and $U_{\alpha}$ are vectors of length $n$ such that $U_{\beta i} = \delta_i - \tau_i t^{\gamma_i}_i$ and $U_{\alpha i} = \delta_i (1 + \gamma_i \log{t_i}) - \tau_i\gamma_i t^{\gamma_i}_i \log{t_i}$; $V_\beta$ and $V_\alpha$ are vectors of lengths $p+1$ and $q+1$ respectively, such that, for $j \ge 0$, $V_{\beta,j+1} = d J_{\lambda_{\beta_j}}(\lvert\beta_j\rvert) / d \beta_j
     = J_{\lambda_{\beta_j}}'(\lvert\beta_j\rvert) d \lvert\beta_j\rvert / d\beta_j$ and $V_{\alpha, j+1} =d J_{\lambda_{\alpha_j}}(\lvert\alpha_j\rvert) / d \alpha_j = J_{\lambda_{\alpha_j}}'({\lvert\alpha_j\rvert}) d \lvert\alpha_j\rvert / d\alpha_j$.
     \par
		Note however, the presence of the absolute value function renders the penalty functions non-differentiable at zero. Various algorithms have been developed to overcome this issue including quadratic programming \citep{tibshirani1996regression}, least angle regression (LARS) \citep{efron2004least}, co-ordinate descent \citep{friedman2007pathwise} and the local quadratic approximation \citep{fan2001variable}. In this paper, we take a different approach, and use an extension of the absolute value function given by {$$a(x) = \sqrt{x^2+\epsilon^2}-\epsilon,$$ where $\lim_{\epsilon \to 0} a(x) = \lvert x \rvert $. This yields a differentiable penalty so that standard gradient-based optimization algorithms can be applied straightforwardly and transparently.} Thus, $a'(x) = x/ \sqrt{\epsilon^2 + x^2}$ (which is an approximation of the signum function) and $a''(x) = {\epsilon^2}/{(\epsilon^2 + x^2)^{3/2}}$. Smaller values of $\epsilon$ bring the approximate penalty closer to the original penalty, but also closer to the penalty being non-differentiable; we have found that fixing $\epsilon=10^{-4}$ generally works well. As we use smooth $J(\cdot)$ functions, and $a(x)$ in place of $|x|$, (\ref{eq:Pen_Lik}) is then smooth in the parameters and can therefore be solved using the Netwon-Raphson algorithm.

{ We denote by $I_\lambda(\theta)$ the matrix of second derivatives of $\ell_\lambda(\theta)$, i.e., $-\nabla_{\theta}\nabla_{\theta}^T\ell_\lambda(\theta)$. Then,
\begin{align*} \label{eq:Imat
I_\lambda(\theta)
= I_0(\theta) +
\begin{pmatrix}
             n\Sigma_{\beta}
             & 0 \\
             0 & n\Sigma_{\alpha}
         \end{pmatrix}
=   \begin{pmatrix}
             X^{T} W_{\beta} X + n\Sigma_{\beta}
             & X^{T} W_{\alpha\beta} Z \\
             Z^{T} W_{\alpha\beta} X & Z^{T} W_{\alpha} Z + n\Sigma_{\alpha}
         \end{pmatrix}
\end{align*}
where $I_0(\theta) = -\nabla_{\theta}\nabla_{\theta}^T\ell(\theta)$ is the usual observed information matrix of the unpenalized likelihood; $\Sigma_{\beta}$ and $\Sigma_{\alpha}$ appear due to the penalties,} and are diagonal matrices of dimension $(p+1)\times(p+1)$ and $(q+1)\times(q+1)$, respectively, such that, for $j \ge 0$, $\Sigma_{\beta,j+1,j+1} = d^2 J_{\lambda_{\beta_j}}(\lvert\beta_j\rvert) / {d \beta_j}^2$ and $\Sigma_{\alpha,j+1,j+1} = d^2 J_{\lambda_{\alpha_j}}(\lvert\alpha_j\rvert) / {d \alpha_j}^2$; and $W_{\beta}$, $W_{\alpha}$, and $W_{\alpha\beta}$ are $n \times n$ diagonal matrices whose $i$th diagonal elements are given by $\tau_i t^{\gamma_i}_{i}$, $\{\tau_i t^{\gamma_i}_i (\gamma_i \log{t_i} + 1 ) -\delta_i\}\gamma_i \log{t_i}$, and $\tau_i \gamma_i t^{\gamma_i}_{i} \log{t_i}$ respectively. Thus, following \citet{ha2014variable}, the resulting system of Newton-Raphson equations, which are iteratively solved for {$\theta_\lambda^{(m+1)} = ({\beta_\lambda^{(m+1)}}^T, {\alpha_\lambda^{(m+1)}}^T)^T$}, can be written compactly as
\begin{equation} \label{eq:NR
			\begin{pmatrix}
             X^{T} W_{\beta}^{(m)} X + n\Sigma_{\beta}^{(m)}
             & X^{T} W_{\alpha\beta}^{(m)} Z \\
             Z^{T} W_{\alpha\beta}^{(m)} X & Z^{T} W_{\alpha}^{(m)} Z + n\Sigma_{\alpha}^{(m)}
         \end{pmatrix}

         {
         \begin{pmatrix}
             \beta_\lambda^{(m+1)} - \beta_\lambda^{(m)}\\
             \alpha_\lambda^{(m+1)} - \alpha_\lambda^{(m)}
         \end{pmatrix}}

         =

         \begin{pmatrix}
             X^T U_{\beta}^{(m)} - nV_{\beta}^{(m)}\\
             Z^T U_{\alpha}^{(m)} - nV_{\alpha}^{(m)}
     \end{pmatrix},
\end{equation}
{ where the various elements super-scripted by $(m)$ depend on $\theta_\lambda^{(m)}$, but this dependence is suppressed for notational convenience; we use unpenalized estimates as the initial values in this iterative procedure, i.e., $\theta^{(0)} = \hat\theta_0$. Having obtained the penalized estimates, $\hat\theta_\lambda$, the covariance can be estimated using the sandwich formula
 		\begin{equation} \label{eq:SE_formula}
 			\hat{cov}(\hat{\theta}_{\lambda}) = \{I_\lambda(\hat{\theta}_\lambda)\}^{-1} I_0(\hat{\theta}_\lambda) \{I_\lambda(\hat{\theta}_\lambda)\}^{-1}
		 \end{equation}}
\citep{fan2001variable, fan2002variable,ha2014variable, park2018penalized}. This formula is known to have good accuracy when the sample size is moderate \citep{fan2001variable, fan2002variable}, and its performance in our MPR setting is investigated in Section \ref{sec:Simulations} through simulation studies.

		\subsection{Tuning Parameter Selection}
		The selection of ``optimal'' tuning parameter(s) is typically done through the use of data-driven criteria such as generalized cross-validation (GCV), Akaike information criterion (AIC) or Bayesian information criterion (BIC). GCV and the AIC are known to be asymptotically ``loss-efficient'' and ``selection inconsistent'' variable selection criteria \citep{shao1997asymptotic, yang2005can, wang2009shrinkage}. \citet{wang2007tuning} provide a formal proof that the shrinkage or tuning parameter selected using GCV may not be able to identify the true model consistently for the SCAD estimator in linear models and partially linear models. Instead, they suggest using the BIC and prove its model selection consistency property. A similar conclusion has been reached by \citet{wang2007unified} for the ALASSO. Hence, due to its widely reported superior empirical performance in variable selection, we use a BIC-type criterion to determine the values of the tuning parameter(s), where

		\begin{equation}
		\text{BIC}(\lambda)= -2 \ell(\hat{\theta}_{\lambda}) + e_{\lambda} \log{n},
		\label{eq:BIC}
		\end{equation}
		$\ell(\hat{\theta}_\lambda)$ is the unpenalized likelihood function defined in (\ref{eq:NonPen_lik}), $n$ is the sample size and {$e_\lambda = \text{tr}[\{I_\lambda(\hat{\theta}_\lambda)\}^{-1} I_0(\hat{\theta}_\lambda)]$} is the effective degrees of freedom \citep{ha2007model}; we define
		\begin{equation}
\lambda^{*} = \underset{\lambda}{\text{argmin}}\ \text{BIC}(\lambda).
		\label{eq:lambdastr}
		\end{equation}
{ Note that, as described in Section \ref{sec:penalizedlike}, $\lambda^*$ will either be one-dimensional (when a common penalty is applied to $\beta$ and $\alpha$) or two-dimensional (when separate penalties are applied).}

		\par
		The simplest method to solve this optimization problem is grid search. While it is straightforward to implement, grid search is known to suffer from the curse of dimensionality, i.e., the number of grid points grows exponentially with the dimension.
		Furthermore, if the grid is too coarse, the minimum may be overlooked. This is especially true in the case of a multi-modal function, such as what we are trying to optimize (see Figure \ref{fig:BIC}). To overcome the issues associated with the grid search algorithm, we consider ``global'' optimization algorithms. In an empirical comparison of various algorithms for continuous global optimization, \citet{mullen2014continuous} found ``\texttt{DEoptim}'' (implemented in \texttt{R}) \citep{mullen2009deoptim} to be among the best. The function implements a differential evolution algorithm, an example of an evolutionary strategy developed by \citet{storn1997differential} (see \citet{mullen2009deoptim} for a detailed overview of the underlying algorithm).

		\begin{figure}[!ht]
			\begin{center}

             \subfloat[]{{\includegraphics[scale=0.6]{BIC_ONE_TUN.pdf}}}\\

             \subfloat[]{{\includegraphics[scale=0.6]{3D_plot_matlab_cropped.pdf}}}
             \caption{The BIC function evaluated at different tuning parameter values for the Weibull MPR model with the LASSO penalty for the lung cancer dataset analysed in Section \ref{Sec:Data}, (a) corresponds to one tuning parameter and (b) corresponds to two.\label{fig:BIC}}
			\end{center}
	 \end{figure}

{\subsection{Variable Selection Algorithm}
The variable selection algorithm described above is summarized in the following bullet points.
\begin{itemize}
\item {\it Initialization.} Set $\theta^{(0)} = \hat\theta_0$ where $\hat\theta_0$ is the vector of unpenalized estimates, i.e., those which minimize $\ell(\theta)$ (defined in (\ref{eq:NonPen_lik})).
\item {\it  Optimization.}
\begin{itemize}
\item {\it  Outer.} Minimize $\text{BIC}(\lambda)$ with respect to $\lambda$ using \texttt{DEoptim}, yielding $\lambda^*$ as defined in (\ref{eq:lambdastr}). Convergence occurs when $|\text{BIC}(\lambda^{(r+1)})-\text{BIC}(\lambda^{(r)})|$ is below a prespecified threshold. (Here, $\lambda^{(r)}$ is the best $\lambda$ value found at step $r$ of the \texttt{DEoptim} algorithm.)
\item {\it Inner.} For a given value of $\lambda$, maximize $\ell_\lambda(\theta)$ by iteratively re-solving the system of equations given in (\ref{eq:NR}) starting from the initial value, $\theta^{(0)}$; this yields $\hat\theta_\lambda$. Convergence occurs when $||\theta_\lambda^{(m+1)}-\theta_\lambda^{(m)}||_\infty$ is below a prespecified threshold. (Here $||y||_\infty = \max_j | y_j |$ is the infinity norm.)
\end{itemize}
\item {\it Output.} The estimates $\hat\theta_{\lambda^*}$ are returned from the above procedure, and the corresponding standard errors are calculated by evaluating (\ref{eq:SE_formula}) at $\hat\theta_{\lambda^*}$.
\end{itemize}
}

\section{Simulation Studies}
\label{sec:Simulations}
\subsection{Setup
	The performance of the proposed variable selection methods is evaluated through simulation studies. The failure time is simulated from a Weibull MPR model with
	$$\log(\tau_i) = x^{T}_{i}(-1.5, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.8, 0.5, 0.0, 0.0)^{T},$$
	$$\log(\gamma_i) = z^{T}_{i}(0.5, 0.4, 0.0, 0.0, 0.0, 0.4, -0.2, 0.0, 0.0, 0.0, 0.0)^{T},$$
	where $x_i = z_i = (1,x_{i1},\dots,x_{i10})^T$ is a vector of correlated variables generated from an AR(1) process with a correlation coefficient $\rho = 0.5$. Each variable is marginally standard normal and the correlation between any two consecutive variables $x_{ij}$ and $x_{ik}$ is given by $\rho^{|j-k|}$. The corresponding censored times were generated from an exponential distribution. This setup was chosen so as to yield realistic survival data, where the true model is sparse and correlation in the covariates exists. Three different sample sizes ($n$ = 100, 500 and 1000) and two censoring proportions ($p_{\text{cen}} = 25\%, 50\%$) are considered. For each scenario, we considered the LASSO, SCAD, and ALASSO penalties with both a single tuning parameter or two tuning parameters (i.e., one for each of the two regression components). Each simulation was repeated 200 times.
	\subsection{Simulation Results}
    The variable selection and estimation procedures described in Sections \ref{sec:Formulation} and \ref{sec:Estimation} are applied to the simulated data and the results are summarized and discussed here. A number of metrics are used to evaluate the performance of the variable selection procedures, namely: the average number of true zero coefficients \emph{correctly} set to zero (C), the average number of true non-zero coefficients \emph{incorrectly} set to zero (IC), and the probability of choosing the true model (PT); for the oracle model, C = 7 and IC = 0. As a measure of prediction accuracy, we also consider the mean squared error (MSE), given by {$\text{MSE}(\hat{\beta}) = (\hat{\beta} - \beta)^{T} V (\hat{\beta} - \beta)$ and $\text{MSE}(\hat{\alpha}) = (\hat{\alpha} - \alpha)^{T} V (\hat{\alpha} - \alpha)$, where $V$, the simulated sample covariance matrix of the covariates, is computed for each simulation replicate \citep{zhang2007adaptive, tibshirani1997lasso}}. These metrics, averaged over simulation replicates for the scenarios with 25\% censoring, are reported in Table \ref{tab:SimResults}. (The results for 50\% censoring, shown in the Appendix, are similar.)
	\par
	As the sample size increases, we see an improvement across all the four metrics, for both the shape and the scale parameters and across all penalties. However, it is evident that the LASSO penalty does not set enough covariates equal to zero (i.e., it selects an overly complex model) irrespective of whether there is one tuning parameter or two. SCAD performs better but over-selects in the shape component, $\alpha$, when there is only one tuning parameter; this is improved by having two separate tuning parameters. The best performance comes from the ALASSO penalty which, for the largest sample size, selects the true scale and shape covariates more than 90\% of the time. Interestingly, the ALASSO performs well even with a single tuning parameter (but it does improve a little with two tuning parameters). Computation times for each of the penalties are given in Table \ref{tab:Comptime} where we see that SCAD is considerably slower than LASSO and ALASSO, while the times for these latter two are comparable. Furthermore, the computation times for the cases with two tuning parameters are 2 - 3 times longer than those with one tuning parameter.
	\par
Besides variable selection, we also consider parameter inference in terms of estimation bias, accuracy of the estimated standard error (SEE) computed using the sandwich formula, (\ref{eq:SE_formula}), and { the empirical coverage probability (CP) of a nominal 95\% confidence interval}. The results for the ALASSO penalty (for the 25\% censoring level) are presented in Table \ref{tab:MSEcheck}. Results for the ALASSO penalty show that the SEE is accurate for moderate sample sizes, but may underestimate the standard error (SE) for smaller samples. In the samples $n = 100$ and $n = 500$, the smaller parameters are overshrunk, i.e., they are biased downwards. For this reason the CPs do not perform well. However, this is not the case for $n = 1000$. In the case of $n = 1000$, {the CPs are close to the nominal $95\%$ level} (although perhaps a little low for the shape coefficients). An improvement can be seen in the case of two tuning parameters. Similar tables for the other penalties (and 50\% censoring) can be found in the Appendix, and, to summarize these additional results: LASSO overshrinks all parameters and the standard errors are underestimated in all cases, whereas the results for SCAD are better, but not as good as the ALASSO; in particular, the coverage for SCAD confidence intervals is much lower than the nominal level even for $n=1000$. { As expected, when the censoring rate is increased from 25\% to 50\%, the variation (i.e., SE, SEE and MSE) of estimates is increased overall across all the three penalties.}
	\begin{table
     \caption{Simulation results: variable selection metrics averaged over 200 simulation replicates.\label{tab:SimResults}}
	 \begin{center}
     \begin{adjustbox}{width=1\textwidth}

		 \begin{threeparttable}
      {\footnotesize \begin{tabular}{cr|cccc|cccc|cccc}
					\hline
                 \multicolumn{2}{c}{}&\multicolumn{12}{c}{One Tuning Parameter}\\
					\hline
					{} & {} & \multicolumn{4}{c|}{LASSO} &\multicolumn{4}{c|}{SCAD} &\multicolumn{4}{c}{ALASSO}\\
					$p_{\text{cen}}$ = 25\% & $n$ & C(0) & IC(7)& PT & MSE &C(7) & IC(0)&PT& MSE &C(7) & IC(0)&PT & MSE\\
					\hline
					{}&100& 5.29 & 0.18 & 0.12 & 0.33
					&6.32 &  0.14 &  0.53 & 0.33
					& 6.17 & 0.05 & 0.38 & 0.19  \\
					Scale ($\beta$)&500 &5.52 & 0.00 & 0.18 & 0.09
					&6.96 & 0.00 & 0.98 & 0.03
					& 6.79 & 0.00 & 0.83 & 0.03\\
					{} & 1000 & 5.83 & 0.00 & 0.31 & 0.05
					&6.99 & 0.00  & 0.99 & 0.01
					& 6.93 & 0.00 &0.94 & 0.01 \\

					\hline
					{}&100& 3.40 &  0.06 & 0.00 & 0.05
					& 4.19 & 0.09 & 0.02 & 0.05
					& 6.25 & 0.16 & 0.40 &  0.04 \\
					Shape ($\alpha$)& 500 & 3.37 &  0.00 & 0.01 & 0.01
					& 5.86  &0.00 & 0.33 & 0.01
					&  6.81 & 0.00 & 0.85 & 0.00 \\
					{} & 1000 & 3.46 & 0.00 & 0.00 & 0.01
					&6.51 & 0.00 & 0.64 & 0.00
					& 6.93 & 0.00 & 0.95 & 0.00 \\
					\hline
					\multicolumn{2}{c}{}&\multicolumn{12}{c}{Two Tuning Parameters}\\

					\hline
					{} & {} & \multicolumn{4}{c|}{LASSO} &\multicolumn{4}{c|}{SCAD}&\multicolumn{4}{c}{ALASSO}\\
					$p_{\text{cen}}$ = 25\% & $n$ & C(7) & IC(0) & PT & MSE &C(7) & IC(0)&PT& MSE &C(7) & IC(0)&PT & MSE\\
					\hline
					{}&100& 4.61 & 0.10 & 0.06 & 0.28
					& 6.26 & 0.10 & 0.55 & 0.28
					& 6.47 & 0.14 & 0.52 & 0.23 \\
					Scale ($\beta$)&500 & 5.40 & 0.00 & 0.20 & 0.08
					& 6.95 & 0.00& 0.97 & 0.02
					& 6.84 & 0.00 &0.88 & 0.03 \\
					{} & 1000 & 5.32 & 0.00 & 0.20 & 0.04 & 6.89& 0.00& 0.96& 0.01 &6.92 & 0.00& 0.92 & 0.01\\
					\hline
					{}&100& 5.27 & 0.25 &0.10 & 0.05
					& 5.15 & 0.20 & 0.11 & 0.06
					&6.30 & 0.25 & 0.44 & 0.05\\
					Shape ($\alpha$)&500 & 5.30 & 0.00 & 0.17 & 0.01
					& 5.98 & 0.00 & 0.36 & 0.01
					&6.93 & 0.00 & 0.93 & 0.00\\
					{} & 1000 & 5.50 & 0.00& 0.26 & 0.01 & 6.51& 0.00 &0.63 & 0.00 & 6.94 &0.00 & 0.94 & 0.00\\
					\hline
\end{tabular}}
{\footnotesize\begin{tablenotes}
\item[]{ {C, average correct zeros; IC, average incorrect zeros; PT, the probability of choosing the true model; MSE, the average mean squared error.}}
\end{tablenotes}}
\end{threeparttable}
\end{adjustbox}
 \end{center}
\end{table}

	\begin{table
     \caption{Average computation time per simulation replicate (in minutes).\label{tab:Comptime}}
     \begin{center}
				{\footnotesize \begin{tabular}{cr|ccc|ccc}
					\hline
					\multicolumn{1}{c}{} & \multicolumn{1}{c}{} &\multicolumn{3}{c}{One Tuning Parameter} &\multicolumn{3}{c}{Two Tuning Parameters}\\
					\hline
					$p_{\text{cen}}$ = 25\% & $n$&LASSO&SCAD&ALASSO&LASSO&SCAD&ALASSO\\
					\hline
					{}&100& 1.5 & 6.1 & 1.9 & 4.6& 11.7 & 4.7\\
					{}&500 & 2.7 & 6.2 & 3.2 & 7.3& 19.9 & 9.0\\
					& 1000 & 4.3 & 11.5 & 5.8 & 13.1 & 37.5 & 17.1 \\
					\hline
     \end{tabular}}
   \end{center}
	\end{table}

	\begin{table
     \caption{ Further simulation results: estimates, standard errors, and confidence intervals.\label{tab:MSEcheck}}
	 \begin{center}
		\begin{adjustbox}{width=1\textwidth}
		\begin{threeparttable}
				{\footnotesize \begin{tabular}{cc|cccc|cccc|cccc}
					\hline
					\multicolumn{14}{c}{ALASSO}\\
					\hline
					\multicolumn{14}{c}{One Tuning Parameter}\\
					\hline
					$p_{\text{cen}}$ = 25\% & & \multicolumn{4}{c|}{$n$ = 100} & \multicolumn{4}{c|}{$n$ = 500} & \multicolumn{4}{c}{$n$ = 1000}\\
					{} & $\theta$ & $\hat{\theta}$ & SE & SEE & CP & $\hat{\theta}$ & SE & SEE & CP & $\hat{\theta}$ & SE & SEE & CP\\
					\hline
					$\beta_0$ & -1.50 & -1.49 & 0.21 & 0.21 & 0.96
					& -1.47 & 0.09 & 0.09 & 0.92
					& -1.48 & 0.06 & 0.06 & 0.95\\
					$\beta_1$ & -1.00 & -0.99 & 0.19 & 0.17 & 0.93
					& -0.98 & 0.07 & 0.07 & 0.93
					&  -0.99 & 0.05 & 0.05 & 0.95\\
					$\beta_7$ & -0.80& -0.76 & 0.19 & 0.15 & 0.87
					& -0.77 & 0.07 & 0.06 & 0.89
					& -0.79 & 0.05 & 0.04 & 0.94\\
					$\beta_8$ & 0.50 & 0.42 & 0.18 & 0.13 & 0.85
					& 0.47 & 0.06 & 0.06 & 0.88
					& 0.48 & 0.04 & 0.04& 0.91 \\
					\hline
					$\alpha_0$ & 0.50 & 0.54 & 0.10 & 0.09&  0.92
					&0.50 &  0.04& 0.04 & 0.95
					& 0.50& 0.03& 0.03& 0.98\\
					$\alpha_1$ & 0.40 & 0.38 & 0.06& 0.06& 0.93
					& 0.40 & 0.02 & 0.02 & 0.93
					& 0.40& 0.02& 0.01& 0.91\\
					$\alpha_5$ & 0.40 & 0.35 & 0.08 & 0.06& 0.77
					& 0.39& 0.03 & 0.02& 0.91
					& 0.39& 0.02& 0.02& 0.91\\
					$\alpha_6$ &-0.20 & -0.15 &0.09 &0.05&0.77
					& -0.18 & 0.03 & 0.02& 0.87
					& -0.19& 0.02& 0.02& 0.91 \\
					\hline
					\multicolumn{14}{c}{Two Tuning Parameters}\\
					\hline
					$p_{\text{cen}}$ = 25\% & & \multicolumn{4}{c|}{$n$ = 100} & \multicolumn{4}{c|}{$n$ = 500} & \multicolumn{4}{c}{$n$ = 1000}\\
					{} & $\theta$ & $\hat{\theta}$ & SE & SEE & CP & $\hat{\theta}$ & SE & SEE & CP & $\hat{\theta}$ & SE & SEE & CP\\
					\hline
					$\beta_0$ & -1.50 & -1.48 & 0.24 & 0.21 & 0.92
					& -1.49 & 0.10 & 0.09 & 0.93
					& -1.50 & 0.06 & 0.06 & 0.97\\
					$\beta_1$ & -1.00 & -0.97 & 0.22 & 0.16 & 0.86
					& -0.99 & 0.08 & 0.07 & 0.92
					&  -1.00 & 0.05 & 0.05 & 0.95\\
					$\beta_7$ & -0.80&  -0.74 & 0.21 & 0.15 & 0.86
					& -0.79 & 0.07 & 0.06 & 0.93
					& -0.79 & 0.05 & 0.04 & 0.92\\
					$\beta_8$ & 0.50 & 0.38 & 0.20 & 0.12 & 0.77
					& 0.47 & 0.07 & 0.06 & 0.87
					& 0.49 & 0.04 & 0.04 & 0.93 \\
					\hline
					$\alpha_0$ & 0.50 & 0.53 & 0.11 & 0.09&  0.89
					&0.50 &  0.05& 0.04 & 0.92
					& 0.50& 0.03& 0.03& 0.97\\
					$\alpha_1$ & 0.40 & 0.37 & 0.07& 0.06& 0.86
					& 0.39 & 0.02 & 0.02 & 0.95
					& 0.40& 0.01& 0.01& 0.95\\
					$\alpha_5$ & 0.40 & 0.35 & 0.10 & 0.06& 0.72
					& 0.39 & 0.03 & 0.02& 0.93
					& 0.40& 0.02& 0.02& 0.95\\
					$\alpha_6$ &-0.20 & -0.14 &0.10 &0.05&0.72
					& -0.19 & 0.03 & 0.02& 0.87
					& -0.19& 0.02& 0.02& 0.93 \\
					\hline
     \end{tabular}}
   {\footnotesize\begin{tablenotes}
	\item[]{ SE, standard deviation of estimates over 200 replications; SEE, average of estimated standard errors over 200 replications; CP, the empirical coverage probability of a nominal 95\% confidence interval.}
	\end{tablenotes}}
\end{threeparttable}
\end{adjustbox}
 \end{center}
\end{table}

 \section{Lung Cancer Data}
 \label{Sec:Data}
 To illustrate the penalized variable selection methods on real data, a lung cancer dataset is considered. This dataset was collected as part of a PhD thesis by \citet{wilkinson1995lungc} (see also \citet{burke2017multi}). This dataset contains all individuals, of all ages, diagnosed with lung cancer in Northern Ireland during the one year period 1 October 1991 to 30 September 1992. Only cases of primary lung cancer were included. The date of diagnosis was taken to be the time origin for an individual and the end point was the earlier of the occurrence of death or the study end date, which was on 30 May 1993. Individuals who were still alive on the study end date were taken to have censored survival times. Individuals who died from another cause or who dropped out of the study were also censored. The final dataset included 855 patients, of which there were 673 deaths and 182 censored times. Besides the survival time and the censoring indicator, a number of other variables were recorded for each of the patients enrolled in the study (reference categories are listed first): age group ($< 40-$, $50-$, $60-$, $70-$, $> 80$), sex (female, male), treatment group (palliative, surgery, chemotherapy, radiotherapy, chemotherapy and radiotherapy), WHO status (normal activity, light work, unable to work, $> 50\%$ walking, bed/chair bound), cancer cell type (squamous cell, small cell, adenocarcinoma, other), serum sodium level ($\ge 136\ mmol/l$, $< 136\ mmol/l$, missing), serum albumen level ($\ge 35\ g/l$, $< 35\ g/l$, missing), metastases (no, yes, unknown), smoking status (non-smoker, current smoker, ex-smoker, missing).

 \subsection{Adequacy of Weibull}

{
Before considering covariates and variable selection, we first carry out an initial check that a baseline Weibull distribution is appropriate for the lung cancer data. The cumulative hazard function for the Weibull model is given by $H(t) = \int_0^t h(u) du = \lambda t^\gamma$, and, hence, $\log H(t) = \log \lambda + \gamma \log t$. Therefore, given an estimate $\hat H(t)$, a plot of $\log \hat H(t)$ against $\log t$ should produce a straight line. This standard Weibull model check is shown in Figure \ref{fig:weibullfit}, and, despite a slight deficiency for very small survival times, it appears that the Weibull model is reasonable.
\begin{figure}[!ht]
   \begin{center}
     \includegraphics[scale=0.5]{Log_CumHaz_Plot.pdf}
     \caption{Weibull model check. Here $\hat H(t)$, along with the 95\% confidence intervals, come from the Kaplan-Meier estimator.\label{fig:weibullfit}}
   \end{center}
\end{figure}
}

 \subsection{Variable Selection Results}
     The variable selection results for the different penalties are summarized in Table \ref{tab:LC_VSsummary}. In line with the results of the simulation study, the LASSO penalty selects the most complex model and the ALASSO penalty selects the least complex. Both ALASSO penalties (one and two tuning parameter cases) are in agreement on the non-importance of sex and smoking status, and although age group is selected in the scale in the case with one tuning parameter, it is not significant. {Interestingly, the two tuning parameter ALASSO selects the same set of covariates as identified by \citet{burke2017multi} using a BIC stepwise procedure (albeit they additionally selected treatment in the shape)}. We also see that, in the two tuning parameters cases, the scale tuning parameter is smaller that of the one tuning parameter case, while the shape tuning parameter is larger. This suggests that the single penalty over-penalizes the scale coefficients and under-penalizes the shape; this is also evident from the scale and shape degrees of freedom. Interestingly, the one tuning parameter ALASSO converges in less than half the time of the two tuning parameter ALASSO, and achieves similar results. We expect this based on our simulation studies, and also expect the results of the two tuning parameter case to be marginally better (albeit it takes longer to converge).
     \par
	Table \ref{tab:LC_VSestA-lasso} displays the estimated coefficients for both ALASSO penalties along with the unpenalized coefficients (we focus the ALASSO due to its superior performance in our simulation studies, but similar tables for LASSO and SCAD can be found in the Appendix). From \citet{burke2017multi}, note that the scale coefficients characterize the overall scale of the hazard (a positive value indicates an increase relative to the reference category), while the shape coefficients characterize its time evolution (a positive value indicates a hazard which increases over time relative to the reference category). We clearly see the similarity of the coefficient values for both the one and two tuning parameter ALASSO penalties, and, furthermore, that the selected variables are broadly in line with those which are statistically significant in the unpenalized model. Focussing on the results of the two tuning parameter case we find that: all treatments (apart from chemotherapy) have a negative scale coefficient suggesting that treatment reduces hazard (relative to palliative care); however, worse WHO status, small cancer cell type, presence of metasteses, and reduced sodium and albumen levels increase the hazard; lastly, sex, age group, and smoking status have no significant effect on the hazard. Since no variable appears in the shape component (i.e., all shape coefficients are set to zero), the selected model is a proportional hazards model, and exponentiating the scale coefficients yields the hazard ratios, e.g., the surgery hazard ratio is $\exp(-0.98) = 0.375$ so that the risk of death is approximately 37.5\% that of a patient receiving palliative care.
     \begin{table}[!ht]
         \caption{Summary of penalized models {(lung cancer dataset)} \label{tab:LC_VSsummary}}
         \begin{center}
         \begin{threeparttable}
         {\footnotesize\begin{tabular}{l|ccc|ccc}
           \hline
           {} & \multicolumn{3}{c|}{One Tuning Parameter}& \multicolumn{3}{c}{Two Tuning Parameters} \\
           {} & LASSO & SCAD & ALASSO & LASSO & SCAD & ALASSO\\
           \hline
           Treatment & $\boldsymbol{\beta}$, {\color{gray} $\alpha$} & \bm{$\beta$}, \bm{$\alpha$} & \bm{$\beta$}, \bm{$\alpha$} &\bm{$\beta$}, {\color{gray}$\alpha$} &\bm{$\beta$}, {\color{gray}$\alpha$} & \bm{$\beta$} \\
           Age group & {\color{gray} $\alpha$}  & {\color{gray}$\alpha$} & {\color{gray} $\beta$} & {\color{gray} $\alpha$} & {\color{gray} $\alpha$} & - \\
           WHO status &\bm{$\beta$}, {\color{gray}$\alpha$} & \bm{$\beta$}, \bm{$\alpha$} & \bm{$\beta$} & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$}, {\color{gray}$\alpha$} & \bm{$\beta$} \\
           Sex & {\color{gray} $\alpha$} & - & - & - & - & -  \\
           Smoking status & {\color{gray} $\alpha$} & {\color{gray}$\alpha$} & - & \bm{$\beta$} & \bm{$\beta$} & - \\
           Cell type  & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$} & \bm{$\beta$} & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$} \\
           Metastases & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$} & \bm{$\beta$} & \bm{$\beta$} & \bm{$\beta$}\\
           Sodium & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$}, {\color{gray} $\alpha$}  & {\color{gray} $\beta$} & \bm{$\beta$} & \bm{$\beta$} & \bm{$\beta$}\\
           Albumen & \bm{$\beta$}, \bm{$\alpha$} &\bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$} & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$}, {\color{gray} $\alpha$} & \bm{$\beta$} \\[15pt]

           Tuning parameter(s) & 0.026 & 0.041 & 0.015 & 0.014 & 0.024 & 0.004 \\
					& & & & 0.080 &  0.074 & 0.045 \\
					Degrees of freedom & 32.5 & 27.1 & 15.5 &
				 25.6 & 24.6 & 15.2\\
					Scale degrees of freedom & 14.2 & 12.0 & 13.2 &
				 18.3 & 17.1 & 14.2 \\
					Shape degrees of freedom & 18.3 & 15.0 & 2.4 &
				 7.4 & 7.6 & 1.0\\
           Computation time (in minutes)& 27.0 & 41.9 & 28.8 & 58.3 & 100.6 & 86.9 \\
           \hline
         \end{tabular}}
           {\footnotesize\begin{tablenotes}
           \item[] $\beta$ = ``selected in scale'', $\alpha$ = ``selected in shape'', and those which are non-significant (at the $5\%$ level) are shown in gray.
           \end{tablenotes}}
           \end{threeparttable}
       \end{center}
	 \end{table}

		\begin{sidewaystable
         \caption{Coefficients estimates and standard errors for the ALASSO penalties {(lung cancer dataset)} \label{tab:LC_VSestA-lasso}}\centering
             \begin{threeparttable}
					\begin{adjustbox}{max width=\linewidth}
					{\footnotesize	\tabcolsep = 12.9pt	\begin{tabular}{ll|ccc|ccc}
						\hline
						{} & {} & \multicolumn{3}{c|}{Scale} & \multicolumn{3}{c}{Shape}\\
						\hline
						\multicolumn{2}{c|}{Covariate} & No Penalty & One Tuning &Two Tuning& No Penalty& One Tuning & Two Tuning\\
						Intercept & & \bf{-3.38 (0.66)} & \bf{-3.12 (0.17)} & \bf{-3.15 (0.17)} &  -0.16 (0.22) & 0.04 (0.03) & 0.04 (0.03)\\
						Treatment & surgery & \bf{-1.69 (0.83)}& \bf{-0.89 (0.21)} & \bf{-0.98 (0.22)} & 0.11 (0.21) & 0.00 (0.00) & 0.00 (0.00) \\
						& chemotherapy &  -0.33 (0.37) & 0.00  (0.00)  & 0.00 (0.00) & -0.03 (0.15) & 0.00 (0.00) & 0.00 (0.00) \\
						& radiotherapy & \bf{-0.85 (0.21)} & -0.16 (0.10) & \bf{-0.21 (0.10)} & \bf{0.22 (0.08)} & 0.00 (0.00) &  0.00 (0.00) \\
						& chemo. \& radio. & \bf{-3.83 (0.98)}& \bf{-2.30 (0.89)} & \bf{-0.63 (0.22)} & \bf{0.77 (0.20)} & \bf{0.51 (0.21)} & 0.00 (0.00)\\
						Age group & $50-$ &  \bf{-0.90 (0.43)} & 0.00 (0.00) & 0.00 (0.00) & \bf{0.39 (0.16)} & 0.00 (0.00) & 0.00 (0.00)\\
						& $60-$ & \bf{-0.94 (0.39)} &0.00 (0.00) & 0.00 (0.00) & \bf{0.40 (0.15)} & 0.00 (0.00) & 0.00 (0.00)\\
						& $70-$ & -0.77 (0.39) &  0.02 (0.08) & 0.00 (0.00)& \bf{0.31 (0.15)} & 0.00 (0.00) & 0.00 (0.00)\\
						& $> 80$ & -0.78 (0.42) & 0.00 (0.00)  & 0.00 (0.00) & 0.31 (0.17) & 0.00 (0.00) & 0.00 (0.00)\\
						WHO status & light work & -0.02 (0.45) & 0.00 (0.00) & 0.00 (0.00) & 0.02 (0.12) & 0.00 (0.00) & 0.00 (0.00) \\
						& unable to work & 0.84 (0.43) & \bf{0.41 (0.10)} & \bf{0.44 (0.10)} &-0.10 (0.13) & 0.00 (0.00) & 0.00 (0.00) \\
						& $> 50\%$ walking & \bf{1.31 (0.44)} & \bf{0.99 (0.11)} & \bf{0.97 (0.11)} & -0.13 (0.14) & 0.00 (0.00) & 0.00 (0.00)\\
						& bed/chair bound & \bf{1.78 (0.50)} & \bf{1.28 (0.28)} & \bf{1.54 (0.25)}&  -0.03 (0.20) & 0.00 (0.00) & 0.00 (0.00) \\
						Sex & male & 0.03 (0.14) & 0.00 (0.00) &0.00 (0.00) & -0.03 (0.05) & 0.00 (0.00) & 0.00 (0.00) \\
						Smoking status & current smoker & 0.10 (0.22) & 0.00 (0.00) & 0.00 (0.00) & 0.15 (0.08) & 0.00 (0.00) & 0.00 (0.00) \\
						& ex-smoker & -0.05 (0.23) & 0.00 (0.00) & 0.00 (0.00) & 0.17 (0.09) & 0.00 (0.00) & 0.00 (0.00) \\
						& missing &  0.29 (0.40) & 0.00 (0.00) & 0.00 (0.00) &0.00 (0.00)  & 0.00 (0.00) & 0.00 (0.00) \\
						Cell type & small cell &  \bf{0.83 (0.26)} &  \bf{0.31 (0.12)} & \bf{0.43 (0.13)} & -0.05 (0.10) & 0.00 (0.00) & 0.00 (0.00) \\
						& adenocarcinoma & 0.28 (0.28) & 0.00 (0.00) & 0.00 (0.00) & 0.03 (0.10) & 0.00 (0.00) & 0.00 (0.00) \\
						& other & 0.32 (0.20) & 0.00 (0.00) & 0.09 (0.09) & -0.04 (0.07) & 0.00 (0.00) & 0.00 (0.00)  \\
						Metastases & yes &  \bf{1.35 (0.28)} & \bf{0.89 (0.12)} & \bf{0.84 (0.12)}& \bf{-0.19 (0.08)} & 0.00 (0.00) & 0.00 (0.00) \\
						& unknown &   \bf{0.83 (0.30)} & \bf{0.53 (0.13)} & \bf{0.41 (0.13)} & -0.14 (0.09) & 0.00 (0.00) & 0.00 (0.00)\\
						Sodium level & $< 136\ mmol/l$ &  \bf{0.33 (0.14)} & 0.14 (0.08) & \bf{0.24 (0.08)}& -0.01 (0.05) & 0.00 (0.00) & 0.00 (0.00)\\
						& missing &  -0.77 (0.45) & 0.00 (0.00) & 0.00 (0.00)&  \bf{0.32 (0.16)} & 0.00 (0.00) & 0.00 (0.00) \\
						Albumen level & $< 35\ g/l$ &  \bf{0.65 (0.16)} & \bf{0.36 (0.09)} &  \bf{0.37 (0.09)} & -0.10 (0.06) & 0.00 (0.00) & 0.00 (0.00) \\
						& missing & \bf{0.59 (0.28)} & 0.00 (0.00) &  \bf{0.27 (0.14)}& 0.09 (0.15) & 0.00 (0.00) & 0.00 (0.00) \\
						\hline
					\end{tabular}}
       \end{adjustbox}
       {\footnotesize\begin{tablenotes}
         \item[] Bold indicates statistically significant at the $5\%$ level.
       \end{tablenotes}}
       \end{threeparttable}
\end{sidewaystable}
\section{Discussion}
\label{sec:conc}
The MPR approach results in flexible models which extend standard models, but the presence of multiple regression components means that variable selection is necessarily more challenging than in standard settings where there is only a single regression component. In this article, we have proposed a penalized variable selection procedure for the simultaneous selection of significant variables in the shape and scale parameters of a Weibull MPR model in the survival analysis setting. The performance of these methods was illustrated using simulation studies and a real data example. While we have considered the Weibull model example in this article, the proposed variable selection procedures can be applied easily to other MPR models.
\par
Given that we model different distributional parameters (a scale and a shape parameter), there is no reason to assume that variable selection can be achieved with a single penalty applied to both regression components; hence, we also investigated the need for a separate tuning parameter for each regression component. We have found that the ALASSO performs very favourably in terms of identifying the true subset of covariates and coverage of calculated confidence intervals. This is true even with a single tuning parameter, however the results are improved when there are two tuning parameters (albeit this is more computationally intensive). On the other hand, SCAD does not perform well in the MPR setting, selecting an overly complex model and with poor confidence interval coverage for shape parameters.

\section*{Acknowledgements}

The first author is funded by the the Irish Research Council. The second author was supported by the Basic Science Research Program
through the National Research Foundation of Korea (NRF) funded by the Ministry of Science \& ICT (No.~NRF-2017R1E1A1A03070747)

\bibliographystyle{apa}


\section{Introduction}
Camera calibration is the estimation of a camera's mapping function between a set of known world points and their measured image coordinates. The parameters that define this mapping are usually divided into two categories: intrinsic and extrinsic parameters.  The intrinsic parameters represent the internal characteristics of the image lens and sensors, while the extrinsic parameters model the relative position and orientation of the camera to the 3D world \cite{Zhang:2000}.

\vspace{0.2cm}

\begin{figure}[htb]
	\begin{center}

		\includegraphics[height=5cm, width=6cm]{figs3/entP.png}
	\end{center}

\vspace{-0.5cm}
	\caption{Illustration of the varying Entrance Pupil (EP). The entrance $c$ coincide with the optical axis $O_c$. The intersection of the incidence rays originated from world points $P^{(1)}, P^{(2)}, P^{(3)}$ intersect the optical axis at off-axis shift $c^{(1)}, c^{(2)}, c^{(3)}$ from the EP.
	}
	\label{fig:osznec}
\end{figure}

\begin{figure}[htb]
		\centering

		\centerline{		\includegraphics[height=7cm, width=6cm]{figs3/projection1.png} }

	\caption{ Image formation model for fisheye cameras under Non-Single Viewpoint (NSVP). $(P^{(1)}, P^{(2)})$ represent the world point projected onto the image plane at $(\tilde{p}^{(1)}, \tilde{p}^{(2)})$ in the presence of lens distortion. Assuming no lens distortion, the projected image points are $(p^{(1)}, p^{(2)})$. These point rays both form incidence angles $(\theta^{(1)}, \theta^{(2)})$ with the optical axis.  }
	\label{fig:fullpattfg}
\end{figure}

Most of the popular conventional imaging system assume that all incidence ray from a world point must pass through the optical axis in order to conform with a Single Viewpoint (SVP) cameras, under thin lens model \cite{Hall:1982}. However, the single viewpoint assumption is not always true since imaging system are made of various optical elements and/or mirrors. In figure \ref{fig:osznec}, we illustrate a Non-Single Viewpoint (NSVP) imaging using a fisheye lens. The locus of all the viewpoints in such system setup form a caustic depicted with thick curve line. The dotted lines represent the incidence ray originated from world points $(P^{(1)}, P^{(2)}, P^{(3)})$ with varying incidence angles from the optical axis. The incidence rays meet at off-axis shift $(c^{(1)}, c^{(2)}, c^{(3)})$ from the Entrance Pupil (EP) denoted as $c$ which is equivalent to the optical center $O_c$. Only the incidence angle that coincides with the  optical axis meets at the EP in the figure illustration. As the incidence angle increase away from the optical axis, the off-axis shift value increase as well.

For the majority of optical lenses with narrow field of view, the varying EP is small and thus can be ignored during the intrinsic camera calibration.  However, this is not the case for special imaging device such as catadioptric and fisheye lenses. A fisheye lens consists of several refractions and reflections in the imaging formation process as a result of large number of optical elements. Therefore, a single viewpoint assumption for fisheye camera calibration in \cite{Kannala06ageneric, fasogbon1:2018} implies an unrealistic constraint on the entrance pupil location, which impacts the parameter's sensitivity and correctness. There are literature that incorporate EP variation, such as \cite{Di:2004, Grossberg:2014, Swaminathan:2001, Kumar:2014} using thin lens and \cite{Kumar:2015} using thick lens for perspective cameras. Also, \cite{Gennery2006} using thin lens for fisheye cameras. Most of the thin lens model either have extreme complexity or more suitable for catadioptric imaging system only.

\vspace{0.5cm}

A new camera model based on varying entrance pupil is proposed for fisheye lens calibration. We introduce EP correction strategy to bring NSVP to a SVP imaging model. In the proposed model, the extrinsic part of the model incorporates the EP parameters that are used to correct the incidence ray from world points, under a non-linear optimization through bundle adjustment procedure. The experimental validation on calibration data shows slightly better reprojection error over two representative state of the art \cite{Kannala06ageneric, Gennery2006} that are under SVP and NSVP categories respectively.

\vspace{0.5cm}

In the next section, we present state of the art image formation model and necessary equations needed to understand the proposed method. The section starts with the equidistant projection formulation \cite{Kannala06ageneric}, followed by the optical lens modeling for fisheye lenses. The proposed image formation model and the suitable calibration procedure is detailed in section 3. In section 4, we perform experiments on calibration image data, and comparison with two representative state of art methods. In the last section, we summarize the paper and explain the future direction of this work.

\section{Image Formation Model}\label{sb:intrinsic}

In the figure \ref{fig:fullpattfg}, $i=1,\ldots, n$ represents the number of world points projected onto the lens.  $O(X,Y,Z), O_c(X_{c},Y_{c},Z_{c})$ represent the world and camera origins. $o(x,y), o_p(u,v)$ represent the image and pixel coordinates origin, $f$ is the camera's focal length and $r$ is the distance between a projected image point $p$ and origin $o$. In addition, $O_c$ correspond to the optical center and entrance pupil $c$ of the camera/lens, and $Z_{c}-$axis describe the camera's optical axis. In the figure, the incidence ray for every object point should pass through the optical center under SVP assumption for thin lens modeling.

\vspace{0.2cm}
\noindent
{\bf Extrinsic :} The camera point $P_c^{(i)}(X_c, Y_c, Z_c)$ is derived from world point $P^{(i)}(X, Y, Z)$ using equation (\ref{fig:p4}). The $3\times 3$ rotation matrix $\mathcal{R}$ and $3\times 1$ translational vector $\mathcal{T}$ are the extrinsic parameters that determine the external pose of a camera to the world coordinates \cite{Hall:1982, Zhang:2000}. Thanks to equation (\ref{fig:p4}), the incidence angle $\theta^{(i)}$ under SVP with the optical axis is defined as \resizebox{.35\hsize}{!}{$\theta^{(i)} = cos^{-1} \frac{Z_c^{(i)}}{\| \overrightarrow{ O_cP_{c}^{(i)}} \|}$}.
\begin{equation}
\begin{split}
\begin{bmatrix}
X_{c}\\
Y_{c}\\
Z_{c}
\end{bmatrix}^{(i)} = [\mathcal{R}|\mathcal{T}]  \begin{bmatrix}
X\\
Y\\
Z\\
1
\end{bmatrix}^{(i)}
\end{split}
\label{fig:p4}
\end{equation}

\vspace{0.2cm}
\noindent
{\bf Intrinsic $\&$ Lens distortion :}
In the presence of optical distortions and using equidistant projection \cite{Gennery2006, Kannala06ageneric, Hall:1982}, distorted image point $\tilde{p}^{(i)}(\tilde{x},\tilde{y})$ is expressed in equation (\ref{kansimp}), where angle $\varphi = tan^{-1} \frac{y}{x}$, radial distance $r(\theta) = \theta + k_1 \theta^3 + k_2 \theta^5  + k_3 \theta^7  + k_4 \theta^9$  \cite{Kannala06ageneric}, and $k_1, k_2, k_3, k_4$ are the lens distortion parameters.

\begin{equation}
\begin{split}
\tilde{x} = r(\theta). cos (\varphi),\\
\tilde{y} = r(\theta). sin (\varphi) \\
\end{split}
\label{kansimp}
\end{equation}

As a final step, image point $\tilde{p}^{(i)}(\tilde{x},\tilde{y})$ is transformed to its pixel coordinate location as : $\tilde{u} = f_x.\tilde{x} + sk.\tilde{y} + u_0$  and $\tilde{v} = f_y.\tilde{x} + v_0$ respectively. The parameter $sk$ is the skew factor, $(f_x, f_y)$ are the focal lengths in pixels, $(u_0$, $v_0)$ are the principal point in pixels \cite{Zhang:2000}.

\section{Proposed Model using Entrance Pupil}
Some fisheye calibration methods \cite{Kannala06ageneric, Scaramuzza4059340, fasogbon1:2018} put a lot of unnecessary pressure on the optical lens model to simplify the distortions introduced as a result of the varying EP. This tend to affect the accuracy of the calibration and proper estimation of the distortion parameters. Most especially, closer object to the camera tends to be the most impacted by the entrance pupil variation. This problem is investigated using a 3D reconstruction of some measured object in an industrial application. As a result, we are motivated in this work to separate the EP formation from the intrinsic modeling of the system. With this approach, EP is then modeled as part of the camera extrinsic parameters.

\vspace{0.2cm}
\noindent
{\bf Extrinsic :}
For the non-single viewpoint in Fig. \ref{fig:fullpattfg}, the incidence angle is calculated as $\theta^{(i)} = cos^{-1} \frac{Z_{c}^{(i)}}{\| \overrightarrow{c^{(i)}P_{c}^{(i)}} \|}$. The EP variation is then modeled as a function of $\theta$ by using higher degree of odd polynomial. This is a reasonable assumption when one assume there is no decentering in the lens system  \cite{Gennery2006, Brown71close-rangecamera}. Our approach ensure that the EP model is separated from image coordinate system during future optimizations, and thus is modeled as part of the extrinsic parameters. As a result, we do not need to worry about decentering effect.

\vspace{0.2cm}
The proposed solution is better illustrated in Figure \ref{fig:onecsols}. We move the off-axis $c^{(i)}$ to optical center $O_c$ in a series of non-linear optimization. This is equivalent to moving point $P^{(i)}$ along the $Z_c-$axis using parameter $E(\theta^{(i)})$ which is $E(\theta)=e_1*\theta^3    +   e_2*\theta^5     +    e_3*\theta^7  +    e_4*\theta^9 $, where $e_1$, $e_2$, $e_3$, $e_4$ are the entrance pupil parameters. Using the parameter $E(\theta)$, our extrinsic model is rewritten as in equation (\ref{err:rep4}).

\begin{equation}
\begin{split}
\begin{bmatrix}
X_c \\
Y_c \\
Z_c\\
\end{bmatrix}^{(i)} =
\begin{bmatrix}
\mathcal{R} & \mathcal{T} \end{bmatrix}\begin{bmatrix}
X \\
Y \\
Z + E(\theta)\\
1 \\
\end{bmatrix}^{(i)}
\end{split}
\label{err:rep4}
\end{equation}

\begin{figure}
	\begin{center}

		\includegraphics[height=5cm, width=6cm]{figs3/offaxis.png}
	\end{center}
	\caption{Our correction for the entrance pupil variation. The world points $(P^{(1)}, P^{(2)})$ are moved in a series of optimizations to an optimal target plane using parameters  $E(\theta^{(1)})$ and $E(\theta^{(2)})$.
	}
	\label{fig:onecsols}
\end{figure}

\vspace{0.2cm}

These complete the proposed extrinsic model to correct the varying EP. The entrance pupil variation most impact the Z-axis, so the proposed correction is significant, while the X and Y-axis corrections are left to be corrected as part of radial distortion in section \ref{sb:intrinsic}. For lenses with large distortions other than radial ones, various models that incorporates tilt, tangential distortions can be easily integrated to this model \cite{Brown71close-rangecamera, fasogbon1:2018}

\subsection{Calibration Procedure}
For the calibration procedure, we follow the popular planar based calibration method in \cite{Zhang:2000}. A planar calibration target with known feature points that is moved at different poses from the camera to be calibrated is used. Let us denote $i=1,\ldots, n$ number of feature points, and $j=1,\ldots, m$ number of poses. A world point and its corresponding distorted image pixel can then be represented as $P^{(i,j)}$ and  $\tilde{p}^{(i,j)}$ respectively.

The calibration procedure consist of three main steps: (i) feature extraction \cite{fasogbon2:2018}, (ii) initial camera parameters estimation \cite{Kannala06ageneric},	(iii) Non-linear refinement of initial camera parameters with the proposed model.

\vspace{0.2cm}
\noindent
{\bf Feature Extraction :} We exert the correspondence between points $P^{(i,j)}$ on the calibration target and the projected image  pixels $\tilde{p}^{(i,j)}(u, v)$ using the method proposed in \cite{fasogbon2:2018}.

\vspace{0.2cm}
\noindent
{\bf Initial camera parameters estimation :} We retrieve some initial camera parameters: intrinsic $(f_x, f_y, u_0, v_0, k_1, k_2)$ and extrinsic $(\mathcal{R}_{(j)} , \mathcal{T}_{(j)})$ for each pose $j$  using method proposed in traditional fisheye calibration \cite{Kannala06ageneric}.

\vspace{0.2cm}
\noindent
{\bf Non-linear optimization with EP:} The remaining parameters that have not been initialized from earlier step are set to zero i.e ($sk=0$, $k_3, k_4=0$, $e_1, \ldots, e_4=0$). Finally, all these parameters are refined in equation (\ref{err:rep}) by minimizing the reprojection error between the true image pixels $\tilde{p}^{(i,j)}$ and the reprojected one $\hat{p}^{(i,j)}$.

\begin{equation}
\begin{split}
min \sum\limits_{j=1}^m\sum\limits_{i=1}^n\|\tilde{p}^{(i,j)} - \hat{p}^{(i,j)}(\hat{\mathcal{K}},\hat{\mathcal{R}}_{(j)},\hat{\mathcal{T}}_{(j)},\hat{e}_1,\ldots,\hat{e}_4, \\ \hat{k}_1, \ldots, \hat{k}_4)\|^2\\
\end{split}
\label{err:rep}
\end{equation}
where $\hat{\mathcal{K}}:\{\hat{f}_x, \hat{f}_y, \hat{sk}, \hat{c}_x, \hat{c}_y\}$ are the estimated intrinsic parameters, $\{ \hat{k}_1,\ldots,\hat{k}_4 \}$ are the estimated radial distortion parameters, and $\{ \hat{\mathcal{R}}_{(j)},\hat{\mathcal{T}}_{(j)}, \hat{e}_1, \ldots, \hat{e}_4 \}$ are the estimated extrinsic parameters.
The optimization is done with bundle adjustment using Levenberg-Marquardt's method \cite{Triggs:1999:BAM:646271.685629}. Please note that the world point is corrected at each iteration to give the estimated $\hat{p}^{(i,j)}$.

\section{Experiments}
The proposed calibration method has been evaluated on both real and synthetic dataset. The experiments have been done with a calibration target in a working distance of between $100 ~to~ 150~mm$ from the camera to be calibrated. This working distance is required to be able to accurately measure the capability of the system in an extreme setup. The calibration target is made up of (11$\times$7) 77 features in total (see figure \ref{fig:oneds}). The acquired image exhibit large optical distortion as a result of the fisheye lens. The feature pixels are the Center of Gravity (CoG) of the projected circular patterns, and have been extracted using the method proposed in \cite{fasogbon2:2018}.

\begin{figure}
	\begin{center}

		\includegraphics[height=4cm, width=6cm]{figs3/test2.jpg}
	\end{center}
	\caption{An example image for pose $j$ of the calibration target. The extracted feature pixels are the center of gravity of the projected circular pattern.
	}
	\label{fig:oneds}
\end{figure}

The name of the fisheye lens camera device used for the calibration test is witheld due to intellectual property issue. However, it is a virtual reality device that is made up of 8 cameras with spatial resolution of 2048 $\times$ 2048. Each fisheye lens has a Field of View (FOV) of 195 $^{\circ}$ and focal length of 3.5mm. We have made several test on all of the 8 cameras but only the result of one principal camera is documented in this paper due to page limitation. The proposed calibration is evaluated with Root Mean Square (RMS) error between the true image pixels and reprojected ones using the estimated camera parameters, over 25 calibration images.

\vspace{0.2cm}
In the remaining part of this section, we present the results of the proposed calibration method and compare the re-projection error with the representative state of art algorithms: (i) Kannala et al.  \cite{Kannala06ageneric} that is based on SVP approach and (ii) Gennery \cite{Gennery2006} that is based on NSVP approach. We have chosen Kannala et al.  \cite{Kannala06ageneric} as the compared fisheye calibration method over Scaramuzza et al. \cite{Scaramuzza4059340} because it provides  better result under our current experimental setup \cite{fasogbon2:2018}. All the evaluated methods are based on thin lens model which is the main focus in this work.

\vspace{0.2cm}
\noindent
{\bf Analysis of real image test:} The result and comparison between the experimented state of art methods and our approach are shown in Table \ref{tab:dfgggadf}. Different imaging conditions have been abbreviated in the caption of the table. The estimate of the principal point for our method is not very far from the physical principal point ($\frac{(u_0, v_0)}{2}$), even with the incorporation of the entrance pupil variation.

By analysing the whole table, one can see that our calibration method provides the lowest mean reprojection error. The reprojection error ({\color{blue} 0.1956 }) is only slightly lower than Gennery method \cite{Gennery2006} but still validates the proposed method.

\begin{table}[htb]
	\begin{center}

		\resizebox{0.85\hsize}{!}{$	\begin{tabular}{|   l | c    || c  |  |c||}
			\hline
			& {\bf Kannala \cite{Kannala06ageneric}} &  {\bf Gennery \cite{Gennery2006}}  &  {\bf Our Method} \\
			\hline     \hline

			$(f_x)$	&	587.9521   & 597.3014 & 591.7301 \\
			$(f_y)$	&593.1520	&   594.0581  & 592.0340 \\  \hline

			$sk$	& 0.2314	&  0.1526  & 0.1978 \\ \hline

			$(u_0)$	& 1020.3841	 & 1013.7982  & 1013.2001 \\
			$(v_0)$	& 1028.0564	&  1024.2808   & 1025.0300 \\ \hline

			$e_1$	& -	& 0.0712  & 0.0851\\
			$e_2$	& -	& -0.1092  & -0.2577\\
			$e_3$	& -	& 0.2917  & 0.3016\\
			$ e_4$	&  - & 0.4861   & 0.5368 \\  \hline

			$k_1$	& 0.0152	&  -0.0096 & 0.0109\\
			$k_2$	& 0.0161	&  0.0028 & -0.0013\\
			$k_3$	& -0.0232	&  -0.0132 & 0.0008\\
			$ k_4$	&  0.0065 &   0.0027 & -0.0004 \\  \hline

			$error$		& 0.2913 &  0.2115 & 0.1956 \\
			$std$	& 0.0891	&  0.0792   & 0.0735\\ \hline

		\end{tabular} $}

\end{center}

\caption{The experimental results and comparism between the proposed calibration method (ours) and state of art methods  under SVP and NSVP.  ($error$- Mean error, $std$ - standard deviation error in pixels). }
\label{tab:dfgggadf}
\end{table}

\vspace{0.5cm}
\noindent
{\bf Analysis of synthetic test:} To demonstrate further the validity of the proposed model, we made synthetic experiment validation using exact setup under the real image test. We use synthetic generated feature points in order to eliminate any bias from circular feature detection under camera calibration. Table  \ref{tab:dfggxsdf} shows the standard deviation of the calibration parameters. In the table, we summarize the distortion parameters (radial and EP), so we only display their mean values. The smallest mean values between the compared methods for each camera parameter is illustrated in {\bf bold} font. From the outcome of the experiments, our method provides the lowest standard deviation of all the compared methods for the focal length and principal point estimation.

\begin{table}[htb]
	\begin{center}

		\resizebox{0.95\hsize}{!}{$	\begin{tabular}{|   l | c    || c  |  |c||}
			\hline
			& {\bf Kannala \cite{Kannala06ageneric}} &  {\bf Gennery \cite{Gennery2006}}  &  {\bf Our Method} \\
			\hline     \hline

			$(f_x, f_y)$	&	0.1145   &  0.9275 & {\bf 0.9081} \\    \hline

			$sk$	& {\bf  0.1982}	& 0.2140   & 0.2047 \\ \hline

			$(u_0, v_0)$	& 0.3023	 & 0.2312  & {\bf 0.2147} \\ 	 \hline

			$e_1,\ldots, e_4$	& -	& 0.0313  & {\bf 0.0305}\\   \hline

			$k_1,\ldots, k_4$	& {\bf 0.0298}	& 0.0311  & 0.0354\\   \hline

		\end{tabular} $}

\end{center}

\caption{Mean standard deviation of the camera parameters estimated using synthetic data feature points. The smallest mean values is illustrated in {\bf bold} font.}
\label{tab:dfggxsdf}
\end{table}

\section{Conclusion}
We propose a new camera model that integrates entrance pupil variation for fisheye camera calibration. We model the entrance pupil variation as part of extrinsic and thus separate it from the intrinsic that includes the conventional lens distortion models. The proposed method gives accurate entrance pupil and other camera parameters estimate when applied on both real and synthetic images.

By taking into account the Entrance Pupil variation, one ensure that the pixels located at the highly distorted extremities of the fisheye lens, and especially image objects that are very close to the lens are corrected. This improves the calibration result and it is expected to greatly improve future 3D reconstruction.

\bibliographystyle{IEEEbib}