Estimation d'un modèle de régression probit avec optim()

Question

Je dois programmer manuellement un modèle de régression probit sans utiliser glm. J'utiliserais optim pour la minimisation directe de log-vraisemblance négative.

J'écrit le code ci-dessous, mais il ne fonctionne pas, ce qui donne une erreur:

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) 

données: https://drive.google.com/file/d/0B06Id6VJyeb5OTFjbHVHUE42THc/view?usp=sharing

Answer 1

sensible à la casse

Y et y sont différents. Vous devez donc utiliser Y et non y dans vos fonctions définies probit.nll et probit.gr.

Ces deux fonctions ne me semblent pas correctes non plus. Le problème le plus évident est l'existence de rnorm. Les suivants sont corrects.

fonction de log-vraisemblance négative

# 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)) 
    } 

fonction gradient

# 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) 
    } 

valeurs initiales des paramètres de lm()

Cela ne ressemble pas à un idée raisonnable. En aucun cas, nous ne devons appliquer une régression linéaire aux données binaires.

Toutefois, en se concentrant uniquement sur l'utilisation de lm, vous avez besoin de + et non de , pour séparer les covariables dans le côté droit de la formule.

---

exemple reproductible

Produisons un jeu de données de jouet

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 

Ils sont très proches les uns des autres!