Solved – Calculating the Log Likelihood of models in glmnet

devianceglmnetlikelihood

glmnet() returns a lambda sequence fitobj$lambda and I would like to calculate the log likelihood of the models (LL_model) defined by the lambda sequence.

The obvious solution is to just take the parameters and calculate the LL of each model manually. However, this not very elegant and very slow. Therefore, I am trying to calculate the LL_model from the deviance-measures that are returned by glmnet.

the glmnet-object gives me:

1) nulldev = 2*(LL_sat – LL_null), where LL_sat is the saturated model and LL_null is the NULL model (one value)

2) dev.ratio = 1 – dev.model/nulldev, where dev.model is the deviance of the model at hand (k values, for k lambda values/models)

3) and glmnet:::deviance.glmnet gives me the dev_model = (1-dev.ratio)nulldev dev.model = 2(LL_sat – LL_model) (k values for, k lambda values/models)

To calculate the LL for each model, I would do the following:

1) Calculate LL_null

2) Solve (1) for LL_sat and calculate LL_sat (one value)

3) solve (3) for LL_model and calculate LL_model (k vector)

Now, my two questions are:

1) How is this NULL model defined? The glmnet-manual says "The NULL model refers to the intercept model." But I am a bit puzzled by the fact that there is only one NULL model and nulldev for the whole lambda sequence. Based on which lambda is this intercept model calculated? I have the feeling I am missing something.

2) Does anybody see an easier way to calculate the LL of each model in the sequence? The final goal is to calculate the (E)BIC of each model. I am surprised by the fact that turns out to be so tedious.

Any help would be greatly appreciated!

Best Answer

Here is one sure way to compute the log-likelihood of logistic (and probit) regressions, no matter how it is estimated. All you need is the dependent dummy and fitted values.

I tested the function on a model with just a constant and it seems to work well. It is NOT that hard to compute a log likelihood -- so just write it down and compute it by hand. Then, just always use the same function everywhere and your computations will be coherent.

I really don't know why logLik is not yet compatible with glmnet. Penalized pseudo-R2 are often used measures of fit, you need it for information criterion and it just seems more convenient to have that value to compute plenty of test statistics. Anyway, here it is for future reference:

 logit_logLik <- function(dummy,fitted_values){
   # Description: Computes log-likelihood of a fitted
   # logit model 

   # Format variables
   y <- as.matrix(dummy)
   p <- as.matrix(fitted_values)
   # Adjust dimensions
   skip <- dim(y)[1] - dim(p)[1]
   y    <- as.matrix(y[-c(1:skip),])

   # Compute log-likelihood
   item <- sapply(1:dim(y)[1], function(i)
                          y[i,]*log(p[i,]) + (1-y[i,])*log(1-p[i,]))
   return(sum(item))
 }
Related Question