受限因变量（limited dependent variable）是指因变量的观测值是连续的，但是受到某种限制，其抽样并非完全随机的，得到的观测值不能完全反应因变量的真实情况。选择性样本（selective sample）是受限因变量的主要形式，其样本观测值是在选择性限制的情况下抽取的。选择性样本理论与方法发展于20世纪70年代，Heckman做了很多基础性的工作，近几年选择性样本模型已经成为微观统计学的主要研究内容之一，并得到了广泛应用，尤其是在劳动经济学和卫生经济学等领域。受限因变量常见的两类数据：截断（truncation）数据和审查（Censoring）数据。受限因变量模型主要包括截断模型（Truncated Model）和审查模型（Censoring Model）两类。这两类模型多应用在调查数据的分析中。
受限因变量模型主要包括截断模型（Truncated Model）和审查模型（Censoring Model）两类。
截断模型问题的提出：
截断数据即不能从全部个体，而只能从一部分个体中随机抽取因变量的样本观测值。比如，我们对农村居民在过去一年中的家庭医疗支出进行随机调查，在所调查的8000户家庭里，其中有2000户在一年内有发生医疗支出，而剩余的6000户没有医疗支出，如果以2000户的样本进行分析，则显然没有医疗支出的6000户被截断了。
审查模型问题的提出：
在过去一年中，您家里人住院的次数有几次？
A.0次
B.1次
C.2次
D.3次
E.4次
F.5次及以上
该问题中，家里人住院的次数应该是一个计数数据，如果以该变量作为因变量可以使用possion回归，但是此处特殊的地方是当住院次数在5次以上时被归并（censor）在5次，我们观察不到6次或6次以上的值，这就是一种审查数据（censored data）。该如何分析这样的数据？
审查模型—婚外性行为的研究
Fair（1978）研究了婚外性行为的问题。该数据是从大约2000份回收的电子问卷中抽取的601份初婚且未离婚的样本对象。有关初始数据前10个观测值如下。
ffairs表示过去几年中婚外性行为的次数，因变量affairs，取值0、1、2、3、7和12六种取值。其中0为过去几年中没有婚外性行为，1、2、3分别为有1、2和3次婚外性行为；当婚外性行为的次数在4-10之间时，affairs定义为7；当婚外性行为频率较高，为每周，每月或每日时，affairs定义为12；gender表示性别虚拟变量，1为男性，0为女性；age表示年龄；Yearsmarried表示结婚年数；children表示有否孩子，有为1，无为0；religiousness表示宗教信仰，5种程度，从1为反对宗教信仰到5为非常相信宗教；education表示教育年限，有9（初中）、12（高中）、14、16、17、18和20（博士或其他）共7种结果；occupation表示职业，根据Hollingshead职业量表分为1-7种职业。rating表示自我对婚姻的评价，1-5个级别，1为非常不高兴，5为非常幸福。该如何建模分析哪些因素会影响婚外性行为的次数？
mdata <- read.csv('http://data.galaxystatistics.com/blog_data/regression/affairs_data_transform.csv',header=T)
head(mdata)## affairs age yearsmarried religiousness occupation rating
## 1 0 37 10.00 3 7 4
## 2 0 27 4.00 4 6 4
## 3 0 32 15.00 1 1 4
## 4 0 57 15.00 5 6 5
## 5 0 22 0.75 2 6 3
## 6 0 32 1.50 2 5 5#描述性统计计算
modeldata <- mdata
summary(modeldata)## affairs age yearsmarried religiousness
## Min. : 0.000 Min. :17.50 Min. : 0.125 Min. :1.000
## 1st Qu.: 0.000 1st Qu.:27.00 1st Qu.: 4.000 1st Qu.:2.000
## Median : 0.000 Median :32.00 Median : 7.000 Median :3.000
## Mean : 1.456 Mean :32.49 Mean : 8.178 Mean :3.116
## 3rd Qu.: 0.000 3rd Qu.:37.00 3rd Qu.:15.000 3rd Qu.:4.000
## Max. :12.000 Max. :57.00 Max. :15.000 Max. :5.000
## occupation rating
## Min. :1.000 Min. :1.000
## 1st Qu.:3.000 1st Qu.:3.000
## Median :5.000 Median :4.000
## Mean :4.195 Mean :3.932
## 3rd Qu.:6.000 3rd Qu.:5.000
## Max. :7.000 Max. :5.000#数据集字段名称转换
renames <- function(x) {
ocol <- ncol(x)
if (ocol==2) {
xnam <- paste0("X")
names(x) <- c('Y', xnam)
}else if(ocol>2) {
xnam <- paste0("X", 1:(ocol-1))
names(x) <- c('Y', xnam)
}
return(x)
}
modeldata <- renames(modeldata)
head(modeldata)## Y X1 X2 X3 X4 X5
## 1 0 37 10.00 3 7 4
## 2 0 27 4.00 4 6 4
## 3 0 32 15.00 1 1 4
## 4 0 57 15.00 5 6 5
## 5 0 22 0.75 2 6 3
## 6 0 32 1.50 2 5 5class(modeldata)## [1] "data.frame"dim(modeldata)## [1] 601 6regFormula <- as.formula(paste("Y ~ ", paste(names(modeldata)[-1], collapse= "+")))
regFormula## Y ~ X1 + X2 + X3 + X4 + X5因变量Y的描述统计
| 取值 | 计数 | 百分比 | 累计数 | 累计百分比 | |
|---|---|---|---|---|---|
| 0 | 451 | 75.04 | 451 | 75.04 | |
| 1 | 34 | 5.66 | 485 | 80.70 | |
| 2 | 17 | 2.83 | 502 | 83.53 | |
| 3 | 19 | 3.16 | 521 | 86.69 | |
| 7 | 42 | 6.99 | 563 | 93.68 | |
| 12 | 38 | 6.32 | 601 | 100.00 | |
| Total | 601 | 100.00 | 601 | 100.00 |
从上表中可以看到，Y中取值为0的样本有451个，占了总样本的75%，而Y大于0的样本只有25%左右。因此，如果用传统的OLS法直接对线性回归模型进行估计，则很难确拟合Y的值。
#####拟合估算：方式方法1#####
###建立Tobit模型可以使用censReg包里的censReg()，使用方法和结果如下。
#install.packages('censReg')
library(censReg)
Tobit.fit <- censReg(regFormula, data=modeldata)
#展示拟合模型的详细结果
summary(Tobit.fit)##
## Call:
## censReg(formula = regFormula, data = modeldata)
##
## Observations:
## Total Left-censored Uncensored Right-censored
## 601 451 150 0
##
## Coefficients:
## Estimate Std. error t value Pr(> t)
## (Intercept) 8.17420 2.74145 2.982 0.00287 **
## X1 -0.17933 0.07909 -2.267 0.02337 *
## X2 0.55414 0.13452 4.119 3.80e-05 ***
## X3 -1.68622 0.40375 -4.176 2.96e-05 ***
## X4 0.32605 0.25442 1.282 0.20001
## X5 -2.28497 0.40783 -5.603 2.11e-08 ***
## logSigma 2.10986 0.06710 31.444 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Newton-Raphson maximisation, 7 iterations
## Return code 1: gradient close to zero
## Log-likelihood: -705.5762 on 7 Df#列出拟合模型的模型参数
coefficients(Tobit.fit)## (Intercept) X1 X2 X3 X4 X5
## 8.1741974 -0.1793326 0.5541418 -1.6862205 0.3260532 -2.2849727
## logSigma
## 2.1098592# predict(Tobit.fit, modeldata[-1])
#####拟合估算：方式方法2#####
###Tobit模型也可以用AER包里的tobit()函数，使用方法和结果如下。
#install.packages('AER')
library(AER)
Tobit.fit2 <- tobit(regFormula, data=modeldata)
#展示拟合模型的详细结果
summary(Tobit.fit2)##
## Call:
## tobit(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = modeldata)
##
## Observations:
## Total Left-censored Uncensored Right-censored
## 601 451 150 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.17420 2.74145 2.982 0.00287 **
## X1 -0.17933 0.07909 -2.267 0.02337 *
## X2 0.55414 0.13452 4.119 3.80e-05 ***
## X3 -1.68622 0.40375 -4.176 2.96e-05 ***
## X4 0.32605 0.25442 1.282 0.20001
## X5 -2.28497 0.40783 -5.603 2.11e-08 ***
## Log(scale) 2.10986 0.06710 31.444 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Scale: 8.247
##
## Gaussian distribution
## Number of Newton-Raphson Iterations: 4
## Log-likelihood: -705.6 on 7 Df
## Wald-statistic: 67.71 on 5 Df, p-value: 3.0718e-13#列出拟合模型的模型参数
coefficients(Tobit.fit2) ## (Intercept) X1 X2 X3 X4 X5
## 8.1741974 -0.1793326 0.5541418 -1.6862205 0.3260532 -2.2849727#predict(Tobit.fit2, modeldata[-1])