From: ronggui <ronggui.huang_at_gmail.com>

Date: Fri 24 Feb 2006 - 21:29:17 EST

Date: Fri 24 Feb 2006 - 21:29:17 EST

I should ues "Not obvious " instead of "invalid", this is my mistake.

2006/2/24, Prof Brian Ripley <ripley@stats.ox.ac.uk>:

> Please note, I told you that the deviance was minus twice log-likelihood

*> unless summ > 0. I had not checked the latter case, where it is not
**> obvious, but I did not say it was invalid.
**>
**> In fact the answer is to be found on p.203 of MASS4 (we do ask people to
**> read the supporting documentation), and this is valid also for summ > 0.
*

I did read the help page in nnet.And MASS4 is not avaivable just now :(

I will check that when I can.

> I will add a comment to the help file.

appreciate!

> On Wed, 22 Feb 2006, ronggui wrote:

*>
**> > 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 not obvious.
**> > {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)
**> >
**> > - - - - - - - - - - - - - - - - - - - - -
**> > The Fitstats are :
**> > G2(6) = 381.351620 ,Prob < 2.22e-16
**> > AIC = 2376.571142 ,ML.R2 = 0.244679
**> > - - - - - - - - - - - - - - - - - - - - -
**> > Cragg.Uhler.R2 = 0.282204
**> > McFadden.R2 = 0.139082
**> > McFadden.Adj.R2 = 0.133247
**> > Count.R2 = 0.596026
**> > Count.adj.R2 = 0.123003
**> >
**> > Note:The maxinum of ML R2 is less than 1
**> > - - - - - - - - - - - - - - - - - - - - -
**> >
**> > ÔÚ 06-2-22£¬ronggui<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
**> >
**> >
**>
**> --
**> 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 UniversityReceived on Fri Feb 24 21:37:22 2006______________________________________________ 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

*
This archive was generated by hypermail 2.1.8
: Sat 25 Feb 2006 - 00:08:47 EST
*