r - 使用 optim() 估计概率回归模型

Question

我需要在不使用 glm 的情况下手动编写概率回归模型。我会使用 optim 直接最小化负对数似然。

我在下面编写了代码，但它不起作用，并给出错误:

cannot coerce type 'closure' to vector of type 'double'

# load data: data provided via the bottom link
Datospregunta2a <- read.dta("problema2_1.dta")
attach(Datospregunta2a)

# model matrix `X` and response `Y`
X <- cbind(1, associate_professor, full_professor, emeritus_professor, other_rank)
Y <- volunteer

# number of regression coefficients
K <- ncol(X)

# initial guess on coefficients
vi <- lm(volunteer ~ associate_professor, full_professor, emeritus_professor, other_rank)$coefficients

# negative log-likelihood
probit.nll <- function (beta) {
  exb <- exp(X%*%beta)
  prob<- rnorm(exb)
  logexb <- log(prob)
  y0 <- (1-y)
  logexb0 <- log(1-prob)
  yt <- t(y)
  y0t <- t(y0)
  -sum(yt%*%logexb + y0t%*%logexb0)
  }

# gradient
probit.gr <- function (beta) {
  grad <- numeric(K)
  exb <- exp(X%*%beta)
  prob <- rnorm(exb)
  for (k in 1:K) grad[k] <- sum(X[,k]*(y - prob))
  return(-grad)
  }

# direct minimization
fit <- optim(vi, probit.nll, gr = probit.gr, method = "BFGS", hessian =  TRUE)

数据:https://drive.google.com/file/d/0B06Id6VJyeb5OTFjbHVHUE42THc/view?usp=sharing

Answer 1

区分大小写

Y 和 y 是不同的。因此，您应该在定义的函数 probit.nll 和 probit.gr 中使用 Y 而不是 y。

这两个函数在我看来也不正确。最明显的问题是 rnorm 的存在。以下是正确的。

负对数似然函数

# requires model matrix `X` and binary response `Y`
probit.nll <- function (beta) {
  # linear predictor
  eta <- X %*% beta
  # probability
  p <- pnorm(eta)
  # negative log-likelihood
  -sum((1 - Y) * log(1 - p) + Y * log(p))
  }

梯度函数

# requires model matrix `X` and binary response `Y`
probit.gr <- function (beta) {
  # linear predictor
  eta <- X %*% beta
  # probability
  p <- pnorm(eta)
  # chain rule
  u <- dnorm(eta) * (Y - p) / (p * (1 - p))
  # gradient
  -crossprod(X, u)
  }

来自lm()的初始参数值

这听起来不像是一个合理的想法。在任何情况下，我们都不应该对二进制数据应用线性回归。

但是，纯粹关注 lm 的使用，您需要 + 而不是 , 来分隔公式右侧的协变量.

---

可重现的示例

让我们生成一个玩具数据集

set.seed(0)
# model matrix
X <- cbind(1, matrix(runif(300, -2, 1), 100))
# coefficients
b <- runif(4) 
# response
Y <- rbinom(100, 1, pnorm(X %*% b))

# `glm` estimate
GLM <- glm(Y ~ X - 1, family = binomial(link = "probit"))

# our own estimation via `optim`
# I am using `b` as initial parameter values (being lazy)
fit <- optim(b, probit.nll, gr = probit.gr, method = "BFGS", hessian = TRUE)

# comparison
unname(coef(GLM))
# 0.62183195  0.38971121  0.06321124  0.44199523

fit$par
# 0.62183540  0.38971287  0.06321318  0.44199659

他们的关系很亲密!