6.8 Tobit model
The dependent variable is partially observed in Tobit models due to sampling schemes, whereas the regressors are completely observed. In particular,
\[\begin{equation} y_i = \begin{Bmatrix} L, & \quad y_i^* < L, \\ y_i^*, & \quad L \leq y_i^* < U, \\ U, & \quad y_i^* \geq U, \end{Bmatrix} \end{equation}\]
where \(y_i^* \stackrel{i.i.d.}{\sim} N(\mathbf{x}_i^{\top} \boldsymbol{\beta}, \sigma^2)\), \(\mathbf{x}_i\) is a \(K\)-dimensional vector of regressors.38
We use conjugate independent priors \(\boldsymbol{\beta} \sim N(\boldsymbol{\beta}_0, \mathbf{B}_0)\) and \(\sigma^2 \sim IG(\alpha_0/2, \delta_0/2)\), and data augmentation using \(\mathbf{y}^*_C\) such that \(y_{C_i}^* \stackrel{i.n.d.}{\thicksim} N(\mathbf{x}_i^{\top} \boldsymbol{\beta}, \sigma^2)\), \(y_{C_i} = \left\{ y_{C_i^L}^* \cup y_{C_i^U}^* \right\}\) are lower and upper censored data. This allows implementing the Gibbs sampling algorithm (S. Chib 1992).
Then,
\[\begin{align} \pi(\boldsymbol{\beta}, \sigma^2, \mathbf{y^*} \mid \mathbf{y}, \mathbf{X}) &\propto \prod_{i=1}^N \left[ \mathbb{1}({y_i = L}) \mathbb{1}({y_{C_i^L}^* < L}) + \mathbb{1}({L \leq y_i < U}) + \mathbb{1}({y_i = U}) \mathbb{1}({y_{C_i^U}^* \geq U}) \right] \\ &\times N(y_i^* \mid \mathbf{x}_i^{\top} \boldsymbol{\beta}, \sigma^2) \times N(\boldsymbol{\beta} \mid \boldsymbol{\beta}_0, \mathbf{B}_0) \times IG(\sigma^2 \mid \alpha_0/2, \delta_0/2) \end{align}\]
The posterior distributions are:
\[\begin{equation} y_{C_i}^* \mid \boldsymbol{\beta}, \sigma^2, \mathbf{y}, \mathbf{X} \sim \begin{Bmatrix} TN_{(-\infty,L)}(\mathbf{x}_i^{\top}\boldsymbol{\beta}, \sigma^2) \ , \ y_i = L \\ TN_{[U, \infty)}(\mathbf{x}_i^{\top}\boldsymbol{\beta}, \sigma^2) \ \ , \ y_i = U \end{Bmatrix} \end{equation}\]
\[\begin{equation} \boldsymbol{\beta} \mid \sigma^2, \mathbf{y}, \mathbf{X} \sim N(\boldsymbol{\beta}_n, \sigma^2 \mathbf{B}_n) \end{equation}\]
\[\begin{equation} \sigma^2 \mid \boldsymbol{\beta}, \mathbf{y}, \mathbf{X} \sim IG(\alpha_n/2, \delta_n/2) \end{equation}\]
where \(\mathbf{B}_n = (\mathbf{B}_0^{-1} + \sigma^{-2} \mathbf{X}^{\top} \mathbf{X})^{-1}\), \(\boldsymbol{\beta}_n = \mathbf{B}_n(\mathbf{B}_0^{-1} \boldsymbol{\beta}_0 + \sigma^{-2} \mathbf{X}^{\top} \mathbf{y}^*)\), \(\alpha_n = \alpha_0 + N\), and \(\delta_n = \delta_0 + (\mathbf{y}^* - \mathbf{X} \boldsymbol{\beta})^{\top}(\mathbf{y}^* - \mathbf{X} \boldsymbol{\beta})\).
Example: The market value of soccer players in Europe continues
We continue with the example of the market value of soccer players from Section 6.1. We specify the same equation but assume that the sample is censored from below, such that we only have information for soccer players whose market value is higher than one million euros. The dependent variable is log(ValueCens), and the left censoring point is 13.82.
The following Algorithm illustrates how to estimate Tobit models in our GUI. Our GUI utilizes the MCMCtobit command from the MCMCpack package.
Algorithm: Tobit models
Select Univariate Models on the top panel
Select Tobit model using the left radio button
Upload the dataset by first selecting whether there is a header in the file, and the kind of separator in the csv file of the dataset (comma, semicolon, or tab). Then, use the Browse button under the Choose File legend. You should see a preview of the dataset
Select MCMC iterations, burn-in, and thinning parameters using the Range sliders
Select dependent and independent variables using the Formula builder table
Click the Build formula button to generate the formula in R syntax. You can modify the formula in the Main equation box using valid arguments of the formula command structure in R
Set the left and right censoring points. To censor above only, specify -Inf in the left censoring box, and to censor below only, specify Inf in the right censoring box
Set the hyperparameters: mean vector, covariance matrix, shape, and scale parameters. This step is not necessary as by default our GUI uses non-informative priors
Click the Go! button
Analyze results
Download posterior chains and diagnostic plots using the Download Posterior Chains and Download Posterior Graphs buttons
We run this application using the same hyperparameters that we set in the example of Section 6.1. All results seem similar to those in the example of linear models. In addition, the posterior chains seem to achieve good diagnostics.
rm(list = ls()); set.seed(010101)
Data <- read.csv("https://raw.githubusercontent.com/besmarter/BSTApp/refs/heads/master/DataApp/1ValueFootballPlayers.csv", sep = ",", header = TRUE, quote = "")
attach(Data)## The following objects are masked from Data (pos = 3):
##
## Age, Age2## The following objects are masked from mydata:
##
## Age, Age2## The following objects are masked from Data (pos = 5):
##
## Age, Age2, Assists, Exp, Exp2, Goals, Goals2, NatTeam, Perf, Perf2,
## Player, Value, ValueCensy <- log(ValueCens)
X <- cbind(1, Perf, Age, Age2, NatTeam, Goals, Exp, Exp2)
k <- dim(X)[2]
N <- dim(X)[1]
# Hyperparameters
d0 <- 0.001; a0 <- 0.001
b0 <- rep(0, k); c0 <- 1000; B0 <- c0*diag(k)
B0i <- solve(B0)
# MCMC parameters
mcmc <- 50000
burnin <- 10000
tot <- mcmc + burnin
thin <- 1
# Posterior distributions using packages: MCMCpack sets the model in terms of the precision matrix
posterior <- MCMCpack::MCMCtobit(y~X-1, b0=b0, B0 = B0i, c0 = a0, d0 = d0, burnin = burnin, mcmc = mcmc, thin = thin, below = 13.82, above = Inf)
summary(coda::mcmc(posterior))##
## Iterations = 1:50000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 50000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X 1.049702 2.639780 1.181e-02 1.676e-02
## XPerf 0.033979 0.004512 2.018e-05 2.288e-05
## XAge 1.024758 0.213319 9.540e-04 1.337e-03
## XAge2 -0.021861 0.004002 1.790e-05 2.546e-05
## XNatTeam 0.847758 0.126113 5.640e-04 6.463e-04
## XGoals 0.010091 0.001649 7.376e-06 7.688e-06
## XExp 0.174196 0.070237 3.141e-04 3.808e-04
## XExp2 -0.005652 0.002977 1.331e-05 1.555e-05
## sigma2 0.982768 0.095944 4.291e-04 6.698e-04
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X -4.192313 -0.713326 1.071600 2.841806 6.1771087
## XPerf 0.025131 0.030948 0.033967 0.036976 0.0428691
## XAge 0.609991 0.880679 1.023709 1.167147 1.4480407
## XAge2 -0.029803 -0.024542 -0.021820 -0.019162 -0.0141026
## XNatTeam 0.603262 0.762408 0.847351 0.932705 1.0950039
## XGoals 0.006877 0.008975 0.010094 0.011193 0.0133502
## XExp 0.037752 0.126722 0.173706 0.221290 0.3127266
## XExp2 -0.011525 -0.007632 -0.005642 -0.003657 0.0001633
## sigma2 0.811752 0.915417 0.976781 1.042729 1.1887701# Gibbs sampling functions
XtX <- t(X)%*%X
PostBeta <- function(Yl, sig2){
Bn <- solve(B0i + sig2^(-1)*XtX)
bn <- Bn%*%(B0i%*%b0 + sig2^(-1)*t(X)%*%Yl)
Beta <- MASS::mvrnorm(1, bn, Bn)
return(Beta)
}
PostYl <- function(Beta, sig2, L, U, i){
Ylmean <- X[i,]%*%Beta
if(y[i] == L){
Yli <- truncnorm::rtruncnorm(1, a = -Inf, b = L, mean = Ylmean, sd = sig2^0.5)
}else{
if(y[i] == U){
Yli <- truncnorm::rtruncnorm(1, a = U, b = Inf, mean = Ylmean, sd = sig2^0.5)
}else{
Yli <- y[i]
}
}
return(Yli)
}
PostSig2 <- function(Beta, Yl){
an <- a0 + length(y)
dn <- d0 + t(Yl - X%*%Beta)%*%(Yl - X%*%Beta)
sig2 <- invgamma::rinvgamma(1, shape = an/2, rate = dn/2)
return(sig2)
}
PostBetas <- matrix(0, mcmc+burnin, k); Beta <- rep(0, k)
PostSigma2 <- rep(0, mcmc+burnin); sig2 <- 1
L <- log(1000000); U <- Inf
# create progress bar in case that you want to see iterations progress
pb <- txtProgressBar(min = 1, max = tot, style = 3)
for(s in 1:tot){
Yl <- sapply(1:N, function(i){PostYl(Beta = Beta, sig2 = sig2, L = L, U = U, i)})
Beta <- PostBeta(Yl = Yl, sig2 = sig2)
sig2 <- PostSig2(Beta = Beta, Yl = Yl)
PostBetas[s,] <- Beta; PostSigma2[s] <- sig2
setTxtProgressBar(pb, s)
}## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%keep <- seq((burnin+1), tot, thin)
PosteriorBetas <- PostBetas[keep,]
colnames(PosteriorBetas) <- c("Intercept", "Perf", "Age", "Age2", "NatTeam", "Goals", "Exp", "Exp2")
summary(coda::mcmc(PosteriorBetas))##
## Iterations = 1:50000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 50000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## Intercept 1.013599 2.639018 1.180e-02 1.680e-02
## Perf 0.033985 0.004520 2.021e-05 2.274e-05
## Age 1.027092 0.213316 9.540e-04 1.341e-03
## Age2 -0.021910 0.003999 1.789e-05 2.555e-05
## NatTeam 0.849940 0.125638 5.619e-04 6.327e-04
## Goals 0.010102 0.001652 7.389e-06 7.747e-06
## Exp 0.175133 0.070469 3.151e-04 3.809e-04
## Exp2 -0.005684 0.002988 1.336e-05 1.561e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## Intercept -4.227386 -0.744950 1.046922 2.811236 6.1216253
## Perf 0.025161 0.030902 0.034004 0.037041 0.0428490
## Age 0.614097 0.881087 1.024736 1.169974 1.4502181
## Age2 -0.029864 -0.024574 -0.021866 -0.019179 -0.0141818
## NatTeam 0.604448 0.765333 0.849641 0.934294 1.0972474
## Goals 0.006839 0.008996 0.010106 0.011223 0.0133239
## Exp 0.036747 0.128063 0.174950 0.222325 0.3149495
## Exp2 -0.011549 -0.007680 -0.005673 -0.003692 0.0002033##
## Iterations = 1:50000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 50000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.9855375 0.0957929 0.0004284 0.0006728
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 0.8154 0.9183 0.9793 1.0468 1.1902References
We can set \(L\) or \(U\) equal to \(-\infty\) or \(\infty\) to model data censored on just one side.↩︎