[R] AIC and logLik for logistic regression in R and S-PLUS

[Ben Bolker]

Leandra Desousa <sousa <at> ims.uaf.edu> writes:

> I am using 'R' version 2.2.1 and 'S-PLUS' version 6.0; and I loaded the 
> MASS library in 'S-PLUS'.
> 
> I am running a logistic regression using glm:
> 
>  >summary(mydata.glm)
> Call:
> glm(formula = COMU ~ MeanPycUpT + MeanPycUpS, family = binomial,data = 
> mydata)

  [snip]

>     Null deviance: 30.316  on 21  degrees of freedom
> Residual deviance: 23.900  on 19  degrees of freedom
> AIC: 29.9
> --------------------------------------------------------------------
> 'R'
> -----------------
>  > AIC(mydata.glm)
> [1] 29.89986
> 
>  > logLik(mydata.glm)
> 'log Lik.' -11.94993 (df=3)
> -----------------
> 
> 'S-PLUS'
> -----------------
>  > AIC(mydata.glm)
> [1] 71.03222
> 
>  > logLik(mydata.glm)
> [1] -31.51611
> -----------------
> 
>
> 1) Which AIC value is the correct one?
> 2) Which log-likelihood value is the correct one?

   AIC and log-likelihood are often defined differently
in software packages -- specifically, additive constants
can be included or excluded as long as they are done consistently,
without affecting inferences from the model.  The absolute
values of AIC and logLik aren't that important; the only thing
that really matters are differences among models.  Have you
tried comparing models within R and within S-PLUS to establish
whether they give the same inferences (I would guess they do)?

> 3) If  'extractAIC' in 'S-PLUS' and all values in 'R' are the correct 
> ones, and the 'AIC' and 'logLik' in 'S-PLUS' values are wrong then:
>    Why 'S-PLUS' cannot retrieve a log-likelihood value from my glm 
> object('mydata.glm'), even though it is using log-likelihood to 
> calculate its residual deviance?

   It's very hard for us to debug S-PLUS!  Some (most?) of
us on the list don't even have S-PLUS installed any more ...
Perhaps you should ask this question on an S-PLUS mailing list,
or of the (paid) S-PLUS technical support team ...
 
   And the obligatory R-list nags:

 * it would help if you upgraded your R version -- R 2.3.0
came out in April 2006, and we're now up to 2.5.1
  * it's the MASS "package", not the MASS "library"

  cheers
    Ben Bolker