How to confine the alpha in scobit model to only positive value when working in R?

Question

The R reference is just for clarity, the question is more of stats theory. I have the following codes to set up the log-likelihood scobit function:

LLscobit <- function(theta, y, X) {

  beta <- theta[1:ncol(X)]
  alpha <- theta[ncol(X)+1]

  mu <- X %*% beta

  alpha <- exp(alpha)

  p <- 1/((1+exp(-mu))^alpha)

  ll <- y * log(p) + (1 - y) * log(1 - p)
  ll <- sum(ll)
  return(ll)
}

Initially my function didn't have the alpha <- exp(alpha) part (for those not knowing R, what it does is it transforms all values of alpha into their exponential values), then when I optim the function it gives off "NaNs produced" message.

startvalsscobit <- c(rep(0, ncol(X)), 1)

optim(par = startvalsscobit, fn = LLscobit, y = y, X = X,
              control = list(fnscale = -1),
              hessian = TRUE,
              method = "BFGS"
              )

> NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

I then figured out that it was because I had to constrain alpha to having only positive value (as defined in the scobit model). Therefore I added that part to take the exponential of alpha to make it always positive and the function worked.

However, I'm still quite skeptical about whether doing that would alter the result and make it wrong, since now the value of alpha that is used to calculate Pr(Y=1) is different? Or it doesn't matter because we only care about the maximum likelihood point, and taking log or taking exponential alters the value of P but not the maximum point? I also tried to transform alpha into absolute value with alpha <- abs(alpha), and it then gave different results for beta, which means one is more correct than the other?

Answer 1

No it won't, remember that your optim model aims to maximize directly the theta but not the alpha which is a transformed parameter of theta itself. And you are right about the transformed alpha not really mattering too much as it only converts the true probability result of the model but does not change the maximum point.

I would suggest though that you start running the optim from the points where the scobit model resembles to a logit model, that is alpha equal to 1. In order to do this, you need to set the last theta to 0 so the exponential of it is 1, therefore:

startvalsscobit <- c(rep(0, ncol(X)), 0)

What I would normally do though is to optim a logit model of the dataset first, and then use the maximum point of the logit model as starting points to run a scobit model and see if the scobit model would provide better results and better flexibility than the logit.

startvalsscobit <- c(logitmodel$par, 0)

I would suggest that you transform theta using exponential rather than absolute value, since abs() only changes the negative values of theta and does not change the positive ones, which makes it somewhat biased to produce a more predictive model. If you need to transform, you need to try to transform every single data point, which could better be done by exp().

Hope this helps.