Example 1: Single Population
Some common examples where this type of problem may arise is when there is a natural truncation imposed by the structure of the data. For instance several applications have an implied zero truncation such as product lifetimes, age, or detection thresholds. To show how to implement the functions within the package, I will demonstrate a simple simulation example.
Assume that we have observations from an underlying truncated normal distribution. In our case we will assume a zero-truncated model by setting a=0. We generate this truncated normal data below and refer to it as y_star.
#loading the package
library(tcensReg)
mu <- 0.5
sigma <- 0.5
a <- 0
#generate random values from the truncated normal distribution using tcensReg function
y_star <- rtnorm(n=1000, mu=mu, sd=sigma, a=a)
#note that the lowerbound will always be non-negative
round(range(y_star), 3)## [1] 0.003 2.063Next, we can imagine a scenario where we have an imprecise measurement of y_star leading to censoring. In our case we assume that values below a limit of detection, nu, are censored. This creates a random variable y.
In the example below we set our limit of detection as nu=0.25.
nu <- 0.25
y <- ifelse(y_star <= nu, nu, y_star)
#calculating the number of censored observations
sum(y == nu)/length(y) ## [1] 0.182We can observe the histogram and density plot for the uncensored data, which shows the zero-truncation. 
We can then compare this to the censored observations below 
We can then estimate the mean, mu, and standard deviation, sigma, using y with the tcensReg package as shown below.
## $theta
## Estimate
## (Intercept) 0.5067501
## log_sigma -0.6774143
##
## $convergence
## [1] 0
##
## $initial_ll
## [1] -691.7561
##
## $final_ll
## [1] -677.5373
##
## $var_cov
## (Intercept) log_sigma
## (Intercept) 0.0008743094 -0.0008562579
## log_sigma -0.0008562579 0.0015877418
##
## $method
## [1] "CG"Note that the this will return parameter estimates, variance-covariance matrix, the number of iterations until convergence, and the initial/final log-likelihood values.
Comparing the values to the truth we see that the estimates are unbiased.
#tcensReg model
output <- tcensReg(y ~ 1, data=dt, a=a, v=nu)
#extracting the point estimates
tcensReg_est <- output$theta
#exponentiating the estimate of log_sigma to obtain sigma
tcensReg_est[2] <- exp(tcensReg_est[2])
#OLS model
lm_output <- lm(y ~ 1, data=dt)
lm_est <- c(coef(lm_output), summary(lm_output)$sigma)
#censored only model, i.e., Tobit model
cens_output <- tcensReg(y ~ 1, data=dt, v=nu) ## Warning: `a` is not specified indicating no truncationcens_est <- cens_output$theta
cens_est[2] <- exp(cens_est[2])
results_df <- data.frame(rbind(c(mu, sigma),
t(tcensReg_est),
lm_est,
t(cens_est)))
names(results_df) <- c("mu", "sigma")
row.names(results_df) <- c("Truth", "tcensReg", "Normal MLE", "Tobit")
results_df$mu_bias <- abs(results_df$mu - mu)
results_df$sigma_bias <- abs(results_df$sigma - sigma)
knitr::kable(results_df, format="markdown", digits=4)| mu | sigma | mu_bias | sigma_bias | |
|---|---|---|---|---|
| Truth | 0.5000 | 0.5000 | 0.0000 | 0.0000 |
| tcensReg | 0.5068 | 0.5079 | 0.0068 | 0.0079 |
| Normal MLE | 0.6746 | 0.3759 | 0.1746 | 0.1241 |
| Tobit | 0.6300 | 0.4424 | 0.1300 | 0.0576 |
Other methods result in significant bias for both mu and sigma.
