In williazo/tcensReg: MLE of a Truncated Normal Distribution with Censored Data
knitr::opts_chunk$set(echo = TRUE)
The goal of this package is to estimate parameters from a linear model when the data comes from a truncated normal distribution with censoring. Maximum likelihood values are returned. There are multiple method available for optimization with the default set as conjugate gradient. This package is also able to return maximum likelihood estimates for truncated only or censored only data similar to truncreg and censReg packages.
Installation
You can install tcensReg from CRAN for the stable release or install the GitHub version for the active development version:
#stable CRAN version
install.packages("tcensReg")
# or ----
#active devel. GitHub version
install.packages("devtools")
devtools::install_github("williazo/tcensReg")
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)
Next, 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)
#collecting the uncensored and censored data together
dt <- data.frame(y_star, y)
We can observe the histogram and density plot for the uncensored data, which shows the zero-truncation.
library(ggplot2)
library(viridis)
ggplot(data=dt, aes(x=y_star, y=..density..))+
geom_histogram(binwidth=0.25, col=viridis(1),fill=viridis(1), alpha=0.6)+
stat_density(bw=0.5, col=viridis(1), fill=viridis(1), alpha=0.3)+
scale_x_continuous(breaks=seq(0, 2, 0.25))+
ylab("Density")+
xlab(expression(Y^"*"))+
theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), panel.grid.major=element_blank(),
panel.background=element_blank(), panel.border=element_rect(fill=NA, color="black"),
axis.text.x=element_text(size=10), legend.title=element_text(size=15),
legend.text=element_text(size=15), legend.key.width=unit(15, units="mm"), strip.text=element_text(size=20),
axis.title=element_text(size=20))
We can then compare this to the censored observations below
ggplot(data=dt, aes(x=y, y=..density..))+
geom_histogram(binwidth=0.25, col=viridis(1, begin=0.5),fill=viridis(1, begin=0.5), alpha=0.6)+
stat_density(bw=0.25, col=viridis(1, begin=0.5), fill=viridis(1, begin=0.5), alpha=0.3)+
scale_x_continuous(breaks=seq(0, 2, 0.25))+
ylab("Density")+
xlab(expression(Y))+
theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), panel.grid.major=element_blank(),
panel.background=element_blank(), panel.border=element_rect(fill=NA, color="black"),
axis.text.x=element_text(size=10), legend.title=element_text(size=15),
legend.text=element_text(size=15), legend.key.width=unit(15, units="mm"), strip.text=element_text(size=20),
axis.title=element_text(size=20))
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)
By 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.
names(t_mod)
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)
cens_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)
Other methods result in significant bias for both mu and sigma.
Example 2: Two Population Model with Separate Variances
As an extension for the single population model above, we can imagine a two independent truncated normal random variables that have common censoring and truncation values but different standard deviations.
We can simulate the underlying truncated normal distributions Y1_star and Y2_star similar to $Y$ above except now we allow them to have separate mean and variances.
For this example we let mu_1=0.5, mu_2=1, sigma_1=0.25, sigma_2=2, and a=0.
mu_1 <- 0.5
mu_2 <- 1
sigma_1 <- 0.25
sigma_2 <- 2
a <- 0
set.seed(032420)
y_1_star <- rtnorm(1000, mu = mu_1, sd = sigma_1, a = a)
y_2_star <- rtnorm(1000, mu = mu_2, sd = sigma_2, a = a)
df <- data.frame(y_star = c(y_1_star, y_2_star),
group = c(rep("Population 1", length(y_1_star)),
rep("Population 2", length(y_2_star))))
Plotting each of these uncensored population densities, we can see the difference in shape based on the underlying parameter selection.
ggplot(data = df, aes(x = y_star, y = ..density.., group = group, fill = group, col = group))+
stat_density(bw = 0.5, alpha = 0.3)+
ylab("Density")+
xlab(expression(Y^"*"))+
geom_vline(xintercept = 0.5, lty = 2, col = viridis(1, begin = 0.25))+
geom_vline(xintercept = 1, lty = 2, col = viridis(1, begin = 0))+
scale_color_manual(name = "", values = viridis(2, begin = 0.25, end = 0))+
scale_fill_manual(name = "", values = viridis(2, begin = 0.25, end = 0))+
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), panel.grid.major = element_blank(),
panel.background = element_blank(), panel.border = element_rect(fill = NA, color = "black"),
axis.text.x = element_text(size = 10), legend.title = element_text(size = 15),
legend.text = element_text(size = 15), legend.key.width = unit(15, units = "mm"), strip.text = element_text(size = 20),
axis.title = element_text(size = 20))
Then censoring each observation at nu, we are left with Y1 and Y2. Again, we let nu=0.25.
nu <- 0.25
df$y <- ifelse(df$y_star<=nu, nu, df$y_star)
We then can fit our model with separate variances for each group using the command tcensReg_sepvar as shown below.
mod_result <- tcensReg_sepvar(y ~ group, a=a, v=nu, group_var="group", method="maxLik", data=df)
mod_result
sepvar_est <- mod_result$theta
mu_1_est <- sepvar_est[1]
mu_2_est <- sum(sepvar_est[1:2])
sigma_1_est <- exp(sepvar_est[3])
sigma_2_est <- exp(sepvar_est[4])
results_df <- data.frame(rbind(c(mu_1, mu_2, sigma_1, sigma_2),
c(mu_1_est, mu_2_est, sigma_1_est, sigma_2_est)))
names(results_df) <- c("mu_1", "mu_2", "sigma_1", "sigma_2")
row.names(results_df) <- c("Truth", "tcensReg")
results_df$mu1_pct_bias <- paste0(round(((results_df$mu_1 - mu_1)/mu_1)*100, 2), "%")
results_df$mu2_pct_bias <- paste0(round(((results_df$mu_2 - mu_2)/mu_2)*100, 2), "%")
results_df$sigma1_pct_bias <- paste0(round(((results_df$sigma_1 - sigma_1)/sigma_1)*100, 2), "%")
results_df$sigma2_pct_bias <- paste0(round(((results_df$sigma_2 - sigma_2)/sigma_2)*100, 2), "%")
knitr::kable(results_df, format="markdown", digits=4)
Performance Comparison: Censored-Only and Truncated-Only (Intercept Only Model)
Note also that the tcensReg can also estimate parameters in the censored-only or truncated-only cases. We show below that by using analytic values in the tcensReg implementation that our method is faster then the alternative estimation procedures while providing better variance estimates. With a small set of covariates and p<<n we can use the Newton Raphson method of optimization, which is computationally fast with few covariates.
library(microbenchmark)
#testing the censored-only regression
library(censReg)
cens <- microbenchmark(tcensReg_method = tcensReg(y ~ 1, data=dt, v=nu, method="Newton"),
censReg_method = censReg(y ~ 1, left=nu, data=dt))
knitr::kable(summary(cens), format="markdown", digits=4)
#point estimates are equivalent
tcensReg_est <- as.numeric(tcensReg(y ~ 1, data=dt, v=nu, method="Newton")$theta)
censReg_est <- as.numeric(coef(censReg(y ~ 1, left=nu, data=dt)))
all.equal(tcensReg_est, censReg_est)
#testing the truncated-only regression
library(truncreg)
trunc <- microbenchmark(
tcensReg_method = tcensReg(y_star ~ 1, data=dt, a=a, method="Newton"),
truncreg_method = truncreg(y_star ~ 1, point=a, data=dt))
knitr::kable(summary(trunc), format="markdown", digits=4)
tcensReg_est <- as.numeric(tcensReg(y_star ~ 1, data=dt, a=a, method="Newton")$theta)
#note truncreg returns sigma not log_sigma so we need to exponentiate our value
tcensReg_est[2] <- exp(tcensReg_est[2])
truncreg_est <- as.numeric(coef(truncreg(y_star ~ 1, point=a, data=dt)))
all.equal(tcensReg_est, truncreg_est)
williazo/tcensReg documentation built on July 24, 2020, 1:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
The goal of this package is to estimate parameters from a linear model when the data comes from a truncated normal distribution with censoring. Maximum likelihood values are returned. There are multiple method available for optimization with the default set as conjugate gradient. This package is also able to return maximum likelihood estimates for truncated only or censored only data similar to truncreg and censReg packages.
Installation
You can install tcensReg from CRAN for the stable release or install the GitHub version for the active development version:
#stable CRAN version install.packages("tcensReg") # or ---- #active devel. GitHub version install.packages("devtools") devtools::install_github("williazo/tcensReg")
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)
Next, 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) #collecting the uncensored and censored data together dt <- data.frame(y_star, y)
We can observe the histogram and density plot for the uncensored data, which shows the zero-truncation.
library(ggplot2) library(viridis) ggplot(data=dt, aes(x=y_star, y=..density..))+ geom_histogram(binwidth=0.25, col=viridis(1),fill=viridis(1), alpha=0.6)+ stat_density(bw=0.5, col=viridis(1), fill=viridis(1), alpha=0.3)+ scale_x_continuous(breaks=seq(0, 2, 0.25))+ ylab("Density")+ xlab(expression(Y^"*"))+ theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), panel.grid.major=element_blank(), panel.background=element_blank(), panel.border=element_rect(fill=NA, color="black"), axis.text.x=element_text(size=10), legend.title=element_text(size=15), legend.text=element_text(size=15), legend.key.width=unit(15, units="mm"), strip.text=element_text(size=20), axis.title=element_text(size=20))
We can then compare this to the censored observations below
ggplot(data=dt, aes(x=y, y=..density..))+ geom_histogram(binwidth=0.25, col=viridis(1, begin=0.5),fill=viridis(1, begin=0.5), alpha=0.6)+ stat_density(bw=0.25, col=viridis(1, begin=0.5), fill=viridis(1, begin=0.5), alpha=0.3)+ scale_x_continuous(breaks=seq(0, 2, 0.25))+ ylab("Density")+ xlab(expression(Y))+ theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), panel.grid.major=element_blank(), panel.background=element_blank(), panel.border=element_rect(fill=NA, color="black"), axis.text.x=element_text(size=10), legend.title=element_text(size=15), legend.text=element_text(size=15), legend.key.width=unit(15, units="mm"), strip.text=element_text(size=20), axis.title=element_text(size=20))
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)
By 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.
names(t_mod)
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) cens_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)
Other methods result in significant bias for both mu and sigma.
Example 2: Two Population Model with Separate Variances
As an extension for the single population model above, we can imagine a two independent truncated normal random variables that have common censoring and truncation values but different standard deviations.
We can simulate the underlying truncated normal distributions Y1_star and Y2_star similar to $Y$ above except now we allow them to have separate mean and variances.
For this example we let mu_1=0.5, mu_2=1, sigma_1=0.25, sigma_2=2, and a=0.
mu_1 <- 0.5 mu_2 <- 1 sigma_1 <- 0.25 sigma_2 <- 2 a <- 0 set.seed(032420) y_1_star <- rtnorm(1000, mu = mu_1, sd = sigma_1, a = a) y_2_star <- rtnorm(1000, mu = mu_2, sd = sigma_2, a = a) df <- data.frame(y_star = c(y_1_star, y_2_star), group = c(rep("Population 1", length(y_1_star)), rep("Population 2", length(y_2_star))))
Plotting each of these uncensored population densities, we can see the difference in shape based on the underlying parameter selection.
ggplot(data = df, aes(x = y_star, y = ..density.., group = group, fill = group, col = group))+ stat_density(bw = 0.5, alpha = 0.3)+ ylab("Density")+ xlab(expression(Y^"*"))+ geom_vline(xintercept = 0.5, lty = 2, col = viridis(1, begin = 0.25))+ geom_vline(xintercept = 1, lty = 2, col = viridis(1, begin = 0))+ scale_color_manual(name = "", values = viridis(2, begin = 0.25, end = 0))+ scale_fill_manual(name = "", values = viridis(2, begin = 0.25, end = 0))+ theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), panel.border = element_rect(fill = NA, color = "black"), axis.text.x = element_text(size = 10), legend.title = element_text(size = 15), legend.text = element_text(size = 15), legend.key.width = unit(15, units = "mm"), strip.text = element_text(size = 20), axis.title = element_text(size = 20))
Then censoring each observation at nu, we are left with Y1 and Y2. Again, we let nu=0.25.
nu <- 0.25 df$y <- ifelse(df$y_star<=nu, nu, df$y_star)
We then can fit our model with separate variances for each group using the command tcensReg_sepvar as shown below.
mod_result <- tcensReg_sepvar(y ~ group, a=a, v=nu, group_var="group", method="maxLik", data=df) mod_result sepvar_est <- mod_result$theta mu_1_est <- sepvar_est[1] mu_2_est <- sum(sepvar_est[1:2]) sigma_1_est <- exp(sepvar_est[3]) sigma_2_est <- exp(sepvar_est[4]) results_df <- data.frame(rbind(c(mu_1, mu_2, sigma_1, sigma_2), c(mu_1_est, mu_2_est, sigma_1_est, sigma_2_est))) names(results_df) <- c("mu_1", "mu_2", "sigma_1", "sigma_2") row.names(results_df) <- c("Truth", "tcensReg") results_df$mu1_pct_bias <- paste0(round(((results_df$mu_1 - mu_1)/mu_1)*100, 2), "%") results_df$mu2_pct_bias <- paste0(round(((results_df$mu_2 - mu_2)/mu_2)*100, 2), "%") results_df$sigma1_pct_bias <- paste0(round(((results_df$sigma_1 - sigma_1)/sigma_1)*100, 2), "%") results_df$sigma2_pct_bias <- paste0(round(((results_df$sigma_2 - sigma_2)/sigma_2)*100, 2), "%") knitr::kable(results_df, format="markdown", digits=4)
Performance Comparison: Censored-Only and Truncated-Only (Intercept Only Model)
Note also that the tcensReg can also estimate parameters in the censored-only or truncated-only cases. We show below that by using analytic values in the tcensReg implementation that our method is faster then the alternative estimation procedures while providing better variance estimates. With a small set of covariates and p<<n we can use the Newton Raphson method of optimization, which is computationally fast with few covariates.
library(microbenchmark) #testing the censored-only regression library(censReg) cens <- microbenchmark(tcensReg_method = tcensReg(y ~ 1, data=dt, v=nu, method="Newton"), censReg_method = censReg(y ~ 1, left=nu, data=dt)) knitr::kable(summary(cens), format="markdown", digits=4) #point estimates are equivalent tcensReg_est <- as.numeric(tcensReg(y ~ 1, data=dt, v=nu, method="Newton")$theta) censReg_est <- as.numeric(coef(censReg(y ~ 1, left=nu, data=dt))) all.equal(tcensReg_est, censReg_est) #testing the truncated-only regression library(truncreg) trunc <- microbenchmark( tcensReg_method = tcensReg(y_star ~ 1, data=dt, a=a, method="Newton"), truncreg_method = truncreg(y_star ~ 1, point=a, data=dt)) knitr::kable(summary(trunc), format="markdown", digits=4) tcensReg_est <- as.numeric(tcensReg(y_star ~ 1, data=dt, a=a, method="Newton")$theta) #note truncreg returns sigma not log_sigma so we need to exponentiate our value tcensReg_est[2] <- exp(tcensReg_est[2]) truncreg_est <- as.numeric(coef(truncreg(y_star ~ 1, point=a, data=dt))) all.equal(tcensReg_est, truncreg_est)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.