heckman 模型的预测函数答案

【问题标题】：Predict function for heckman modelheckman 模型的预测函数
【发布时间】：2012-12-22 19:19:21
【问题描述】：

我使用 sampleSelection 包中的示例

## Greene( 2003 ): example 22.8, page 786
    data( Mroz87 )
    Mroz87$kids  <- ( Mroz87$kids5 + Mroz87$kids618 > 0 )
    # Two-step estimation
    test1  = heckit( lfp ~ age + I( age^2 ) + faminc + kids + educ,
                     wage ~ exper + I( exper^2 ) + educ + city, Mroz87 ) 
    # ML estimation
    test2 =  selection( lfp ~ age + I( age^2 ) + faminc + kids + educ,
                        wage ~ exper + I( exper^2 ) + educ + city, Mroz87 ) 
    pr2 <- predict(test2,Mroz87)
    pr1 <- predict(test1,Mroz87)

我的问题是预测功能不起作用。我收到此错误：

    Error in UseMethod("predict") : 
  no applicable method for 'predict' applied to an object of class "c('selection', 'maxLik', 'maxim', 'list')"

predict 函数适用于许多模型，所以我想知道为什么赫克曼回归模型会出错。

-----------更新----------- 我取得了一些进展，但我仍然需要你的帮助。我建立了一个原始的 heckman 模型进行比较：

data( Mroz87 )
Mroz87$kids  <- ( Mroz87$kids5 + Mroz87$kids618 > 0 )
test1  = heckit( lfp ~ age + I( age^2 ) + faminc + kids + educ,
                 wage ~ exper + I( exper^2 ) + educ + city, Mroz87[1:600,] )

之后我开始自己构建它。 Heckman 模型需要一个选择方程：

zi* = wi γ + ui
where zi =1 if zi* >0  and zi = 0 if zi* <=0
after you calculate yi = xi*beta +ei ONLY for the cases where zi*>0

我先建立概率模型：

library(MASS)
#probit1 = probit(lfp ~ age + I( age^2 ) + faminc + kids + educ, Mroz87, x = TRUE, print.level = print.level -      1, iterlim = 30)
myprobit <- glm(lfp ~ age + I( age^2 ) + faminc + kids + educ, family = binomial(link = "probit"), 
                data = Mroz87[1:600,])
summary(myprobit)

模型与 heckit 命令完全相同。

然后我建立一个lm模型：

#get predictions for the variables (the data is not needed but I specify it anyway)
selectvar <- predict(myprobit,data = Mroz87[1:600,])
#bind the prediction to the table (I build a new one in my case)
newdata = cbind(Mroz87[1:600,],selectvar)
#Build an lm model for the subset where zi>0
lm1 = lm(wage ~ exper + I( exper^2 ) + educ + city , newdata, subset = selectvar > 0)
summary(lm1)

我现在的问题是 lm 模型与 heckit 创建的模型相比并不多。我不知道为什么。有什么想法吗？

【问题讨论】：

没有predict.selection函数。
感谢 DWin。有什么办法可以让我自己做吗？也许将概率用于选择模型，然后将其与 lm 函数一起使用？然后预测应该起作用。所以基本上我会使用 probit 来生成模型，然后对这个模型使用 predict 并在创建 lm 模型时将其作为变量应用。这行得通吗？
预测选择对象中包含的lm对象是否有意义，即predict(test2$twoStep$lm)？
它是这样工作的，但是当我使用模型来预测新数据时它仍然是一个问题。例如，表中的前 10 行是我们的测试样本。如果我要运行 predict(test2$twoStep$lm,Mroz87[1:10,])，那么我会收到警告 Warning message: 'newdata' had 10 rows but variable(s) found have 753 rows ，它会返回 753 行而不是我请求的 10 行。

标签： r

【解决方案1】：

实施

这里是 predict.selection 函数的实现——它产生 4 种不同类型的预测（解释 here）：

library(Formula)
library(sampleSelection)
predict.selection = function(objSelection, dfPred, 
                                 type = c('link', 'prob', 'cond', 'uncond')) {

  # construct the Formula object
  tempS = evalq(objSelection$call$selection)
  tempO = evalq(objSelection$call$outcome)

  FormHeck = as.Formula(paste0(tempO[2], '|', tempS[2], '~', tempO[3], '|', tempS[3]))

  # regressor matrix for the selection equation
  mXSelection = model.matrix(FormHeck, data = dfPred, rhs = 2)

  # regressor matrix for the outcome equation
  mXOutcome = model.matrix(FormHeck, data = dfPred, rhs = 1)

  # indices of the various parameters in selectionObject$estimate
  vIndexBetaS = objSelection$param$index$betaS
  vIndexBetaO = objSelection$param$index$betaO
  vIndexErr = objSelection$param$index$errTerms

  # get the estimates
  vBetaS = objSelection$estimate[vIndexBetaS]
  vBetaO = objSelection$estimate[vIndexBetaO]

  dLambda = objSelection$estimate[vIndexErr['rho']]*
    objSelection$estimate[vIndexErr['sigma']]

  # depending on the type of prediction requested, return
  # TODO allow the return of multiple prediction types
  pred = switch(type, 
         link = mXSelection %*% vBetaS,
         prob = pnorm(mXSelection %*% vBetaS),
         uncond = mXOutcome %*% vBetaO,
         cond = mXOutcome %*% vBetaO + 
           dnorm(temp <- mXSelection %*% vBetaS)/pnorm(temp) * dLambda)
  return(pred)
}

测试

假设您使用 MLE 估计以下 Heckman 样本选择模型：

data(Mroz87)

# define a new variable
Mroz87$kids  = (Mroz87$kids5 + Mroz87$kids618 > 0)

# create the estimation sample
Mroz87Est = Mroz87[1:600, ]

# create the hold out sample
Mroz87Holdout = Mroz87[601:nrow(Mroz87), ]

# estimate the model using MLE
heckML =  selection(selection = lfp ~ age + I(age^2) + faminc + kids + educ,
                    outcome = wage ~ exper + I(exper^2) + educ + city, data = Mroz87Est) 
summary(heckML)

不同类型的预测计算如下：

vProb = predict(objSelection = heckML, dfPred = Mroz87Holdout, type = 'prob')
vLink = predict(objSelection = heckML, dfPred = Mroz87Holdout, type = 'link')
vCond = predict(objSelection = heckML, dfPred = Mroz87Holdout, type = 'cond')
vUncond = predict(objSelection = heckML, dfPred = Mroz87Holdout, type = 'uncond')

您可以在产生这些输出的平台上验证这些计算，例如Stata。

【讨论】：