# Variable Selection in Doubly Truncated Regression

Ming Zheng, Chanjuan Lin, and Wen Yu
Department of Statistics, School of Management, Fudan University

###### Abstract

Doubly truncated data arise in many areas such as astronomy, econometrics, and
medical studies. For the regression analysis with doubly truncated response
variables, the existence of double truncation may bring bias for estimation as
well as affect variable selection. We propose a simultaneous estimation and
variable selection procedure for the doubly truncated regression, allowing a
diverging number of regression parameters. To remove the bias introduced by
the double truncation, a Mann-Whitney-type loss function is used. The adaptive
LASSO penalty is then added into the loss function to achieve simultaneous
estimation and variable selection. An iterative algorithm is designed to
optimize the resulting objective function. We establish the consistency and
the asymptotic normality of the proposed estimator. The oracle property of the
proposed selection procedure is also obtained. Some simulation studies are
conducted to show the finite sample performance of the proposed approach. We
also apply the method to analyze a real astronomical data.

Keywords: Adaptive LASSO; Double truncation; Diverging number of parameters;
Least absolute deviation; Oracle property; Variable selection.

Accepted by SCIENTIA SINICA Mathematica (in Chinese)

##  1 Introduction

Truncated data arise in astronomy, econometrics, and survival analysis, etc.
In truncation, only those subjects fall with an interval can be observed along
with the interval. The subjects fall out of their respective intervals are not
known to exist and consequently, have no chance to be observed. When the
truncation interval is unbounded from above, the truncation is called left-
truncation; when the interval is unbounded from below, it is called right-
truncation. Double truncation occurs when the truncation interval is bounded
in both sides.

Truncation may bring systematic bias to statistical analysis if it is not
dealt with appropriately. For nonparametric approach of distribution
estimation, Turnbull (1976) developed a general algorithm for finding the
nonparametric maximum likelihood estimator for arbitrarily grouped, censored
or truncated data. Lyden-Bell (1971) obtained a similar estimator for singly
truncated data. The counting process techniques are then adopted by Wang et
al. (1986), Keiding and Gill (1990), Lai and Ying (1991a), etc. For regression
analysis with single truncation, one can refer to Bhattacharya et al. (1983)
for extended Mann-Whitey approach, Tsui et al. (1988) for iterative bias
adjustment, Tsai (1990) for Kendall’s tau correlation, Lai and Ying (1991b)
for rank-based estimation, and so on. Some more recent developments include
Greene (2012) for econometrics and Kim et al. (2013) and Liu et al. (2016) for
general biased sampling.

Double truncation is technically more challenging to deal with compared with
single truncation, so fewer results have been found in existing literature.
For distribution estimation, Efron and Petrosian (1999) and Shen (2010)
developed nonparametric maximum likelihood estimation. Moreira and Alvarev
(2012) proposed kernel-type density estimation approach. For two-sample
problem, Bilker and Wang (1996) and Shen (2013a) extended the Mann-Whitney
test. For regression ananlsis with doubly truncated responses, Moreia and
Alvarev (2016) proposed kernel type approach for low dimensional covariates.
Shen (2013b) considered estimating a class of semiparametric transformation
models. More recently, Ying et al. (2020) proposed an extended Mann-Whitney-
type loss function for parameter estimation in linear regression model with
fixed covariates dimension. A Mann-Whitney-type rank estimator is then defined
as its minimizer.

Besides the parameter estimation, double truncation may also bring difficulty
for variable selection in regression analysis. To the best of our knowledge,
no literature has studied variable selection procedure for linear regression
models with doubly truncated responses. Based on the extended Mann-Whitney-
type loss function proposed by Ying et al. (2020), we adopt the adaptive LASSO
penalty (Zou, 2006) to develop a simultaneous estimation and variable
selection procedure for the doubly truncated linear models. The proposed
procedure can not only deal with fixed-dimensional covariates, but also be
applicable for the diverging model size where the dimension of covariates goes
to infinity with the sample size. We show that the adaptive LASSO penalty
leads to the oracle properties of the selection procedure. Meanwhile, an
iterative algorithm is designed for minimizing the proposed objective
function. In each iteration, a least absolute deviation (LAD) optimization is
involved and can be solved efficiently through standard LAD algorithm. A
modified BIC approach is used for selecting the tuning parameter. We also
develop a random weighting approach for estimating the standard errors of the
non-zero regression parameters. Numerical studies are conducted to illustrate
the finite sample behaviors of the proposed approach.

The rest of the paper is organized as follows. The next section introduces
necessary notation and the estimation procedure of doubly truncated regression
models. In Section 3, we give out the main results of the variable selection
procedure, including the basic idea of the selection method, the algorithm,
the oracle properties of the selection procedure, and the approach to
determine the tuning parameter. In Section 4, the results of simulation
studies are presented to show the finite sample behavior of the proposed
approach and a real dataset is analyzed to illustrate the approach. Some
concluding remarks are given in Section 5. All the technical details are
summarized in the Appendix.

##  2 Notation, model specification, and estimation

First some notations are introduced. Let $\tilde{Y}$ be the response variable
and $\tilde{X}$ be a $p$-dimensional covariate vector. We consider the linear
location-shift model

$\displaystyle\tilde{Y}=\beta^{\top}\tilde{X}+\tilde{\varepsilon},$ (1)

where $\beta$ is the $p$-dimensional regression coefficient vector and
$\tilde{\varepsilon}$ is the error term independent of $\tilde{X}$, is assumed
to hold. The response $\tilde{Y}$ is subject to double truncation. Let
$\tilde{L}$ and $\tilde{R}$ be the left and right truncation variables,
respectively. $(\tilde{Y},\tilde{L},\tilde{R},\tilde{X})$ can be observed if
and only if $\tilde{Y}$ falls into the interval $(\tilde{L},\tilde{R})$. We
also assume that $\tilde{Y}$ is conditionally independent of
$(\tilde{L},\tilde{R})$ given $\tilde{X}$. This is a common assumption for
truncation data analysis. Under model (1), this assumption is equivalent to
assume that $\tilde{\varepsilon}$ is independent of $(\tilde{L},\tilde{R})$.
The distribution and density function of $\tilde{\varepsilon}$ is denoted by
$F$ and $f$, respectively.

Let $(\tilde{Y}_{i},\tilde{L}_{i},\tilde{R}_{i},\tilde{X}_{i})$,
$i=1,\ldots,\tilde{n}$, be independent and identically distributed (i.i.d.)
copies of $(\tilde{Y},\tilde{L},\tilde{R},\tilde{X})$. As we mentioned,
$\tilde{Y}_{i}$ is observed if and only if
$\tilde{L}_{i}<\tilde{Y}_{i}<\tilde{R}_{i}$, together with
$(\tilde{L}_{i},\tilde{R}_{i},\tilde{X}_{i})$. Denote the number of the
observed data points by $n$, that is,
$n=\\#\\{i:\tilde{L}_{i}<\tilde{Y}_{i}<\tilde{R}_{i}\\}$. We use
$(Y_{i},L_{i},R_{i},X_{i})$, $i=1,\ldots,n$ to denote the observed data
points. Moreover, for $i=1,\ldots,n$, let
$e_{i}(\beta)=Y_{i}-\beta^{\top}X_{i}$,
$L_{i}(\beta)=L_{i}-\beta^{\top}X_{i}$, and
$R_{i}(\beta)=R_{i}-\beta^{\top}X_{i}$.

To correct the bias brought by the double truncation, Ying et al. (2020)
proposed a modified Mann-Whitney-type rank estimation equation

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}I\\{L_{j}(\beta)<e_{i}(\beta)<R_{j}(\beta),L_{i}(\beta)<e_{j}(\beta)<R_{i}(\beta)\\}(X_{i}-X_{j})\mbox{sgn}\\{e_{i}(\beta)-e_{j}(\beta)\\}=0,$
(2)

where $I\\{\cdot\\}$ is the indicator function and $\mbox{sgn}\\{\cdot\\}$ is
the sign function. By introducing the indicators, the estimation function
takes the summation of all the so-called comparable pairs and becomes the
unbiased when $\beta$ takes the true value. Meanwhile, solving (2) is
equivalent to minimizing the following loss function

$\displaystyle
l_{n}(\beta)=\sum_{i=1}^{n}\sum_{j=1}^{n}\left|\left[(e_{i}(\beta)-e_{j}(\beta))\wedge(R_{j}-Y_{j})\wedge(Y_{i}-L_{i})\right]\vee(L_{j}-Y_{j})\vee(Y_{i}-R_{i})\right|,$
(3)

where $\wedge$ and $\vee$ are the minimum and maximum operator, respectively.
The minimizer of $l_{n}(\beta)$, denoted by $\hat{\beta}$, is the proposed
estimator of $\beta$ by Ying et al. (2020). They showed that the estimator is
consistent and asymptotically normal under some regularity conditions.

##  3 Main results of variable selection

###  3.1 The method of variable selection

We consider variable selection for model (1) based on the observed data
points. The loss function (3) motivates a penalized approach for variable
selection. Instead of restricting to fixed dimension of $X_{i}$, we allow the
dimension to go to infinity with $n$, that is, $p$ depends on $n$, denoted by
$p_{n}$. Then all the quantities that are functions of the covariates depend
on $n$. For the simplicity of notation, we shall suppress the subscript $n$
when there is no confusion. Meanwhile, $n/\tilde{n}$ is assumed to converge to
a positive constant less than 1 as $\tilde{n}$ goes to infinity.

Because of the $L_{1}$ form of the loss function (3), the LASSO-type penalties
are naturally considered. Specifically, the objective function with usual
LASSO penalty is given by

$\displaystyle l_{n}(\beta)+\lambda_{n}\sum_{j=1}^{p_{n}}|\beta_{j}|,$

where $\lambda_{n}$ is a (data-dependent) tuning parameter. There exist two
main concerns if we want to use the above LASSO penalty. The first one is
about the computation. One can see that $l_{n}(\beta)$ is not convex in
$\beta$, so the optimization is no longer convex. However, we note that
$l_{n}(\beta)$ can be represented as

$\displaystyle
l_{n}(\beta)=\sum_{i=1}^{n}\sum_{j=1}^{n}I\\{L_{j}(\beta)<e_{i}(\beta)<R_{j}(\beta),L_{i}(\beta)<e_{j}(\beta)<R_{i}(\beta)\\}|e_{i}(\beta)-e_{j}(\beta)|.$

The non-convexity comes from the indicator function, while
$|e_{i}(\beta)-e_{j}(\beta)|$ is convex in $\beta$. Thus, one can develop an
iterative algorithm where the parameter in the indicator function is kept
fixed in each iteration step and the optimization is implemented for the rest
part of the objective function. Some more details of the computation are
discussed in Section 3.2 and Appendix A.1.

The second one is that it is well known that the LASSO estimator does not have
the unbiasedness property and there are scenarios where the LASSO selection is
not consistent. To overcome these drawbacks, Zou (2006) proposed adaptive
LASSO penalty, in which proper weights are introduced for penalizing different
coefficients in the lasso penalty. The weights used by Zou (2006) are defined
based on root-$n$-consistent estimator of $\beta$. He showed that with a
proper choice of the tuning parameter $\lambda_{n}$, the adaptive LASSO
possesses the oracle properties when the dimension of $\beta$ is fixed.

Trying to keep the oracle property, we also consider the adaptive LASSO
approach for variable selection. Specifically, the objective function with
adaptive LASSO is given by

$\displaystyle l_{n}^{\mbox{\tiny
AL}}(\beta)=l_{n}(\beta)+\lambda_{n}\sum_{j=1}^{p_{n}}\hat{w}_{j}|\beta_{j}|,$
(4)

where $\hat{w}_{j}$ is the adaptive weight of $\beta_{j}$. Originally, Zou
(2006) proposed to use $\hat{w}_{j}=1/|\tilde{\beta}_{j}|^{\gamma}$, where
$\tilde{\beta}_{j}$ is a root-$n$-consistent estimator of $\beta_{j}$ and
$\gamma$ is a predetermined positive constant. In our case, apparently, even
when $p_{n}$ diverges with $n$, we can still get $\hat{\beta}$ by minimizing
$l_{n}(\beta)$ as long as $p_{n}$ does not grow too fast. Based on
$\hat{\beta}$, we define $\hat{w}_{j}=1/|\hat{\beta}_{j}|^{\gamma}$ in (4),
where $\hat{\beta}_{j}$ stands for the $j$-th component of $\hat{\beta}$. The
proposed adaptive LASSO estimator, conceptually, is defined to be the
minimizer of $l_{n}^{\mbox{\tiny AL}}(\beta)$.

###  3.2 The algorithm

As we mention in Section 3.1, $l_{n}(\beta)$ is not convex, neither is
$l_{n}^{\mbox{\tiny AL}}(\beta)$, so we turn to an iterative procedure to get
the proposed estimator. Specifically, we modify $l_{n}^{\mbox{\tiny
AL}}(\beta)$ to

$\displaystyle l_{n}^{\mbox{\tiny
AL}}(\beta,b)=\sum_{i=1}^{n}\sum_{j=1}^{n}I\\{L_{j}(b)<e_{i}(b)<R_{j}(b),L_{i}(b)<e_{j}(b)<R_{i}(b)\\}|e_{i}(\beta)-e_{j}(\beta)|+\lambda_{n}\sum_{j=1}^{p_{n}}\hat{w}_{j}|\beta_{j}|.$

Let $\hat{\beta}^{(0)}$ be an initial estimate; for example, one may use
$\hat{\beta}$ as the start. Then, in the $k$-th iteration step, let
$\hat{\beta}^{(k)}=\mbox{argmin}_{\beta}l_{n}^{\mbox{\tiny
AL}}(\beta,\hat{\beta}^{(k-1)})$ for $k\geqslant 1$. If $\hat{\beta}^{(k)}$
converges to a limit as the number of $k\to\infty$, then we use the limit as
the proposed estimator.

The main reason that we use the above iterative algorithm is because in the
$k$th iteration, the objective function can be written as

$\displaystyle l_{n}^{\mbox{\tiny
AL}}(\beta,\hat{\beta}^{(k-1)})=\sum_{i,j\in{\cal
S}_{k}}|e_{i}(\beta)-e_{j}(\beta)|+\lambda_{n}\sum_{j=1}^{p_{n}}\hat{w}_{j}|\beta_{j}|,$
(5)

where ${\cal
S}_{k}=\\{(i,j)|L_{j}(\hat{\beta}^{(k-1)})<e_{i}(\hat{\beta}^{(k-1)})<R_{j}(\hat{\beta}^{(k-1)}),L_{i}(\hat{\beta}^{(k-1)})<e_{j}(\hat{\beta}^{(k-1)})<R_{i}(\hat{\beta}^{(k-1)})\\}$.
Thus, the optimization problem becomes a linear programming problem which can
be solved efficiently. Moreover, we find that after some algebraic operations,
the optimization can be transformed into a least absolute deviation (LAD)
problem which can be solved efficiently through standard LAD procedure, such
as the quantreg package in R. More details on the transformation are provided
in the Appendix A.1.

###  3.3 Large sample properties

Next we establish the oracle property of the proposed adaptive lasso approach
under some regularity conditions. The oracle property requires the consistency
of $\hat{\beta}$. Although Ying et al. (2020) proved root-$n$-consistency of
$\hat{\beta}$ under suitable conditions, their result is only confined to
fixed dimension. When $p_{n}$ diverges with $n$, the consistency of
$\hat{\beta}$ needs further investigation. Let $\beta_{0}$ denote the true
value of $\beta$. The following proposition, proved in the Appendix A.2, gives
out the consistency of $\hat{\beta}$ with the diverging dimension.

######  Proposition 1.

Assume that conditions C1 to C6 listed in the Appendix A.2 hold. If $p_{n}\log
n/\sqrt{n}=o(1)$, then $\|\hat{\beta}-\beta_{0}\|=O_{p}((p_{n}/n)^{1/2})$.

We then show the existence of an adaptive lasso estimator and it is consistent
with the convergence rate of $O_{p}((p_{n}/n)^{1/2})$. The following theorem,
proved in the Appendix A.2, states the results.

######  Theorem 1.

Assume that the conditions C1 to C6 listed in the Appendix A.2 hold. If
$\sqrt{n}\lambda_{n}=o(1)$ and $p_{n}\log n/\sqrt{n}=o(1)$, then with
probability tending to 1 there exists a local minimizer of $l_{n}^{\mbox{\tiny
AL}}(\beta)$, denoted by $\hat{\beta}^{\mbox{\tiny AL}}$, satisfying that
$\|\hat{\beta}^{\mbox{\tiny AL}}-\beta_{0}\|=O_{p}((p_{n}/n)^{1/2})$.

To present the following oracle properties, we need some more notation. Let
$\xi_{ij}(\beta)=I\\{L_{j}(\beta)<e_{i}(\beta)<R_{j}(\beta),L_{i}(\beta)<e_{j}(\beta)<R_{i}(\beta)\\}(X_{i}-X_{j})\mbox{sgn}\\{e_{i}(\beta)-e_{j}(\beta)\\}$,
$h_{ij}(\beta)=\left|\left[(e_{i}(\beta)-e_{j}(\beta))\wedge(R_{j}-Y_{j})\wedge(Y_{i}-L_{i})\right]\vee(L_{j}-Y_{j})\vee(Y_{i}-R_{i})\right|$,
and $\bar{l}(\beta)=\mathsf{E}(h_{ij}(\beta))$, where the expectation is taken
with respect to the distribution of the observed data. We define two
$p_{n}\times p_{n}$ matrices
$V=\mathsf{E}(\xi_{ij}(\beta_{0})\xi_{ik}^{\top}(\beta_{0}))$ and
$A=\partial^{2}\bar{l}(\beta)/\partial\beta\partial\beta^{\top}|_{\beta=\beta_{0}}$.
We further partition $\beta_{0}$ into $\beta_{0,1}$ and $\beta_{0,2}$, where
$\beta_{0,1}$ stands for the nonzero component of $\beta_{0}$ and
$\beta_{0,2}$ for the zero component. The dimension of $\beta_{0,1}$ is
denoted by $p_{n,1}$. Let $\hat{\beta}^{\mbox{\tiny AL}}_{1}$ and
$\hat{\beta}^{\mbox{\tiny AL}}_{2}$ be the adaptive lasso estimator for
$\beta_{0,1}$ and $\beta_{0,2}$, respectively. Moreover, let $V_{11}$ denote
the $p_{n,1}\times p_{n,1}$ components of $V$ corresponding to $\beta_{0,1}$
and $A_{11}$ denote the $p_{n,1}\times p_{n,1}$ components of $A$.

######  Theorem 2.

Assume that the conditions C1 to C6 listed in the Appendix A.2 hold. If
$\sqrt{n}\lambda_{n}=o(1)$ and $p_{n}\log n/\sqrt{n}=o(1)$, then with
probability tending to 1 the adaptive lasso estimator
$\hat{\beta}^{\mbox{\tiny AL}}$ has the following properties: i)
$\hat{\beta}^{\mbox{\tiny AL}}_{2}=0$; ii) and for any nonzero $p_{n,1}\times
1$ constant vector $u_{n}$ with $\|u_{n}\|=1$,
$\sqrt{n}u_{n}^{\top}(\hat{\beta}^{\mbox{\tiny
AL}}_{1}-\beta_{0,1})/[2(u_{n}^{\top}A_{11}^{-1}V_{11}Z_{11}^{-1}u_{n})^{-1/2}]$
converges in distribution to $N(0,1)$.

We give out the proof of Theorem 2 in the Appendix A.2. By choosing a suitable
$u_{n}$, it can be seen from Theorem 2 that the estimator of each individual
component of $\beta_{0,1}$ is asymptotically normal. The results can be used
to make inference on individual coefficients. The estimation of the variance
is discussed in Section 3.4.

###  3.4 The selection of tuning parameter and variance estimation

Typical methods of selecting tuning parameter include information criteria
methods such as AIC, BIC, and data-driven procedures such as $K$-fold cross-
validation. Because of the computational issue, we consider information
criteria approaches here. It is well known that BIC usually leads to the
consistency of selection for fixed dimension scenarios. However, as we allow
the dimension of the parameter to diverge, the standard BIC does not guarantee
the selection consistency. Thus, we consider the following modified BIC:

$\displaystyle\mbox{BIC}(\lambda)=l_{n}(\hat{\beta}^{\mbox{\tiny
AL}}_{\lambda})+\frac{\log n}{n}C_{n}\mbox{df}_{\lambda},$

where $\hat{\beta}^{\mbox{\tiny AL}}_{\lambda}$ represents the adaptive lasso
estimate obtained at the tuning parameter with level $\lambda$, $C_{n}$ is a
sequence going to infinity with $(\log n/n)C_{n}\to 0$, and
$\mbox{df}_{\lambda}$ stands for the size of the non-zero component of the
estimate at level $\lambda$. We find that $C_{n}=\log\log p_{n}$ works well in
this paper. The tuning parameter we choose is given by

$\displaystyle\hat{\lambda}=\mbox{argmin}_{\lambda}\mbox{BIC}(\lambda).$

The proposed modified BIC has the consistency of selection. Let ${\cal
M}^{\ast}$ and ${\cal M}_{\hat{\lambda}}$ be the true model and the selected
model, respectively. We have the following result on selection consistency.

######  Theorem 3.

Assume that the conditions C1 to C6 listed in the Appendix A.2 hold. If
$\sqrt{n}\lambda_{n}=o(1)$, $\lambda_{n}n/p_{n}\to\infty$, $p_{n}\log
n/\sqrt{n}=o(1)$, and $\sqrt{n/\\{C_{n}p_{n}\log
n\\}}\times\liminf_{n\to\infty}\\{\min_{j\in{\cal
M}^{\ast}}|\beta_{0,j}|\\}\to\infty$, then $\mathsf{P}({\cal
M}_{\hat{\lambda}}={\cal M}^{\ast})\to 1$ as $n\to\infty$.

Finally, we also need to give out the variance estimation of the proposed
estimator for the nonzero component of $\beta_{0}$. Note that the asymptotic
variance-covariance matrix of $\hat{\beta}^{\mbox{\tiny AL}}_{1}$ involves the
density function of $\tilde{\varepsilon}$ and is complicated to estimate
directly. We turn to use a random weighting method to obtain the variance
estimate. To be specific, we defined the perturbed penalized estimator,
denoted by $\hat{\beta}^{\ast}$, to be the minimizer of the following
perturbed version of the objective function $l^{\mbox{\tiny AL}}_{n}(\beta)$:

$\displaystyle
l_{n}^{\ast}(\beta)=\sum_{i=1}^{n}\sum_{j=1}^{n}(W_{i}+W_{j})h_{ij}(\beta)+\hat{\lambda}\sum_{j=1}^{p_{n}}\hat{w}_{j}|\beta_{j}|,$

where $W_{i}$, $i=1,\ldots,n$ are i.i.d. non-negative and bounded random
variables with mean 0.5 and variance 1, independent of the observed data. Let
$\hat{\beta}^{\ast}_{1}$ be component of $\hat{\beta}^{\ast}$ corresponding to
$\hat{\beta}^{\mbox{\tiny AL}}_{1}$, i.e., the non-zero part of
$\hat{\beta}^{\mbox{\tiny AL}}$. It can be shown that under suitable
conditions, given the observed data, the conditional asymptotic distribution
of $\hat{\beta}^{\ast}_{1}$ is the same as the asymptotic distribution of
$\hat{\beta}^{\mbox{\tiny AL}}_{1}$. Thus, by generating repeatedly the
sequence of $W_{i}$, $i=1,\ldots,n$, we can obtain a large number of
replications of $\hat{\beta}^{\ast}_{1}$ and then approximate the variance of
$\hat{\beta}^{\mbox{\tiny AL}}_{1}$.

##  4 Numerical results

###  4.1 Simulation studies

We conduct some simulation studies to illustrate the finite sample performance
of the proposed approach. We set $p_{n}=[7\tilde{n}^{1/5}]$, where $\tilde{n}$
is the size of the full sample, as defined in Section 2. The size of the non-
zero components is set to be $p_{n,1}=[p_{n}/3]$. We consider two scenarios
with $\tilde{n}=300$ and $\tilde{n}=500$, respectively. In the first scenario,
$p_{n}=21$ and $p_{n,1}=7$. We set
$\beta_{1,0}=(3.12,2.20,-0.86,0.92,-2.49,1.95,-1.32)^{\top}$. For the
covariates $\tilde{X}=(\tilde{X}_{1},\ldots,\tilde{X}_{21})$, we independently
generate $\tilde{X}_{1}$ from Binomial(0.25) (the Binomial distribution with
success probability 0.25), $\tilde{X}_{2}$ from Binomial(0.8), $\tilde{X}_{8}$
from Uniform(0,2) (the uniform distribution from 0 to 2), $\tilde{X}_{9}$ from
Binomial(0.5), and $\tilde{X}_{10}$ from Uniform(-2,0). All the rest covariats
are generated from the standard normal distribution, with the correlation
between $\tilde{X}_{i}$ and $\tilde{X}_{j}$ being $0.3^{|i-j|}$. In the second
scenario, $p_{n}=24$ and $p_{n,1}=8$. We set
$\beta_{1,0}=(3.12,2.20,-0.86,0.92,-2.49,1.95,-1.32,-2.13)^{\top}$. We
generate $\tilde{X}_{1}$ from Binomial(0.25), $\tilde{X}_{2}$ from
Binomial(0.8), $\tilde{X}_{9}$ from Uniform(0,2), $\tilde{X}_{10}$ from
Binomial(0.5), $\tilde{X}_{11}$ from Uniform(-2,0), and the rest from the
standard normal distribution, with the correlation between $\tilde{X}_{i}$ and
$\tilde{X}_{j}$ being $0.3^{|i-j|}$. For the error distribution in the linear
model, we consider the standard normal ($N(0,1)$) distribution and the extreme
value distribution with location parameter 0 and scale parameter 1
($EV(0,1)$). The left truncation variable $\tilde{L}$ is generated from
Uniform($a$,$b^{\top}X$) and the right one $\tilde{R}=\tilde{L}+c$, where $a$,
$b$, and $c$ are constants to yield truncation percentage of 30% and 40%. The
percentages of left and right truncation are almost the same.

Under each scenario, 1000 replications are carried out. For variable
selection, we define several criteria to assess the performance. The first one
is the model error (ME), being defined as
ME$=(\tilde{\beta}-\beta_{0})^{\top}\mathsf{E}(X_{i}X_{i}^{\top})(\tilde{\beta}-\beta_{0})$,
where $\tilde{\beta}$ is an estimator of $\beta_{0}$ derived from a specific
approach. The median and median absolute deviation (MAD) of the 1000 MEs are
recorded to evaluate the prediction performance of different procedures. The
second one is the average number of total estimated zero coefficients (TN),
that is, the average number of zero estimates obtained in the 1000
replications, including correctly estimated zero number (CN) and incorrectly
estimated zero number (IN). For CN, the closer it is to $p_{n}-p_{n,1}$ the
better a procedure performs, while for IN, the closer it is to 0 the better
the procedure performs. The third one is the ratio of the correctly identified
models (RCM), that is, the proportion of the replications that can identify
the true model correctly among the 1000 replications. The closer it is to
100%, the better a procedure performs. We consider three procedures, including
the proposed selection procedure (Proposed), the variable selection procedure
with adaptive lasso penalty without taking truncation into consideration
(Naive), and the oracle procedure where the correct subset of covariates is
used to fit the model (Oracle). The results are summarized in the following
tables.

[Insert Table 1 here]

[Insert Table 2 here]

From the tables, we see that the proposed method outperforms the naive
approach significantly in terms of the ME. This is expectable since the naive
approach yields biased estimators for the regression coefficients. For the TN
and RCM, the proposed method also has better results than the naive approach,
while the extend of the improvement is not very large. The TN and RCM are
criteria for variable selection performance. Although selecting important
variables seems to be less affected by the double truncation, correcting the
bias induced by double truncation still helps in increasing the precision of
the variable selection procedure.

###  4.2 A real data example

We analyze an astronomical data from the 46420 quasar catalog produced by the
Sloan Digital Sky Survey (SDSS) team. In this data, the main purpose is to
predict the red shift of a quasar by using some related features. Thus, the
linear model is used and the dependent variable is the red shift of the
quasars. Compared with the complete SDSS catalog, the data set we use here
omits some technical columns and therefore contains 23 covariates. Some more
details about the data set can be found in
https://astrostatistics.psu.edu/datasets.

It might not be that straightforward to see why there exists double truncation
in the red shift data. This mainly comes from the nature of the observation
feature of the astronomical telescope. Because of the aperture of the
astronomical telescope, the red shift of a quasar one can observe is limited
from the above. For the current data, we consider the maximum value of the red
shift (5.4135) as the right truncation bound of the dependent variable. That
is, the red shift larger this value is right truncated. For the low red shift,
the measurement might be unreliable. The values should be omitted if the value
of the red shift is lower its estimated error. The estimated error is then
treated as the left truncation bound of the red shift observation. Therefore,
the red shift value of the quasar is considered to be doubly truncated. Among
the 23 covariates, some are the measurement errors of the brightness which are
determined by the SDSS team from knowledge of the observing conditions,
detector background, and other technical considerations. These covariates are
not included in the model. We use 13 main covariates (shown in Table 3) to fit
the linear model. Due to the computation capacity limitation, we do not use
the full data to fit the model. A random sample with sample size 929 is drawn
to illustrate the proposed method. The sampled data is split into the training
part and testing part according to the proportion of 7:3. We use the proposed
variable selection procedure on the training set to get a list of ”important”
covariates and then used the selected model to do prediction on the testing
set. The results of the variable selection are summarized in Table 4.

[Insert Table 3 here]

[Insert Table 4 here]

The proposed variable selection procedure gives out 5 ”important” covariates
with non-zero coefficients. We use the selected covariates to do the
prediction on the testing set. In Figure 1, we draw the line plot to compare
the predicted values and the true response values. We see that the predicted
values (shown by the dotted line) are quite close to the true values (shown by
the solid line). The selected covariates give out reasonable predictions.

[Insert Figure 1 here]

##  5 Concluding remarks

We combine the Mann-Whitney-type loss function and the adaptive LASSO penalty
to develop a simultaneous estimation and variable selection procedure for
double truncated linear models, which allows the number of regression
parameters to grow to infinity with the sample size. To overcome the non-
convexity of the objective function, an iterative algorithm is designed. In
each iteration, the objective function to be optimized can be transformed into
a LAD problem. The oracle property of the proposed variable selection
procedure is proved. We use the modified BIC to select the tuning parameter
and adopt random weighting approach for variance estimation. Simulation
studies show the reasonable finite sample performances of the proposed
approach.

One of the main contributions of the current approach is to allow the
dimension of the covariates to go to infinity with the sample size. However,
by using the penalty based method, the growing rate of the covariates
dimension can not be very large. For the ultra-high dimension situations, some
screening approach is needed first. This is one of the interested directions
that can be explored in future study.

##  Appendix A Appendix

###  A.1 More details on the algorithm

Here we give out some more details on the algorithm discussed in Section 3.2.
Define $\check{Y}=Y_{i}-Y_{j}$ and $\check{X}=X_{i}-X_{j}$ for $(i,j)\in{\cal
S}_{k}$. Let $m$ be an appropriate index of ${\cal S}_{k}$. Thus
$\sum_{i,j\in{\cal
S}_{k}}|e_{i}(\beta)-e_{j}(\beta)|=\sum_{m=1}^{N}|\check{Y}_{m}-\beta^{\top}\check{X}_{m}|$.
Further define the matrices

$\begin{matrix}\check{\mathbf{Y}}=\begin{pmatrix}Y_{i_{1}}-Y_{j_{1}}\\\
\vdots\\\ Y_{i_{N}}-Y_{j_{N}}\end{pmatrix},&\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\check{\mathbf{X}}(\beta)=\begin{pmatrix}\beta^{\top}(X_{i_{1}}-X_{j_{1}})\\\
\vdots\\\ \beta^{\top}(X_{i_{N}}-X_{j_{N}})\end{pmatrix}\end{matrix}$

and the augmented matrices

$\begin{matrix}\check{\mathbf{Y}}_{\mbox{\tiny
aug}}=\begin{pmatrix}\check{\mathbf{Y}}/n(n-1)\\\
\mathbf{0}_{p_{n}}\end{pmatrix},&\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\check{\mathbf{X}}_{\tiny{\mbox{aug}}}(\beta,\lambda_{n})=\begin{pmatrix}\check{\mathbf{X}}(\beta)/n(n-1)\\\
\mathbf{diag}(\lambda_{n}\hat{w})\end{pmatrix}\end{matrix},$

where $\mathbf{0}_{p}$ is a $p_{n}$-vector of zeros and
$\mathbf{diag}(\lambda_{n}\hat{w})$ is a $p_{n}\times p_{n}$ diagonal matrix
with the elements
$(\lambda_{n}\hat{w}_{1},\cdots,\lambda_{n}\hat{w}_{p_{n}})$. Based on these
notations, to minimize the objective function in (5) can be formulated to
minimize $\|\check{\mathbf{Y}}_{\mbox{\tiny
aug}}-\check{\mathbf{X}}_{\tiny{\mbox{aug}}}(\beta,\lambda_{n})\|_{1}$ with
respect to $\beta$. This is a linear progamming problem and can be solved
quite efficiently through standard unpenalized LAD procedure, such as quantreg
package in R. The specific estimation procedure can be summarized as follows.

(i) Get an initial estimate value, for instance, the estimate value obtained
by ignoring the double truncation and the regularization.

(ii) Iteratively solve the unpenalized optimization problem until convergence
to get $\hat{\beta}$.

(iii) Use $\hat{\beta}$ to construct the objective function
$l_{n}^{\mbox{\tiny AL}}(\beta)$ and the corresponding modified one
$l_{n}^{\mbox{\tiny AL}}(\beta,b)$.

(iv) Use $\hat{\beta}^{(0)}$ as the initial value. Iteratively solve the
optimization problem $\min_{\beta}\|\check{\mathbf{Y}}_{\mbox{\tiny
aug}}-\check{\mathbf{X}}_{\tiny{\mbox{aug}}}(\beta,\lambda_{n})\|_{1}$ by
using linear programming until convergence to get $\hat{\beta}^{\mbox{\tiny
AL}}$.

###  A.2 Proof of the large sample properties

To obtain the large sample properties presented in Section 3.3, some
regularity conditions are needed. Let $Z=(Y,R,L,X)$ be the population of the
observed data. Let $\cal N$ be an $o(1)$ neighborhood of $\beta_{0}$ and $S$
be the sample space of $Z$. We use the notation $h(Z_{i},Z_{j};\beta)$ to
represent $h_{ij}(\beta)$, where $Z_{i}=(Y_{i},R_{i},L_{i},X_{i})$,
$i=1,\ldots,n$. Define $\tau(z,\beta)=\mathsf{E}[h(z,Z;\beta)]$. We use
$\nabla_{k}$ to represent the $k$-th order derivative of a function of $\beta$
with respect to $\beta$. We assume the following conditions hold.

C1. The density function $f$ is bounded and has a bounded and continuous
derivative.

C2. The parameter space $\Theta$ is compact and the true parameter value
$\beta_{0}$ is an interior point of $\Theta$.

C3. Each component of the covariate vector has a bounded second moment.

C4. For all $z\in S$, $\nabla_{2}\tau(z;\beta)$ exist on $\mathcal{N}$ and
continuous at $\beta_{0}$.

C5. Let $I(z)=\sup_{\alpha\in
S_{p}}\alpha^{\top}\nabla_{2}\tau(z;\beta_{0})\alpha$ with $S_{p}=\\{\alpha\in
R^{p_{n}}:\|\alpha\|=1\\}$. $\mathsf{E}[I^{2}(Z)]<\infty$.

C6. There exist positive constants $C_{1},C_{2},C_{3}$, and $C_{4}$ such that
for all $n$,

$0<C_{1}<\lambda_{min}(A)\leqslant\lambda_{max}(A)<C_{2}<\infty,0<C_{3}<\lambda_{min}(V)\leqslant\lambda_{max}(V)<C_{4}<\infty,$
(6)

where $\lambda_{min}(\cdot)$ and $\lambda_{max}(\cdot)$ are the largest and
smallest eigenvalue of a matrix, respectively.

C1 and C2 are mild conditions. C2 is used to prove the consistency of
$\hat{\beta}$, which is crucial to guarantee the consistency of the proposed
estimator $\hat{\beta}^{\tiny\mbox{AL}}$. C3 is sufficient to meet the
condition 5 in Wang et al. (2019) so that some results on high dimensional
Hoeffding’s decomposition can used to obtain the quadratic approximation of
$l_{n}(\beta)-l_{n}(\beta_{0})$. C4 is for the Taylor’s expansion of
$\tau(z;\beta)$ at $\beta_{0}$. C5 is used to bound the empirical term in the
Hoeffding’s decomposition of $l_{n}(\beta)-l_{n}(\beta_{0})$. C6 is used to
confine the Hessian matrix of the quadratic approximation of the objective
function and ensure that the covariance matrix of the score function is
positively definite with uniformly bounded eigenvalues for all $n$. This
provides justification for the component-wise asymptotic normality of
$\hat{\beta}^{\tiny\mbox{AL}}$. Similar conditions can be found in some
related literature, e.g., Peng and Fan (2004), Cai et al. (2005), Cho and Qu
(2013), and Ni et al. (2016).

We then give out the proof of the proposition and the theorems.

Proof of Proposition 1: For a U-process of order 2 with the kernel function
$h\in{\cal H}$, where $\mathcal{H}=\\{h(\cdot;\beta):\beta\in\Theta\\}$ is a
measurable class of symmetric functions on $S^{2}\equiv S\otimes S$ with a
non-negative envelope. Denote a set as
$\mathbf{h}_{\omega}(\cdot;\beta)=\\{h(Z_{1}(\omega),Z_{2}(\omega);\beta),\cdots,h(Z_{n-1}(\omega),Z_{n}(\omega);\beta)\\}$
with length $n(n-1)$ and a class
$\mathcal{H}_{\omega}=\\{\mathbf{h}_{\omega}(\cdot;\beta):\beta\in\Theta\\}$.
According to the Example 1 in Wang et al. (2019) and Lemma 4.4 in Pollard
(1990), we have that the pseudo-dimension of $\mathcal{H}_{\omega}$ is at most
$100p_{n}$ for every $\omega$. If $p_{n}\log n/\sqrt{n}=o(1)$, it is
sufficient to guarantee the condition on the pseudo-dimension in Proposition 1
of Wang et al. (2019). Hence, we have the following quadratic approximation

$l_{n}(\beta)-l_{n}(\beta_{0})=\frac{1}{2}(\beta-\beta_{0})^{\top}A(\beta-\beta_{0})+(\beta-\beta_{0})^{\top}\hat{U}_{n}+o_{p}(\|\beta-\beta_{0}\|^{2})+o_{p}\left(\frac{\sqrt{p_{n}}\|\beta-\beta_{0}\|}{\sqrt{n}}\right),$

where $\hat{U}_{n}=2\sum_{i=1}^{n}\xi_{ij|i}/n$ with
$\xi_{ij|i}=\mathsf{E}[\xi_{ij}(\beta_{0})|Z_{i}]$. By constraining $p_{n}$,
the pseudo-dimension is at the rate of $o(\sqrt{n}/\log(n))$. Therefore, all
the conditions required listed in Theorem 3 in Wang et al. (2019) are
satisfied. Then the consistency of $\hat{\beta}$ can be obtained.∎

To prove Theorem 1, a lemma is needed.

######  Lemma A.1.

Assume that the conditions C1 to C6 hold. Then $\forall b\in R^{p_{n}}$ with
$\|b\|_{2}=1$, $\sqrt{n}(b^{\top}Vb)^{-1/2}b^{\top}\hat{U}_{n}/2$ converges in
distribution to $N(0,1)$.

Proof: It is not difficult to see that

$\frac{\sqrt{n}}{2}(b^{\top}Vb)^{-1/2}b^{\top}\hat{U}_{n}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(b^{\top}Vb)^{-1/2}b^{\top}\xi_{ij|i}=:I_{1}.$

Note that $\mathsf{E}(\xi_{ij|i})=E(\xi(Z_{i},Z_{j},\beta_{0}))=0$ and
$Var(I_{1})=1$. By C3 and C4, $b^{\top}Vb$ is bounded. By the Cauchy-Schwarz
inequality, we have that

$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\mathsf{E}\left|(b^{\top}Vb)^{-1/2}\cdot
b^{\top}\xi_{ij|i}\right|^{4}$ $\displaystyle=$ $\displaystyle
O(1)n^{-2}\sum_{i=1}^{n}\mathsf{E}\left|b^{\top}\xi_{ij|i}\right|^{4}\leqslant
O(1)n^{-2}\sum_{i=1}^{n}\|\xi_{ij|i}\|_{2}^{4}$ $\displaystyle=$
$\displaystyle O_{p}(n^{-2}\cdot n\cdot
p_{n}^{2})=O_{p}(\frac{p_{n}^{2}}{n})=o_{p}(1).$

Thus, the Lyapunov condition holds for $I_{1}$. By Lyapunov central limit
theorem, $I_{1}$ converges in distribution to $N(0,1)$.∎

Proof of Theorem 1: Let $\alpha_{n}=(p_{n}/n)^{1/2}$. It can be shown that for
any $\epsilon>0$ and any constant vector $u_{n}$ with $\|u_{n}\|=C$, there
exists a large enough $C$ such that

$\mathsf{P}\left\\{\inf_{\|u_{n}\|=C}l^{\mbox{\tiny
AL}}_{n}(\beta_{0}+\alpha_{n}u_{n})>l^{\mbox{\tiny
AL}}_{n}(\beta_{0})\right\\}\geqslant 1-\epsilon.$

This implies the existence of a local minimizer $\hat{\beta}^{\mbox{\tiny
AL}}$ such that $\|\hat{\beta}^{\mbox{\tiny
AL}}-\beta_{0}\|=O_{p}(\alpha_{n})$. Then we have that

$\displaystyle l^{\mbox{\tiny
AL}}_{n}(\beta_{0}+\alpha_{n}u_{n})-l^{\mbox{\tiny AL}}_{n}(\beta_{0})$
$\displaystyle=$ $\displaystyle
l_{n}(\beta_{0}+\alpha_{n}u_{n})-l_{n}(\beta_{0})+\sum_{j=1}^{p_{n}}\\{P_{\lambda_{n}j}(|\beta_{0j}+\alpha_{n}u_{nj}|)-P_{\lambda_{n}j}(|\beta_{0j}|)\\}$
$\displaystyle\geqslant$
$\displaystyle\\{l_{n}(\beta_{0}+\alpha_{n}u_{n})-l_{n}(\beta_{0})\\}+\sum_{j=1}^{p_{n,1}}\\{P_{\lambda_{n}j}(|\beta_{0j}+\alpha_{n}u_{nj}|)-P_{\lambda_{n}j}(|\beta_{0j}|)\\}$
$\displaystyle=:$ $\displaystyle U_{1}-U_{2}.$

We write that

$U_{1}=\alpha_{n}u_{n}^{\top}\hat{U}_{n}+\frac{1}{2}\alpha^{2}_{n}u_{n}^{\top}Au_{n}+(\|u_{n}\|^{2}+\|u_{n}\|)o_{p}(\alpha^{2})=:U_{11}+U_{12}+U_{13},$

where

$|U_{11}|=|\alpha_{n}u_{n}^{\top}\hat{U}_{n}|\leqslant\alpha_{n}\|u_{n}\|\|\hat{U}_{n}\|=\alpha_{n}\|u_{n}\|O_{p}\left(\sqrt{\frac{p_{n}}{n}}\right)=\|u_{n}\|O_{p}(\alpha_{n}^{2})$

and

$|U_{12}|=|\frac{1}{2}\alpha^{2}_{n}u_{n}^{\top}Au_{n}|\geqslant\frac{1}{2}\alpha^{2}_{n}\|u_{n}\|^{2}\lambda_{min}(A)>\alpha^{2}_{n}\|u_{n}\|^{2}C_{1}/2.$

Therefore, for large enough $\|u_{n}\|$, $|U_{12}|$ dominates $|U_{11}|$ and
$|U_{13}|$ in $U_{1}$. Meanwhile,

$\displaystyle|U_{2}|$ $\displaystyle=$
$\displaystyle\left|\sum_{j=1}^{p_{n,1}}\frac{\lambda_{n}}{|\hat{\beta}_{j}|}|\beta_{0j}+\alpha_{n}u_{nj}|-\frac{\lambda_{n}}{|\hat{\beta}_{j}|}|\beta_{0j}|\right|\leqslant\alpha_{n}\lambda_{n}\left|\sum_{j=1}^{p_{n,1}}\frac{|u_{nj}|}{|\hat{\beta}_{j}|}\right|$
$\displaystyle=$
$\displaystyle\alpha_{n}\lambda_{n}\sqrt{p_{n,1}}O_{p}(1)\|u_{n}\|\|u_{n}\|o_{p}(\alpha^{2}_{n}).$

For large enough $C=\|u_{n}\|$, $|U_{12}|$ dominates $|U_{2}|$. Since $A$ is
positive definite, then $U_{12}>0$. Consequently, $l^{\mbox{\tiny
AL}}_{n}(\beta_{0}+\alpha_{n}u_{n})-l^{\mbox{\tiny AL}}_{n}(\beta_{0})>0$ with
probability tending to 1 as $n\to\infty$ for large enough $C$.∎

To prove Theorem 2, another lemma is needed.

######  Lemma A.2.

Assume that the conditions C1 to C6 hold. If $\sqrt{n}\lambda_{n}=o(1)$,
$\lambda_{n}n/p_{n}\to\infty$, and $p_{n}\log n/\sqrt{n}=o(1)$, then with
probability tending to 1, for any given $\beta_{1}$ satisfying
$\|\beta_{1}-\beta_{0,1}\|=O(p_{n}^{1/2}n^{-1/2})$ and any constant $C$,
$l^{\mbox{\tiny
AL}}_{n}\\{(\beta_{1}^{\top},0^{\top})^{\top}\\}=\min_{\|\beta_{2}\|\leq
Cp_{n}^{1/2}n^{-1/2}}l^{\mbox{\tiny
AL}}_{n}\\{(\beta_{1}^{\top},\beta_{2}^{\top})^{\top}\\}$.

Proof: It is sufficient to show that with probability tending to 1, for any
$\beta_{1}$ satisfying $\|\beta_{1}-\beta_{0,1}\|=O(p_{n}^{1/2}n^{-1/2})$ and
$\|\beta_{2}\|\leqslant Cp_{n}^{1/2}n^{-1/2}$, $\partial l^{\mbox{\tiny
AL}}_{n}(\beta)/\partial\beta_{j}$, $\beta_{j}$ have the same signs for
$j=(p_{n,1}+1),\cdots,p_{n}$. We have that

$\displaystyle\frac{\partial l^{\mbox{\tiny
AL}}_{n}(\beta_{n})}{\partial\beta_{j}}$ $\displaystyle=$
$\displaystyle\hat{U}_{n}(\beta_{0})_{j}+\sum_{k=1}^{p_{n}}A_{jk}(\beta_{k}-\beta_{0k})+o_{p}\left(\|\beta_{n}-\beta_{0}\|+\sqrt{\frac{p_{n}}{n}}\right)+P^{\prime}_{\lambda_{n}j}(|\beta_{j}|)sign(\beta_{j})$
$\displaystyle=:$ $\displaystyle V_{1}+V_{2}+V_{3}+V_{4}.$

It can be shown that

$V_{1}=O_{p}(\frac{1}{\sqrt{n}})=o_{p}\left(\sqrt{\frac{p_{n}}{n}}\right),$
$|V_{2}|=\left|\sum_{k=1}^{p_{n}}A_{jk}(\beta_{k}-\beta_{0k})\right|\leqslant\|\beta-\beta_{0}\|\left(\sum_{k=1}^{p_{n}}A_{jk}^{2}\right)^{\frac{1}{2}}=O_{p}\left(\sqrt{\frac{p_{n}}{n}}\right)O(1)=O_{p}\left(\sqrt{\frac{p_{n}}{n}}\right),$

and $|V_{3}|=o_{p}(\sqrt{p_{n}/n})$. Thus,
$V_{1}+V2+V_{3}=O_{p}(\sqrt{p_{n}/n})$. Following this, we have that

$\displaystyle\frac{\partial l^{AL}_{n}(\beta_{n})}{\partial\beta_{j}}$
$\displaystyle=$ $\displaystyle
P^{\prime}_{\lambda_{n}j}(|\beta_{j}|)sign(\beta_{j})+O_{p}\left(\sqrt{\frac{p_{n}}{n}}\right)=\frac{\lambda_{n}}{|\hat{\beta}_{nj}|}sign(\beta_{j})+O_{p}\left(\sqrt{\frac{p_{n}}{n}}\right)$
$\displaystyle=$
$\displaystyle\frac{\lambda_{n}}{|\hat{\beta}_{nj}|}\left\\{sign(\beta_{j})+O_{p}\left(\frac{p_{n}}{n\lambda_{n}}\right)\right\\}.$

Since $\lambda_{n}n/p_{n}\to\infty$, $sign(\beta_{j})$ dominates
$O_{p}(p_{n}/(n\lambda_{n}))$. Therefore, $\partial l^{\mbox{\tiny
AL}}_{n}(\beta_{n})/\partial\beta_{j}$ and $\beta_{j}$ have the same signs for
$j=(p_{n,1}+1),\cdots,p_{n}$ with probability tending to 1 as $n\to\infty$.∎

Proof of Theorem 2: The part i) follows directly from Lemma A.2. To prove ii),
we first denote that
$B_{n}=\\{P^{\prime}_{\lambda_{n1}}(|\beta_{01}|)sign(\beta_{01}),\cdots,P^{\prime}_{\lambda_{np_{n,1}}}(|\beta_{0k}|)sign(\beta_{0p_{n,1}})\\}^{\top}$.
It is easy to see that $B_{n}$ converges to 0 as $n\to\infty$. Since
$\hat{\beta}^{\mbox{\tiny AL}}$ is the minimizer of $l^{\mbox{\tiny
AL}}_{n}(\beta)$, $\partial l^{\mbox{\tiny AL}}_{n}(\hat{\beta}^{\mbox{\tiny
AL}})/\partial\beta_{1}=0$. Combining $\hat{\beta}^{\mbox{\tiny AL}}_{2}=0$,
we get that

$\hat{U}_{1}+A_{11}(\hat{\beta}^{\mbox{\tiny
AL}}_{1}-\beta_{0,1})+o_{p}(\|\beta_{n}-\beta_{0}\|+\sqrt{\frac{p_{n}}{n}})+B_{n}=0,$
(7)

where $\hat{U}_{1}$ consists of the first $p_{n,1}$ components of
$\hat{U}_{n}$. Rearrange (7) to get

$\sqrt{n}(\hat{\beta}^{\mbox{\tiny
AL}}_{1}-\beta_{0,1})=-\sqrt{n}A_{11}^{-1}\hat{U}_{1}-\sqrt{n}A_{11}^{-1}B_{n}+o_{p}(\sqrt{p_{n}}).$

By the condition that $\sqrt{n}\lambda_{n}=o(1)$, we have that
$\sqrt{n}A_{11}^{-1}B_{n}\to 0$. By C1 and C3, it can be seen that
$\mathsf{E}|\nabla_{1j}\tau(Z_{1}^{\prime}\beta_{0})|^{2}<\infty$ for
$j=1,\cdots,p_{n}$. Thus, $\|A_{11}^{-1/2}\hat{U}_{1}\|=O_{p}(\sqrt{p_{n}/n})$
and

$\sqrt{n}(\hat{\beta}^{\mbox{\tiny
AL}}_{1}-\beta_{0,1})=\sqrt{n}A_{11}^{-1}\hat{U}_{1}(1+o_{p}(1)).$

Since for any nonzero $p_{n,1}\times 1$ constant vector $u_{n}$ with
$\|u_{n}\|=1$,
$\sqrt{n}u_{n}^{\top}A_{11}^{-1}\hat{U}_{1}/[2\left[u_{n}^{\top}A_{11}^{-1}V_{11}A_{11}^{-1}u_{n}\right]^{-1/2}]$
converges in distribution to $N(0,1)$. Thus, by Slutsky’s theorem, we conclude
that $\sqrt{n}u_{n}^{\top}(\hat{\beta}^{\mbox{\tiny
AL}}_{1}-\beta_{0,1})/[2\left[u_{n}^{\top}A_{11}^{-1}V_{11}A_{11}^{-1}u_{n}\right]^{-1/2}]$
converges in distribution to $N(0,1)$.∎

Proof of Theorem 3: To state the proof more clearly, some more notation are
needed. The underlying true model is denoted by
$\mathcal{M}^{\ast}=\\{j:\beta_{0j}\neq 0,j=1,\ldots,p_{n}\\}$. If we use
$\hat{\beta}_{\lambda}$ to represent the estimator at the tuning parameter
value $\lambda$, then the corresponding model is denoted by
$\mathcal{M}_{\lambda}=\\{j:\hat{\beta}_{\lambda,j}\neq
0,j=1,\ldots,p_{n}\\}$. For any model $\mathcal{M}=\\{j_{1},\ldots,j_{m}\\}$
with size $m$, the corresponding estimated coefficients are
$\hat{\beta}_{\mathcal{M}}=(\beta_{j_{1}},\ldots,\beta_{j_{m}})$. Note that
$\hat{\beta}_{\mathcal{M}}$ and $\hat{\beta}_{\lambda}$ are two different
terms. We need to distinguish them carefully. Next, we divide the $R^{p}$
space into the following three disjoint regions: i)
$\Omega_{\minus}=\\{\lambda>0:\mathcal{M}_{\lambda}\not\supset\mathcal{M}^{\ast}\\}$,
ii) $\Omega_{0}=\\{\lambda>0:\mathcal{M}_{\lambda}=\mathcal{M}^{\ast}\\}$, and
iii)
$\Omega_{+}=\\{\lambda>0:\mathcal{M}_{\lambda}\supset\mathcal{M}^{\ast},\mathcal{M}_{\lambda}\neq\mathcal{M}^{\ast}\\}$.
This division divides the range of $\lambda$ into 3 categories, corresponding
to the under-fitted model, the correctly fitted model and the over-fitted
model. For the under-fitted model, there exists $j\in\mathcal{M}^{\ast}$ and
$j\notin\mathcal{M}_{\lambda}$. The over-fitted model refers to the model such
that for all $j\in\mathcal{M}^{\ast}$, $j\in\mathcal{M}_{\lambda}$, but there
exists $k\in\mathcal{M}_{\lambda}$ and $k\notin\mathcal{M}^{\ast}$. We prove
the theorem by considering the two parts: the under-fitted model and the over-
fitted model.

For the under-fitted model part, when $\lambda\in\Omega_{\minus}$, the length
of the vector $\hat{\beta}^{\mbox{\tiny AL}}_{S}$ is inconsistent with
$\beta_{0}$. In order to make $\hat{\beta}^{\mbox{\tiny AL}}_{S}$ and
$\beta_{0}$ comparable, let us slightly modify the definition of
$\hat{\beta}_{S}$ here. In this part, we define for any model $S$,

$\hat{\beta}_{S}=\mbox{argmin}_{\beta\in R_{k_{n}}:\beta_{j}=0,\forall i\notin
S}l_{n}(\beta).$

By definition,

$\displaystyle\inf_{\lambda\in\Omega_{\minus}}\mbox{BIC}_{\lambda}$
$\displaystyle=$
$\displaystyle\inf_{\lambda\in\Omega_{\minus}}l_{n}(\hat{\beta}_{\lambda})+\frac{\log
n}{n}C_{n}\cdot\mbox{df}_{\lambda}\geqslant\inf_{\lambda\in\Omega_{\minus}}l_{n}(\hat{\beta}_{\lambda})$
$\displaystyle=$
$\displaystyle\inf_{\lambda\in\Omega_{\minus}}l_{n}(\beta_{0})+\frac{1}{2}(\hat{\beta}_{\lambda}^{\ast}-\beta_{0})^{\top}A(\hat{\beta}_{\lambda}^{\ast}-\beta_{0})\geqslant
l_{n}(\beta_{0})+\frac{1}{2}\lambda_{min}(A)\inf_{\lambda\in\Omega_{\minus}}\|\hat{\beta}_{\lambda}^{\ast}-\beta_{0}\|^{2}$
$\displaystyle\geqslant$ $\displaystyle
l_{n}(\beta_{0})+\frac{1}{2}\lambda_{min}(A)\inf_{S\not\supset\mathcal{M}^{\ast}}\|\hat{\beta}_{S}-\beta_{0}\|^{2}\geqslant
l_{n}(\beta_{0})+\frac{1}{2}\lambda_{min}(A)\min_{j\in\mathcal{M}_{*}}(\beta_{0,j}^{2}).$

By C6, we can obtain that

$\inf_{\lambda\in\Omega_{\minus}}\mbox{BIC}_{\lambda}-l_{n}(\beta_{0})\geqslant\frac{1}{2}C_{1}\min_{j\in\mathcal{M}^{\ast}}(\beta_{0,j}^{2}).$
(8)

By the condition $\sqrt{n/\\{C_{n}p_{n}\log
n\\}}\times\liminf_{n\to\infty}\\{\min_{j\in S_{T}}|\beta_{0,j}|\\}\to\infty$,
the term on the right side of (8) is $O(C_{n}p_{n}\log n/n)$, which is
$O(\log\log p_{n}\times p_{n}\log n/n)$ as we choose $C_{n}=\log\log p_{n}$
here. Similarly, we derive that

$\displaystyle\mbox{BIC}_{\hat{\lambda}}$ $\displaystyle=$ $\displaystyle
l_{n}(\hat{\beta})+\frac{\log n}{n}C_{n}\cdot\mbox{df}\leqslant
l_{n}(\hat{\beta})+\frac{p_{n}\log n\log\log p_{n}}{n}$ $\displaystyle=$
$\displaystyle
l_{n}(\beta_{0})+\frac{1}{2}(\hat{\beta}-\beta_{0})^{\top}A(\hat{\beta}-\beta_{0})+\frac{p_{n}\log
n\log\log p_{n}}{n}$ $\displaystyle\leqslant$ $\displaystyle
l_{n}(\beta_{0})+\frac{1}{2}\lambda_{max}(A)\|\hat{\beta}-\beta_{0}\|^{2}+\frac{p_{n}\log
n\log\log p_{n}}{n}$ $\displaystyle=$ $\displaystyle
l_{n}(\beta_{0})+O_{p}\left(\dfrac{p_{n}}{n}\right)+\frac{p_{n}\log n\log\log
p_{n}}{n}$ $\displaystyle=$ $\displaystyle
l_{n}(\beta_{0})+o_{p}\left(\frac{p_{n}\log n\log\log
p_{n}}{n}\right)+\frac{p_{n}\log n\log\log p_{n}}{n}$

and

$\mbox{BIC}_{\hat{\lambda}}-l_{n}(\beta_{0})=O_{p}\left(\frac{p_{n}\log
n\log\log p_{n}}{n}\right).$

Combining the two results, we can obtain that

$\mathsf{P}\left(\inf_{\lambda\in\Omega_{\minus}}\mbox{BIC}_{\lambda}>\mbox{BIC}_{\hat{\lambda}}\right)\to
1.$ (9)

For the over-fitted model part, for any $\lambda\in\Omega_{+}$,

$\displaystyle\inf_{\lambda\in\Omega_{+}}\left(\mbox{BIC}_{\lambda}-\mbox{BIC}_{\hat{\lambda}}\right)$
$\displaystyle=$
$\displaystyle\inf_{\lambda\in\Omega_{+}}l_{n}(\hat{\beta}_{\lambda})-l_{n}(\hat{\beta})+(\mbox{df}_{\lambda}-\mbox{df}_{\hat{\lambda}})\frac{\log
n}{n}C_{n}$ $\displaystyle=$
$\displaystyle\inf_{\lambda\in\Omega_{+}}\frac{1}{2}(\hat{\beta}_{\lambda}-\beta_{0})^{\top}A(\hat{\beta}_{\lambda}-\beta_{0})-\frac{1}{2}(\hat{\beta}-\beta_{0})^{\top}A(\hat{\beta}-\beta_{0})$
$\displaystyle+(\mbox{df}_{\lambda}-\mbox{df}_{\hat{\lambda}})\frac{\log
n}{n}C_{n}$ $\displaystyle\geqslant$
$\displaystyle-\frac{1}{2}\lambda_{max}(A)\|\hat{\beta}-\beta_{0}\|^{2}+(\mbox{df}_{\lambda}-\mbox{df}_{\hat{\lambda}})\frac{\log
n}{n}C_{n}$ $\displaystyle\geqslant$
$\displaystyle-\frac{1}{2}C_{2}\|\hat{\beta}-\beta_{0}\|^{2}+\frac{\log
n}{n}C_{n}.$

The last line holds since for any $\lambda\in\Omega_{+}$, it must satisfy
$\mbox{df}_{\lambda}-\mbox{df}_{\hat{\lambda}}\geqslant 1$. Furthermore, as
$\|\hat{\beta}^{\mbox{\tiny AL}}-\beta_{0}\|^{2}=O_{p}(p_{n}/n)$, so

$\inf_{\lambda\in\Omega_{+}}\left(\mbox{BIC}_{\lambda}-\mbox{BIC}_{\hat{\lambda}}\right)\geqslant(1+o_{p}(1))\frac{p_{n}\log
n\log\log p_{n}}{n}>0.$

This means that

$\mathsf{P}\left(\inf_{\lambda\in\Omega_{+}}\mbox{BIC}_{\lambda}>\mbox{BIC}_{\hat{\lambda}}\right)\to
1.$ (10)

Combining (9) and ( 10) we can conclude that:

$\mathsf{P}(\mathcal{M}_{\hat{\lambda}}=\mathcal{M}^{\ast})\to 1$

as $n\to\infty$.∎

Acknowledgment

The research of Ming Zheng was supported by the National Natural Science
Foundation of China Grants (11771095). The research of Wen Yu was supported by
the National Natural Science Foundation of China Grants (12071088).

## References

  1. [1]

Bhattacharya, P. K., Chernoff, H., and Yang, S. S. (1983), Nonparametric
estimation of the slope of a truncated regression, The Annals of Statistics,
11, 505-514.

  2. [2]

Bilker, W., and Wang, M.-C. (1996), Generalized Wilcoxon statistics in
semiparametric truncation models, Biometrics, 52, 10-20.

  3. [3]

Cai, J., Fan, J., Li, r., and Zhou, H. (2005). Variable selection for
multivariate failure time data. Biometrika, 92, 303?16.

  4. [4]

Cho, H. and Qu, A. (2013), Model selection for correlated data with diverging
number of parameters. Statistica Sinica, 23, 901?27.

  5. [5]

Efron, B., and Petrosian, V. (1999), Nonparametric methods for doubly
truncated data, Journal of the American Statistical Association, 94, 824-834.

  6. [6]

Greene, W. H. (2012), Econometric Analysis (7th Ed.), Prentice Hall, Upper
Saddle River, NJ.

  7. [7]

Keiding, N., and Gill, R. D. (1990), Random truncation models and Markov
processes, The Annals of Statistics, 18, 582-602.

  8. [8]

Kim, J.P., Lu, W., Sit, T. and Ying, Z. (2013), A unified approach to
semiparametric transformation models under generalized biased sampling
schemes, Journal of the American Statistical Association, 108, 217-227.

  9. [9]

Lai, T. L., and Ying, Z. (1991a), Estimating a distribution function with
truncated and censored data, The Annals of Statistics, 19, 417-442.

  10. [10]

Lai, T. L., and Ying, Z. (1991b), Rank regression methods for left-truncated
and right-censored data, The Annals of Statistics, 19, 531-556.

  11. [11]

Liu, H., Ning, J., Qin, J. and Shen, Y. (2016), Semiparametric maximum
likelihood inference for truncated or biased-sampling data”, Statistica
Sinica, 26, 1087-1115.

  12. [12]

Lynden-Bell, D. (1971), A method of allowing for known observational selection
in small samples applied to 3CR quasars, Monthly Notices of the Royal
Astronomical Society, 155, 95-118.

  13. [13]

Moreira, C. and Una-Alvarez, J. (2012), Kernel density estimation with doubly
truncated data, Electronic Journal of Statistics, 6, 501-521.

  14. [14]

Moreira, C. and Una-Alvarez, J. (2016), Nonparametric regression with doubly
truncated data, Computational Statistics and Data Analysis, 93, 294-307.

  15. [15]

Ni, A., Cai, J., and Zeng, D. (2016), Variable selection for case-cohort
studies with failure time outcome, Biometrika, 103, 547-562.

  16. [16]

Peng, H. and Fan, J. (2004). Nonconcave penalized likelihood with a diverging
number of parameters, The Annals of Statistics, 32, 928-961.

  17. [17]

Pollard, D. (1990), it Empirical Processes: Theory and Applications,
IMS,Hayward, CA.

  18. [18]

Shen, P.-S., (2010), Nonparametric analysis of doubly truncated data, Annals
of the Institute of Statistical Mathematics, 62, 835-853.

  19. [19]

Shen, P.-S. (2013a), A class of rank-based tests for doubly-truncated data,
TEST, 22, 83-102.

  20. [20]

Shen, P.-S. (2013b), Regression analysis of interval censored data and doubly
truncated data. Computational Statistics, 28, 581-596.

  21. [21]

Tsai, W.-Y. (1990), Testing the independence of truncation time and failure
time, Biometrika, 77, 167-177.

  22. [22]

Tsui, K.-L., Jewell, N. P., and Wu, C. F. J. (1988), A nonparametric approach
to the truncated regression problem, Journal of the American Statistical
Association, 83, 785-792.

  23. [23]

Turnbull, B. W. (1976), The empirical distribution function with arbitrarily
grouped, censored and truncated data, Journal of the Royal Statistical
Society, Ser. B, 38, 290–295.

  24. [24]

Wang, M.-C., Jewell, N. P., and Tsai, W.-Y. (1986), Asymptotic properties of
the product limit estimate under random truncation, The Annals of Statistics,
14, 1597-1605.

  25. [25]

Wang, Z., Lin, Y., Liu, W., and Shao, Q. (2019), U-processes with increasing
dimensional parameters, working paper.

  26. [26]

Ying, Z., Yu, W., Zhao, Z., and Zheng, M. (2020), Regression analysis of
doubly truncated data, Journal of the American Statistical Association, 115,
810-821.

  27. [27]

Zou, H. (2006), The adaptive lasso and its oracle properties. Journal of the
American Statistical Association, 101, 1418-1429.

Table 1: Simulation results for variable selection when the random error follows the standard normal distribution. $\tilde{n}$ | Truncation | Method | ME | TN | RCM
---|---|---|---|---|---
percentage | Median | MAD | CN | IN | (%)
300 | 0.3 | Proposed | 0.220 | 0.137 | 13.921 | 0.000 | 92.836
Naive | 0.487 | 0.242 | 13.908 | 0.000 | 91.364
Oracle | 0.109 | 0.076 | 14.000 | 0.000 | 100.000
0.4 | Proposed | 0.385 | 0.232 | 13.862 | 0.000 | 88.679
Naive | 1.048 | 0.423 | 13.821 | 0.000 | 84.608
Oracle | 0.158 | 0.112 | 14.000 | 0.000 | 100.000
500 | 0.3 | Proposed | 0.133 | 0.073 | 15.995 | 0.000 | 99.503
Naive | 0.289 | 0.135 | 15.966 | 0.000 | 96.723
Oracle | 0.062 | 0.043 | 16.000 | 0.000 | 100.000
0.4 | Proposed | 0.217 | 0.120 | 15.977 | 0.000 | 98.014
Naive | 0.611 | 0.229 | 15.938 | 0.000 | 94.048
Oracle | 0.084 | 0.057 | 16.000 | 0.000 | 100.000
Table 2: Simulation results for variable selection when the random error follows the extreme minimum value distribution. $\tilde{n}$ | Truncation | Method | ME | TN | RCM
---|---|---|---|---|---
percentage | Median | MAD | CN | IN | (%)
300 | 0.3 | Proposed | 0.226 | 0.162 | 13.85 | 0.00 | 86.905
Naive | 0.780 | 0.353 | 13.82 | 0.00 | 83.730
Oracle | 0.129 | 0.098 | 14.00 | 0.00 | 100.000
0.4 | Proposed | 0.380 | 0.293 | 13.81 | 0.00 | 83.765
Naive | 1.232 | 0.565 | 13.72 | 0.00 | 75.896
Oracle | 0.184 | 0.132 | 14.00 | 0.00 | 100.000
500 | 0.3 | Proposed | 0.118 | 0.087 | 15.98 | 0.00 | 98.214
Naive | 0.489 | 0.208 | 15.96 | 0.00 | 95.933
Oracle | 0.075 | 0.048 | 16.00 | 0.00 | 100.000
0.4 | Proposed | 0.197 | 0.146 | 15.90 | 0.00 | 91.144
Naive | 0.989 | 0.369 | 15.88 | 0.00 | 89.154
Oracle | 0.107 | 0.072 | 16.00 | 0.00 | 100.000
Table 3: Description of the 13 covariates to fit the linear model for the SDSS data Name | Description
---|---
R.A. | Right Ascension
Dec. | Declination
u_mag | Brightness in the u (ultraviolet) band in magnitudes.
g_mag | Brightness in the g (green) band
r_mag | Brightness in the r (red) band
i_mag | Brightness in the i (more red) band
z_mag | Brightness in the z (even more red) band
Radio | Brightness in the radio band
X-ray | Brightness in the X-ray band
J_mag | Brightness in the near-infrared J band
H_mag | Brightness in the near-infrared H band
K_mag | Brightness in the near-infrared K band
M_i | The absolute magnitude in the i band.

Table 4: Estamtion and variabel selection results for the SDSS data

Covariate |  EST |  SEE
---|---|---
1 | R.A. | 0 | -
2 | Dec. | 0 | -
3 | u_mag | 0.340 | 0.024
4 | g_mag | -0.320 | 0.046
5 | r_mag | 0 | -
6 | i_mag | 0.500 | 0.043
7 | z_mag | 0 | -
8 | Radio | 0 | -
9 | X.ray | 0 | -
10 | J_mag | 0.017 | 0.007
11 | H_mag | 0 | -
12 | K_mag | 0 | -
13 | M_i | -0.803 | 0.017

EST: estimation of coefficients,
SEE: estimation of standard error

Figure 1: The predicted and original red shift values for the testing set


# Different intermediate water cluster with distinct nucleation dynamics among
mono layer ice nucleation

Yuheng Zhao  Yi Qin Gao<EMAIL_ADDRESS>Institute of Theoretical and
Computational Chemistry, College of Chemistry and Molecular Engineering,
Peking National Laboratory for Molecular Science, Peking University, Beijing
100871, China

###### Abstract

Recent first-principle calculations unveiled a distinctive dynamic behavior in
water molecule rotation during the melting process of highly confined water,
indicating a notable time-scale separation in diffusion. In this short paper,
we conducted molecular dynamics (MD) simulations to explore the rotation
dynamics during the mono-layer ice nucleation process to investigate the
possible intermediate states characterized by the differences in rotation of
water molecules. Our study reveals two types of ice clusters with similar ice
geometric structure but possess distinctly different rotational behaviors. In
terms of molecular rotation, one type cluster is ice like (ILC) and can be
regarded as small ice nuclei while the other is supercooled liquid water like
(SCC). We found distinct nucleation pathways, thermodynamic properties, and
phase transition dynamics to associate with these intermediate clusters, which
yielded an unexpectedly complex picture of mono-layer ice nucleation.

##  1 INTRODUCTION

Ice nucleation process plays a pivotal role in studies of numerous fields and
has been a subject of study for over a century. Gaining a comprehensive
understanding of this process can significantly advance our knowledge of the
formation and growth of cloud [1, 2, 3], the condensation on the
polyelectrolytes brush or silicates surface [4, 5, 6] and the self-assembly of
proteins process [7, 8, 9]. Both experimental [10, 11] and simulation-based
studies [12, 13, 14, 15] have contributed to uncover important properties
involved in ice nucleation [16, 17, 18], prompting the development of several
successful models for ice-water phase transitions [17, 19, 20, 21, 22, 23].
However, fundamental questions regarding ice nucleation still remain elusive
and in some cases controversial [24, 25].

As research continues to delve into the bulk water-ice phase transition, the
dimensionality effect, which significantly alters the thermodynamic properties
of water [26], has also attracted the interests of researchers [27, 28]. The
confined environment, in comparison to bulk water systems, leads to fewer
complex theoretical models [29] but results in a rich array of anomalous
phenomena [30, 31]. In the past decades the research field of confined liquid
water has observed many important efforts [32, 33, 34, 35, 36] which provided
significant new insights into the deep supercooled bulk water phase like ’no
man’s land’ [37], a term widely used in describing the phase transition
mechanism in the inaccessible regions of bulk water [38].

Mono layer water exhibits quasi-2D properties with novel features, such as
higher rotation but lower transition entropy compared to bulk water which
leads to a significantly slower diffusion on the plane [39]. In addition,
novel phases showing abnormal water dynamics have been identified. For
example, Kapil et al. observed a significant time-scale separation between
diffusion and rotation dynamics during melting. The rotation of water is
highly sensitive to increasing temperature, showing a continuous phase
transition in rotation speed [40]. In this study, we examine how rotation
dynamics contributes to the various properties of single water molecules as
well as water clusters.

Due to technological challenges associated with experimental approaches,
achieving in-situ observation of confined water phase transition behaviors
remains difficult. As an alternative, Molecular Dynamics (MD) simulations
offer a practical solution to reproduce confined water-ice phase transitions,
ensuring a reasonably realistic representation of dynamic behaviors with
atomic details. Molecule simulation dynamics (MD) method can be used to study
the local behavior of water and is also efficient enough to simulate water-ice
phase transition in confined environment. In our previous investigations, we
utilized the Molecular Dynamics (MD) simulations to explore the mono-layer
water-ice phase transition. Employing the Classical Nucleation Theory (CNT)
model, we successfully characterized the complete nucleation process [41]. We
observed that the cluster properties were predominantly dictated by their
size, and the comprehensive nucleation process could be satisfactorily
characterized by utilizing maximum cluster sizes as the order parameter. In
this study, we pay attention to dynamic details of 2D ice formation and in
particular, and examine how the intermediate clusters affect the pathways of
this phase transition.

##  2 SIMULATION DETAILS

In this study, we focus on the rotation dynamics in the process of mono layer
water nucleation. Specifically, we trace the evolution of rotational behavior
of individual water molecules within ice nuclei clusters. Keeping most
simulation parameters the same as in our previous study [41], in this study we
saved coordinates every 1000 steps, corresponding to a time interval of 1 ps,
to record the detailed rotational dynamics. In the interest of simulation
efficiency and to scrutinize temperature dependency, we initiated 10 parallel
simulations at each of the three temperatures: 275 K, 278 K, and 280 K. To
introduce rotation dynamics behavior as an order parameter for characterizing
nucleation progress, we compute the average angular speed ($\omega$) over 4 ps
intervals, which allows an effective discrimination between molecules in the
homogeneous ice and those in a liquid water state.

$\omega=\frac{\sum\theta}{t}$ (1)

##  3 RESULTS AND DISCUSSION

In ice formation, the evolution of $\omega$ follows a continuous change rather
than a first-order-like phase transition. We found that the structure-based
order parameter and $\omega$ align most closely when a cutoff of 0.9 is used
for $\omega$ to distinguish between solid-like water and liquid-like water
(SFigure.1). Additionally, relatively broad distributions of $\omega$ are
found to exist within individual clusters, indicating heterogeneous nature of
the rotational dynamics. It is thus a potentially useful parameter for
investigating detailed water properties within individual ice nuclei. Upon
scrutinizing trajectories obtained under various simulation temperatures, we
observed the consistent existence of two distinct types of ice nuclei
(Figure.1(b)/(c)) throughout the studied temperature range and even above the
melting point (SFigure.2). Both types of clusters can spontaneously form
during the ice cluster generation process. They both have the general square-
lattice-like structure of the typical mono-layer ice. Although they have
nearly the same sizes, the two types of clusters show significant differences
in rotation dynamics. In the first type of clusters, the majority of water
molecules exhibit water-like rotational properties, even they mostly have a
near perfect square-like ice structure. We call these clusters of fast
rotational dynamics ’Super Cooled-liquid cluster’ (SCC). In contrast, the
order-parameters of the second type of clusters occupy mainly the solid-like
region of the $\omega$ order parameter (see SFigure.3) and will be called
’Ice-like cluster’ (ILC).

Figure 1:  (a) The evolution of maximum cluster sizes and the rotation
dynamics order parameter $\omega$ within a typical trajectory with target
temperature of 278 K. Each dot represents a frame with the two order parameter
values. The color represents the time in simulation. We find that there are
two different clusters(I/II) which have similar sizes but different rotation
dynamics. (b) A snapshot for ’Ice-like cluster’(ILC). (c) A snapshot for
’Super Cooled-liquid cluster’(SCC).

We then proceed to investigate the nucleation processes for clusters
exhibiting different rotational properties, and found that clusters of
rotational differences also exhibit different growth dynamics. Two different
nucleation pathways were observed: (i) For ILC, upon formation, a significant
reduction in the angular speed of individual molecules occurs, coinciding with
its size growth. A complete transition to an ice-like cluster appears when the
cluster size reaches about 400, subsequently triggering a fast growth of the
ice cluster. (Figure.2a) (ii) For SCC, the cluster size continuously
increases, accompanied by a gradual reduction in the proportion of water-like
molecules, culminating in complete freezing. (Figure.2b). These two types of
clusters can both serve as seed clusters and eventually lead to the complete
freezing of the simulated system. In addition, we calculated the average
waiting time for seed clusters to undergo total freezing. For SCC, this
process takes several nanoseconds, while for ILC, it is less than 10
ps(SFigure.4). These findings thus reveal a pronounced difference in
nucleation between the clusters of different rotational dynamics and suggest
that the nucleation process is more appropriately described by the Ostwald’s
rule [42] than by CNT [43].

Figure 2: The two distinct trajectories with target temperatures of 275 K,
representing two different nucleation pathways initiated by (a) ILC, and (b)
SCC.

To understand the detailed phase transition mechanism, an effective sampling
of the two types of clusters is imperative. Due to the metastable nature of
both types of clusters, these simulations were performed at a temperature of
278 K that allows us to effectively model both ice nucleation and reverse
dissolution. In these simulations we successfully obtained trajectories that
encompass the progression of the generation, stabilization, and dissolution of
the two distinct types of clusters. As seen in Figure.3, a significant energy
difference exists between SCC and ILC. On average, the potential energy for
each water molecule in ILC is approximately 0.21 kcal/mol more stable than
SCC. This observation holds true across various parallel trajectories and
those obtained under different simulation temperatures. The observed lower
energy of the ILC compared to the SCC is consistent with the enhanced
stability associated with slower rotation dynamics in the former. However,
despite the great difference in the potential energy for the two types of
clusters, a small difference is observed between their populations at 278 K.
Their life times (See SI for method) in an overall 660 ns pre-freezing
trajectories are 112 and 94 ns, respectively. Such a result suggests that the
Gibbs free energies of ILC and SCC are similar at this temperature. This near-
degeneracy of the two rotationally different clusters suggests that CNT, with
the nuclei size as the order parameter, can be a reasonable model in
describing the nucleation process [43], although clusters of different
properties and the transition between them (see below) are involved in the ice
nucleation.

Figure 3: Potential energy surface obtained from one typical trajectory at 278
K. Although no bulk nucleation occurs in this trajectory, many transitions
occur between the liquid water, and both SCC and ILC. ILC is marked by I and
SCC, II.

Next, we investigated how the stability of the two types of clusters responds
to temperature. We selected typical structures of both SCC and ILC of similar
sizes at both 275 K ($\sim$ 450 water molecules) and 278 K ($\sim$ 250 water
molecules) and then performed trajectory shooting initiated from these
structures but at different temperatures to study their evolution as a
function of temperature. These shooting simulation results (Figure.4) reveal
that the stability of SCC decreases with increasing temperature, while ILCs
appear to have less temperature sensitivity and can stably exist for a long
period time once they are formed. This observation aligns with the finding
that ILC is of lower potential energy than SCC and is thus energetically more
stable. Besides thermodynamics stability, we also examined how temperature
affects the kinetic properties of these clusters. In particular, we examined
the temperature dependence of the “first-generation” cluster formation during
the simulations. Although the SCC are characterized by highly rotating water
molecules and are energetically less stable structure, they are kinetically
more easily formed, since they were seen to dominate first-generation clusters
especially in the high temperature simulations. At 275 K, the SCC: ILC ratio
is 7 to 3, and at the higher temperature 278 K all observed first-generation
clusters are SCC. This temperature dependence indicates a higher entropy of
SCC than ILC.

Figure 4:  The distribution of cluster size after shooting simulations
performed at (a) 278 K, and (b) 275 K. The red line marks the size of the
clusters at the beginning of the shooting trajectory. The orange and blue
histogram represents the size distribution for trajectory shooting initiated
from the SCC and ILC, respectively.

Next we examined the different nucleation pathways arising from the two
different seed clusters. Interestingly, the pathway through the SCC diminishes
at high temperatures while the continuously growing pathway appears to be less
sensitive to temperature. At high temperatures, although the SCC dominates the
first-generation cluster, they do not lead to efficient complete freezing. In
contrast, the nucleation pathway involving ILC remains effective in terms of
ice formation as temperature increases, positioning them as dominant seed
clusters at high temperatures. At 275 K, the seed cluster ratio for SCC: ILC
is 6:3, while at 278K, this ratio becomes 1:2. The observed disparity in
continuous growth capabilities in SCC and ILC is similar to the growth
capabilities for super cooled water and ice nuclei which again support our
conclusion that SCC is close to a supercooled cluster while ILC represents a
perfect ice-like nuclei.

Another intriguing discovery through these simulation studies is that the
melting process leads to a crossover between two independent nucleation
pathways, which involves a direct transition between SCC and ILC(Figure.5). In
our examination of cluster dynamics, we have observed that these transitions
are irreversible processes, exclusively leading SCC to transform into ILC, but
not the other way around. Moreover, these pathways overlap with the melting of
SCC, in that this transition between SCC and ILC always accompany with the
partial melting of the former. The results thus reveal the existence of two
pathways for SCC melting: i) direct dissolution into liquid water (SFigure.5);
ii) partial melting, leaving behind a cluster that retains a more ice-like
structure, the ILC. In addition, we found such a direct transition between the
SCC and ILC is relatively rare and only occurs at the higher temperature (278
K) to large SCCs (sizes large than 400). At low temperatures, the cluster
tends to continuously grow (Figure4(b)). Such a partial melting process has
also been found in 3D ice formation[44], representing a common relaxation
pathway for unstable cluster during the nucleation process.

Figure 5:  An example of SCC to ILC transformation observed from a simulation
at 278K.

The large difference between two different nucleation pathways prompts us to
also examine the reverse melting process for the bulk monolayer ice. We thus
performed melting simulations initiated from a fully frozen ice state, under a
target temperature set at 290 K to simulate the melting processes of a mono-
layer ice structure. It was found that the homogeneous ice which were grown
from ILC tend to pause for a short period of time at an intermediate state
during its melting, while the melting from ice formed through SCC nucleation
does not show such a behavior (SFigure.6). This observation of a clear
hysteresis is somehow beyond our expectation but suggests that the melting
dynamics of ice formed through water-ice transition simulation is dependent of
the freezing path. Not only the nucleation dynamics between SCC and ILC is
different, its final products also show slightly different properties and take
a long time to reach real equilibrium.

##  4 SUMMARY

In this study, we performed simulations on 2-D ice formation. From the
analysis of simulation trajectories, we find that the liquid water can
spontaneously assemble into clusters of ice-like structures. Clusters
exhibiting distinct rotation dynamics properties were found to form through
these spontaneous processes. The ILC, with their more pronounced ice-like
features, serve as seed clusters leading to a homogeneous frozen state better
characterized by a 1st order phase transition. In contrast, the SCC as seeds
result in gradual growth into bulk ice with a more continuous phase transition
behavior. Therefore, there exist two distinct nucleation pathways in 2-D ice
formation. Compared to the SCC, the ILC are energetically more stable but have
lower entropy . As a result, the SCC and ILC are nearly degenerate in terms of
Gibbs free energy when the temperature is slightly lower than the melting
point. The relative population of seed clusters of different rotation dynamics
in a single bulk nucleation event is thus temperature-dependent. At lower
temperatures, the one involving SCC represents the dominant nucleation
pathway. With an increase in temperature, this pathway becomes less prominent
due to the lower stability of the SCC, which is of a supercooled water nature.
Furthermore, a unidirectional transition pathway exists for the SCC to
transform to the ILC, and the probability of this transition sharply decreases
with the decrease of the size of SCC, which can melt before transition. The
origin of the ILC is thus not solely from spontaneous formation in liquid
water but they can also arise from the partial melting of larger (SCC)
clusters. The differences within the clusters introduce subtle difference
between ice formed from ILC and that from SCC through simulation: the former
but not the latter was found to involve an intermediate state in melting. This
observation suggests that although complete freezing can be observed in the
simulation through both pathways, the ice formed through simulation may take a
relatively long time to reach equilibrium and subtle difference exists in
these solids, although their geometry appear similar.

Finally we note that the highly unstable nature of supercooled clusters,
coupled with their less structured features, presents challenges to obtain
enough statistics on supercooled structures, even in extensive simulations.
Although the current study has provided more details to the 2D water to ice
phase transitions and associated thermodynamics properties, fundamental
questions such as the origin of supercooled cluster and the detailed dynamics
of hydrogen bonding[45] within intermediate cluster remain under-explored. We
hope that an improved sampling strategy will be able to further solve this
problem.

## References

  * Hoose et al. 2010 Hoose, C.; Kristjánsson, J. E.; Burrows, S. M. How important is biological ice nucleation in clouds on a global scale? _Environmental Research Letters_ 2010, _5_ , 024009
  * Murray et al. 2012 Murray, B.; O’sullivan, D.; Atkinson, J.; Webb, M. J. C. S. R. Ice nucleation by particles immersed in supercooled cloud droplets. 2012, _41_ , 6519–6554
  * Zhao et al. 2019 Zhao, B.; Wang, Y.; Gu, Y.; Liou, K.-N.; Jiang, J. H.; Fan, J.; Liu, X.; Huang, L.; Yung, Y. L. Ice nucleation by aerosols from anthropogenic pollution. _Nature Geoscience_ 2019, _12_ , 602–607
  * He et al. 2016 He, Z.; Xie, W. J.; Liu, Z.; Liu, G.; Wang, Z.; Gao, Y. Q.; Wang, J. Tuning ice nucleation with counterions on polyelectrolyte brush surfaces. _Science Advances_ 2016, _2_ , e1600345
  * Jin et al. 2020 Jin, Y.; Wu, C.; Yang, Y.; Wu, J.; He, Z.; Wang, J. Inhibiting Condensation Freezing on Patterned Polyelectrolyte Coatings. _ACS Nano_ 2020, _14_ , 5000–5007
  * Kumar et al. 2018 Kumar, A.; Marcolli, C.; Luo, B.; Peter, T. Ice nucleation activity of silicates and aluminosilicates in pure water and aqueous solutions – Part 1: The K-feldspar microcline. _Atmos. Chem. Phys._ 2018, _18_ , 7057–7079
  * Sear 2007 Sear, R. P. Nucleation: theory and applications to protein solutions and colloidal suspensions. _Journal of Physics: Condensed Matter_ 2007, _19_ , 033101
  * Linse et al. 2007 Linse, S.; Cabaleiro-Lago, C.; Xue, W.-F.; Lynch, I.; Lindman, S.; Thulin, E.; Radford, S. E.; Dawson, K. A. Nucleation of protein fibrillation by nanoparticles. _Proceedings of the National Academy of Sciences_ 2007, _104_ , 8691
  * Sleutel and Van Driessche 2018 Sleutel, M.; Van Driessche, A. E. S. Nucleation of protein crystals – a nanoscopic perspective. _Nanoscale_ 2018, _10_ , 12256–12267
  * Kringle et al. 2020 Kringle, L.; Thornley, W. A.; Kay, B. D.; Kimmel, G. A. Reversible structural transformations in supercooled liquid water from 135 to 245 K. _Science_ 2020, _369_ , 1490
  * Peng et al. 2018 Peng, J.; Guo, J.; Hapala, P.; Cao, D.; Ma, R.; Cheng, B.; Xu, L.; Ondráček, M.; Jelínek, P.; Wang, E.; Jiang, Y. Weakly perturbative imaging of interfacial water with submolecular resolution by atomic force microscopy. _Nature Communications_ 2018, _9_ , 122
  * Cai et al. 2019 Cai, X.; Xie, W. J.; Yang, Y.; Long, Z.; Zhang, J.; Qiao, Z.; Yang, L.; Gao, Y. Q. Structure of water confined between two parallel graphene plates. _The Journal of Chemical Physics_ 2019, _150_ , 124703
  * Tanaka 1996 Tanaka, H. Phase behaviors of supercooled water: Reconciling a critical point of amorphous ices with spinodal instability. _The Journal of Chemical Physics_ 1996, _105_ , 5099–5111
  * Nagata et al. 2015 Nagata, Y.; Usui, K.; Bonn, M. Molecular Mechanism of Water Evaporation. _Physical Review Letters_ 2015, _115_ , 236102
  * Sosso et al. 2016 Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.; Michaelides, A. Crystal Nucleation in Liquids: Open Questions and Future Challenges in Molecular Dynamics Simulations. _Chemical Reviews_ 2016, _116_ , 7078–7116
  * Cheng et al. 2018 Cheng, B.; Dellago, C.; Ceriotti, M. Theoretical prediction of the homogeneous ice nucleation rate: disentangling thermodynamics and kinetics. _Physical Chemistry Chemical Physics_ 2018, _20_ , 28732–28740
  * Li et al. 2011 Li, T.; Donadio, D.; Russo, G.; Galli, G. Homogeneous ice nucleation from supercooled water. _Physical Chemistry Chemical Physics_ 2011, _13_ , 19807–19813
  * Kaufmann et al. 2017 Kaufmann, L.; Marcolli, C.; Luo, B.; Peter, T. Refreeze experiments with water droplets containing different types of ice nuclei interpreted by classical nucleation theory. _Atmos. Chem. Phys._ 2017, _17_ , 3525–3552
  * Zobrist et al. 2008 Zobrist, B.; Marcolli, C.; Peter, T.; Koop, T. Heterogeneous Ice Nucleation in Aqueous Solutions: the Role of Water Activity. _The Journal of Physical Chemistry A_ 2008, _112_ , 3965–3975
  * Lupi et al. 2016 Lupi, L.; Peters, B.; Molinero, V. Pre-ordering of interfacial water in the pathway of heterogeneous ice nucleation does not lead to a two-step crystallization mechanism. _The Journal of Chemical Physics_ 2016, _145_ , 211910
  * Christenson 2013 Christenson, H. K. Two-step crystal nucleation via capillary condensation. _CrystEngComm_ 2013, _15_ , 2030–2039
  * Shi and Tanaka 2018 Shi, R.; Tanaka, H. Microscopic structural descriptor of liquid water. _The Journal of Chemical Physics_ 2018, _148_ , 124503
  * Tanaka 2010 Tanaka, H. Bond orientational ordering in a metastable supercooled liquid: a shadow of crystallization and liquid–liquid transition. _Journal of Statistical Mechanics: Theory and Experiment_ 2010, _2010_ , P12001
  * Bartels-Rausch et al. 2012 Bartels-Rausch, T. et al. Ice structures, patterns, and processes: A view across the icefields. _Reviews of Modern Physics_ 2012, _84_ , 885–944
  * Amann-Winkel et al. 2016 Amann-Winkel, K.; Böhmer, R.; Fujara, F.; Gainaru, C.; Geil, B.; Loerting, T. Colloquium: Water’s controversial glass transitions. _Reviews of Modern Physics_ 2016, _88_ , 011002
  * Chakraborty et al. 2017 Chakraborty, S.; Kumar, H.; Dasgupta, C.; Maiti, P. K. Confined water: structure, dynamics, and thermodynamics. _Accounts of chemical research_ 2017, _50_ , 2139–2146
  * Bellissent-Funel 2001 Bellissent-Funel, M. Structure of confined water. _Journal of Physics: Condensed Matter_ 2001, _13_ , 9165
  * Gilijamse et al. 2005 Gilijamse, J.; Lock, A.; Bakker, H. Dynamics of confined water molecules. _Proceedings of the National Academy of Sciences_ 2005, _102_ , 3202–3207
  * Dubrovskii 2014 Dubrovskii, V. G. _Nucleation theory and growth of nanostructures_ ; Springer, 2014
  * Fumagalli et al. 2018 Fumagalli, L.; Esfandiar, A.; Fabregas, R.; Hu, S.; Ares, P.; Janardanan, A.; Yang, Q.; Radha, B.; Taniguchi, T.; Watanabe, K., et al. Anomalously low dielectric constant of confined water. _Science_ 2018, _360_ , 1339–1342
  * Zangi and Mark 2004 Zangi, R.; Mark, A. E. Electrofreezing of confined water. _The Journal of chemical physics_ 2004, _120_ , 7123–7130
  * Gao et al. 2018 Gao, Z.; Giovambattista, N.; Sahin, O. Phase Diagram of Water Confined by Graphene. _Scientific Reports_ 2018, _8_ , 6228
  * Seeley and Seidler 2001 Seeley, L. H.; Seidler, G. T. Two-Dimensional Nucleation of Ice from Supercooled Water. _Physical Review Letters_ 2001, _87_ , 055702
  * Corsetti et al. 2016 Corsetti, F.; Matthews, P.; Artacho, E. Structural and configurational properties of nanoconfined monolayer ice from first principles. _Scientific Reports_ 2016, _6_ , 18651
  * Algara-Siller et al. 2015 Algara-Siller, G.; Lehtinen, O.; Wang, F. C.; Nair, R. R.; Kaiser, U.; Wu, H. A.; Geim, A. K.; Grigorieva, I. V. Square ice in graphene nanocapillaries. _Nature_ 2015, _519_ , 443–445
  * Nomura et al. 2017 Nomura, K.; Kaneko, T.; Bai, J.; Francisco, J. S.; Yasuoka, K.; Zeng, X. C. Evidence of low-density and high-density liquid phases and isochore end point for water confined to carbon nanotube. _Proceedings of the National Academy of Sciences_ 2017, _114_ , 4066
  * Mishima and Stanley 1998 Mishima, O.; Stanley, H. E. The relationship between liquid, supercooled and glassy water. _Nature_ 1998, _396_ , 329–335
  * Cerveny et al. 2016 Cerveny, S.; Mallamace, F.; Swenson, J.; Vogel, M.; Xu, L. Confined water as model of supercooled water. _Chemical reviews_ 2016, _116_ , 7608–7625
  * Cai et al. 2019 Cai, X.; Xie, W. J.; Yang, Y.; Long, Z.; Zhang, J.; Qiao, Z.; Yang, L.; Gao, Y. Q. Structure of water confined between two parallel graphene plates. _The Journal of Chemical Physics_ 2019, _150_ , 124703
  * Kapil et al. 2022 Kapil, V.; Schran, C.; Zen, A.; Chen, J.; Pickard, C. J.; Michaelides, A. The first-principles phase diagram of monolayer nanoconfined water. _Nature_ 2022, _609_ , 512–516
  * Qiao et al. 2019 Qiao, Z.; Zhao, Y.; Gao, Y. Q. Ice Nucleation of Confined Monolayer Water Conforms to Classical Nucleation Theory. _The Journal of Physical Chemistry Letters_ 2019, _10_ , 3115–3121
  * Van Santen 1984 Van Santen, R. The Ostwald step rule. _The Journal of Physical Chemistry_ 1984, _88_ , 5768–5769
  * Kalikmanov 2013 Kalikmanov, V. I. In _Nucleation Theory_ ; Kalikmanov, V. I., Ed.; Springer Netherlands: Dordrecht, 2013; pp 17–41
  * Li et al. 2023 Li, M.; Huang, Y.; Xia, Y.; Chen, D.; Yang, L.; Yang, Y. I.; Gao, Y. Q. Seek Homogeneous Critical Ice Nucleus at Moderate Supercooling. 2023
  * Yun et al. 2022 Yun, Y.; Khaliullin, R. Z.; Jung, Y. Correlated local fluctuations in the hydrogen bond network of liquid water. _Journal of the American Chemical Society_ 2022, _144_ , 13127–13136