## 3 Estimation Procedure Based on StEM  
### The Models and Notations  
In this section we describe the models and methods in a general form to illustrate their applicability to other applications. Let \(y_{ij}\), with \(j=1,2,\ldots,n_{i}\), represent the longitudinal measurements for subject \(i=1,2,\ldots,n\), taken at time \(t_{ij}\). We consider a general ZIE-based nonlinear random change point model:  
\[\begin{split}& y_{ij}=g((t_{ij}-\tau_{i})^{-},\mathbf{a}_{i})+h((t_{ij}- \tau_{i})^{+},\mathbf{b}_{i})+e_{ij},\\ &\mathbf{a}_{i}\sim N(\mathbf{\alpha},A),\quad\mathbf{b}_{i}\sim N(\mathbf{\beta },B),\quad\tau_{i}\sim ZIE(\pi,\lambda),\quad e_{ij}|\mathbf{a}_{i},\mathbf{b}_{i},\tau_ {i}\sim N(0,\sigma_{e}^{2}).\end{split} \tag{2}\]  
Here, \(g(\cdot)\) and \(h(\cdot)\) are known nonlinear functions, and \(\mathbf{\alpha}\) and \(\mathbf{\beta}\) are vectors of population parameters. The random effects \(\mathbf{a}_{i}\) and \(\mathbf{b}_{i}\) follow normal distribution with mean \(\mathbf{\alpha}\) and \(\mathbf{\beta}\), and covariance matrices \(A\) and \(B\), respectively. The random change point \(\tau_{i}\) is assumed to follow a zero-inflated exponential distribution with parameters \(\pi\) and \(\lambda\). The within-individual variance is denoted as \(\sigma_{e}^{2}\). The function \(x(\cdot)^{-}\) and \(x(\cdot)^{+}\) are defined as before. We assume that \(\mathbf{a}_{i}\), \(\mathbf{b}_{i}\) are independent and both are independent of \(\tau_{i}\), which is introduced externally.  
To estimate and infer the model (2) for left-censored data, we employ a likelihood-based estimation procedure using the observed data \(\{(q_{ij},c_{ij}),i=1,\ldots,n,j=1,\ldots,n_{i}\}\) where \(c_{ij}\) is the censoring indicator such that \(y_{ij}\) is observed if \(c_{ij}=0\) and \(y_{ij}\) is left-censored if \(c_{ij}=1\). That is, \(y_{ij}=q_{ij}\) if \(c_{ij}=0\) and \(y_{ij}\leq q_{ij}\) (\(\equiv d\)) if \(c_{ij}=1\), where \(d\) is the detection limit. Extension to the \"doubly-censored\" case, in which the response may either be left-censored or right-censored, is straightforward.  
Let \(\mathbf{\theta}=(\mathbf{\alpha},\mathbf{\beta},A,B,\sigma_{e}^{2},\pi,\lambda)\) denote the collection of all unknown parameters, and let \(f(\cdot)\) be a generic density function, with \(f(X|Y)\) denoting the conditional density of \(X\) given \(Y\). The observed data likelihood is given by the following:  
\[\begin{split} L(\mathbf{\theta})=\prod_{i=1}^{n}\Big{\{}\int\int \int\Big{[}\prod_{j=1}^{n_{i}}f(y_{ij}|\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i};\mathbf{ \theta})^{1-c_{ij}}& F(d|\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i};\mathbf{ \theta})^{c_{ij}}\Big{]}\\ &\times f(\mathbf{a}_{i})f(\mathbf{b}_{i})f(\tau_{i})d\mathbf{a}_{i}d\mathbf{b}_{ i}d\tau_{i}\Big{\}},\end{split} \tag{3}\]  
where \(F(d|\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i};\mathbf{\theta})=\int_{-\infty}^{d}f(y_{ij}| \mathbf{a}_{i},\mathbf{b}_{i},\tau_{i};\mathbf{\theta})dy_{ij}\). Directly maximizing the likelihood (3) is challenging due to the presence of mixed-type distributions, nonlinear models, and nested integrals. The numerical methods, e.g., Gauss-Hermite Quadrature [16], can be prohibitively intensive for the computation. We therefore resort to EM algorithm-based methods. By treating \(\mathbf{y}_{cen,i}\), the censored component of \(\mathbf{y}_{i}\), and the unobserved random effects \(\mathbf{a}_{i}\), \(\mathbf{b}_{i}\), and \(\tau_{i}\) as \"missing data\", we have the \"complete data\" \(\big{\{}(y_{i},\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i}),i=1,\ldots,n\big{\}}\). The complete-data log-likelihood function for individual \(i\) is expressed as  
\[l_{e}^{(l)}(\mathbf{\theta})=\log f(\mathbf{y}_{i}|\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i};\mathbf{ \theta})+\log f(\mathbf{a}_{i};A)+\log f(\mathbf{b}_{i};B)+\log f(\tau_{i};\mathbf{\theta}). \tag{4}\]