0

I am currently trying to write a code in R that estimates an ordered probit model via the maximum likelihood estimation. The dependent variable has three categories. There are several independent variables that are continuous as well as binomial. This is code:

# Define the likelihood function for the ordered probit model
ordered_probit_likelihood <- function(params, X, y, n_cats) {
  beta <- params[1:(length(params) - n_cats + 1)]
  thresholds <- c(-Inf, params[(length(params) - n_cats + 2):length(params)], Inf)
  
  p_lower <- pnorm((thresholds[1:(length(thresholds) - 1)] - X %*% beta))
  p_upper <- pnorm((thresholds[2:length(thresholds)] - X %*% beta))
  
  probs <- p_upper - p_lower
  ll <- sum(log(probs[y == 1])) + sum(log(1 - probs[y == 0]))
  
  return(-ll)
}

# Custom optimization function
ordered_probit_fit <- function(X, y, n_cats) {
  n_params <- ncol(X) + n_cats - 1
  init_params <- runif(n_params)
  
  fit <- optim(
    par = init_params,
    fn = ordered_probit_likelihood,
    method = "BFGS",
    control = list(maxit = 1000),
    X = X,
    y = y,
    n_cats = n_cats
  )
  
  return(list(coef = fit$par, logLik = -fit$value))
}

# Fit the model
model <- ordered_probit_fit(X, y, n_cats)

And this is a data generating process:

# Simulate data
set.seed(123)

n <- 500
n_cats <- 3

x <- cbind(1, rnorm(n * 2), rbinom(n, size = 1, prob = 0.5))
x <- as.matrix(x)
Xbeta <- c(-1, 1, 0.5)
thresholds <- c(-1, 2)
y_star <- x %*% Xbeta + rnorm(n)
y <- findInterval(y_star, vec = thresholds) + 1

This creates the following error: Error in optim(par = init_params, fn = ordered_probit_likelihood, method = "BFGS", : initial value in 'vmmin' is not finite

As soon as I change the size for rbinom to 2, x <- cbind(1, rnorm(n * 2), rbinom(n, **size = 2**, prob = 0.5)), this error does not come up anymore.

Has anyone the idea, why this happens? Thank you for your help.

Anna Mogilevskaja
  • 1
  • 3
  • When you hit this kind of error, it can be useful to step through your objective function at the initial parameter values to see what's going on, as in [this example](https://stackoverflow.com/questions/76020993/problem-specifying-arguments-when-using-optim-for-mle/76023536#76023536) ... – Ben Bolker Apr 18 '23 at 21:40

0 Answers0