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.001 2.189Next, 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.174We 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.
#loading the package
library(tcensReg)
t_mod <- tcensReg(y ~ 1, data=dt, a=0, v=0.25)
summary(t_mod)##
## Call:
## tcensReg(formula = y ~ 1, a = 0, v = 0.25, data = dt)
##
## Assumed Distribution:
## Truncated Normal with Censoring
##
## Count of Observations:
## Total Observations Censored Observations
## 1000 174
##
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.4902 0.0305 16.0765
## sigma 0.5092 0.0139 36.6421
##
## Log Likelihood: -657.218
## Information Criterion: AIC=1318.4359 BIC=1328.2514
## Optimization Method: CG
## Psuedo R2: 0 method - mckelvey_zavoinaBy default the coefficients are returned along with log likelihood and other fit criterion statistics. Note that the Pseudo R2 in the case of an intercept only model is exactly equal to zero.
## [1] "theta" "convergence" "initial_ll"
## [4] "final_ll" "var_cov" "method"
## [7] "info_criteria" "model_matrix" "call"
## [10] "n_count" "latent_assumption"Note that the this object contains parameter estimates theta, convergence criterion code, initial/final log-likelihood values, variance-covariance matrix, method of optimization, information criterion, design matrix used from the model, formula call, count of total/censored observations, and latent distributional assumption.
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 <- coef(output) #this returns sigma rather than log sigma
#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 <- coef(cens_output)
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 <- results_df$mu - mu
results_df$sigma_bias <- 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.4902 | 0.5092 | -0.0098 | 0.0092 |
| Normal MLE | 0.6643 | 0.3751 | 0.1643 | -0.1249 |
| Tobit | 0.6219 | 0.4375 | 0.1219 | -0.0625 |
Other methods result in significant bias for both mu and sigma.
