5
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ For 1), you're missing that the intercept is unregularized, so $\lambda$ has no effect on the calculated intercept. $\endgroup$ – Matthew Drury Jul 5 '15 at 16:13
  • $\begingroup$ thank you! do you have a reference? or is this somehow incredibly obvious? $\endgroup$ – jmb Jul 5 '15 at 16:30
  • 1
    $\begingroup$ It's not obvious, I learned it from this paper: core.ac.uk/download/pdf/6287975.pdf $\endgroup$ – Matthew Drury Jul 5 '15 at 20:26
  • $\begingroup$ well, that would have been the obvious reference ;) thanks for your help! $\endgroup$ – jmb Jul 6 '15 at 7:33
  • $\begingroup$ Jmb do you have a code to compute the BIC for each model using glmnet procedure ? It would be of great help ! $\endgroup$ – user130069 Sep 5 '16 at 12:10
1
$\begingroup$

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))
 }
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.