# Zero-inflated Smoothing Spline (ZISS) Models for Individual-level Single-
cell Temporal Data

Yifu Tang
Tsinghua University
<EMAIL_ADDRESS>Yi Zhang
University of Arizona
<EMAIL_ADDRESS>Yue Wang
University of Colorado
<EMAIL_ADDRESS>Jingyi Zhang*
Tsinghua University
<EMAIL_ADDRESS>Xiaoxiao Sun*
University of Arizona
<EMAIL_ADDRESS>
###### Abstract

Recent advancements in single-cell RNA-sequencing (scRNA-seq) have enhanced
our understanding of cell heterogeneity at a high resolution. With the ability
to sequence over 10,000 cells per hour, researchers can collect large scRNA-
seq datasets for different participants, offering an opportunity to study the
temporal progression of individual-level single-cell data. However, the
presence of excessive zeros, a common issue in scRNA-seq, significantly
impacts regression/association analysis, potentially leading to biased
estimates in downstream analysis. Addressing these challenges, we introduce
the Zero Inflated Smoothing Spline (ZISS) method, specifically designed to
model single-cell temporal data. The ZISS method encompasses two components
for modeling gene expression patterns over time and handling excessive zeros.
Our approach employs the smoothing spline ANOVA model, providing robust
estimates of mean functions and zero probabilities for irregularly observed
single-cell temporal data compared to existing methods in our simulation
studies and real data analysis.

Keywords: Functional data analysis; Individual-level; Single-cell RNA-seq;
Pseudotime

* To whom correspondence should be addressed.

##  1 Introduction

Recent advances in single-cell RNA-sequencing (scRNA-seq) techniques hold
great promise in improving our understanding of cell-to-cell heterogeneity and
cell differentiation at a high resolution [1, 2, 3]. With the aid of high-
throughput sequencing methods, researchers can process over 10,000 cells per
hour, enabling the collection of large scRNA-seq datasets from multiple
subjects with different disease statuses [4]. In the realm of single-cell data
analysis, trajectory inference methods stand out as pivotal tools. Their main
goal is to order individual cells along a trajectory, often referred to as
pseudotime, based on their gene expression profiles, resulting in the
generation of single-cell temporal data at the subject level. [5, 6, 7]. This
approach opens up fresh avenues for studying biological questions at the high-
resolution level, including association studies between the phenotypes of
participants and cellular processes over time. These investigations play a
crucial role in unraveling the complexities of disease progression and are
thus of utmost significance in scientific research [8, 9].

Due to the low capturing and sequencing efficiency in scRNA-seq, there is a
common dropout issue, referring to the presence of excessive zero counts in
the single-cell data [10, 11, 12]. This zero-inflation phenomenon notably
impacts data analysis, potentially skewing the inferred patterns of true
associations of interest. In addition, for the single-cell temporal data, zero
proportions are changing over pseudotime. Gene expression data are also
correlated across pseudotime. Disregarding these features of single-cell
temporal data during data analysis can lead to estimation bias, as
demonstrated in the following motivation example.

######  Motivation Example.

In this motivation example, we analyzed the scRNA-seq data obtained from a
subject with acute kidney failure [13]. Pseudotime estimation was obtained for
the data [14, 15]. Our objective is to estimate the smooth gene expression
curve over time from the zero-inflated noisy data. The resultant curve plays a
crucial role in downstream analysis of scRNA-seq data, such as differential
analysis. In Figure 1 a), we show the curves fitted for a subject with acute
kidney failure using the proposed method (ZISS) and the original smoothing
spline method applied to the data without zeros (NZSS) [16]. The original
smoothing spline method did not converge when applied to the original data due
to the excessive zeros. Hence, it is crucial to address the issue of excessive
zeros in modeling single-cell temporal data. In Figure 1 b), we present the
data distribution across pseudotime for the same gene discussed in a). This
illustration shows that the proportions of zeros vary throughout the
pseudotime, and the number of observations also differs at each pseudotime
point. Moreover, the proportions of zeros exhibit correlation over pseudotime.
These factors must be taken into account when estimating the probabilities of
zeros.

Figure 1: The motivation example. a). The black dots represent gene expression
levels along pseudotime for different cells from a gene (Ensembl ID:
ENSG00000137673). The red curve represents the estimated expression profiles
using the proposed method (ZISS), while the blue curve represents the
estimated curve obtained by applying the original smoothing spline method to
the data without zeros (NZSS). b). The violin plots depict gene expression
across pseudotime, which has been divided into four distinct intervals for
analysis.

In analyzing single-cell temporal data, the methods developed must tackle two
significant challenges, as illustrated in the motivation example. First,
addressing the issue of excessive zeros is imperative to avoid numerical
errors and estimation bias. Second, considering the correlations between cells
and their corresponding gene expressions over time, it becomes essential to
incorporate smoothness assumptions when estimating the temporal patterns of
gene expression levels from zero-inflated data. To tackle the first challenge,
a commonly employed approach involves constructing a two-component mixture
model. This model specifies the distribution of a random variable by using a
probabilistic mixture of zero and a regular component, the distribution of
which typically belongs to the exponential family. Several approaches have
been developed for the analysis of non-temporal scRNA-seq data [17, 18, 19].
These methods effectively tackle the challenges posed by excessive zeros using
zero-inflated Poisson and/or negative binomial models. However, these methods
are only applicable to regular scRNA-seq data and are not suitable for the
analysis of single-cell temporal data. To surmount the second challenge,
semiparametric models were considered to analyze zero-inflated temporal data
[20, 21]. These models involve modeling the non-zero-inflated response using
flexible semiparametric predictors, while the probability of zero-inflation is
constrained to an unknown constant. However, in practical applications, the
assumption that the zero-inflation probability remains constant over time may
not hold. To address this, zero-inflated generalized additive models (ZIGAMs)
were developed to extend the models to more general cases [22]. ZIGAM employs
smooth functions to model the mean pattern of the non-zero-inflated response
along with time and the non-zero-inflation probability, then estimates the two
smooth functions through two GAMs. Building on ZIGAM, the Constraint ZIGAM, or
the COZIGAM, [23, 24, 25] was introduced to establish a connection between the
mean function of the response and the non-zero-inflation function in cases
where both processes may be influenced by some common factors. To estimate the
unknown link functions, ZIGAM and COZIGAM utilize the regression spline
method. However, these methods may encounter challenges when pseudotime points
are irregularly spaced or records are unevenly distributed along pseudotime,
as the regression spline estimation is highly sensitive to the choice of knots
and can prove to be unreliable in the presence of varying levels of sparsity.
Therefore, novel approaches are needed to analyze the single-cell temporal
data. Additionally, the single-cell temporal data exhibits varying numbers of
observations over time, posing numerical challenges for estimation using
regression spline methods. Our approach, utilizing smoothing spline ANOVA
models, effectively addresses this issue.

In this paper, we introduce the ZISS method for the analysis of single-cell
temporal data. Our approach is specifically designed to address the challenges
posed by excessive zeros in scRNA-seq data, stemming from both biological and
technical factors. While biological zeros represent genuine absence of gene
expression, technical zeros arise due to limited genetic materials in
individual cells. Our proposed model comprises two components: one for
modeling the underlying gene expression patterns over time and another for
accommodating the excessive zeros. We utilize the smoothing spline ANOVA model
to estimate the link functions, offering greater robustness across various
sparsity levels by mitigating the dependence on knot selection. Furthermore,
our method yields more precise estimates of link functions and zero
probabilities compared to existing approaches.

The remainder of the article is structured as follows. Section 2 presents the
model setup for the proposed ZISS model, along with an overview of the model
estimation process. In Sections 3 and 4, we present the results of simulation
studies and real data analysis, respectively. Section 5 provides additional
remarks and conclusions. Notations and derivation are provided in the Appendix
for reference.

##  2 Model Formulation and Estimation

In this section, we present a novel Zero Inflated Smoothing Spline (ZISS)
method and provide estimation details of this approach. A two-component
mixture model will be employed to define the data distribution. One component
will be used to model the true signals based on the smoothing spline ANOVA
models, whereas the other component will be utilized for modeling the
excessive zeros. To adaptively manage the fluctuating zero probabilities over
time, we suggest the utilization of B-spline basis functions to model the
excessive zero probability functions. The estimation of the model parameters
will be conducted using the Expectation–Maximization (EM) algorithm. The
Newton-Raphson algorithm is implemented to estimate the zero probability
functions.

###  2.1 Zero Inflated Smoothing Spline (ZISS)

As demonstrated in the motivation example (see Figure 1), the temporal
dependence should be considered not only in the underlying non-zero inflated
gene expression levels, but also in the probability of excessive zeros when
analyzing single-cell temporal data. To formulate the temporal dependence, we
denote the underlying non-zero inflated gene expression level as $\mu(t)$, and
the probability of zero-inflation, which is also referred to as the dropout
rate, as $p(t)$. For each gene, we further assume the gene expression data are
observed at $N$ pseudotime points, $t_{min}<t_{1}<t_{2}<\cdots<t_{N}<t_{max}$.
At each pseudotime point $t_{i}$, we observe $M_{i}$ cells, which are always
recorded as the counting responses $y_{i,1},\cdots,y_{i,M_{i}}\in\mathbb{N}$,
where $\mathbb{N}$ denotes the set of non-negative integers. We then consider
the zero-inflated model as follows,

$y_{i,j}\sim\begin{cases}\textbf{Poisson}(\mu(t_{i})),\quad\textit{if}~{}g_{i,j}=1,\\\
0,\quad\textit{if}~{}g_{i,j}=0.\end{cases}$ (1)

where $g_{i,j}\sim Ber(p(t_{i}))$ is a Bernoulli variable indicating zero-
inflation. The zero probability function $p(t_{i})$ is used to model the
excessive zeros. When $g_{i,j}=1$, $y_{i,j}$ follows a Poisson distribution
with a mean of $\mu(t_{i})$ at pseudotime $t_{i}$, and is zero otherwise.

To address the irregularly observed single-cell temporal data, we propose
utilizing the smoothing spline ANOVA model [16] to estimate the mean function
$\mu$. This model is knot-free, making it suitable for irregularly observed
data. In particular, we consider the mean function in a reproducing kernel
Hilbert space, $\mathcal{H}=\left\\{\mu:J(\mu)<\infty\right\\},$ where
$J(\mu)$ is a quadratic functional [16, Section 2.4.5], measuring the
roughness of functions. For example, for a twice differentiable function
$\eta:[0,1]\rightarrow\mathbb{R}$, $J(\cdot)$ can be chosen as
$J(\eta)=\int_{0}^{1}\left(\frac{d^{2}\eta}{dx^{2}}\right)^{2}dx.$ We
decompose the space $\mathcal{H}$ into the null (kernel) space of $J(\cdot)$
and its orthogonal complement $\mathcal{H}_{J}$, i.e.
$\mathcal{H}=\mathcal{N}_{J}\bigoplus\mathcal{H}_{J}$, where the null space of
$J$ is defined as $\mathcal{N}_{J}=\left\\{\mu:J(\mu)=0\right\\}.$ Based on
the Representer Theorem [16], the mean function $\mu(t)$ has the following
form,

$\mu(t)=\sum_{\nu=1}^{m}d_{\nu}\phi_{\nu}(t)+\sum_{i=1}^{n}c_{i}R_{J}(t_{i},t)=\boldsymbol{\phi}^{\top}\boldsymbol{d}+\boldsymbol{\xi}^{\top}\boldsymbol{c},$
(2)

where $\left\\{\phi_{\nu}\right\\}_{\nu=1}^{m}$ are basis functions in
$\mathcal{N}_{J}$, $R_{J}(\cdot,\cdot)$ is the kernel function of
$\mathcal{H}_{J}$, $\boldsymbol{\xi}$ and $\boldsymbol{\phi}$ are the vectors
of the corresponding basis functions with coefficient vectors $\boldsymbol{c}$
and $\boldsymbol{d}$. For the Gaussian-type responses, the penalized least
squares functional in a reproducing kernel Hilbert space $\mathcal{H}$ can be
written as $L(\mu)+(\lambda/2)J(\mu)$, where $L(\mu)$ measures the goodness-
of-fit between the responses and the function $\mu$ (e.g., least squares) and
$\lambda$ is the smoothing parameter to control the smoothness of the function
estimate. We then plug the equation (2) into the penalized least squares
functional to obtain the estimates of coefficients. More details about the
estimation for Gaussian-type responses can be found in [16, Chapter 3]. In
line with the Gaussian cases, we will apply the penalized likelihood method to
estimate the mean function $\mu$ in the Poisson processes.

In the motivation example, the zero probabilities exhibit temporal dependency.
Therefore, it is reasonable to assume that the probability is a smooth
function. We assume that $p(t)$ is a piece-wise smooth function, with a
summation equal to one over the temporal domain. Specifically, we assume that

$p(t)=\frac{1}{1+\exp(\sum_{i=1}^{m}\alpha_{i}b_{i}(t))},$

where $b_{i}(t)$ are the B-spline basis functions over the pseudotime interval
and $\alpha_{1},\cdots,\alpha_{m}\in\mathbb{R}$ are the corresponding
coefficients of these basis functions. This nonparametric form offers a
flexible and smooth estimation of zero probabilities over time. In addition,
it simplify the overall estimation process.

###  2.2 Estimation

We utilize the maximum penalized likelihood method to estimate the mean
function $\mu(t)$ and the probability function $p(t)$. In particular, given
observations $t_{i},y_{i,j}$, the latent variables $g_{i,j}$, and the
roughness functional $J(\cdot)$, we minimize the following penalized negative
log-likelihood functional to obtain the estimates of $\mu(t)$ and $p(t)$,

$\displaystyle l(p(\cdot),\mu(\cdot);g_{i,j},y_{i,j})=$
$\displaystyle-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(g_{i,j}\log
p(t_{i})+(1-g_{i,j})\log(1-p(t_{i}))\right)$ (3)
$\displaystyle-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(g_{i,j}y_{i,j}\log\mu(t_{i})-g_{i,j}\mu(t_{i})\right)+\frac{\lambda}{2}J(\mu)$
$\displaystyle+\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}g_{i,j}\log(y_{i,j}!),$

where the smoothing parameter $\lambda$ controls the trade-off between the
goodness-of-fit of the model and the roughness of $\mu(t)$. As the latent
variables are not observed, we utilize the Expectation-Maximization (EM)
algorithm to perform estimation [26]. In the Expectation (E) step, we compute
the conditional expectation of the latent variable given the data and the
current estimates of the mean and probability functions. In the Maximization
(M) step, we maximize the penalized likelihood functional to obtain the
estimates for the mean and probability functions. In the estimation of
$\mu(t)$ and $p(t)$, the last term of the penalized negative likelihood, as
shown in (3), is not involved. Therefore, it is excluded during the EM
process. As a result, the M step is further divided into two sub-problems,
both optimized simultaneously. For the first problem, we employ the Newton-
Raphson method to obtain the estimate of the probability function. For the
second problem, we utilize the penalized likelihood method to obtain estimates
of the mean function. In the following, we will present the EM steps in
estimation. The derivation details can be found in Appendix.

#### Expectation (E) step.

Denoting the fitted mean and probability functions as $\hat{\mu}_{k}$ and
$\hat{p}_{k}$ at the $k$-th step (If $k=0$, $\hat{\mu}_{k}$ and $\hat{p}_{k}$
represent the initial values), we calculate the conditional expectation of the
latent variable given the current estimates and observed data, as
$\mathbb{E}\left[g_{i,j}|y_{i,j},\hat{\mu}_{k},\hat{p}_{k}\right]$ via the
following proposition.

######  Proposition 1 (Computing conditional expectation in the E step).

Given the observed data $\left\\{t_{i},y_{i,j}\right\\}_{1\leq i\leq N,1\leq
j\leq M_{i}}$ and the fitted curves $\hat{\mu}_{k}$ and $\hat{p}_{k}$ at the
$k$-th ($k\geq 0$) step, we have the conditional expectation of latent
variable $g_{i,j}$,

$\mathbb{E}\left[g_{i,j}|y_{i,j}=0,\hat{\mu}_{k},\hat{p}_{k}\right]=\frac{e^{-\hat{\mu}_{k}(t_{i})}\hat{p}_{k}(t_{i})}{e^{-\hat{\mu}_{k}(t_{i})}\hat{p}_{k}(t_{i})+1-\hat{p}_{k}(t_{i})},$
(4)

and

$\mathbb{E}\left[g_{i,j}|y_{i,j}>0,\hat{\mu}_{k},\hat{p}_{k}\right]=1.$

Denoting the obtained conditional expectation by
$q_{i,j}=\mathbb{E}\left[g_{i,j}|y_{i,j},\hat{\mu}_{k},\hat{p}_{k}\right]$, we
use the conditional expectation to replace the latent variables in the log-
likelihood functional $l$ in the M step.

#### Maximization (M) step.

In the M step, we maximize the penalized log-likelihood functional given the
conditional expectation in the previous E step. Equivalently, we minimize the
negative penalized log-likelihood functional (3) via the following
proposition.

######  Proposition 2 (Estimating $\mu$ and $p$ in the M Step).

Given
$q_{i,j}=\mathbb{E}\left[g_{i,j}|y_{i,j},\hat{\mu}_{k},\hat{p}_{k}\right]$
obtained in the E Step, we perform the following two optimization problems in
the M step at the same time,

$(P_{1})~{}~{}\min_{p\in\mathcal{P}}-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(q_{i,j}\log
p(t_{i})+(1-q_{i,j})\log(1-p(t_{i}))\right),$

and

$(P_{2})~{}~{}\min_{\mu\in\mathcal{H}}-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(q_{i,j}y_{i,j}\log\mu(t_{i})-q_{i,j}\mu(t_{i})\right)+\frac{\lambda}{2}J(\mu),$

where for given B-spline functions $b_{1}(t),\cdots,b_{m}(t)$ defined in the
pseudo-time interval,

$\mathcal{P}=\left\\{\frac{1}{1+\exp(\sum_{l=1}^{m}\alpha_{l}b_{l}(t))}:\alpha_{1},\cdots,\alpha_{m}\in\mathbb{R}\right\\},$

and $\mathcal{H}$ denotes the reproducing kernel Hilbert space wherein we
confine the estimates of $\mu$.

For the zero probability function, we have
$\log\frac{p(t)}{1-p(t)}=-\sum_{i=1}^{m}\alpha_{i}b_{i}(t)$. Thus in Problem
$(P_{1})$, the equation takes the form,

$q_{i,j}\log
p(t)+(1-q_{i,j})\log(1-p(t))=-q_{i,j}\sum_{i=1}^{m}\alpha_{i}b_{i}(t)+\log(1-p(t)).$

Algorithm 1 Zero inflated smoothing spline (ZISS)

1: Input: Smoothing parameter $\lambda$, single-cell data
$\left\\{(t_{i},y_{i,j})\right\\}_{1\leq i\leq N,1\leq j\leq M_{i}}$, initial
values $\hat{p}^{(0)},\hat{\mu}^{(0)}$, converge criterion (threshold gap
$\epsilon$)

2: for $m$ in $1:$max-iteration step do

3: (E Step) Compute the conditional expectation of $g_{i,n}$ to obtain
$q_{i,j}^{(m)}$ from equation (4).

4: (M Step) Minimize the objective functional to obtain $\hat{p}^{(m+1)}$ and
$\hat{\mu}^{(m+1)}$ in $(P_{1})$ and $(P_{2})$

5: if converge criterion then

6: Break

7: end if

8: end for

9: if GCV in the final step then

10: Use the smoothing parameter $\lambda_{0}$ selected by GCV to fit final
$\hat{\mu}$ and $\hat{p}$

11: end if

12: return $\hat{p},\hat{\mu}$

This formulation streamlines the minimization process. The second-order
derivative of the first term with respect to $\alpha$ becomes zero, markedly
reducing computational complexity in Newton’s iteration. Moreover, it is
essential that the probability function $p(t)$ remains within the [0,1] range.
The chosen form of $p$ naturally ensures this constraint. In contrast, using a
traditional B-spline function representation, such as
$p(t)=\sum_{i}\alpha^{\prime}b_{i}(t)$, could result in values exceeding 1.
For each iteration, assuming the total number of samples is $n=N*M_{i}$, the
theoretically optimal convergence rate for $\lambda$ is $O(n^{-\frac{2}{9}})$
[16, Section 9]. In our algorithm, we initially set
$\lambda=10n^{-\frac{2}{9}}$ for the smoothing parameter. The smoothing
parameter $\lambda$ is chosen by the generalized cross-validation (GCV) in the
final step [27, 28]. The iterative process is monitored for convergence by
comparing two successive estimations, $\hat{\mu}_{k}$ and $\hat{\mu}_{k+1}$.
Convergence is assumed when the norm of their difference satisfies the
condition $\|\hat{\mu}_{k+1}-\hat{\mu}_{k}\|\leq\epsilon,$ where $\epsilon$ is
a predefined threshold and $\|\cdot\|$ is the $l^{2}$-norm. The proposed
method is detailed in Algorithm 1.

##  3 Simulation Studies

This simulation study is designed to demonstrate the finite-sample performance
of the proposed ZISS algorithm, specifically in terms of mean curve fitting
and probability curve estimation. We compared ZISS with the following existing
methods: (1) the direct smoothing spline (DSS) method, which fits smoothing
spline ANOVA to all data points including zeros, (2) the non-zero fit (NZSS)
method, which fits smoothing spline ANOVA to all non-zero data points, and (3)
the ZIGAM method [22], fitting all data points.

###  3.1 Simulation Settings

We designed two simulation settings over the pseudotime interval $[0,1]$ to
compare different methods. Our first experiment is designed to ensure that
$p_{1}$ and $\mu_{1}$ are symmetric with respect to $t=0.5$, satisfying the
condition which $p$ and $\mu$ are related, while the second experiment is
designed with uncorrelated $p_{2}$ and $\mu_{2}$. In the second scenario, the
data exhibit a higher incidence of zeros, making it challenging to accurately
estimate the mean and zero probability functions.

For the first setting, the mean function is set as

$\mu_{1}(t)=2\sin(9t)+2.5,$

and the zero probability function is chosen as

$p_{1}(t)=\frac{1}{1+e^{-0.5(t-0.5)^{2}+1}},$

where the pseudotime points $t_{i},i=1,2,\cdots,N=41$ are uniformly
distributed within the time interval $[0,1]$. In this setting, $p_{1}$ and
$\mu_{1}$ are both symmetric with respect to the central node $0.5$. There is
a link function [22, 24] $g$ to show the association between $\mu$ and $p$,
such that $g(p)=\mu$,

$\mu=2\sin\left(9\sqrt{2-2\log\left(\frac{1}{p}-1\right)}+4.5\right)+2.5.$

At each pseudotime point $t_{i}$, we independently generated $M_{i}=M=80$
samples, $y_{1},y_{2},\cdots,y_{M}$ (corresponding to the respective
pseudotime), using the zero-inflated Poisson model proposed in (1).

In the second experiment setting, we generated the data with the mean function
$\mu(\cdot)$ according to

$\mu_{2}(t)=\frac{8}{\sqrt{2\pi}}e^{-10(t-0.2)^{2}}+\frac{6}{\sqrt{2\pi}}e^{-100(t-0.7)^{2}},$

and the true zero probability function

$p_{2}(t)=\frac{1}{4}\sin(6t)+0.5.$

In the second scenario, the mean function values within a specific region
approximate zero, making them hard to distinguish from the true zeros in
Poisson-distributed data. This presents a challenge for zero-inflated models
that lack a flexible framework for the estimation of zero probability
functions. Furthermore, the mean and zero probability functions are
uncorrelated in this setting. We summarize the two settings in Table 1.

Table 1: Summary of the two simulation settings Curves | Setting 1 | Setting 2
---|---|---
$\mu(t)$ | $2\sin(9t)+2.5$ | $\frac{8}{\sqrt{2\pi}}e^{-10(t-0.2)^{2}}+\frac{6}{\sqrt{2\pi}}e^{-100(t-0.7)^{2}}$
$p(t)$ | $\frac{1}{1+e^{-0.5(t-0.5)^{2}+1}}$ | $\frac{1}{4}\sin(6t)+0.5.$
Correlation | Link function exists | Link function does not exist
Zero Inflation | Low | High

Figure 2: Comparison of the mean curves estimated by different methods. In the
figure, the $+$ symbol represents the observed data at each pseudotime point.
Due to the possibility of multiple observations having identical values at a
given pseudotime point, they might overlap in the figure. The black line
illustrates the true mean curve $\mu_{i}$, while gray, orange, blue, and red
lines represent the fitting results of DSS, NZSS, ZIGAM and ZISS,
respectively. Below, the comparison between the estimated probability curve
$\hat{p}$ and the true probability function $p$ is displayed.

To illustrate the numerical differences between these two simulation settings,
we applied the methods for comparison to the simulated data from a single
replication. In Figure 2, panels (a) and (c) depict the fitting results of
Setting 1, while panels (b) and (d) present those of Setting 2. In Setting 1,
ZISS and ZIGAM demonstrate comparable performances with accurate fitting
patterns, attributed to the symmetry of true response curve $\mu$ and
probability curve $p$. Conversely, DSS fitting underperforms, consistently
falling below the true curve due to the inflated number of zeros in this
setting. The NZSS method is effective when $\mu(t)$ is significantly different
from 0 but struggles when $\mu(t)\approx 0$. In Setting 2, the ZISS method
offers the only estimate that aligns well with the true mean curve. ZIGAM
appears over-smoothed, and DSS fitting remains below the true value, similar
to its performance in Setting 1. The NZSS curve tends to overestimate the true
curve. This setting’s distinctive outcomes arise for several reasons: the zero
inflation is predominantly due to zeros from two processes (i.e., Possion and
zero processes) are not easily separable.

To demonstrate the superiority of the ZISS method in both simulation settings,
we conducted and reported the results from one hundred replicates. In each
replicate, we generated new data according to the same mean and probability
functions, and then fitted the data using the methods for comparison. We
calculated the mean square error (MSE) between the estimated $\hat{\mu}$ and
the true $\mu$ at each replicate and displayed the distributions of the MSEs
for each method with their corresponding standard deviation (Std). The results
of these repeated experiments are shown in Table 2. The results affirm the
excellent performance of the ZISS method, characterized by the lowest MSEs
coupled with the smallest Stds.

Methods | ZISS (Proposed) | NZSS | ZIGAM | DSS
---|---|---|---|---
MSE(Setting 1) | 0.033 | 0.185 | 0.033 | 5.404
Std (Setting 1) | 0.012 | 0.033 | 0.012 | 0.165
MSE(Setting 2) | 0.169 | 0.550 | 0.368 | 0.971
Std (Setting 2) | 0.009 | 0.0198 | 0.051 | 0.052
Table 2: Comparison of mean square error (MSE) and standard deviation (Std)
across different methods in 100 repeated experiments

###  3.2 Study on the Impact of Zero-inflation

In single-cell data, dropout (manifested as zero-inflation) and over-
dispersion are two prevalent phenomena. In this section, we explore the impact
of zero-inflation and over-dispersion, focusing on how various levels of zero-
inflation and over-dispersion affect the performance of the proposed
algorithm. To examine the increase in zero-inflation and over-dispersion, we
utilized the simulation study as outlined in Ridout et al. (2001) [29]. Such a
study serves as an effective method for analyzing the consequences of
gradually diminishing signal intensity, which corresponds to an increase in
zero-inflation. For the mean function $\mu$, we generated data that follows a
negative binomial distribution with a mean of $\mu$ and a variance of
$\mu(1+\sigma^{2}\mu)$. Here $\sigma^{2}=a$ is designated as the over-
dispersion parameter. For each scenario and each over-dispersion parameter, we
conducted 100 repetitions of the experiment and recorded the mean square error
(MSE) for each method. To study the effect of reduced zero-inflation, we
assessed the performance of the methods for comparison as the mean curve $\mu$
increasingly deviates from 0. For each specified up-moving parameter $h$, we
conducted the experiments under both simulation settings, using data generated
from $\mu+h$. These experiments were carried out with a progressively
increasing average value of the mean curve.

Figure 3: Comparison of various methods across different over-dispersion and
up-moving parameters in Setting 1 (panels a and c) and Setting 2 (panels b and
d).

Figure 3 summarizes the results of over-dispersion and up-moving experiments
for the two settings outlined in Section 3.1. Panels (a) and (b) illustrate
the over-dispersion cases for settings $\mu_{1}$ and $\mu_{2}$, respectively.
In addition, panels (c) and (d) show the MSE results of the methods for
comparison for the up-moving experiments. For Setting 1, the results indicate
that the ZIGAM and the ZISS methods perform similarly, whether the over-
dispersion parameter or $h$ increases. Conversely, in Setting 2, the ZISS
method consistently demonstrates superior performance in both over-dispersion
and up-moving scenarios. This suggests that the proposed method is
particularly effective in detecting and handling excessive zeros and over-
dispersion within the data.

Additionally, the boxplots reveal some notable trends. The performance of
NZSS, a method that removes all zeros, tends to improve as the signal strength
increases. However, this trend does not apply to the DSS method. Several
factors contribute to these observations. Firstly, when the mean curve $\mu$
moves up, the separation between the zero and Poisson processes becomes more
distinct. Consequently, in both settings, the performance of NZSS improves
with an increase in the parameter $h$. Secondly, for the DSS fitting, when the
mean function $\mu$ is modified to $\mu+h$, DSS directly fits the mean of
$\mu+h$ and 0. This results in a larger discrepancy between the DSS fitted
mean and the true mean curves as $h$ increases, leading to a rise in the MSE
for DSS fitting. Lastly, the methods ZIGAM, NZSS, and ZISS tend to deteriorate
as the over-dispersion parameter increases. This decline in performance can be
attributed to the fact that increased over-dispersion leads to a larger
variance in the data, causing more signals to be mistakenly identified as part
of the zero process in the zero-inflated model.

##  4 Real Data Analysis

In this study, we conducted an analysis of scRNA-seq data obtained from the
human kidney cell atlas. The researchers in this kidney disease study employed
single-cell assays on diseased kidneys, providing valuable insights into the
mechanisms underlying the increased risk of chronic kidney disease and kidney
failure resulting from acute kidney injury events [13]. We obtained the raw
count matrix of integrated single-nucleus and scRNA-seq data of the adult
human kidney, which was generated using the 10x Genomics platform. The raw
data consists of 33,920 genes and 304,652 cells. To focus on the epithelial
cell of the proximal tubule in scRNA-seq data, the data is reduced to 33,920
genes and 24,359 cells. To ensure data quality, cells with fewer than 500
reads and more than 50% of reads mapping to the mitochondrial genome were
identified as low quality and excluded from further analysis. Our focus in
this study is on patients with acute kidney failure (n=12) and chronic kidney
disease (n=15), resulting in a total of 19,507 cells from 27 patients being
used for data analysis. Additionally, we utilized 122 kidney-disease related
genes specifically identified in this kidney study[13]. Considering that the
data was collected from multiple laboratories, we applied the Harmony method
to address any batch effects and ensure accurate analysis [30]. We applied the
cluster-based minimum spanning tree (cMST) method to estimate the pseudotime
for each cell [14, 15]. Once the pseudotime was acquired, we partitioned the
entire time interval (0-20,000) into 150 sub-intervals. Thus, cells were
assigned to one of these sub-intervals based on their estimated pseudotime. In
total, we obtained 150 pseudotime points for each gene of a participant. In
this data analysis, we conducted a comparison between the proposed method and
ZIGAM under various conditions. Our findings demonstrate that the proposed
method consistently performs well across various levels of sparsity (i.e.,
proportion of zeros) and different numbers of observations.

###  4.1 The performance of the proposed method remains consistent across
various levels of sparsity.

In this kidney study, a key research objective was to investigate differential
gene expression patterns between participants with acute kidney failure and
chronic kidney disease. To address this question, we selected two important
genes (Ensembl IDs: ENSG00000137673 and ENSG00000118785) associated with
chronic kidney disease and applied both the proposed method and ZIGAM to the
data. The gene expression levels in the participant with acute kidney failure
show a high expression across the pseudotime from 10,000 to 20,000. In
contrast, the gene expression levels in the participant with chronic kidney
disease exhibit a high expression throughout the pseudotime from 0 to 20,000.
In Figure 4 (a), the gene expression data show high sparsity levels across the
pseudotime range of 0 to 5,000. The proposed method demonstrates superior
performance, capturing the underlying pattern accurately, while the ZIGAM
method exhibits an undersmoothed fit. When the sparsity levels are low (e.g.,
pseudotime from 10,000 to 20,000), both methods perform comparably in
capturing the gene expression dynamics.

Figure 4: The fitted curves for two subjects with acute kidney failure and
chronic kidney disease using two different methods: ZISS (a) and ZIGAM (b).
The curves represent the expression profiles of genes identified by their
Ensembl IDs, with ENSG00000137673 used for panel (a) and ENSG00000118785 used
for panel (b).

###  4.2 The performance of the proposed method remains consistent for
different numbers of observations.

We utilized the IDEAS method [31], a differential gene expression detection
approach for individual-level single-cell data, to analyze the dataset,
resulting in the identification of 13 differentially expressed genes out of
the initial set of 122 genes. Among these genes, we selected the one with the
highest total read count for further analysis. In Figure 5, we display the
fitted curves generated by the proposed method and ZIGAM for this particular
gene in four participants. As shown in Figure 5 (a) and (b), when a larger
amount of data is available, both the proposed method and ZIGAM exhibit
satisfactory performance. However, when the number of observations is limited,
the proposed method outperforms ZIGAM, which tends to suffer from issues
related to under-smoothing, see Figure 5 (c) and (d). The issue may stem from
selecting an inappropriate smoothing parameter using GCV in ZIGAM due to the
limited sample size.

Figure 5: The proposed method is robust to various data levels in the data of
four participants.

##  5 Conclusion

In conclusion, our comprehensive analysis and experimental results strongly
advocate for the effectiveness of the Zero-Inflated Smoothing Spline (ZISS)
method in handling datasets characterized by zero inflation and over-
dispersion. The ZISS method has demonstrated a robust capability in accurately
identifying and interpreting signals within such complex datasets,
particularly excelling in scenarios where traditional models struggle.

In summary, the ZISS method exhibits exceptional performance in scenarios
characterized by high zero-inflation rates, adeptly differentiating between
true zeros and technical zeros in single-cell datasets. Additionally, this
method demonstrates resilience in environments with various levels of over-
dispersion, consistently maintaining both accuracy and reliability. Moreover,
while specifically designed for single-cell temporal data, the ZISS method is
versatile enough to be applied effectively in other domains featuring temporal
and/or spatial data with excessive zeros.

Appendix

##  Appendix A List of Notations

Notations |  Explanations
---|---
$t_{i}$ |  the $i$th pseudotime point
$M_{i}$ |  the number of observations at the pseudotime point $t_{i}$
$y_{i,j}$ |  the $j$-th observation at the $i$th pseudotime point
$g_{i,j}$ |  the latent variable indicating zero processes
$p(t_{i})$ |  the probability function, denoting the probability from the Poisson process at the psutotime point $t_{i}$.
$\alpha_{i}$ |  the coefficients of B-spline basis functions
$b_{i}(\cdot)$ |  the B-spline basis functions
$J(\cdot)$ |  the penalized functional to measure the smoothness of functions, quadratic semi-norm in reproducing kernel Hilbert space
$\boldsymbol{\alpha}$ |  a vector of coefficients $(\alpha_{1},\cdots,\alpha_{m})$
$\mu$ |  a smooth function denoting the mean of the Poisson process
$\lambda$ |  smoothing parameter to control the smoothness of mean function estimate
$\hat{\mu}$ |  the estimate of $\mu$
$\hat{p}$ |  the estimate of $p$

##  Appendix B Derivation Details

Given the single-cell temporal data $y_{i,j}$, the likelihood function $L$ is
given as

$L(\mu,p,y)=\prod_{i=1}^{N}\prod_{j=1}^{M_{i}}\left(p(t_{i})\frac{\mu(t_{i})^{y_{i,j}}}{y_{i,j}!}\right)^{g_{i,j}}\left(1-p(t_{i})\right)^{1-g_{i,j}}.$

We add a penalty on the smoothness of $\mu$ to the logarithm of $L$. Then the
negative penalized log-likelihood function is

$\displaystyle l(\mu,p,y)=$
$\displaystyle-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(g_{i,j}\log
p(t_{i})+(1-g_{i,j})\log(1-p(t_{i}))\right)$ (5)
$\displaystyle-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(g_{i,j}y_{i,j}\log\mu(t_{i})-g_{i,j}\mu(t_{i})\right)+\frac{\lambda}{2}J(\mu)+Const,$

where $J(\cdot)$ is the smoothness penalty functional and $\lambda$ is the
smoothing parameter. We apply Expectation-Maximum (EM) algorithm to estimate
the functions $\mu$ and $p$, where $\mu$ is estimated by the smoothing spline
ANOVA model and $p$ is obtained by using the $B$-spline basis expansion.

#### The Expectation (E) step

We first show the proof of Proposition 1.

Proof of Proposition 1. Given $\hat{\mu}^{(m)}$ and $\hat{p}^{(m)}$ from the
$m$th iteration, we compute the expectation of each latent variable $g_{i,j}$
given observed data $y_{i,j}$, denoted by $q_{i,j}^{(m)}$. Using the Bayesian
rule, we have

$\displaystyle\mathbb{E}\left[g_{i,j}|y_{i,j}=0,\hat{\mu},\hat{p}\right]=\mathbb{P}\left[g_{i,j}=1|y_{i,j}=0,\hat{\mu},\hat{p}\right]$
(6)
$\displaystyle=\frac{\mathbb{P}\left[g_{i,j}=1,y_{i,j}=0|\hat{\mu},\hat{p}\right]}{\mathbb{P}\left[y_{i,j}=0|\hat{\mu},\hat{p}\right]}$
$\displaystyle=\frac{\mathbb{P}\left[y_{i,j}=0|g_{i,j}=1,\hat{\mu},\hat{p}\right]\mathbb{P}[g_{i,j}=1|\hat{\mu},\hat{p}]}{\mathbb{P}\left[y_{i,j}=0|g_{i,j}=1,\hat{\mu},\hat{p}\right]\mathbb{P}[g_{i,j}=1|\hat{\mu},\hat{p}]+\mathbb{P}\left[y_{i,j}=0|g_{i,j}=0,\hat{\mu},\hat{p}\right]\mathbb{P}[g_{i,j}=0|\hat{\mu},\hat{p}]}$
$\displaystyle=\frac{e^{-\hat{\mu}(t_{i})}\hat{p}(t_{i})}{e^{-\hat{\mu}(t_{i})}\hat{p}(t_{i})+1-\hat{p}(t_{i})}$

Otherwise if $y_{i,j}\neq 0$, this implies that $g_{i,j}$ has to be $1$, as a
result,

$\mathbb{E}[g_{i,j}|y_{i,j}\neq 0,\hat{\mu},\hat{p}]=1.$

Thus, we complete the proof. $\square$

#### The Maximization (M) step

Given conditional expectation of $g_{i,j}$’s, i.e., $q_{i,j}$, we minimize the
negative log-likelihood function $l$. Given that $l$ can be represented as the
sum of individual terms for $\mu$ and $p$, it is sufficient to minimize these
two terms separately. Proposition 2 implies that we only need to separately
optimize problems $P_{1}$ and $P_{2}$. For $P_{1}$, we have

$(P_{1})~{}~{}\min_{p\in\mathcal{P}}-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(q_{i,j}\log
p(t_{i})+(1-q_{i,j})\log(1-p(t_{i}))\right),$

where $\mathcal{P}$ is some function space we limit the estimate of $p$ into.
We estimate the probability function by assuming $p(t)$ has the form

$p(t)=\frac{1}{1+exp(\sum_{i=1}^{m}\alpha_{i}b_{i}(t))},$

where $b_{i}(t)$’s are the B-spline basis functions and
$\alpha_{1},\cdots,\alpha_{m}$ are the corresponding coefficients. In this
case,
$\log\left(\frac{p(t)}{1-p(t)}\right)=-\sum_{i=1}^{m}\alpha_{i}b_{i}(t)$. The
advantages of adopting this particular form are detailed in Section 2.2. Then
in the M step, we need to estimate the parameters
$\alpha_{1},\cdots,\alpha_{m}\in\mathbb{R}$ to obtain the estimate of $p$. In
particular, we aim to maximize

$\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}q_{i,j}\log
p(t_{i})+(1-q_{i,j})\log(1-p(t_{i})).$

With the form of $p(t)$, we only need to maximize

$F(\alpha_{1},\cdots,\alpha_{m})=-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}q_{i,j}\sum_{k=1}^{m}\alpha_{k}b_{k}(t_{i})+\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\log\frac{1}{1+exp(-\sum_{k=1}^{m}\alpha_{k}b_{k}(t_{i}))}.$

This is a multi-variable function of
$\boldsymbol{\alpha}=(\alpha_{1},\cdots,\alpha_{m})$, so we use the Newton’s
method [32] to find the maximum of $F(\alpha_{1},\cdots,\alpha_{m})$ with the
following updating rule,

$\boldsymbol{\alpha}\leftarrow\boldsymbol{\alpha}-(\nabla_{\boldsymbol{\alpha}}^{2}F(\boldsymbol{\alpha}))^{-1}\nabla_{\boldsymbol{\alpha}}F(\boldsymbol{\alpha}).$

The derivatives on $\alpha_{k}$,

$\frac{\partial
F}{\partial\alpha_{k}}=-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}q_{i,j}b_{k}(t_{i})+\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\frac{b_{k}(t_{i})}{1+\exp(\sum_{l}\alpha_{l}b_{l}(t_{i}))},$

constructs the vector $\nabla_{\boldsymbol{\alpha}}F(\boldsymbol{\alpha})$,
while the second-order derivatives with respect to $\alpha_{k}$ and
$\alpha_{l}$,

$\frac{\partial^{2}F}{\partial\alpha_{k}\partial\alpha_{l}}=\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\frac{b_{k}(t_{i})b_{l}(t_{i})\exp(\sum_{r=1}^{m}\alpha_{r}b_{r}(t_{i})))}{(1+\exp(\sum_{r=1}^{m}\alpha_{r}b_{r}(t_{i})))^{2}},$

constructs the Hessian of $F$ with respect to $\boldsymbol{\alpha},$
$\nabla^{2}_{\boldsymbol{\alpha}}F(\boldsymbol{\alpha})$.

We use the penalized likelihood method to solve the problem $P_{2}$,

$(P_{2})~{}~{}\min_{\mu}-\sum_{i=1}^{N}\sum_{j=1}^{M_{i}}\left(q_{i,j}y_{i,j}\log\mu(t_{i})-q_{i,j}\mu(t_{i})\right)+\frac{\lambda}{2}J(\mu).$

Based on the Representer Theorem, we incorporate the form of $\mu(t)$ as given
in equation (2) into the objective function above and apply Newton’s method
for optimization. Readers may see [16, Section 5] for more details of
estimation.

## References

  * [1] Angela R Wu, Norma F Neff, Tomer Kalisky, Piero Dalerba, Barbara Treutlein, Michael E Rothenberg, Francis M Mburu, Gary L Mantalas, Sopheak Sim, Michael F Clarke, and Stephen R Quake.  Quantitative assessment of single-cell RNA-sequencing methods.  Nature Methods, 11(1):41, 2014.
  * [2] Florian Buettner, Kedar N Natarajan, F Paolo Casale, Valentina Proserpio, Antonio Scialdone, Fabian J Theis, Sarah A Teichmann, John C Marioni, and Oliver Stegle.  Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells.  Nature Biotechnology, 33(2):155–160, 2015.
  * [3] Efthymia Papalexi and Rahul Satija.  Single-cell RNA sequencing to explore immune cell heterogeneity.  Nature Reviews Immunology, 18(1):35, 2018.
  * [4] Aviv Regev, Sarah A Teichmann, Eric S Lander, Ido Amit, Christophe Benoist, Ewan Birney, Bernd Bodenmiller, Peter Campbell, Piero Carninci, Menna Clatworthy, Hans Clevers, Bart Deplancke, Ian Dunham, James Eberwine, Roland Eils, Wolfgang Enard, Andrew Farmer, Lars Fugger, Berthold Göttgens, Nir Hacohen, Muzlifah Haniffa, Martin Hemberg, Seung Kim, Paul Klenerman, Arnold Kriegstein, Ed Lein, Sten Linnarsson, Emma Lundberg, Joakim Lundeberg, Partha Majumder, John C Marioni, Miriam Merad, Musa Mhlanga, Martijn Nawijn, Mihai Netea, Garry Nolan, Dana Pe’er, Anthony Phillipakis, Chris P Ponting, Stephen Quake, Wolf Reik, Orit Rozenblatt-Rosen, Joshua Sanes, Rahul Satija, Ton N Schumacher, Alex Shalek, Ehud Shapiro, Padmanee Sharma, Jay W Shin, Oliver Stegle, Michael Stratton, Michael J T Stubbington, Fabian J Theis, Matthias Uhlen, Alexander van Oudenaarden, Allon Wagner, Fiona Watt, Jonathan Weissman, Barbara Wold, Ramnik Xavier, Nir Yosef, and Human Cell Atlas Meeting Participants.  The human cell atlas.  Elife, 6, 2017.
  * [5] Wouter Saelens, Robrecht Cannoodt, Helena Todorov, and Yvan Saeys.  A comparison of single-cell trajectory inference methods.  Nature Biotechnology, 37(5):547–554, 2019.
  * [6] Robrecht Cannoodt, Wouter Saelens, and Yvan Saeys.  Computational methods for trajectory inference from single-cell transcriptomics.  European Journal of Immunology, 46(11):2496–2506, 2016.
  * [7] Zehua Liu, Huazhe Lou, Kaikun Xie, Hao Wang, Ning Chen, Oscar M Aparicio, Michael Q Zhang, Rui Jiang, and Ting Chen.  Reconstructing cell cycle pseudo time-series via single-cell transcriptome data.  Nature Communications, 8(1):1–9, 2017.
  * [8] Amos Tanay and Aviv Regev.  Scaling single-cell genomics from phenomenology to mechanism.  Nature, 541(7637):331–338, 2017.
  * [9] Martin Etzrodt, Max Endele, and Timm Schroeder.  Quantitative single-cell approaches to stem cell research.  Cell Stem Cell, 15(5):546–558, 2014.
  * [10] Xi Chen, J Christopher Love, Nicholas E Navin, Lior Pachter, Michael JT Stubbington, Valentine Svensson, Jonathan V Sweedler, and Sarah A Teichmann.  Single-cell analysis at the threshold.  Nature Viotechnology, 34(11):1111–1118, 2016.
  * [11] Wei Vivian Li and Jingyi Jessica Li.  An accurate and robust imputation method scImpute for single-cell RNA-seq data.  Nature Communications, 9(1):1–9, 2018.
  * [12] Peng Qiu.  Embracing the dropouts in single-cell RNA-seq analysis.  Nature Communications, 11(1):1–9, 2020.
  * [13] Blue B Lake, Rajasree Menon, Seth Winfree, Qiwen Hu, Ricardo Melo Ferreira, Kian Kalhor, Daria Barwinska, Edgar A Otto, Michael Ferkowicz, Dinh Diep, et al.  An atlas of healthy and injured cell states and niches in the human kidney.  bioRxiv, pages 2021–07, 2021.
  * [14] Zhicheng Ji and Hongkai Ji.  Tscan: Pseudo-time reconstruction and evaluation in single-cell rna-seq analysis.  Nucleic acids research, 44(13):e117–e117, 2016.
  * [15] Wenpin Hou, Zhicheng Ji, Zeyu Chen, E John Wherry, Stephanie C Hicks, and Hongkai Ji.  A statistical framework for differential pseudotime analysis with multiple single-cell rna-seq samples.  bioRxiv, 2021.
  * [16] Chong Gu and Chong Gu.  Smoothing spline ANOVA models, volume 297.  Springer, 2013.
  * [17] Ruochen Jiang, Tianyi Sun, Dongyuan Song, and Jingyi Jessica Li.  Statistics or biology: the zero-inflation controversy about scrna-seq data.  Genome biology, 23(1):1–24, 2022.
  * [18] Rui Fang, BD Wagner, JK Harris, and SA Fillon.  Zero-inflated negative binomial mixed model: an application to two microbial organisms important in oesophagitis.  Epidemiology & Infection, 144(11):2447–2455, 2016.
  * [19] Dominique Lord, Simon P Washington, and John N Ivan.  Poisson, poisson-gamma and zero-inflated regression models of motor vehicle crashes: balancing statistical fit and theory.  Accident Analysis & Prevention, 37(1):35–46, 2005.
  * [20] Ludwig Fahrmeir and Leyre Osuna Echavarría.  Structured additive regression for overdispersed and zero-inflated count data.  Applied Stochastic Models in Business and Industry, 22(4):351–369, 2006.
  * [21] Chen Xue-Dong.  Bayesian analysis of semiparametric mixed-effects models for zero-inflated count data.  Communications in Statistics—Theory and Methods, 38(11):1815–1833, 2009.
  * [22] Simon C Barry and Alan H Welsh.  Generalized additive modelling and zero inflated count data.  Ecological Modelling, 157(2-3):179–188, 2002.
  * [23] Hai Liu and Kung-Sik Chan.  Generalized additive models for zero-inflated data with partial constraints.  Scandinavian Journal of Statistics, 38(4):650–665, 2011.
  * [24] Hai Liu and Kung-Sik Chan.  Constrained generalized additive model with zero-inflated data.  Technical report, Citeseer, 2008.
  * [25] Hai Liu and Kung-Sik Chan.  Introducing cozigam: An r package for unconstrained and constrained zero-inflated generalized additive model analysis.  Journal of Statistical Software, 35:1–26, 2010.
  * [26] Todd K Moon.  The expectation-maximization algorithm.  IEEE Signal processing magazine, 13(6):47–60, 1996.
  * [27] Peter Craven and Grace Wahba.  Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation.  Numerische Mathematik, 31(4):377–403, 1978.
  * [28] Chong Gu and Grace Wahba.  Minimizing GCV/GML scores with multiple smoothing parameters via the Newton method.  SIAM Journal on Scientific and Statistical Computing, 12(2):383–398, 1991.
  * [29] Martin Ridout, John Hinde, and Clarice GB Demétrio.  A score test for testing a zero-inflated poisson regression model against zero-inflated negative binomial alternatives.  Biometrics, 57(1):219–223, 2001.
  * [30] Ilya Korsunsky, Nghia Millard, Jean Fan, Kamil Slowikowski, Fan Zhang, Kevin Wei, Yuriy Baglaenko, Michael Brenner, Po-ru Loh, and Soumya Raychaudhuri.  Fast, sensitive and accurate integration of single-cell data with harmony.  Nature methods, 16(12):1289–1296, 2019.
  * [31] Mengqi Zhang, Si Liu, Zhen Miao, Fang Han, Raphael Gottardo, and Wei Sun.  IDEAS: individual level differential expression analysis for single-cell RNA-seq data.  Genome Biology, 23(1):1–17, 2022.
  * [32] Joseph Frederick Traub and H Woźniakowski.  Convergence and complexity of newton iteration for operator equations.  Journal of the ACM (JACM), 26(2):250–258, 1979.