Guillaume Gerber

Yohann Le Faou
email: yohann.lefaou@courriel.upmc.fr

Olivier Lopez
email: olivier.lopez@sorbonne-universite.fr

Michael Trupin

The impact of churn on client value in health insurance, evaluation using a random forest under various censoring mechanisms

The impact of churn on client value in health insurance, evaluation using a random forest under various censoring mechanisms

Introduction

In insurance brokerage, an important problem is to evaluate what is known as "prospect value".

Roughly speaking, this represents the rentability of a potential policyholder. More effort could then be dedicated to attract a prospect with probable high profitability. A first approach to evaluate this prospect value is to predict how long a given current customer will keep their contract.

To investigate this, brokers now have access to large databases of potentially relevant information. In the present paper, our aim is to discuss an extension of the random forest algorithm which takes into account specific details of this framework. Among these, a crucial issue is to deal with right-censoring, which is a common problem when dealing with duration variables.

Due to the temporal phenomenon that we study, a significant number of observations are incomplete and require particular care to counterbalance the corresponding lack of information. The method we propose consists of an appropriate weighting of the observations to compensate for censoring. We discuss different approaches and compare them using simulation and real data analysis.

The random forest algorithm was first proposed in [START_REF] Breiman | Random forests[END_REF]. The underlying aim is to estimate E[φ(T )|X] with the help of bootstrap data augmentation and aggregation of regression trees. In our case, T is a right-censored variable, X a vector of covariates, and φ a known gain function, namely the value of a customer remaining a time T in the portfolio. In the presence of censoring [START_REF] Hothorn | Survival ensembles[END_REF] proposed a generalization extending the random forest algorithm. Using the inverse probability of censoring weights (IPCW) approach, described in [START_REF] Van Der Laan | Unified methods for censored longitudinal data and causality[END_REF], they introduce weights to compensate for censoring. Additionally, [START_REF] Steingrimsson | Doubly robust survival trees[END_REF][START_REF] Steingrimsson | Censoring unbiased regression trees and ensembles[END_REF] combined doubly robust estimators with survival trees. The novelty of our approach is to compute IPCW in situations where censoring depends on the covariates, as it was suggested in [START_REF] Molinaro | Tree-based multivariate regression and density estimation with right-censored data[END_REF], and we will discuss its performance compared to other weighting strategies. The IPCW technique is a very general one: once the weights have been determined, any regression technique can be generalized to this framework -see for example [START_REF] Koul | Regression analysis with randomly rightcensored data[END_REF] for linear regression, [START_REF] Goldberg | Support vector regression for right censored data[END_REF] for support vector machines, and [START_REF] Lopez | Tree-based censored regression with applications in insurance[END_REF] for classification and regression trees (CART). The difficulty in defining appropriate weights is in determining which identifiability assumption is reasonable given the situation. The weights may be considerably different if censoring is allowed to depend on covariates X (see e.g., [START_REF] Lopez | Single index regression models in the presence of censoring depending on the covariates[END_REF]), or not (see e.g., [START_REF] Stute | Kaplan-meier integrals[END_REF]), and, as we will show on real data, this has important consequences on the results of the procedure.

Another modification of random forests to the presence of censoring is the random sur-vival forest (RSF) of [START_REF] Ishwaran | Random survival forests[END_REF] (see also [START_REF] Ishwaran | Random survival forests for r[END_REF]). In this procedure, a log-rank test is used to split the observations at each computational step of the the regression trees. This type of procedure for growing a survival tree was previously proposed in [START_REF] Ciampi | Stratification by stepwise regression, correspondence analysis and recursive partition: a comparison of three methods of analysis for survival data with covariates[END_REF] and [START_REF] Segal | Regression trees for censored data[END_REF]. It was also studied in [START_REF] Leblanc | Survival trees by goodness of split[END_REF] and in [START_REF] Hothorn | Unbiased recursive partitioning: A conditional inference framework[END_REF] in the context of conditional inference trees. Moreover, [START_REF] Zhu | Recursively imputed survival trees[END_REF] have studied the impact of recursively fitted RSF models on prediction quality. [START_REF] Leblanc | Relative risk trees for censored survival data[END_REF] introduced the relative risk tree (RRT) algorithm, where a proportional hazard likelihood is used as the criterion at the splitting step. Forests of relative risk trees were investigated in [START_REF] Hothorn | Bagging survival trees[END_REF] and [START_REF] Ishwaran | Relative risk forests for exercise heart rate recovery as a predictor of mortality[END_REF]. The extension to left-truncation of this approach was considered in [START_REF] Fu | Survival trees for left-truncated and right-censored data, with application to time-varying covariate data[END_REF]. Additional splitting criteria examined in the literature include the exponential log-likelihood in [START_REF] Davis | Exponential survival trees[END_REF], and splits determined through analysis of residuals of the Cox model [START_REF] Ahn | Tree-structured proportional hazards regression modeling[END_REF]). More recently, [START_REF] Zhu | Tree-based methods for survival analysis and high-dimensional data[END_REF] and [START_REF] Zhu | Reinforcement learning trees[END_REF] proposed the reinforcement learning trees (RLT) algorithm which consists in fitting an embedded model at each step of the tree construction to improve the selection of the splitting variable, while [START_REF] Li | Censored quantile regression forests[END_REF] explored censored quantile regression using random forests.

While the theoretical properties of tree and forest-based learning algorithms are not fully understood yet, the consistency of RSF was first investigated in Ishwaran & Kogalur (2010) under the assumption that X takes a finite number of values. These properties were further studied in Cui, Zhu, & Kosorok (2017) and [START_REF] Cui | Consistency of survival tree and forest models: splitting bias and correction[END_REF] which provided a theoretical framework to consider the consistency of survival forests, and established consistency under specific conditions that include random splitting rules and splitting rules with marginal signal checking. [START_REF] Cui | Consistency of survival tree and forest models: splitting bias and correction[END_REF] also underlined the problem of non optimal split selection for usual survival tree methods. The method we study in this article does not suffer from such problem since it requires less assumptions. Finally, a review of the literature on survival trees can be found in [START_REF] Bou-Hamad | A review of survival trees[END_REF].

The rest of the paper is organized as follows. In Section 2, we present the specific details of the types of observations we examine here, and describe the relevant random forest algorithm we have adapted to censoring. A simulation study is performed in Section 3, while the realworld behavior of our procedure is illustrated in Section 4, where we present an application of this weighted random forest to model the churn behavior of health insurance policy holders.

Further results as well as explanations about the choice of the random forest parameters are provided in supplementary material.

2 Description of the method

The survival regression setting

To study the termination risk, we introduce the lifetime of a contract, denoted T . This duration is not directly observed due to the presence of right-censoring (which is a classical issue in survival analysis). Instead of observing T , we observe a pair of variables (Y, δ) defined as

Y = min(T, C), δ = 1 T ≤C .
The need for the censoring variable C is due to the fact that some contracts may not have been terminated at the end of the observation period. Additionally, a vector of covariates X ∈ X ⊂ R p is available, in order to identify certain characteristics that may influence the termination time. Observations correspond to i.i.d. replicates (Y i , δ i , X i ) 1≤i≤n of this set-up. For the sake of simplicity, we consider the case where T and C are continuous random variables. For any continuous random variable U, we will denote by S U (t) = P(U > t) its survival function.

Our aim is to estimate the function f

(x) = E[φ(T )|X = x].
In our application, φ is a pricing function that associates a unitary profit with a customer remaining a time T with the portfolio in question.

Since T is not directly observed because of censoring, classical estimators of f (random forests, as developed below, or any other regression estimator) are biased if no attempt is made to correct for this phenomenon. The procedure we propose is based on the inverse probability of censoring weighting (IPCW, see e.g., [START_REF] Van Der Laan | Unified methods for censored longitudinal data and causality[END_REF]) which, through an appropriate weighting of the observations, aims to suppress this bias. This general method is described in Section 2.2 below.

Inverse probability of censoring weighting 2.2.1 Introduction

IPCW is a general principle which consists of correcting the bias introduced by censoring.

Given a continuous random variable V (possibly multidimensional), the key task is to determine which weights should be put on an uncensored observation in order to retrieve the distribution of the variable of interest.

Proposition 1 (The IPCW principle) Let γ =        1 if V is not censored 0 if V is censored and V = γV.
Let p(v) = P(γ = 1|V = v) and assume that ∀v, p(v) > 0.

Then for any function ψ,

E W • ψ(V ) = E ψ(V) , with W = γ p(V )
.

This result states that given a random variable V that is not always observable, it is still possible to estimate its distribution based on the observations (V i , γ i ) i=1,..,n . This is done by attributing weights γ i p(V i ) to the observations, with p(v) the probability of V being non-censored, given V = v.

In our setting, we can apply the IPCW routine to the vector (T, X). This results in the equality:

E W • ψ(Y, X) = E ψ(T, X) , (1)
where

W = δ p(Y,X) and p(t, x) = P(δ = 1|T = t, X = x) = P(t ≤ C|T = t, X = x).
In the survival setting, it is generally impossible to infer values for p(t, x) since -as it is well known -the identifiability of the model requires assumptions on the dependence between T and C that cannot be statistically verified (see Section 4.1 in [START_REF] Lagakos | General right censoring and its impact on the analysis of survival data[END_REF]). However, with assumptions on the dependence between T and C, it is possible to compute p(t, x).

Let H1 and H2 denote the following hypotheses:

H1 : P(t ≤ C|X = x, T = t) = S C (t), H2 : P(t ≤ C|X = x, T = t) = S C (t|X = x), with S C (•|X = x) the conditional distribution function of C given X = x. Sufficient conditions
for these hypotheses to be satisfied are, respectively, (T, X) C (H1) and T C conditionally on X (H2).

Clearly H2 is more general than H1, but requires estimation of a conditional survival function, which is more tricky. Hence, if the strongest assumption H1 is reasonable, it may be interesting to use it instead of H2, in order to facilitate computation of the appropriate weights.

The obvious drawback is that, if the censoring is dependent on the covariates, H1 does not hold.

In what follows, we consider these two assumptions as separate cases, but both are handled simultaneously. Depending on the hypothesis we make, we note W = δ/S C (Y) (H1) or W = δ/S C (Y|X) (H2), so that in both cases equation (1) holds. We also note for all i ∈ {1, .., n},

W i = δ i /S C (Y i ) (or W i = δ i /S C (Y i |X = X i ))
the exact IPCW weights. However, the function S C (resp. S C (•|X)) is unknown and we have to estimate it in order to compute the weights W.

Computation of the weights

Computation of the IPCW weights requires us to estimate the survival function of the censoring variable. In many applications, the roles of the variables T and C are symmetric, and the censoring variable needs to be studied using survival models. In our work, we have used three strategies to estimate S C (resp. S C (•|X)):

• S C estimated with the Kaplan-Meier (KM) estimator [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF]):

Ŝ C (t) = Y i ≤t 1 -(1 -δ i )/ n j=1 1 Y j ≥Y i .
• S C (•|X) estimated with the Cox model [START_REF] Cox | Regression models and life-tables[END_REF]). This model assumes that the hazard rate has the form: λ C (t|X = x) = λ 0 (t) • e t β•x with β a vector of coefficients, and for some random variable U, λ U = -S U /S U .

• S C (•|X) estimated with the random survival forest (RSF) algorithm [START_REF] Ishwaran | Random survival forests[END_REF]).

An alternative idea for estimating S C (•|X) would be to use a kernel estimator such as the conditional Kaplan-Meier estimator of [START_REF] Beran | Nonparametric regression with randomly censored survival data[END_REF] and [START_REF] Dabrowska | Uniform consistency of the kernel conditional kaplan-meier estimate[END_REF]. Nevertheless, this would rely on kernel smoothing whose behavior deteriorates when the dimension of X increases (i.e. p > 3) as noted in [START_REF] Lopez | Single index regression models in the presence of censoring depending on the covariates[END_REF]. This dimensional constraint is the reason for imposing a stronger assumption on the conditional distribution of the censoring as a feasible compromise.

Given estimators Ŝ

C (resp. Ŝ C (•|X)) of S C (resp. S C (•|X)), the weights ( Ŵi ) i=1,..n are com- puted using the formula Ŵi = δ i / Ŝ C (Y i ) (resp. Ŵi = δ i / Ŝ C (Y i |X i ))
. Therefore, each method used to estimate S C leads to different IPCW weights, then to different estimators of the joint distribution of (T, X). In the following section, the weights ( Ŵi ) i=1,..n refer to any collection of weights computed with one of the three methods.

A weighted random forest algorithm for the regression of right-censored data

Random forest is an ensemble learning algorithm which consists of an aggregation of elemental regression trees (base learners), see e.g., [START_REF] Breiman | Random forests[END_REF] and [START_REF] Biau | A random forest guided tour[END_REF]. A regression tree (see [START_REF] Breiman | Classification and regression trees[END_REF]) produces a partition of a dataset using successive binary splits based on values of the covariates. In each subset obtained (called a leaf ), the observations correspond to a homogeneous group. The key behind group formation is the choice of the splitting criterion (e.g., least squares) which aims to reduce heterogeneity at each step.

The random forest algorithm is a combination of two methods: a classification and regression tree (CART) algorithm to compute each tree, and a bagging algorithm (bagging means bootstrap aggregating, see [START_REF] Breiman | Bagging predictors[END_REF]) to introduce randomness into partition formation by building bootstrap samples. Compared to a single regression tree approach, tree aggregation in random forests stabilizes results by making them less sensitive to changes in the database.

The rest of this section is devoted to adapting the random forest algorithm to the presence of right censoring. In Section 2.3.1, we propose a modification of the splitting criterion to work in our framework. Computation of predictions given a tree is shown in Section 2.3.2. An overall description of the algorithm is given in Section 2.3.3, and we discuss its parameters in Section 2.3.4.

Split selection

In the building of a binary tree, the main interest is in the way splits are determined. Here, we describe the split selection procedure used at each node of the tree, for the growing of one tree of the forest. This tree is built using a bootstrap sample of the initial data.

Let D n be a list of indices in {1, . . . , n} which represents a bootstrap sample of size n of the initial data, drawn uniformly with replacement. Given a set B ⊂ X, let n D n ,w (B) =

i∈D n W i 1 X i ∈B be the weighted number of observations of D n that belong to B, and φD n ,w

B = 1/n D n ,w (B) i∈D n W i φ(Y i )1 X i ∈B the weighted mean of φ(Y i ) for the observations within B (we set φD n ,w B = +∞ if n D n ,w (B) = 0).
Let A ⊂ X be a node of a tree, j ∈ {1, . . . , p} a variable, and u real number. Define:

L w (u, j, A, D n ) = 1 n D n ,w (A) • i∈D n W i • φ(Y i ) -φD n ,w A l 1 X ( j) i ≤u -φD n ,w A r 1 X ( j) i >u 2 • 1 X i ∈A ,
where

A l = A ∩ {X ( j) ≤ u} and A r = A ∩ {X ( j) > u}.
We also define n D n , ŵ(B), φD n , ŵ B and L ŵ(u, j, A, D n ) as the same quantities with the weights W i replaced by the weights Ŵi .

The binary split chosen at node

A is A = A * l ∪ A * r , with A * l = A∩{X ( j * ) ≤ u * }, A * r = A∩{X ( j * ) > u * },
and (u * , j * ) given by:

(u * , j * ) = argmin j ∈ {1, .., p} u ∈ X ( j) i : i ∈ D n s.t. X i ∈ A L ŵ(u, j, A, D n ).
For a heuristic understanding of the algorithm, consider a fixed set A. We abusively assume that A is independent from the data, which is not correct in practice since A is selected from the data. Under this assumption, we can study the asymptotics of L w (u, j, A, D n ) since the observations (Y i , X i ) i∈D n which belong to A are i.i.d. and follow the conditional distribution

L ((Y, X)|X ∈ A).
Until the end of this section, we assume that A is independent from the data.

With A z standing for either A l or A r , we have from the law of large numbers and equation

(1),

φD n ,w A z -→ n→+∞ E W • φ(Y)|X ∈ A z = E φ(T )|X ∈ A z , and
1 n D n ,w (A z ) • i∈D n X i ∈A z W i • φ(Y i ) -φD n ,w A z 2 -→ n→+∞ E W • φ(Y) -E[W • φ(Y)|X ∈ A z ] 2 X ∈ A z = E φ(T ) -E[φ(T )|X ∈ A z ] 2 X ∈ A z .
Then, by splitting the sum in L w (u, j, A, D n ) between A l and A r , we get:

L w (u, j, A, D n ) -→ n→+∞ L ∞ (u, j, A), with L ∞ (u, j, A) = P(X ∈ A l ) P(X ∈ A) • E (φ(T ) -E[φ(T )|X ∈ A l ]) 2 X ∈ A l + P(X ∈ A r ) P(X ∈ A) • E (φ(T ) -E φ(T )|X ∈ A r ) 2 X ∈ A r .
Therefore we observe that, at the cost of replacing Ŵi by W i , the selected split is asymptotically the one which minimizes the sum of the within variances of the two child nodes, weighted by their relative importance.

It is also important to see that given D n and the splitting criteria L, the growing of a tree is deterministic. Therefore we can define a function T (D, L) which outputs a tree given a sample and a splitting criterion.

Tree prediction

The growing of a binary tree T results in a partition (X k ) k=1,..,K of the input space. In order to predict the mean value of the target variable, a natural estimator is

m ŵ(x) = K k=1 φD n , ŵ X k 1 x∈X k . Under H1/H2 we have for the denoised version m w , m w (x) -→ n→+∞ K k=1 E φ(T )|X ∈ X k 1 x∈X k ,
which shows that m ŵ approximates a piecewise constant estimator of f . However, though the growing step is done using a least-square criterion -which is welladapted to mean regression -the prediction made in the terminal leaf may not necessarily be a sample mean (e.g., it may be a quantile -see [START_REF] Meinshausen | Quantile regression forests[END_REF]). We can formalize this by defining a function M(T , D, m) which outputs a predictor given a binary tree, a sample, and a type of terminal leaf-based estimator. In the following section, we propose a second terminal leaf estimator for the mean within a leaf.

Both for tree prediction and split selection, note that the method we propose is a generalization of the classical CART regression algorithm. Indeed, setting the weights W i all equal to 1 results in Breiman's original algorithm. We note L and m the split criteria and the terminal leaf estimator obtained when setting the weights W i all equal to 1.

Description of the algorithms

The tools we developed in previous sections can be used in different random forest algorithms.

We have studied three survival weighted random forests (swRF) in this article: swRF1, swRF2 and swRF3. These algorithms differ in the way weights Ŵi are introduced.

In swRF1 the weights Ŵi are taken into account in the bootstrap part of the random forest.

Based on the whole training set, we build Ŝ C (or Ŝ C (•|X)) and deduce ( Ŵi ) i=1,...,n . Then, all bootstrap samples of the forest are drawn with replacement and with sampling probabilities proportional to the weights ( Ŵi ) i=1,...,n . A CART tree is then grown on each bootstrap sample, using the classical algorithm (i.e. with L and m which denote the splitting criterion and the terminal leaf estimator with uniform weights).

In algorithms swRF2 and swRF3, bootstrap samples are drawn uniformly with replacement, as it is done in the classical random forest algorithm, and the weights are computed on the bootstrap samples. This means that the weights of a given observation may vary from tree to tree.

The only difference between swRF2 and swRF3 is the terminal leaf estimator used. In swRF2, m ŵ is used as described in Section 2.3.2, whereas in swRF3 we use

m ŵloc (x) = K k=1 φD n , ŵloc X k 1 x∈X k .
In fact, m ŵloc is similar to m ŵ except that the weights used in each terminal leaf are computed using a Kaplan-Meier estimator inside the leaf.

Pseudocode of the three algorithms we study is shown in Algorithms 1 and 2 below.

Algorithm 1: swRF1 Input : Data : (Y i , δ i , X i ) i=1,..,n , M > 0 : number of trees Output: Ensemble predictor swRF1 1 Compute weights ( Ŵi ) i=1,..,n : Ŵi = δ i / Ŝ C (Y i ) (or Ŵi = δ i / Ŝ C (Y i |X i )) 2 for j = 1, . . . ,
: i ∈ D n, j ) : ŴD n, j i = δ i / Ŝ D n, j C (Y i ) or ŴD n, j i = δ i / Ŝ D n, j C (Y i |X i ) 4 Build T j = T (D n, j , L ŵ) 5 mj = M(T j , D n, j , m ŵ) if mode = 1 ( mj = M(T j , D n, j , m ŵloc ) if mode = 2) 6 end return : mswRF(2,3) = 1 M M i=1 mj
Remark 1 For each of swRF1, swRF2 or swRF3, we can build three models which correspond to the three ways to estimate the weights (W i ) i=1,...,n developed in Section 2.2.2.

A practical limitation of the weighted random forest algorithm is that for a large noncensored observation Y i , the corresponding weight tends to be very big (Cui, Zhu, & Kosorok (2017)). This is even more so in the conditional case if X is multidimensional, since the estimated survival function Ŝ C (•|X = x) may decrease quickly for some values of x. Overly-large weights may give too much importance to a single observation, both in terms of split selection and terminal leaf estimation.

We propose two ways to deal with this problem. The first is to threshold the weights Ŵi so that the maximum ratio between two nonzero weights does not exceed some value. Let r max be the maximum ratio allowed between two nonzero IPCW weights. The choice of the value for r max corresponds to a typical bias-variance trade-off. Indeed, thresholding the weights introduces bias in the swRF procedure since large non-censored observations get smaller weights, and thus predicted values tend to be lower, but on the other hand this lowers the estimator's variance, which is high when some observations have very large weights.

The second strategy is to use swRF to model the time T = min(T, T max ) and not the original time T . This is natural because we are working on the estimation of the expected value for T , and the problem of estimating the average of a right-censored time is generally hard to solve; indeed, there is high variance in the estimation of the tail of the distribution, which causes fluctuations in the estimated mean.

In the following, we will combine the two strategies to achieve good accuracy using swRF.

Parameters of the random forest algorithm

In our earlier presentation of random forests, we made no mention of the parameters which need to be specified to build them. Since there is a large number of parameters, and since parameters may differ according to the random forest implementation in question, we only mention the most common ones -those which have been shown to have the strongest impact on the resulting model.

The main parameters in a random forest are summarized in Tab. 1. The parameters minlea f and maxdepth impact tree size in the forest, and we discuss their influence on the performances of the algorithms in supplementary material. As for the parameter mtry, it is especially important in high dimensional settings.

Parameter Description minlea f minimum value for the number of observations that should be present in a leaf; a split that would result in a leaf of less than minlea f observations cannot be selected maxdepth length of the longest downward path from the root node to a leaf of the tree; a tree only consisting of the root node has depth 0. mtry number of candidate variables that are randomly sampled at each node when choosing the best split

Tab. 1: Parameters of the random forest algorithm.

Assessing the quality of a model's fit

In the applications which follow, our strategy to evaluate a fitted model relies on train-test approaches: let D tr and D te be train and test sets of indices. However, due to the right censoring, the usual accuracy criterion used in regression settings cannot be computed on the test set.

In fact, since the problems of bias on the test and train samples are similar, we can use for model validation the same sample fitting method with IPCW that we use for model training.

Let ( Ŵi ) i∈D te be the estimated IPCW and f the predictor function of a fitted model. Then, the In our simulations, the covariates X are distributed as marginal uniform distributions on [-1, 1] linked through a Gaussian copula with covariance given by an AR(1) covariance matrix with coefficient ρ, i.e., a matrix K ρ such that K ρ i, j = ρ |i-j| . The coefficient ρ is random and uniformly distributed on [0, 0.6].

quantity 1/n te • i∈D te Ŵi • (φ(Y i ) -f (X i )) 2 is
We investigate three experimental settings which correspond to three different models for L (T |X), the distribution of T given X. In the following description of the three cases, λ T and k T are constants that will be specified later, and Weibull(λ, k) refers to the distribution with density

g(u) = k/λ • (u/λ) k-1 e -(u/λ) k 1 u≥0 : • Case 1 (Weibull): T |(X, β T ) ∼ Weibull(λ T e -t β T •X , k T ), with β T ∼ N(µ = 0, σ 2 = 0.04 • I p )
and β T X.

• Case 2 (Independent mixture of Weibulls): T |(X, G, (β T, j ) j=1,...,4 ) ∼ Weibull(λ T e -t β T,G •X , k T ) with:

-G ∼ Uni f {1, 2, 3, 4} and G X, (β T, j ) j=1,...,4 , -(β T, j ) j=1,...,4 i.i.d. with ∀ j ∈ 1, 4 , β T, j ∼ N(µ = 0, σ 2 = 0.09 • I p ) and (β T, j ) j=1,...,4 X.

• Case 3 (Covariate-dependent mixture of Weibulls, p ≥ 2):

T |(X, G, (β T, j ) j=1,...,4 ) ∼ Weibull(λ T e -t β T,G •X , k T ) with: -G = 1 X 1 ≥0,X 2 ≥0 + 2 • 1 X 1 ≥0,X 2 <0 + 3 • 1 X 1 <0,X 2 ≥0 + 4 • 1 X 1 <0,X 2 <0 , -(β T, j ) j=1,..,4 i.i.d. with ∀ j ∈ 1, 4 , β T, j ∼ N(µ = 0, σ 2 = 0.04 • I p ) and (β T, j ) j=1,...,4 X.
The nature of L (C|X) is the same in every setting:

C|X ∼ Weibull(λ C e -t β C •X , k C ) with λ C
and β C fitted to satisfy certain conditions (see Algorithm 3).

In each setting, we test the algorithm for two different functions φ: φ(t) = t and φ(t) = log(t + 1). We also assess the sensitivity of the results to the censoring rate of the simulated data (q = 0.1, 0.3 or 0.5) and to the strength of the dependence between C and X (measured in terms of percentage of explained variance: R2 C = 0.05 or 0.1).

The process of data simulation is detailed in Algorithm 3.

Algorithm 3: Simulated data generation

Input : n : number of simulated observations p : dimension of X q : proportion of censored observations R2 C : proportion of explained variance for C|X k C : shape parameter of C (λ T , k T ) : scale and shape parameter for T T max : threshold for T Output: Dataset of generated data

1 Draw ρ uniformly on [0, 0.6] 2 Simulate (X i ) i∈1,...,n i.i.d. with X ∼ GaussianCopula(K ρ ) and marginals Uni f [-1, 1] 3 Generate (T i ) i=1,...,n from case 1, 2 or 3 4 Draw β C 0 ∼ N(µ = 0, σ 2 = 0.04 • I p ), let β C = η • β C 0 ,
and calibrate η so that the empirical estimate of the explained variance R2 C satisfies R2 C ≈ R2 C 5 Calibrate λ C so that the empirical proportion of censored observations q satisfies q ≈ q 6 Generate (C i ) i=1,...,n with ∀i,

C i ∼ Weibull(λ C e -t β C •X i , k C ) 7 Check on the simulated (C i ) i=1,..,n that R2 C ≈ R2 C and q ≈ q 8 Build ∀i = 1, . . . , n, Y i = min(T i , C i ), δ i = 1 T i ≤C i 9 Build ∀i = 1, . . . , n, T i = min(T i , T max ), Y i = min(T i , C i ), δ i = 1 T i ≤C i (cf. Section 2.2.2.
Remark: q at step 5 is derived from T and not T ) return : Dataset with columns X, T , C, Y, δ, T , Y , δ

We take n = 2000, p = 6, λ T = 100, k T = 1, and k C = 0.8 as parameters for the data simulation. The parameter T max is set to 1000.

Different models

In our experiments, we compare the performance of sixteen models. First, we have three swRF1 models, associated with , Cox (swRF12) and random survival forest (swRF13) weights, respectively.

For swRF2 and swRF3, many weight computations are necessary since weights are computed on the bootstrap samples. The use of RSF to compute the weights is thus computationally intensive, so for these two models we only use KM (swRF21 and swRF31) and Cox (swRF22 and swRF32) weights.

We also test five other algorithms to compare with swRF. They all follow the same idea: first, fit a model to estimate S T (•|X) (which gives an estimator Ŝ T (•|X)), and second, integrate

Ŝ T (•|X) to get an estimator of E[φ(T )|X = x]: f (x) = -φ • d Ŝ T (•|X = x).
Five different models are used to estimate S T (•|X) at the first step : the Cox model, the RSF algorithm, the forest of relative risk trees (RRT) algorithm of [START_REF] Ishwaran | Relative risk forests for exercise heart rate recovery as a predictor of mortality[END_REF], and two versions of the reinforcement learning trees (RLT) algorithm described in the Section 4 of Zhu (2013) (one which uses reinforcement learning and one which does not). This gives five benchmarks for swRF that we call respectively RSFr, Cr, RRTr, RLTr (with reinforcement) and nRLTr (without reinforcement).

Since here we are working with simulated data, it is interesting to consider, as a baseline indicator, the performance of the random forest algorithm as if there was no censoring in the data. This random forest is trained on the non-censored data (X i , T i ) i=1,..n , and the associated model is denoted RF. Similarly, it is interesting to consider the exact IPCW weights W i , as described in Section 2.2.1, for the use of swRF. The exact IPCW weights correspond to theoretical weights computed using the true conditional survival function of the censoring S C (•|X) which is known in our simulated data application. This gives three additional comparison models called swRF14, swRF24 and swRF34.

So that the different random forest models involved are comparable, all of the random forests in the simulations are set with the same values for the three parameters we evoked in Section 2.3.4. The parameter mtry is set to 6 so that mtry = p and there is no randomness involved in split selection. We also take maxdepth = 4 and minlea f = 50, as justified in supplementary material.

The swRF models are trained on the threshold data (Y , δ , X) with Y = min(T , C), setting r max = 50 (see the definitions of T and r max in the Remark 1 of Section 2.3.3). RF is trained on the data (T , X), whereas Cr, RS Fr, RRT r, RLT r and nRLT r are trained on the data (Y, δ, X). For the latter models, we then compute the prediction at the point

x with f (x) = -φ•d Ŝ T (•|X = x), where Ŝ T (t|X) = Ŝ T (t|X) • 1 t≤T max .

Results and analysis

Performance of the models

To compare the models, we generated 100 simulated datasets of 2000 observations with Algorithm 3. For each iteration, we train all models on the same 1000 observations, and evaluate them on the 1000 remaining ones. Model accuracy is measured in terms of the mean squared error (MSE), which can can be computed in the simulated data context since we have access to all of the φ(T i ) values. The means of the MSE over the 100 i.i.d. replicates of the simulation process are given in Fig. 1. To keep the figure visually understandable, we chose to represent in the main text only the results for the most illuminating models. Nevertheless, we emphasize that the following analysis is consistent with the complete results, obtained for the sixteen models, that are given in supplementary material.

The results on simulated data help us to learn more about the various models. An initial observation is that the results are quite affected by the φ function being considered, and less by changes in the distribution of L (T |X). Of course, RF is generally the best model, except when the constrained form of the Cox model is well-suited to the problem, whereby Cr may outperform RF. Recall that the RF algorithm should only be seen as a benchmark that cannot be used in practice, since it relies on the complete (uncensored) observations. For q = 0.1, the performance of the swRF models is very close to that of RF, which is natural since swRF is equivalent to RF in the non-censored case (i.e., q = 0).

We can also see that swRF models are more sensitive to an increase in censoring rate than the comparison methods which model directly S T (•|X), namely: Cr, RS Fr, RRT r and RLT r.

The swRF models perform well for q = 0.3 and 0.1, but the MSE is much larger when q = 0.5.

For φ(t) = t, RS Fr and Cr perform very well overall, but for φ(t) = log(t + 1), the swRF models are usually more accurate, especially when q = 0.1 and q = 0.3, or the data is simulated as a covariate-dependent mixture (Case 3).

We can also compare the swRF models with each other. We first remark that the swRF scores are organized in terms of the groups of swRF (swRF1 and swRF3). While for φ(t) = t and q = 0.5 swRF3 models achieve high MSE due to particular iterations where swRF3 does not work well, the results for φ(t) = log(t + 1) are very different for swRF1 and swRF3. The difference is mostly due to the terminal leaf estimator used (this is confirmed by the results of swRF2 in supplementary material). Indeed, we can see that the results of swRF3 are closer to the results of RS Fr, RRT r and RLT r, and the terminal leaf estimators used in these models are almost the same; swRF3 relies on a within leaf nonparametric estimation of the distribution of T using a KM estimator (i.e. m ŵloc ), whereas the other models rely on a Nelson-Allen estimator [START_REF] Aalen | Nonparametric inference for a family of counting processes[END_REF]) as explained in [START_REF] Ishwaran | Random survival forests for r[END_REF]. Also, the type of weights we consider has a second-order impact on the results, which becomes more significant with larger censoring rates. However, it is interesting to note that comparing the results of swRF11 and swRF13, we see that conditional weights computed with RSF tend to give better results than the Kaplan-Meier ones. By comparing the results of swRF32 and swRF34, we observe that For each setting, the mean of the MSE values over 100 i.i.d. replicates of the simulation process is shown. The censoring rate q is equal to 0.1, 0.3, or 0.5, while the percentage of explained variance of C given X: R2 C , is set to 0.05 or 0.1.

Cox IPCW and exact IPCW lead to very close MSE. This shows that Cox IPCW (but also RSF IPCW as we will see in the next section) are good proxies for exact IPCW.

Correlations between weighted MSE (IPCW) and MSE

The results shown in Fig. 1 are estimated MSE which do not suffer from the effects of censoring.

Since such estimators are not available for real data studies, in practice we rely -as explained in Section 2.4 -on weighted estimators (involving KM, Cox, RSF, or uniform weights) when comparing models. Therefore, it is of interest to compute, in the simulated data case, correlations between the non-censored MSE and weighted approximations of it. We added to the results the correlations with the exact (theoretical) weights which are available in the simulated data case.

The average values of these correlations are displayed in Tab. 2. These results show that RSF and Cox weights better replicate the non-censored MSE than the KM ones, and therefore offer better comparison criteria for model selection. Let us note that RSF and Cox weights are even more accurate than exact weights. Such type of phenomena have already been observed in other censored regression models. For example, in the case of linear regression, [START_REF] Koul | Regression analysis with randomly rightcensored data[END_REF] noticed that the asymptotic variance of their slope estimator was smaller using KM weights rather than the true distribution function of the censoring (see their Remark 4.5 p. 1280). As a possible interpretation, as pointed by [START_REF] Koul | Regression analysis with randomly rightcensored data[END_REF], one may claim that the information contained in the censored observations would be completely lost if we were using exact weights, while estimation of the distribution of the censoring relies on the censored observations. On the other hand, the criterion based on uniform weighting of the non-censored observations has little correlation with the non-censored criterion, so we should not rely on it. Mean Spearman correlations between the non-censored and weighted (IPCW) MSE, as a function of q and R2 C .

For each of the six cases, we have 600 calculations of the MSE and weighted MSE (2 choices for φ × 3 choices for L(T |X) × 100 iterations), and we calculate correlations between different fitting criteria.

4 Real data application

Modeling the churn behavior of policy holders

The problem of predictive modeling has gained interest in recent years in the insurance community, and optimization of underwriting is one relevant application [START_REF] Frees | Predictive modeling applications in actuarial science[END_REF]). For insurance contract brokers, forecasting customer value represents an opportunity to improve margins in a competitive environment [START_REF] Cummins | The economics of insurance intermediaries[END_REF], [START_REF] Maas | How insurance brokers create value-a functional approach[END_REF]). Moreover, churn modeling is important for the estimation of customer value, as described in [START_REF] Verhoef | Predicting customer potential value an application in the insurance industry[END_REF]. Here, we focus on applying the methodology developed in Section 2 to build a predictive model of the impact of churn on prospect value (as defined in Section 1) as predicted by insurance brokers.

A broker's approach to estimate the prospect values is given by the following multiplicative formula, where the hats indicate we are referring to estimated quantities:

value = p sub • pr • f ew • f ch ,
with p sub the probability of subscription, pr the predicted premium, f ew the early withdrawal factor, and f ch the churn factor. We are only interested in modeling f ch here. In this factorization, the churn factor only depends on the termination time of the contract via a function φ ch as shown in Fig. 2.

The data that we study here is the customer base of an insurance broker from 1 st October 2009 to 31 st July 2016, and we focus in particular on complementary health insurance contracts. Data are available about the effective dates of the contracts, current states of contracts (active, terminated), and termination dates of contracts (if terminated). Before the underwriting, prospective information is given as 6 variables: age (8 levels starting from 18 years old), gender (2 levels: female, male), number of children insured (5 levels: 0, 1, 2, 3, ≥ 4), social security regime (7 levels: agricultural employee, employee, retired, retired self-employed, selfemployed, student, unemployed), level of insurance (3 levels: low, medium, high), geographical zone (4 levels: Ile-de-France region, north, southern, other). Other client characteristics are available in the database but are only known after subscription occurs, hence it is impossible to use them to evaluate prospects.

We use this information to predict churn factors, which correspond to the variable φ(T ).

This variable is censored for a contract which is still active on January 15 th , 2016. The censoring Here, φ ch is the broker's commission, expressed per unit of annual premium. It divides into two parts: a withholding part during the first year -which reaches 50% of the annual premium after one year -and a subsequent part, starting after one year, equal to 10% of annual premiums (taking into account annual resets of the premium).

variable C is the age of the contract, i.e., the duration between the date of effect of the contract and January 15 th , 2016. The 6 variables above constitute the covariate vector X.

Additional details about the experiments

The same models as in the simulated data application (see Section 3.1.2) are compared in this real data application, excluding RF and the swRF models based on the exact weights (swRF14, swRF24, and swRF34), which cannot be used in the real data case. Also, T max is set at 1465 days; this makes sense from the point of view of this application since we are then estimating the expected return of a prospect within the four years following subscription. All swRF, RS Fr, RRT r and RLT r models are set with the parameters mtry = 6, maxdepth = 5, and minlea f = 100. For swRF models, r max is set to 50.

We test the various models with four different φ functions: φ = φ ch (Fig. 2), φ(t) = log(t+1), φ(t) = t, and φ(t) = 1 t>380 . For the latter function, the choice of the threshold 380 is related to the willingness to consider churns that occur in the first contract year. Most of these churns occur before 365 days, but a significant number is reported a few days later, explaining the margin of 15 days that we take (T max was set to 381 in this case since, from 381 days on information on φ(T ) is not censored). 2.8 RS Fr 5.9 Cr 6.9 RRT r 6.9 RLT r 3.7

Tab. 3: Mean ranks of the models; φ = φ ch , with RSF weights.

Conclusion

In this paper, we have considered a class of weighted random forest algorithms, where the weight put on each observation is designed to compensate for censoring. A classical issue in censored regression models is the identifiability assumption that defines the dependence structure between the censoring and the variables involved in the model. Therefore, we have distinguished between two situations, namely the case where censoring is independent, and the case where this variable is allowed to depend on the duration variable T through the covariates.

The latter case is the more general of the two, but leads to difficulties in computing appropriate weights, due to the strong impact of the dimension of the covariate vector. For this reason, we have proposed a compromise by modeling the conditional censoring distribution using a Cox model or a RSF model. We showed on simulated data that our method is competitive with other statistical methods in terms of accuracy. Also, this approach appeared to give the most accurate results in our application to modeling the commission received by an insurance broker, which is a non linear function of the time at which an insurance policyholder surrenders their contract.

Moreover, it gives more accurate results than competing approaches such as random survival forest, relative risk forest, and reinforcement learning trees, in the setting we have considered.

Finally, the weighting strategies we have proposed easily generalize to other regression-based approaches such as, to give one example, quantile regression.

As future work, it would be interesting to investigate the performances of the weighted random forest method in a setting where the dimension of the covariate vector is higher, which is an important application case [START_REF] Zhu | Tree-based methods for survival analysis and high-dimensional data[END_REF], Ishwaran et al. (2010)). Indeed, in this situation the estimation of the conditional IPCW is hard and may require to adapt the algorithm. A suitable method might be to sample, for each tree of the forest, a subset of covariates to take into account in the computation of the conditional IPCW. Also, doubly robust survival trees studied in [START_REF] Steingrimsson | Doubly robust survival trees[END_REF] seems to be a great research direction. Promising results are already presented in [START_REF] Steingrimsson | Censoring unbiased regression trees and ensembles[END_REF] using a relative risk tree to estimate S C (•|X), and we might wonder if we could observe the same phenomenon in the results as we obtain in our work regarding the influence of the estimation of S C (•|X): i.e. an improvement when using a conditional estimator instead of the KM estimator, and a benefit of using ensemble models such as RSF. Finally, theoretical work may be of interest as well; since the weighted random forest generalizes the (non-censored) random forest for regression, it may be possible to transpose consistency results obtained in the non-censored case (e.g. in [START_REF] Scornet | Consistency of random forests[END_REF]) to the weighted case.

M do 3

 3 Draw D n, j : sample n observations w.r.t. weights ( Ŵi ) i=1,..,n , with replacement 4 Build T j = T (D n, j , L)5 mj = M(T j , D n, j , m) swRF2 & swRF3 Input : Data : (Y i , δ i , X i ) i=1,..,n , M > 0 : number of trees, mode ∈ {1, 2} : type of terminal leaf estimator Output: Ensemble predictor swRF2 (or swRF3) 1 for j = 1, . . . , M do 2 Draw D n, j : sample n observations uniformly with replacement 3

  an estimator of the mean squared error (MSE) of the model f (with n te = Card(D te ) the number of test observations). This is the method suggested in[START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF] and the one we adopt in this article. One disadvantage is that again this approach raises the question of the choice of weights to consider for IPCW; in the same way as for training, we choose here to compute the validation error corresponding to three types of weight: KM, Cox and RSF. Therefore, we use not one but three criteria to compare model performance. Computation of the IPCW weights is performed separately on the training and test sets.The problem of model validation in the presence of censoring has received considerable attention in the literature, and other performance measures exists. The C-index[START_REF] Harrell | Evaluating the yield of medical tests[END_REF]), which generalizes the Kendall tau for right-censored data and is related to the area under the ROC curve[START_REF] Heagerty | Survival model predictive accuracy and roc curves[END_REF]), is very popular, so we also compute it when comparing models.

Fig. 1 :

 1 Fig. 1: Results on simulated data.

Fig. 2 :

 2 Fig. 2: The churn factor as a function of termination time.

Fig. 3 :

 3 Fig. 3: Results on real data.

Fig. 4 :

 4 Fig. 4: Boxplot of the performances (MSE) of the models; φ = φ ch , with RSF weights.

Acknowledgements

We would like to thank the company Santiane, who provided the data that served for the experiments, and was always available to answer our questions.Also, we greatly thank the two referees for their valuable comments and their insightful suggestions.

Software

All analyses were performed with the R package sword (github.com/YohannLeFaou/sword) which was developed for the purposes of this article. The code used for producing the results is available at the address: github.com/YohannLeFaou/impact-churn-health-insurance.

Results and analysis

As in the simulated data case, we use a train-test approach (5000 observations each) for model testing, with 100 repetitions. Each training and testing set are drawn uniformly without replacement from the original dataset in such a way that the training and test sets do not overlap. The results in terms of MSE computed with KM, Cox and RSF weights, and on the C-index, are shown in Fig. 3. We give the results for the same models as in Section 3.2.1 except for RF and swRF34 which are not usable anymore with real data, and for swRF22 that we choose to incorporate into the figure. The results for the other models are given in supplementary material.

We first observe that the results for φ = φ ch , φ(t) = log(t + 1), and φ(t) = t look very similar, especially for φ = φ ch and φ(t) = log(t +1). Hence, the strong influence of the function φ that we noticed for the simulated data does not carry over to the real data. Moreover, there is a tendency in the results that swRF models fitted using a certain type of weights perform very well with the MSE computed using the same type of weights. This is especially clear for the Cox and KM weights (swRF11 and swRF22 models), indicating that the type of weights considered has a real impact on the estimated distribution of (T, X).

We choose to consider the RSF weighted estimation of the MSE as the reference criteria for comparing the models. Indeed, we saw in Section 3.2.2 that Cox and RSF weights give approximations of the MSE that are the most correlated with the MSE computed on non-censored data. In addition, there is a risk of overfitting the type of weights we use, and since RSF weights are only involved in the model swRF13, they constitute a good choice. In Fig. 4 and Tab. 3, we show the performance of each model, calculated with RSF weights and for φ = φ ch . Fig. 4 indicates the model swRF32 achieves on average the lowest error, followed by RLT r, an ordering confirmed by the table of rank statistics.

The C-index's statistic [START_REF] Harrell | Evaluating the yield of medical tests[END_REF]) corresponds to the proportion of ordered pairs in the test set which are well-ordered by the model f , i.e. :

.

We can see in Fig. 3 that RS Fr, Cr, RRT r and RLT r achieve the best C-index out of the set of models. This illustrates that it is important to consider a quadratic error criterion rather than a rank one when the goal is to estimate a mean value. It is the case in our situation since we are interested in getting the best prediction for the churn factor f ch = φ ch (T ). The C-index values nevertheless show the benefit of using conditional weights within the swRF algorithm in order to get models that rank the observations well.


M Melosso

L Bizzocchi
email: bizzocchi@mpe.mpg.de

O Sipilä

B M Giuliano

L Dore

F Tamassia

M.-A Martin-Drumel

O Pirali

E Redaelli

P Caselli

First detection of NHD and ND2 in the interstellar medium

Keywords: Astrochemistry -ISM, molecules -Line, identification -ISM, abundances -Submillimeter, ISM

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Despite the low cosmic abundance of deuterium (D/H = 1.5 × 10 -5 , Linsky 2003), deuterated molecules can be found in the interstellar medium (ISM) with enormously large amounts and are routinely used as proxies for the cold environments typical of the early stages of the star formation process. It is well known that deuterium fractionation processes, i.e., the deuterium enrichment in a given molecular species, can occur at very high rate under some specific physical conditions including, above all, a low gas temperature (5-20 K) and strong CO depletion (Caselli & Ceccarelli 2012). Inheriting such conditions from their pre-stellar core phase, Class 0 protostars are the astronomical objects which exhibit the largest molecular deuteration found so far [START_REF] Ceccarelli | Protostars and Planets VI[END_REF]). In the last few decades, several deuterium-bearing species (including doubly-and triplydeuterated forms) have been detected in these sources, e.g., water D 2 O (Butner et al. 2007), formaldehyde D 2 CO [START_REF] Ceccarelli | Protostars and Planets V[END_REF], hydrogen sulfide D 2 S (Vastel et al. 2003), methanol CD 3 OH (Parise et al. 2004), and ammonia ND 3 (Roueff et al. 2005;van der Tak et al. 2002).

The Class 0 protostar IRAS 16293-2422, located in the nearby ρ-Ophiuchi star-forming region (at a distance of ∼120 pc), is considered a specimen for super-deuteration phenomena [START_REF] Ceccarelli | Protostars and Planets V[END_REF]. Indeed, besides an exceptionally rich chemistry in its inner and warmer regions (van Dishoeck et al. 1995;Jørgensen et al. 2012;Fayolle et al. 2017), this source harbors a large number of interstellar deuterated species in the outer and cold envelope [START_REF] Loinard | [END_REF]). Among them, the family of neutral nitrogen hydrides (NH, NH 2 , and NH 3 ) is well-represented, given that all the deuterated forms of ammonia (NH 2 D, NHD 2 , and ND 3 ) and the imidogen radical ND (Bacmann et al. 2010(Bacmann et al. , 2016) ) have been detected so far. However, deuterium enrichment has never been observed for the amidogen radical NH 2 , neither in IRAS 16293-2422 nor in any other astronomical sources.

Several factors may be imputed to the current lack of observation of interstellar NHD and ND 2 . Until recently, laboratory data were limited to the frequency region below 500 GHz (Kanada et al. 1991;Kobayashi et al. 1997). At low temperatures, neither species possess any bright transitions in this portion of the spectrum, with b-type fundamental transitions lying high into the submillimeter-wave domain (around 770 GHz for NHD and around 527 GHz for ND 2 ). Furthermore, the amidogen radical is a textbook floppy molecule, for which frequency extrapolation from low frequency measurements is extremely unreliable. Besides the deficiency of laboratory spectral data, groundbased astronomical observations of these light species are hindered by the atmosphere opacity at submillimeter-wavelength. Therefore, even with high-frequency facilities (such as ALMA or SOFIA), the observation of light hydrides remains a challenging task. On-board observations using the Herschel Space Observatory mission, during which numerous projects attempted at gaining more insights into the chemistry of the ISM, offer a unique chance to overcome the latter issue. Also, now that recent laboratory studies of nitrogen hydrides in their various deuterated forms have been reported up to the terahertz domain (Melosso et al. 2017(Melosso et al. , 2019a(Melosso et al. ,b, 2020;;Bizzocchi et al. 2020), confident searches for these species can be undertaken in Herschel archival surveys.

In this paper, we report the first detection of NHD and ND 2 in the interstellar medium, thanks to HIFI1 observations carried out as part of the CHESS2 key program (PI: C. Ceccarelli). Both species were observed in absorption against the continuum of the low-luminosity protostar IRAS 16293-2422. From the modeling of their hyperfine structures, column densities of NHD and ND 2 have been determined and the deuterium fractionation of amidogen radical evaluated.

Observations

The observations analysed in the present paper have been retrieved from the HSA3 . They were collected as part of CHESS key program, whose observational details have been given extensively in Ceccarelli et al. (2010). Briefly, the Class 0 protostar IRAS 16293-2422 was observed in March 2010 with the HIFI spectrometer (de Graauw et al. 2010;Roelfsema et al. 2012) for a total of 50 h. The targeted coordinates were α 2000 = 16 h 23 m 22.75 s , δ 2000 = -24 • 28 34.2 . Seven different HIFI bands were employed to perform a spectral survey covering the frequency range between 480 and 1902 GHz almost continuously. The observations were carried out in double-sideband (DSB) using the wide band spectrometer (WBS), whose resolution is ∼ 1.1 MHz.

For the present study we have employed bands 1a (488-552 GHz), 2b (724-792 GHz), and 4a (967-1 042 GHz), corresponding to the IDs number 1342191499, 1342192332, and 1342191619, respectively. HSA provides level 2.5 data products obtained with HIFI pipeline of the Herschel interactive processing environment (HIPE, Ott 2010). These data generally offer a scientific quality sufficient to perform the line analysis. However, given the weakness of the NHD and ND 2 spectral signals, we have also gone through the final processing manually. Starting from level 2 products, we have inspected all the data, manually removed bad scans, and exported all the "good" 1 GHz chunks to CLASS/GILDAS4 format.

Side-band deconvolution has been performed with the minimization algorithm of [START_REF] Comito | [END_REF] implemented into CLASS90, and the resulting spectra have been crosschecked with the original level 2.5 products. In a few instances, the deconvolution procedure produced a small increase of the noise level compared to the averaged level 2 scans. In such cases, we directly use the DSB spectra after verifying that no interfering emission was present in the image band. At the later stage, intensive baseline subtraction on H-and V-polarisation spectra has been applied in the vicinity of each target feature; then, the two set of data have been re-sampled and averaged to obtain the final spectrum.

Analysis

The presence of amidogen radical in IRAS 16293-2422 has already been established by Hily-Blant et al. (2010), during the early phase of the CHESS key program. The results obtained in the previous study for the parent species NH 2 have been used to test the validity of our data processing and the subsequent spectral analysis. Figure 1 shows the strongest J = 3 2 -1 2 component of the fundamental N Ka, Kc = 1 1,1 -0 0,0 transition of o-NH 2 around 952 GHz. Here, most of the HFS components are resolved and the intensity of the absorption lines is well above the 3σ noise level threshold. The resulting model of the full HFS (plotted with the solid blue line in Fig. 1) has been computed using a custom Python3 code which exploits a similar approach to that of the HFS method implemented in CLASS. Differently from the latter, our code uses the column density (N) and the excitation temperature (T ex ) as fit parameters, whereas the total opacity of the transition (τ) is regarded as a derived quantity. To this aim, the routine makes use of the physical line parameters and implements the computation of the rotational partition function for any given T ex . The set of the fit parameters also includes the line systemic velocity ( 0 ) in the local standard of rest, and the full-width half maximum of the lines (FWHM, ∆ ). Briefly, for a given input the total line opacity is first calculated through:

τ = ln 2 16π 3 c 3 A ul g u N ν 3 Q rot (T ex )∆ exp -E low k B T ex 1 -exp -hν 0 k B T ex , (1)
where ν 0 is the rest frequency of the (unsplit) rotational transition, c is the speed of light, Q rot (T ex ) is the rotational partition function at T ex , g u is the degeneracy of the upper level, A ul is the Einstein coefficient for spontaneous emission, E low is the energy of the lower state, h is the Planck constant, and k B is the Boltzmann constant. Then, assuming the local thermodynamic equilibrium (LTE) approximation and Gaussian profiles, the opacity of the i-th HFS components is computed as τ i = τ R i , where R i is the corresponding normalised relative intensity. The i-th component opacity at each velocity channel is then given by:

τ i ( ) = τ i exp -4 ln 2 -0,i ∆ 2 , (2)
where 0,i is the corresponding velocity position expressed as an offset with respect to the systemic velocity 0 . Finally, the continuum-subtracted antenna temperature at a given velocity channel is modeled as

T ant ( ) = [J ν (T ex ) -J ν (T CMB ) -T C ] 1 -e -i τ i ( ) , (3)
where J ν (T ) is the Rayleigh-Jeans radiation temperature, T C is the temperature of the continuum and T CMB is the cosmic background temperature. These data are compared to the continuumsubtracted observed antenna temperatures, T A -T C , in a nonlinear least-squares fashion to determine the best fit values of N, T ex , ∆ , and 0 . Whenever possible, all four parameters have been adjusted during the analysis. This sometimes produced strong correlation between them, specially when fitting weak features (e.g., those of ND 2 ), for which the relative intensities of the observed HFS components are poorly constrained because of the low signal-to-noise ratio (S/N). In these cases, suitable assumption for T ex and ∆ had to be used and the corresponding parameters were kept fixed in the least-squares fit. The single-sideband continuum temperature T C , inserted in the radiative transfer equality of Eq. ( 3), has been estimated using the linear formula given by Hily-Blant et al. ( 2010), T C /[K] = 1.10 ν -0.42, with ν expressed in THz.

The spectroscopic data for NH 2 , NHD, and ND 2 have been taken from the most recent laboratory studies (Martin-Drumel et al. 2014;Bizzocchi et al. 2020;Melosso et al. 2017), and the Q rot values used in Eq. ( 1) have been computed at any given temperature by cubic interpolation on a grid of finely spaced entries spanning the 2.7-19 K interval. These were obtained by direct summation on all rotational levels through the SPCAT spectral tool (Pickett 1991), with the ortho and para species of the symmetric isotopologues NH 2 and ND 2 treated separately. A brief explanation of the energy levels structure of NHD and ND 2 is given in Appendix A, while the complete list of HFS components used for the analysis of each transition is reported in Appendix C. The partition function values of both species computed at temperatures between 2.725 and 300 K are given in Appendix B.

Results

The line parameters obtained from the analysis of the observed transitions of NH 2 , NHD, and ND 2 are collected in Table 1. In the first two rows, the results obtained for the finestructure line of o-NH 2 detected previously by Hily-Blant et al.

(2010) are reported. Our fitted N and T ex values compare very well with those obtained from the CLASS HFS method employed in that study. Our weighted average of the column density is (5.4 ± 0.4) × 10 13 cm -2 , which is consistent within the uncertainty with the value (4.4 ± 0.7) × 10 13 cm -2 obtained previously. The excitation temperatures determined for the two components are 1σ coincident (8.44 ± 0.12 K the average), and are also close to the ones obtained by Hily-Blant et al. (2010) (8.5 K for the J = 3 2 -1 2 and 9.5 K for J = 1 2 -1 2 ). In the same spectral survey, we detect for the first time two groups of hyperfine transitions belonging to the singlydeuterated form of amidogen radical, NHD. Figure 2 shows the two fine-structure components of the fundamental N Ka,Kc = 1 1,1 ← 0 0,0 transition of NHD, both observed in absorption around 770 and 776 GHz in the 2b HIFI band. Since the J = 3 2 -1 2 fine-structure line is detected at ca. 10σ level, we could adjust all the four parameters in the model. The J = 1 2 -1 2 component is about twice weaker and its detection is only at ca. 5σ level. In this case, the noise prevented to adjust T ex and ∆ , which were then fixed at the values obtained for the strongest component of the fine-structure doublet. For NHD, we found an excitation temperature of 7.4 ± 0.1 K, slightly lower than that of NH 2 . This small difference is likely to be produced by the bigger Herschel beam size in the 2b band (∼ 27 ) compared to that in the 4a band (∼ 22 ). Indeed, a larger region is sampled for the NHD lines, and the contribution of the inner

Species

Transition Notes. Numbers in parentheses refer to 1σ uncertainties expressed in units of the last quoted digit. (a) Rest frequency corresponding to the hypothetically unsplit fine-structure transition. (b) Half-power beam width of the observations. (c) A conservative uncertainty of 0.25 MHz produced by the HIPE pipeline is summed in quadrature. (d) Derived through Eq. ( 1).

Frequency a HPBW b T C T ex LSR ∆ c N τ d N K a ,Kc , J [GHz] [arcsec] [K] [K] [km s -1 ] [km s -1 ] [10 13 cm -2 ] o-NH 2 1 1,1 , 3 2 -0 0,
-30 -20 -10 warmer gas is expected to be slightly smaller. The average value of the column density is N = (4.7 ± 0.7) × 10 13 cm -2 . Compared to the parent isotopologue, the doubly deuterated form of amidogen has a denser spectrum in the sub-millimeter region, thus both the N Ka, Kc = 1 1,1 ← 0 0,0 and N Ka, Kc = 2 1,2 ← 1 0,1 rotational lines fall within the coverage of the CHESS spec-tral survey. A clear detection has been obtained in the 2b band for both J = 5 2 -3 2 and J = 3 2 -1 2 fine-structure components of the N Ka, Kc = 2 1,2 ← 1 0,1 transitions of the para species. They are shown in Fig. 3. The N Ka, Kc = 1 1,1 ← 0 0,0 , J = 3 2 -1 2 of the ortho species, illustrated in Fig. 4, is located at 527 GHz. In this spectral region, covered by the HIFI 1a band, the continuum emission of the source is only 200 mK; hence, the resulting absorption feature is rather faint and the detection is only marginal, reaching at most 5σ. The other fine-structure component (intrinsically less intense) is just below the noise level and not detected. Given the weakness of all the three components observed for ND 2 , their analysis had to be carried out adopting suitable assumptions for T ex and ∆ , as these parameters could not be set free in the least-squares fits. Given the similarity of the beam size (∼ 27 ), for the p-ND 2 N Ka, Kc = 2 1,2 ← 1 0,1 lines we have assumed the same excitation temperature determined for NHD. A smaller value is expected to hold for the o-ND 2 N Ka, Kc = 1 1,1 ← 0 0,0 line, as the Herschel beam HPBW at this frequency is as broad as 40 , thus sampling a comparatively large volume of cold gas. For this transition we have fixed T ex = 4.5 K, a value determined for ND in the same source (Bacmann et al. 2010) through the analysis of its N = 1 -0 fine-structure doublet centered at 534 GHz. In fact, this value allows for a satisfactory fit of the absorption o-ND 2 feature, whereas for higher excitation temperatures (e.g., already at 5.5 K) the LTE model predicts the line to be observable in emission. Constraints for the FWHM line width ∆ have been derived from the values determined for NH 2 after deconvolution from the WBS spectrometer resolution (1.1 MHz). The column densities derived for p-ND 2 using two components are consistent within 2σ and the weighted average value is N = (0.241 ± 0.014) × 10 13 cm -2 . The observed ortho-to-para ratio is 2.7 ± 0.4. By considering (i) the energy difference of 10.77 cm -1 (E/k ≈ 15.5 K) between the lowest para and ortho rotational levels, and (ii) the spin statistical weights, this ratio corresponds to an equilibrium temperature of 9 ± 1 K. Due to the different beam sizes used to observe the N Ka, Kc = 1 1,1 ← 0 0,0 , and 2 1,2 ← 1 0,1 transitions, the p-ND 2 is preferentially sampled in the inner warmer region; therefore, this value is likely to be slightly overestimated.

0 10 20 30 40 v LSR / km s -1 -0.4 -0.3 -0.2 -0.1 0.0 0.1 T A -T C / K
For the same reason, the evaluation of the abundance ratios between amidogen isotopologues requires some careful considerations. From the averaged N values we obtained [NHD]/[NH 2 ] = 0.87 ± 0.14, hinting to a very high deuteration level for amidogen, lying in the 70-100% range. We consider this results as robust. The HPBW of the observations differs by only 25%, thus the source area sampled by the antenna is similar for both NHD and NH 2 lines. Furthermore, the energy difference between ortho and para species of NH 2 is 21.11 cm -1 (E/k B ≈ 30.4 K), and the contribution of the unobserved p-NH 2 is expected to be less than 5% at 10 K. Finally, since the used transitions have their lower level located at "zero" energy, the product of two exponential terms in Eq. ( 1) tends to 1 for small hν 0 /k B T ex , and the relation between N and τ depends only mildly on the temperature. Hence, the retrieved column density values are not much affected by the inaccuracies in the corresponding T ex determination.

The derivation of sound values of the isotopic ratios involving ND 2 is less straightforward, as various approaches can be adopted to estimate the overall ND 2 column density from the available observational results. By assuming that the two spin species are in thermal equilibrium, total N values of 0.70 × 10 13 cm -2 and 0.98 × 10 13 cm -2 can be obtained from the o-ND 2 and p-ND 2 results, respectively. A direct sum of the ortho and para column density with no assumptions gives instead 0.90 × 10 13 cm -2 . Based on the resulting 30% discrepancies, we conservatively quote N = (0.9 ± 0.3) × 10 13 cm -2 as the total column density of ND 2 . This yields [ND 2 ]/[NHD] = 0.19 ± 0.07 and [ND 2 ]/[NH 2 ] = 0.17 ± 0.06, that is, deuteration levels in the 10-30% range. The amidogen D/H ratio measured in the low-mass protostar IRAS16293 is very high and is comparable to the one derived for the related species imidogen et al. 1995). The ND 2 species is very abundant, yielding a remarkably high deuteration ratio, even higher than the one found for methanol ([CHD 2 OH]/[CH 3 OH]= 6%, Parise et al. 2004), a species for which the triply deuterated species was also detected in this source.

Chemical simulation

A steady-state model of cold and CO-depleted clouds was published by Roueff et al. (2005), but their predictions on deuterium enhancement, i.e., 10% for NHD and below 1% for ND 2 , are significantly underestimated with respect to our findings, which are 87% and 18% for NHD and ND 2 , respectively.

In an effort to reproduce the observed tendencies of the D/H and spin-state abundance ratios, we have carried out new chemical simulations of (deuterated) amidogen using a model of IRAS16293-2422. We followed the approach described in Brünken et al. (2014) and Harju et al. (2017), where more details can be found. Here, we provide only a brief description of the simulations. For the source model we adopted the one used in the above references, where the physical model from Crimier et al. (2010), which represents the protostellar core, was extended with a layer of gas representing the molecular cloud that the protostellar system lies in (with n(H 2 ) = 10 4 cm -3 and T dust = T gas = 10 K). The source model is illustrated in Fig. 5. We considered a two-stage model, in which the chemical evolution was first calculated over time t 1 in physical conditions corresponding to the core extension (stage 1), and then the protostellar core model including the low-density extension was run over time t 2 (stage 2) using the initial chemical abundances obtained from the stage 1. We explored a range of t 1 and t 2 values to search for the best fit to the observed column densities.

The simulation results depend on how the chemical reactions related to (deuterated) amidogen formation are assumed to proceed. Here we used our gas-grain chemical code (Sipilä et al. 2019b) and ran a parameter-space exploration, testing the effects of multilayer (three-phase) ice chemistry (Sipilä et al. 2016) and direct proton/deuteron hop (Sipilä et al. 2019a) as opposed to full scrambling in proton-donation reactions (Sipilä et al. 2015) on the various column densities. As an example of the results, Fig. 5 shows a comparison of the beam-convolved modelled and observed column densities and column density ratios using the multi-layer ice model and full scrambling for proton-donation reactions. The chemical model reproduces very well the column density of o-NH 2 and the ND 2 /NHD ratio, but fails to reproduce the NHD/NH 2 ratio as well as the NHD and (ortho and para) ND 2 column densities. Other parameter combinations reproduce instead very well the NHD and (ortho and para) ND 2 column densities, but strongly overpredict the o-NH 2 column density. None of the parameter combinations we explored lead to the observed level for the NHD/NH 2 ratio. The chemical model predicts that most of the NHD column originates near the edge of the protostellar core, where the temperature in the source model (∼14 K at the very edge, rising inward; cf. Fig. 5) inhibits NHD production with respect to NH 2 . Indeed, a NHD/NH 2 ratio of unity or above can be produced by the chemical model if the gas is sufficiently dense (n(H 2 ) 10 5 cm -3 ) and cold (T ∼ 10 K), suggesting that our two-phase approach to the core modelling is inadequate for the modelling of amidogen chemistry, and that a substantial amount of time (of the order of 10 6 yr) is required in the pre-stellar phase for the NHD/NH 2 ratio to increase to the observed level. We stress that in the present work, we did not carry out radiative transfer simulations of the amidogen absorption lines, and therefore the comparison of the modelled column densities against the observed ones does not take into account any excitation effects. We will present an in-depth comparison of the model versus the observations in an upcoming work, where we will also carry out line simulations.

Discussion and conclusions

This work presents the very first detection of both deuterated forms of amidogen radical in the ISM. From the observations, column densities and deuterium fractionation values have been derived allowing for a direct comparison with the predictions obtained from chemical networks. We have run a parameter-space exploration using our gas-grain chemical code, implementing a two-stage approach to the dynamical evolution of the source. The observed abundances of all NH 2 , NHD, and ND 2 isotopic variants could not be simultaneously reproduced with a single set of modelling parameters, thus indicating that the two-stage approach to the core modelling is inadequate for simulating amidogen isotopic chemistry. More insights on this topic would come from additional observations of deuterated amidogen isotopologues. They are, however, partially hampered by the atmospheric opacity. In fact, The most intense transition of NHD, located at 770 GHz, is not accessible with the ground-based facilities because of nearby strong water absorption lines. Hence, NHD observations can rely mainly on the atype N Ka, Kc = 1 0, 1 ← 0 0, 0 transition around 413 GHz. Luckily, spectral coverage at this frequency is supported by the ALMA 8 Band and the new nFLASH receiver of APEX. As for ND 2 , the transitions detected in this work fall in spectral windows accessible either with the incoming band 1 of the 4GREAT receiver on board SOFIA (490-635 GHz) or ALMA 10 Band (787-950 GHz).

The conclusions of the present work can be summarized as follows:

1. Thanks to the recently obtained laboratory data, two deuterated forms of amidogen radical, namely NHD and ND 2 , have been identified for the first time in the ISM. The analysis of hyperfine resolved transitions, observed in absorption against the continuum of the Class 0 protostar IRAS 16293-2422, allowed the accurate determination of column densities and excitation temperatures. 2. The deuterium fractionation of amidogen results to be 70-100% for NHD and 10-30% for ND 2 . Our derived values differ significantly, being larger by about one order of magnitude, from those predicted by earlier astrochemical models. 3. The low temperature and the high degree of deuteration observed suggest that amidogen, similarly to other nitrogencontaining species, is less depleted and subsists longer in the gas-phase, thus leading to high abundances of NHD and ND 2 . 4. The amidogen D/H and spin-state abundance ratios in IRAS 16293-2422 cannot be reproduced with our gas-grain chemical model unless the pre-stellar stage lasts long, of the order of one million years.

Fig. 1 .

 1 Fig. 1. Spectrum of the J = 3 2 -1 2 component of the N Ka,Kc = 1 1,1 -0 0,0 o-NH 2 transition around 952 GHz observed towards IRAS 16293-2422. The filled gray histograms show the observed spectra, while the blue trace is our HFS fit. The red sticks represent the position of the HFS components assuming LTE intensities.

Fig. 2 .

 2 Fig. 2. Spectra of the fine-components of the N Ka,Kc = 1 1,1 ← 0 0,0 NHD transition observed towards IRAS 16293-2422. Upper panel: J = 3 2 -1 2 component at ca. 770 GHz. Lower panel: J = 1 2 -1 2 component at ca. 776 GHz. The filled gray histograms show the observed spectra, while the blue trace is our HFS fit. The red sticks represent the position of the HFS components assuming LTE intensities.

Fig. 3 .

 3 Fig. 3. Spectra of the fine-structure components of the N Ka,Kc = 2 1,2 ← 1 0,1 p-ND 2 transition observed towards IRAS 16293-2422. Upper panel: J = 5 2 -3 2 component at ca. 785 GHz. Lower panel: J = 3 2 -1 2 component at ca. 787 GHz. The filled grey histograms show the observed spectra, while the blue trace is our HFS fit. The red sticks represent the position of the HFS components assuming LTE intensities.

Fig. 4 .

 4 Fig. 4. Spectrum of the J = 3 2 -1 2 component of the N Ka,Kc = 1 1,1 -0 0,0 o-ND 2 transition around 527 GHz observed towards IRAS 16293-2422. The filled gray histograms show the observed spectra, while the blue trace is our HFS fit. The red sticks represent the position of the HFS components assuming LTE intensities.

  ([ND]/[NH]= 30 -70%, Bacmann et al. 2010), while ammonia shows a lower level ([NH 2 D]/[NH 3 ]∼ 10%, van Dishoeck

Fig. 5 .

 5 Fig. 5. Left panel: Density (red) and temperature (blue) structure in the IRAS16293 source model in stages 1 and 2 (see Sect. 5 for details). Middle panel: Modelled beam-convolved molecular column densities (solid lines, labeled in the plot) as functions of time in stage 2. Dashed lines represent the observed values, omitting error bars for clarity. Right panel: Modelled D/H ratios (solid lines, labeled in the plot) as functions of time in stage 2. Here, the NH 2 and ND 2 abundances have been summed over the ortho and para forms. Dashed lines represent the observed values, for which the error bars have again been omitted.

Table 1 .

 1 Summary of the transitions observed for NH 2 , NHD, and ND 2 towards IRAS 16293-2422.

  TableC.5. List of HFS components used to reproduce the J = 3 2 ← 1 2 , N Ka, Kc = 2 1, 2 ← 1 0, 1 transition of ND 2 .

	1 2 , N Ka, Kc = 1 1, 1 ← 0 0, 0 transition of NHD. Table C.3. List of HFS components used to reproduce the J = 3 2 ← 1 2 , N Ka, Kc = 1 1, 1 ← 0 0, 0 transition of ND 2 .	Table C.2. List of HFS components used to reproduce the J = 1 2 ← 1 2 , N Ka, Kc = 1 1, 1 ← 0 0, 0 transition of NHD. Table C.4. List of HFS components used to reproduce the J = 5 2 ← 3 2 , N Ka, Kc = 2 1, 2 ← 1 0, 1 transition of ND 2 .
	Frequency Frequency Frequency	LGINT LGINT LGINT	Up. state Up. state Up. state		Lo. state Lo. state Lo. state		Frequency Frequency	LGINT LGINT	Up. state Up. state		Lo. state Lo. state
	(MHz) (MHz) (MHz)	(nm 2 MHz) F 1 (nm 2 MHz) F 1 (nm 2 MHz) F 1	F 2 F I D, tot F I D, tot F	F 1 F 1 F 1	F 2 F I D, tot I D, tot	F F	(MHz) (MHz)	(nm 2 MHz) F 1 (nm 2 MHz) F 1	F 2 F I D, tot F	F 1 F 1	F 2 F I D, tot	F
	770684.033 527133.0435 786826.0110	-3.5757 -4.2749 -3.4286	0.5 0.5 1.5	0 1 2	1 0.5 1.5 1.5 0.5 0.5	1	1 2	1	0.5 1.5	775951.7799 784901.6100	-3.3743 -3.3624	1.5 2.5	2 2	2 3.5 2.5 1.5	1	2	1	3.5
	770697.605 527137.3007 786828.1584	-3.4229 -4.2731 -3.2317	1.5 0.5 1.5	1 1 2	2 0.5 1.5 1.5 1.5 0.5	1	1 2	1	1.5 1.5	775954.2373 784901.8686	-3.0381 -3.3381	1.5 2.5	2 0	3 2.5 2.5 1.5	1	0	2	2.5
	770701.070 527138.6811 786832.1370	-3.0725 -4.2681 -3.2842	1.5 0.5 1.5	1 1 2	2 1.5 1.5 1.5 2.5 0.5	1	1 2	2	2.5 1.5	775983.7269 784903.2122	-3.2975 -3.2352	1.5 2.5	1 2	2 4.5 2.5 1.5	2	2	3	4.5
	770703.764 527145.1770 786832.2338	-3.5532 -4.1674 -3.0963	1.5 1.5 2.5	1 1 2	0 2.5 1.5 1.5 3.5 1.5	1	1 2	1	1.5 3.5	776004.0628 784906.4925	-3.1888 -3.4213	0.5 1.5	1 2	2 2.5 1.5 1.5	1	2	2	2.5
	770706.007 527147.6220 786834.3398	-3.2813 -3.9877 -2.9932	1.5 1.5 2.5	1 1 2	1 1.5 1.5 1.5 2.5 1.5	1	1 2	2	0.5 2.5	776004.6459 784906.9759	-3.3020 -3.2482	1.5 1.5	2 2	1 3.5 1.5 1.5	2	2	1	3.5
	770713.045 527151.0397 786834.9400	-3.3739 -3.8630 -3.0538	1.5 1.5 2.5	2 1 2	1 0.5 1.5 1.5 1.5 1.5	2	1 2	1	0.5 1.5	776009.0444 784907.7212	-3.1157 -3.3377	1.5 1.5	2 0	2 1.5 1.5 1.5	2	0	2	1.5
	770714.635 527151.8792 786835.0841	-3.2184 -3.5628 -3.2113	1.5 1.5 2.5	2 1 2	2 1.5 1.5 1.5 0.5 1.5	2	1 2	2	1.5 0.5	776009.3879 784912.0082	-3.3118 -3.4046	1.5 1.5	1 2	2 0.5 1.5 0.5	1	2	1	1.5
	770717.184 527155.0073 786836.6556	-3.0809 -3.2897 -3.1210	1.5 1.5 1.5	2 1 2	3 2.5 1.5 1.5 2.5 0.5	2	1 2	3	2.5 2.5	776009.5205 784913.9865	-3.5730 -3.4905	0.5 1.5	0 2	1 1.5 1.5 1.5	1	2	1	2.5
	770729.939 527155.2969 786836.7849	-3.1592 -3.8423 -3.0547	0.5 1.5 1.5	1 1 0	2 0.5 1.5 0.5 1.5 0.5	1	1 0	2	1.5 0.5	776012.9847 784920.0864	-3.3589 -2.8823	0.5 1.5	0 2	1 2.5 0.5 1.5	1	2	2	1.5
	770731.407 527161.7095 786836.9956	-2.6209 -3.8657 -3.2843	2.5 1.5 2.5	2 1 2	3 1.5 1.5 1.5 1.5 1.5	1	1 2	2	2.5 0.5	776013.1471 784924.5270	-3.2111 -2.4925	1.5 1.5	1 2	2 3.5 0.5 0.5	1	2	2	2.5
	770732.303 527175.2302 786836.9995	-3.3683 -3.3744 -2.6241	0.5 0.5 2.5	1 1 0	1 1.5 0.5 0.5 2.5 1.5	1	1 0	1	1.5 1.5	776013.9374 784925.9631	-3.5115 -2.8414	1.5 1.5	1 0	1 1.5 0.5 0.5	1	0	2	0.5
	770733.192 527179.4875 786838.1735	-3.5106 -3.2692 -2.9018	0.5 0.5 2.5	1 1 2	1 1.5 0.5 0.5 2.5 1.5	1	1 2	0	0.5 1.5	776016.0039 784927.5805	-2.8268 -2.9946	1.5 1.5	2 2	3 1.5 0.5 1.5	2	2	3	1.5
	770733.891 527182.9457 786838.3241	-2.9068 -2.9117 -2.3974	2.5 2.5 2.5	2 1 2	2 2.5 1.5 1.5 4.5 1.5	1	1 2	1	1.5 3.5	776022.8395 784929.8205	-3.5879 -2.9028	1.5 2.5	1 2	2 2.5 1.5 0.5	0	2	1	1.5
	770735.685 527183.5235 786838.9615	-3.4570 -3.1343 -2.6275	0.5 2.5 2.5	1 1 2	0 1.5 1.5 0.5 3.5 1.5	1	1 2	1	0.5 2.5	776023.6298 784929.8382	-3.3538 -2.6114	1.5 2.5	1 2	1 3.5 1.5 0.5	0	2	1	2.5
	770735.961 527183.6802 786842.7431	-3.2356 -3.4747 -2.7050	2.5 0.5 1.5	2 1 2	1 0.5 0.5 1.5 3.5 0.5	1	1 2	0	1.5 2.5	776042.6716 784931.7629	-3.3983 -2.3779	1.5 2.5	2 2	2 4.5 1.5 0.5	1	2	1	3.5
	770737.355 527185.6745 786864.1688 770737.593 527187.7807 786868.6874 770739.631 527187.9375 786872.8893 770743.306 527191.5565 786873.4617 770744.024 527192.7760 786877.9803 770745.755 527198.2587 786887.3254	-3.1949 -2.6845 -3.4603 -3.2680 -3.5532 -3.0107 -2.8049 -4.1230 -3.1436 -2.7096 -3.0134 -2.8638 -2.5392 -3.5051 -3.0524 -3.3952 -3.5357 -3.4508	2.5 2.5 0.5 2.5 2.5 0.5 0.5 0.5 0.5 2.5 1.5 0.5 2.5 2.5 0.5 0.5 1.5 1.5	2 1 2 2 1 2 1 1 0 3 1 2 3 1 2 1 1 2	2 3.5 1.5 1.5 1.5 0.5 1 1.5 1.5 1.5 1.5 0.5 2 0.5 0.5 0.5 0.5 0.5 2 1.5 2.5 0.5 2.5 0.5 3 1.5 2.5 1.5 2.5 0.5 1 0.5 1.5 0.5 2.5 1.5	1 1 0 2 2 0	1 2 1 2 1 0 1 2 1 2 1 2	2 1 1 1 2 1	2.5 1.5 1.5 2.5 0.5 0.5 1.5 1.5 2.5 2.5 1.5 3.5	5. ALMA represents the best candidate to perform future ob-776044.2348 -3.6053 1.5 2 1 0.5 1 0 784932.0249 -3.2870 2.5 2 1.5 1.5 2 0.5 servations of NHD and ND 2 , not only in IRAS 16293-2422, 776045.4241 -3.1291 1.5 2 3 0.5 1 2 784932.2158 -3.3227 1.5 2 0.5 0.5 2 1.5 but also in other sources where deuterated ammonia has been 776046.4309 -3.5307 1.5 2 2 0.5 1 2 784932.7588 -2.6164 2.5 0 2.5 1.5 0 1.5 detected. More observations would enable to make a step to-776052.6624 -3.5309 0.5 1 0 1.5 2 1 784934.3875 -3.0821 2.5 2 1.5 1.5 2 1.5 wards a thorough understanding of ammonia formation and 776057.7495 -3.1714 0.5 1 1 1.5 2 2 784934.6498 -3.0523 1.5 2 2.5 0.5 2 2.5 deuteration in the ISM. 776065.8293 -2.9268 0.5 1 2 1.5 2 3 784934.6762 -3.2289 2.5 2 0.5 1.5 2 0.5
	770746.127 527201.6764 786890.0744	-2.3691 -4.3765 -3.3125	2.5 1.5 1.5	3 1 2	4 0.5 0.5 1.5 1.5 1.5	2	1 2	3	1.5 2.5	776104.1715 784936.4341	-3.4331 -3.0168	0.5 2.5	0 2	1 2.5 1.5 0.5	1	2	2	2.5
	770746.604 527202.5160 786891.7607	-2.6141 -3.6385 -3.4671	1.5 1.5 1.5	2 1 2	3 1.5 0.5 0.5 0.5 1.5	1	1 2	2	0.5 1.5	784937.6522	-2.5509	3.5	2	3.5 2.5		2	2.5
	770748.263 527205.9337 786893.4129	-2.9477 -3.7120 -2.9898	1.5 1.5 1.5	2 1 2	2 0.5 0.5 0.5 3.5 1.5	1	1 2	1	0.5 3.5	784937.6570	-2.7232	3.5	2	2.5 2.5		2	1.5
	770748.379 786893.6763	-3.4022 -3.1369	2.5 1.5	3 0	2 1.5 1.5 1.5	2	0	2	1.5	784938.3512	-2.9047	3.5	2	1.5 2.5		2	0.5
	770751.746 786894.0531	-3.4616 -3.4158	1.5 1.5	2 2	1 2.5 1.5 0.5	1	2	1	2.5	784938.5498	-2.3902	3.5	2	4.5 2.5		2	3.5
	770751.990 786897.8867	-3.3905 -3.4000	2.5 1.5	3 2	3 2.5 1.5 1.5	2	2	3	1.5	784940.0965	-2.4157	3.5	0	3.5 2.5		0	2.5
	770752.022	-3.4349	1.5	2	2	0.5	1		2		784940.4444	-3.1259	2.5	2	3.5 1.5		2	3.5
	770752.634	-3.3344	1.5	2	1	0.5	1		0		784940.4531	-2.2396	3.5	2	5.5 2.5		2	4.5
	770762.836	-3.6026	1.5	1	2	1.5	2		3		784941.3117	-3.2175	3.5	2	1.5 2.5		2	1.5
	770778.684	-3.4439	0.5	0	1	0.5	1		2		784942.8337	-3.0461	3.5	2	2.5 2.5		2	2.5
	770792.256	-3.4614	1.5	1	2	0.5	1		2		784945.3078	-3.0176	3.5	2	3.5 2.5		2	3.5
	770793.173	-3.2079	2.5	2	3	1.5	2		3		784948.8334	-3.1526	3.5	2	4.5 2.5		2	4.5
	770801.949	-3.5909	1.5	1	2	0.5	0		1
	770806.886	-3.6006	1.5	1	1	0.5	0		1

Heterodyne Instrument for the Far-Infrared.

Chemical Herschel surveys of star forming regions.

ESA Herschel Science Archive http://archives.esac.esa. int/hsa/whsa.

See GILDAS home page at the URL: http://www.iram.fr/ IRAMFR/GILDAS.

Acknowledgements. This study was supported by Bologna University (RFO funds) and by MIUR (Project PRIN 2015: STARS in the CAOS, Grant Number 2015F59J3R). The work at SOLEIL was supported by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A.1: Singly-deuterated amidogen radical

The manifold of rotational energy levels is very complex in NHD and its spectrum is characterized by many interactions. First, NHD is a light asymmetric rotor far from the prolate limit (κ = -0.66) whose electric dipole moment lies in the ab principal symmetry plane, with components µ a = 0.67 D and µ b = 1.69 D (Brown et al. 1979). Because of the large values of the rotational constants, the room temperature spectrum is quite sparse and peaks in the far-infrared (FIR) region.

NHD is an open-shell molecule with a X 2 A ground electronic state. Due to the presence of an unpaired electron, the electronic spin (S = 1 /2) couples with the rotational angular momentum N, thus splitting each rotational level N K a , K c into two fine-structure sub-levels. The fine-structure levels have either J = N + 1 /2 or J = N -1 /2 quantum numbers, with the exception of the N = 0 level, for which only the J = 1 /2 level exists. The total angular momentum J further couples with each nuclear spin (I N and I D = 1, I H = 1 /2), giving rise to the socalled hyperfine structure (HFS). The angular momentum coupling scheme adopted takes into account the magnitude of the interaction (from strongest to weakest):

A more detailed accounts of the NHD spectroscopy can be found in Bizzocchi et al. (2020).

Appendix A.2: Doubly-deuterated amidogen radical

The higher symmetry of ND 2 with respect to NHD requires additional considerations in the description of the energy levels manifold. ND 2 belongs to the C 2v symmetry point group and a πrotation along the b axis exchanges the two identical deuterium nuclei. Since deuterium is a boson, the total wave-function of ND 2 must be symmetric upon the permutation of the two particles. Given that the electronic ground state is antisymmetric ( X 2 B 1 ) and the vibrational ground state is symmetric, deuterium nuclear spin functions (I D, tot = I D 1 + I D 2 ) can only combine with rotational levels of opposite symmetry. That is, symmetric rotational levels (K a + K c = even) combine with antisymmetric deuterium spin functions, with I D, tot = 1 (ortho states), while antisymmetric rotational levels (K a + K c = odd) combine with the spin functions I D, tot = 0, 2 (para states). Only transitions within the states are allowed by the permanent electric dipole moment of ND 2 (µ b = 1.82(5) D Brown et al. 1979). At low temperature, given the energy difference of 15.5 K between the ortho and para states, o-ND 2 and p-ND 2 should be considered as different species.

The network of fine and hyperfine interactions in ND 2 is analogous to that described earlier for NHD, while the coupling scheme of angular momentum changes as follows:

A detailed account of the spectroscopic features of ND 2 can be found in Melosso et al. (2017).

Appendix B: Partition function values

The spectroscopic parameters from Bizzocchi et al. (2020) and Melosso et al. (2017) have been input in the SPCAT subroutine (Pickett 1991) to evaluate the rotational partition function of NHD and ND 2 , respectively. All the energy levels up to J = 20 and K a = 15 have been considered. The temperature-dependence of Q rot has been computed analytically for both species at temperatures between 2.725 and 300 K. These values are given in Table B.1.

Appendix C: List of the used HFS components

Tables C.1-C.5 list the rest frequencies of the observed hyperfine components of NHD and ND 2 . The LTE intensities (computed at 300 K and labeled LGINT in the tables) and the HFS quantum numbers are also given. Spectral predictions are based on the spectroscopic studies reported in Melosso et al. (2017) and in Bizzocchi et al. (2020).