In clustered survival analysis applications, researchers frequently fit frailty models using parametric and nonparametric approaches to obtain the estimates for the parameters associated with the survival model covariates and heterogeneity (frailty). Availability of the off- the-shelve implementations and freely available R software packages makes it convenient for the practitioners to fit these complicated models easily. Even though there has been a couple of studies assessing the stability of the older packages (e.g., survival, coxme) under a variety of scenarios, some of the newer implementations (e.g., frailtySurv, JM and parfm) have not gone through similar rigorous assessment. It is worth evaluating these new software implementations, and comparing them with the older packages. In the current work, via simulations, we will examine the estimates from all of these popularly used software implementations under a variety of scenarios when the corresponding assumptions related to the baseline hazard and frailty distributions are misspecified. Additionally, true heterogeneity parameter, censoring patterns and number of clusters were varied in the simulations to assess respective impacts on the estimates. From these simulations, we observed that when there is a large number of clusters and mild censoring, Cox PH frailty models fitted using a newer semiparametric estimation technique (from the frailtySurv package) produced regression and heterogeneity parameter estimates that were associated with unusually large bias and variability. On the other hand, when the true heterogeneity parameter is substantially large, the Cox PH frailty models fitted using the coxme package were often producing highly variable estimates of the heterogeneity parameter. The simulation findings then guided our choice of appropriate frailty model in the context of determining the birth interval dynamics in Bangladesh.