The main focus of Cai and Betensky (2003) was for the analysis of interval censored data under the Cox model. They proposed to weakly parameterize the baseline hazard using a piecewise linear spline and maximize the likelihood function with a PQL approximation. Theoretical and numerical results demonstrated the validity of the proposed point and interval estimators for the general case. For the special case with right censored data, Cai and Betensky (2003) suggested the use of Poisson mixed model to approximate the likelihood for computational ease. However, as pointed out by Bove and Held (2013), the Poisson approximation given in Cai and Betensky (2003) was incorrect. Bove and Held (2013) suggested to create pseudo observations yij = I(i = j)δi and fit a Poisson mixed model to {yij : i = 1, …, n, j = 1, …, i}. This approach will indeed lead to a valid approximation for the point estimator. On the other hand, since their new proposal involves a pseudo data with n(n + 1)/2 observations, there might be some computational burden in fitting such a Poisson mixed model even with moderate sample sizes such as n = 500. Furthermore, one may not directly obtain correct standard error estimates from the fitting since the yij’s are not independent observations. An alternative strategy to obtain the maximizer of l(θ;σb)=δT(X−Tθ−+ZTβ)−∑i=1nΛi(ti)−12σb2bTb−K2σb2 is to consider an iterative procedure based on ∑i=1nΛi(ti)≈∑i=1n∑j=1iqijλ0(Tj)exp(ZiTβ)=∑i=1n∑j=1iexp{η0(Tj)}qijexp(ZiTβ) (1) using the same notation as Bove and Held, where θ− = (α0, α1, bT)T and X− = [1, Ti, (Ti − κ)+]i=1,…,n. Specifically, obtain an initial estimate of β, β(0), say as the standard maximum partial likelihood estimator. Then one may iterate via the following steps starting from m = 1: Let qij(m*)=qijexp(ZiTβ^(m−1)) and q¯(m*)=[∑i=jnqij(m*)]j=1,…,n. Maximize lp(m*)(θ−;σb)≈δTX−Tθ−−1Texp{X−Tθ−+log(q¯(m*))}−12σb2bTb−K2σb2 with respect to {θ−, σb} to obtain {θ^−(m),σ^b(m)} and the corresponding η^0(m)(Tj). Let qij(m†)=qijexp{η^0(m)(Tj)} and q¯(m†)=[∑j=1iqij(m†)]i=1,…,n. Maximize lp(m†)(β)≈δTZTβ−1Texp{ZTβ+log(q¯(m†))} with respect to β to obtain β(m). Let m = m + 1 and go back to Step (1) until convergence. The maximizations can be achieved by fitting a Poisson mixed model with offset in Step (1) and a standard Poisson model with offset in Step (2).