## 5 Simulation Studies  
The performance of the StEM algorithm was evaluated through simulation. Here, we design two simulation studies aimed at assessing and comparing their estimation accuracy under different specifications for the post-ART segment, inspired by the real data analysis.  
We model irregular viral load recording times since HIV diagnosis using a progressive state-transition model, assuming a first-order Markov process. This means that the length of time between viral records depends on the previous recording time. To generate a stochastic measurement time \((t_{i1},t_{i2},\ldots,t_{in_{i}})\), we use parameters obtained from fitting the model to the actual viral load test dates in HCCD. Specifically, we assume that the viral load measurement time \(T\) follows an exponential distribution with parameter \(\xi>0\). Given the previous recording time \(u>0\), the next recording time, conditioned on \(u\), is given by \(T_{|u}=-\frac{1}{\xi}\log(X)+u\), where \(X\sim\text{Uniform}(0,1)\). For this simulation, we set \(\xi=1.6\), which is estimated from the real data. The recording time is terminated when the time exceeds 3 years to simulate the real data. For left-censoring, in addition to simulating the 60% as in real data, we also run the model without any censoring to assess the StEM algorithm under such an ideal scenario.  
To efficiently manage computational resources, each simulated dataset included 1000 individuals. In order to thoroughly evaluate the accuracy and bias of the estimated parameter values, we conducted 200 simulation runs for each scenario. Across these simulation replicates, we calculated both the mean squared error (MSE) and the bias (Bias) by comparing the estimated parameter values with the true values: \(\text{Bias}=\frac{1}{\xi}\sum_{s=1}^{S}(\hat{\theta}^{(s)}-\theta^{\text{true}}),\text{MSE}=\frac{1}{5}\sum_{s=1}^{S}100\times\frac{\sqrt{(\hat{\theta}^{(s)} -\theta^{\text{true}})^{2}+\xi\text{SE}(\hat{\theta}^{(s)})}}{|\theta^{\text{ true}}|}\). Here, \(S\) is the number of replicates, and \(\hat{\theta}^{(s)}\) is the estimate from simulation \(s\).  
Tables 2 and 3 display the simulation results for studies 1 and 2, respectively. We see that there is no major performance difference between M1 and M2. For each model, the fixed effects are estimated better than the dispersion parameters in general. The largest MSE occurred in the estimation of \(\sigma_{A}^{2}\) under M1 when 60% of viral loads are left-censored. For M2, we see \(\sigma_{\text{B11}}^{2}\) is estimated with the largest bias, while \(\beta_{4}\), the second-phase decay rate, is estimated with the biggest MSE. Such sub-optimal performances are likely due to the sparsity of the observation frequency where insufficient data are available to provide the best estimation.  
\begin{table}
\begin{tabular}{c c c c c c c c c c c c} \hline \hline \(\theta\) & \(\alpha\) & \(\beta_{1}\) & \(\beta_{2}\) & \(\beta_{3}\) & \(\sigma_{A}^{2}\) & \(\sigma_{\text{B11}}^{2}\) & \(\sigma_{\text{B22}}^{2}\) & \(\sigma_{\text{B33}}^{2}\) & \(\sigma_{e}^{2}\) & \(\pi\) & \(\lambda\) \\
**True** & **0.43** & **10.51** & **3.92** & **1.69** & **0.49** & **4.61** & **0.64** & **1.92** & **0.10** & **0.31** & **1.11** \\ \hline \multicolumn{10}{c}{0\% left-censored} \\ Est & 0.36 & 10.40 & 3.92 & 1.67 & 0.57 & 5.27 & 0.65 & 1.98 & 0.11 & 0.32 & 1.03 \\ MSE & 0.40 & 0.02 & 0.02 & 0.08 & 0.40 & 0.20 & 0.25 & 0.14 & 0.29 & 0.31 & 0.16 \\ Bias & \(-0.7\) & \(-0.11\) & 0.00 & \(-0.02\) & 0.08 & 0.66 & 0.01 & 0.06 & 0.03 & 0.01 & \(-0.11\) \\ \multicolumn{10}{c}{60\% left-censored} \\ Est & 0.39 & 10.38 & 4.02 & 2.01 & 1.20 & 6.59 & 0.76 & 1.34 & 0.14 & 0.34 & 0.99 \\ MSE & 0.49 & 0.03 & 0.04 & 0.26 & 1.50 & 0.45 & 0.35 & 0.52 & 0.43 & 0.31 & 0.17 \\ Bias & \(-0.04\) & \(-0.13\) & 0.10 & 0.32 & 0.71 & 1.98 & 0.12 & \(-0.58\) & 0.05 & \(-0.02\) & \(-0.09\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Simulation results on the performance of StEM algorithm when the post-ART trajectory is modeled with one-compartment model (M1).