Re: [R] does multinomial logistic model from multinom (nnet) has logLik?

From: ronggui <ronggui.huang_at_gmail.com>
Date: Thu 23 Feb 2006 - 01:50:51 EST

Here is a function for calculating the measures of fit for multinomial logistic model (using nnet::multinom).If anything wrong ,I hope experts point it out.Thank you.

fitstat <- function(object) {
#thanks Ripley, B. D. for telling how to get the LogLik and when is invalid. {if (!is.null(object$call$summ) && !identical(object$call$summ,0))

   stop("when 'summ' argument is not zero,can NOT get Loglik") } object.base <- update(object,.~1,trace=FALSE) dev.base <- deviance(object.base) ; L.base <- - dev.base/2 dev.full <- deviance(object) ; L.full <- - dev.full/2 G2 <- dev.base - dev.full
df <- object$edf - object.base$edf
LR.test.p <- pchisq(G2,df,lower=F)

aic <- object$AIC

n<-dim(object$residuals)[1]

#get the predict value to cal count R2
pre <- predict(object,type="class")
y <- eval.parent(object$call$data)[,as.character(object$call$formula[[2]])] if (!identical(length(y),length(pre))) stop("Length not matched.") tab <- table(y,pre)
if (!identical(dim(tab)[1],dim(tab)[2])) stop("pred and y have diff nlevels") ad <- max(rowSums(tab))#max of row sum

#cal R2
ML.R2 <- 1-exp(-G2/n)
McFadden.R2 <- 1-(L.full/L.base)
McFadden.Adj.R2 <- 1-((L.full-mod$edf)/L.base)

Cragg.Uhler.R2 <- ML.R2/(1-exp(2*L.base/n))
Count.R2 <- sum(diag(tab))/sum(tab)
Count.adj.R2 <- (sum(diag(tab))-ad)/(sum(tab)-ad)

#get the result
res<-list(LR=G2,df=df,LR.test.p =LR.test.p ,aic=aic,ML.R2=ML.R2,Cragg.Uhler.R2=Cragg.Uhler.R2,McFadden.R2 =McFadden.R2 ,McFadden.Adj.R2=McFadden.Adj.R2,Count.R2=Count.R2,Count.adj.R2=Count.adj.R2)

#print the result
cat("\n",

    paste(rep("-",21)),
    "\n The Fitstats are : \n",
    sprintf("G2(%d) = %f",df,G2),
    " ,Prob ",format.pval(LR.test.p),
    "\n",sprintf("AIC = %f",aic),
    sprintf(",ML.R2 = %f \n",ML.R2),
    paste(rep("-",21)),"\n",

    sprintf("Cragg.Uhler.R2  = %f \n",Cragg.Uhler.R2),
    sprintf("McFadden.R2     = %f \n",McFadden.R2),
    sprintf("McFadden.Adj.R2 = %f \n",McFadden.Adj.R2),
    sprintf("Count.R2        = %f \n",Count.R2),
    sprintf("Count.adj.R2    = %f \n",Count.adj.R2),
    "\n Note:The maxinum of ML R2 is less than 1 \n",     paste(rep("-",21)),"\n")
invisible(res)
}

#example
require(nnet)
data(mexico,package="Zelig")
mod <- multinom(vote88 ~ pristr + othcok + othsocok,mexico) summary(mod,cor=F)
fitstat(mod)

#reference:
#J. SCOTT LONG and JEREMY FREESE,REGRESSION MODELS FOR CATEGORICAL DEPENDENT VARIABLES USING STATA.
> fitstat(mod)

 Note:The maxinum of ML R2 is less than 1

06-2-22ronggui<ronggui.huang@gmail.com> д
> So it's valid to get logLik (deviance/-2) when the summ argument is unused?

>

> Thank you.
>

> 2006/2/22, Prof Brian Ripley <ripley@stats.ox.ac.uk>:
> > On Wed, 22 Feb 2006, ronggui wrote:
> >
> > > I want to get the logLik to calculate McFadden.R2 ,ML.R2 and
> > > Cragg.Uhler.R2, but the value from multinom does not have logLik.So my
> > > quetion is : is logLik meaningful to multinomial logistic model from
> > > multinom?If it does, how can I get it?
> >
> > From the help page:
> >
> > Value:
> >
> > A 'nnet' object with additional components:
> >
> > deviance: the residual deviance.
> >
> > So it has a residual deviance. That is -2 log Lik in many cases (but not
> > if the argument 'summ' is used)
> >
> > > Thank you!
> > >
> > > ps: I konw VGAM has function to get the multinomial logistic model
> > > with logLik, but I prefer use the function from "official" R
> > > packages .
> > >
> > > --
> > > ronggui
> > > Deparment of Sociology
> > > Fudan University
> >
> > --
> > Brian D. Ripley, ripley@stats.ox.ac.uk
> > Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> > University of Oxford, Tel: +44 1865 272861 (self)
> > 1 South Parks Road, +44 1865 272866 (PA)
> > Oxford OX1 3TG, UK Fax: +44 1865 272595
> >
>
>

> --
> ٹ
> Deparment of Sociology
> Fudan University
>
--
ronggui
Deparment of Sociology
Fudan University



______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Received on Thu Feb 23 03:20:05 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:39 EST