regressCens computes estimates of regression parameters under the censored (Tobit) model
The model is estimated by Maximum Likelihood (ML) assuming a Gaussian (normal) distribution of the error term. The maximization of the likelihood function is done by function fminunc of the optimization toolbox.
Example of right censoring.out
=regressCens(y
,
X
,
Name, Value
)
In the example of Kleiber and Zeileis (2008, p. 142), the number of a person's extramarital sexual inter-courses ("affairs") in the past year is regressed on the person's age, number of years married, religiousness, occupation, and won rating of the marriage. The dependent variable is left-censored at zero and not right-censored. Hence this is a standard Tobit model which can be estimated by the following lines
load affairs.mat X=affairs(:,["age" "yearsmarried" "religiousness" "occupation" "rating"]); y=affairs(:,"affairs"); out=regressCens(y,X,"dispresults",true);
Observations: Total LeftCensored Uncensored RightCensored _____ ____________ __________ _____________ 601 451 150 0 Coefficients betaout sebetaout tout pval ________ _________ _______ __________ (Intercept) 8.176 2.7425 2.9812 0.0029883 age -0.17937 0.079165 -2.2657 0.023828 yearsmarried 0.55418 0.13463 4.1162 4.3959e-05 religiousness -1.6864 0.40385 -4.1757 3.4159e-05 occupation 0.32603 0.25444 1.2814 0.20056 rating -2.2851 0.40801 -5.6006 3.2682e-08 sigma 8.2473 0.55409 14.884 8.1345e-43 Number of observations: 601, Error degrees of freedom:594 Log-likelihood: -705.5762
When left=-Inf and right=0 indicates that there is no left-censoring but there is a right censoring at 0. The same model as above but with the negative number of affairs as teh dependent variable can be estimated by
load affairs.mat X=affairs(:,["age" "yearsmarried" "religiousness" "occupation" "rating"]); y=affairs(:,"affairs").*(-1); out=regressCens(y,X,"dispresults",true,"left",-Inf,"right",0); % This estimation returns beta parameters that have the opposite sign of % the beta parameters in the original model, but the estimate of sigma does % not change.
Observations: Total LeftCensored Uncensored RightCensored _____ ____________ __________ _____________ 601 0 150 451 Coefficients betaout sebetaout tout pval ________ _________ _______ __________ (Intercept) -8.176 2.7425 -2.9812 0.0029883 age 0.17937 0.079165 2.2657 0.023828 yearsmarried -0.55418 0.13463 -4.1162 4.3959e-05 religiousness 1.6864 0.40385 4.1757 3.4159e-05 occupation -0.32603 0.25444 -1.2814 0.20056 rating 2.2851 0.40801 5.6006 3.2682e-08 sigma 8.2473 0.55409 14.884 8.1345e-43 Number of observations: 601, Error degrees of freedom:594 Log-likelihood: -705.5762
When left=-Inf and right=0 indicates that there is no left-censoring but there is a right censoring at 0. The same model as above but with the negative number of affairs as teh dependent variable can be estimated by
load affairs.mat X=affairs(:,["age" "yearsmarried" "religiousness" "occupation" "rating"]); y=affairs(:,"affairs").*(-1); out=regressCens(y,X,"dispresults",true,"left",-Inf,"right",0); % This estimation returns beta parameters that have the opposite sign of % the beta parameters in the original model, but the estimate of sigma does % not change.
Observations: Total LeftCensored Uncensored RightCensored _____ ____________ __________ _____________ 601 0 150 451 Coefficients betaout sebetaout tout pval ________ _________ _______ __________ (Intercept) -8.176 2.7425 -2.9812 0.0029883 age 0.17937 0.079165 2.2657 0.023828 yearsmarried -0.55418 0.13463 -4.1162 4.3959e-05 religiousness 1.6864 0.40385 4.1757 3.4159e-05 occupation -0.32603 0.25444 -1.2814 0.20056 rating 2.2851 0.40801 5.6006 3.2682e-08 sigma 8.2473 0.55409 14.884 8.1345e-43 Number of observations: 601, Error degrees of freedom:594 Log-likelihood: -705.5762
This example is taken from https://stats.oarc.ucla.edu/r/dae/tobit-models/ Consider the situation in which we have a measure of academic aptitude (scaled 200-800) which we want to model using reading and math test scores, as well as, the type of program the student is enrolled in (academic, general, or vocational). The problem here is that students who answer all questions on the academic aptitude test correctly receive a score of 800, even though it is likely that these students are not “truly” equal in aptitude. The same is true of students who answer all of the questions incorrectly. All such students would have a score of 200, although they may not all be of equal aptitude.
XX=readtable("https://stats.idre.ucla.edu/stat/data/tobit.csv","ReadRowNames",true); XX.prog=categorical(XX.prog); % The dataset contains 200 observations. The academic aptitude variable is % "apt", the reading and math test scores are read and math respectively. The % variable prog is the type of program the student is in, it is a % categorical (nominal) variable that takes on three values, academic % general, and vocational (prog = 3). % The scatterplot matrix is shown below spmplot(XX(:,[1 2 4])); % Now let’s look at the data descriptively. Note that in this dataset, the % lowest value of apt is 352. That is, no students received a score of 200 % (the lowest score possible), meaning that even though censoring from % below was possible, it does not occur in the dataset. summary(XX) % Define y and X y=XX(:,"apt"); X=XX(:,["read", "math" "prog"]); % Call regressCens out=regressCens(y,X,'right',800,'left',-Inf,'dispresults',true); % Show the plot of fitted vs residuals fitted=out.X*out.Beta(1:end-1,1); residuals=(y{:,1}-fitted)/out.Beta(end,1); figure scatter(fitted,residuals) xlabel('Fitted values (yhat)') ylabel('Residuals')
Variables: read: 200×1 double Values: Min 28 Median 50 Max 76 math: 200×1 double Values: Min 33 Median 52 Max 75 prog: 200×1 categorical Values: academic 45 general 105 vocational 50 apt: 200×1 double Values: Min 352 Median 633 Max 800 Observations: Total LeftCensored Uncensored RightCensored _____ ____________ __________ _____________ 200 0 183 17 Coefficients betaout sebetaout tout pval _______ _________ _______ __________ (Intercept) 209.57 32.774 6.3943 1.1763e-09 read 2.6979 0.61879 4.36 2.1086e-05 math 5.9145 0.7099 8.3314 1.4325e-14 prog_general -12.715 12.406 -1.0249 0.3067 prog_vocational -46.144 13.724 -3.3623 0.00093125 sigma 65.677 3.4827 18.858 9.5605e-46 Number of observations: 200, Error degrees of freedom:194 Log-likelihood: -1041.0629
y
— Response variable.
Vector.Response variable, specified as a vector of length n, where n is the number of observations. Each entry in y is the response for the corresponding row of X.
Data Types: array or table
X
— Predictor variables in the regression equation.
Matrix.Matrix of explanatory variables (also called 'regressors') of dimension n x (p-1) where p denotes the number of explanatory variables including the intercept. Rows of X represent observations, and columns represent variables. By default, there is a constant term in the model, unless you explicitly remove it using input option intercept, so do not include a column of 1s in X. Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows) with missing or infinite values will automatically be excluded from the computations.
Data Types: array or table
Specify optional comma-separated pairs of Name,Value
arguments.
Name
is the argument name and Value
is the corresponding value. Name
must appear
inside single quotes (' '
).
You can specify several name and value pair arguments in any order as
Name1,Value1,...,NameN,ValueN
.
'bsb',[3 5 20:30]
, 'dispresults',true
, 'left',1
, 'initialbeta',[3 8]
, 'intercept',false
, 'nocheck',false
, 'right',800
,'optmin.Display','off'
bsb
—units forming subset.vector.m x 1 vector.
The default value of bsb is 1:length(y), that is all units are used to compute parameter estimates.
Example: 'bsb',[3 5 20:30]
Data Types: double
dispresults
—Display results of final fit.boolean.If dispresults is true, labels of coefficients, estimated coefficients, standard errors, tstat and p-values are shown on the screen in a fully formatted way. The default value of dispresults is false.
Example: 'dispresults',true
Data Types: logical
left
—left limit for the censored dependent variable.scalar.If set to -Inf, the dependent variable is assumed to be not left-censored; default value of left is zero (classical Tobit model).
Example: 'left',1
Data Types: double
initialbeta
—initial estimate of beta.vector.p x 1 vector. If initialbeta is not supplied (default) standard least squares is used to find initial estimate of beta
Example: 'initialbeta',[3 8]
Data Types: double
intercept
—Indicator for constant term.true (default) | false.Indicator for the constant term (intercept) in the fit, specified as the comma-separated pair consisting of 'Intercept' and either true to include or false to remove the constant term from the model.
Example: 'intercept',false
Data Types: boolean
nocheck
—Check input arguments.boolean.If nocheck is equal to true no check is performed on supplied structure model
Example: 'nocheck',false
Data Types: logical
right
—right limit for the censored dependent variable.scalar.If set to Inf, the dependent variable is assumed to be not right-censored; default value of left is Inf (classical Tobit model).
Example: 'right',800
Data Types: double
optmin
—It contains the options dealing with the
maximization algorithm.structure.Use optimset to set these options.
Notice that the maximization algorithm which is used is fminunc if the optimization toolbox is present else is fminsearch.
Example: 'optmin.Display','off'
Data Types: double
out
— description
StructureStructure which contains the following fields
Value | Description |
---|---|
Beta |
p-by-3 matrix containing: 1st col = Estimates of regression coefficients; 2nd col = Standard errors of the estimates of regression coefficients; 3rd col = t-tests of the estimates of regression coefficients; |
LogL |
scalar. Value of the maximized log likelihood (ignoring constants) |
Exflag |
Reason fminunc stopped. Integer. out.Exflag is equal to 1 if the maximization procedure did not produce warnings or the warning was different from "ILL Conditiioned Jacobian". For any other warning which is produced (for example, "Overparameterized", "IterationLimitExceeded", 'MATLAB:rankDeficientMatrix") out.Exflag is equal to -1; |
X |
design matrix which has been used. This output is present just in X is a table containing categorical variables. |
In the standard Tobit model (Tobin 1958), we have a dependent variable $y$ that is left-censored at zero: \begin{eqnarray} y_i^* & = & x_i^{\prime} \beta+\varepsilon_i \\ y_i = & 0 & \text { if } y_i^* \leq 0 \\ y_i = & y_i^* & \text { if } y_i^*>0 \end{eqnarray} Here the subscript $i=1, \ldots, n$ indicates the observation, $y_i^*$ is an unobserved ("latent") variable, $x_i$ is a vector of explanatory variables, $\beta$ is a vector of unknown parameters, and $\varepsilon_i$ is a disturbance term.
The censored regression model is a generalisation of the standard Tobit model.
The dependent variable can be either left-censored, right-censored, or both left-censored and right-censored, where the lower and/or upper limit of the dependent variable can be any number: \begin{eqnarray} y_i^*=x_i^{\prime} \beta+\varepsilon_i \\ y_i & = & a \qquad \text { if } y_i^* \leq a \\ & = & y_i^* \qquad \text { if } a<y_i^*<b \\ & = & b \qquad \text { if } y_i^* \geq b \end{eqnarray} Here $a$ is the lower limit and $b$ is the upper limit of the dependent variable. If $a=-\infty$ or $b=\infty$, the dependent variable is not left-censored or right-censored, respectively.
Censored regression models (including the standard Tobit model) are usually estimated by the Maximum Likelihood (ML) method. Assuming that the disturbance term $\varepsilon$ follows a normal distribution with mean 0 and variance $\sigma^2$, the log-likelihood function is \begin{aligned} \log L=\sum_{i=1}^N & \left[ I_i^a \log \Phi\left(\frac{a-x_i^{\prime} \beta}{\sigma}\right)+I_i^b \log \Phi\left(\frac{x_i^{\prime} \beta-b}{\sigma}\right) \right. \\ & \left.+\left(1-I_i^a-I_i^b\right)\left(\log \phi\left(\frac{y_i-x_i^{\prime} \beta}{\sigma}\right)-\log \sigma\right)\right], \end{aligned} where $\phi( \cdot)$ and $\Phi( \cdot )$ denote the probability density function and the cumulative distribution function, respectively, of the standard normal distribution, and $I_i^a$ and $I_i^b$ are indicator functions with \begin{eqnarray} I_i^a & = & 1 \text { if } y_i=a \\ & = & 0 \text { if } y_i>a \\ I_i^b & = & 1 \text { if } y_i=b \\ & = & 0 \text { if } y_i<b \\ \end{eqnarray} In this file the censored log likelihood above is maximized with respect to the parameter vector (\beta' \sigma) using routine fminunc of the optimization toolbox
Greene, W.H. (2008), "Econometric Analysis, Sixth Edition", Prentice Hall, pp. 871-875.
Henningsen, A. (2012), Estimating Censored Regression Models in R using the censReg Package, [https://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf
Kleiber C., Zeileis A. (2008), "Applied Econometrics with R", Springer, New York.
Tobin, J. (1958), Estimation of Relationships for Limited Dependent Variables, "Econometrica", 26, pp. 24-36.