[R] multinom and contrasts

John Fox jfox at mcmaster.ca
Thu Apr 14 23:26:51 CEST 2005


Dear array,

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of array chip
> Sent: Thursday, April 14, 2005 3:35 PM
> To: John Fox
> Cc: R-help at stat.math.ethz.ch
> Subject: RE: [R] multinom and contrasts
> 
> Dear John,
> 
> My dataset has a response variable with 6 levels, and
> 12 independent variables, 10 of them are continuous variable, 
> one is a categorical variable with 2 levels, the other one is 
> also a categorical variable with 4 levels. total 206 
> observations. I attached my dataset with the email "sample.txt".
> 

You are, therefore, fitting a complex model with 5*(1 + 10 + 1 + 3) = 75
parameters to a data set with 206 observations.
 
> library(MASS)
> library(nnet)
> 
> sample<-read.table("sample.txt",sep='\t',header=T,row.names=1)
> wts2<-sample$wts
> sample<-sample[,-4]
> 
> options(contrasts=c('contr.helmert','contr.poly'))
> obj1<-multinom(class~.,sample,weights=wts2,maxit=1000)
> options(contrasts=c('contr.treatment','contr.poly'))
> obj2<-multinom(class~.,sample,weights=wts2,maxit=1000)
> 
> predict(obj1,type='probs')[1:5,]
> predict(obj2,type='probs')[1:5,]
> 
> Interestingly, if I change the values of the variable "bkgd" 
> for 2 observations (from "a", to "f"), then I can get 
> convergence with helmert contrast, but still not converged 
> with treatment contrast:
> 
> sample$bkgd[201]<-'f'
> sample$bkgd[205]<-'f'
> 
> options(contrasts=c('contr.helmert','contr.poly'))
> obj1<-multinom(class~.,sample,weights=wts2,maxit=1000)

If you look at the covariance matrix of the coefficients, you'll see that
there are still problems with the fit:

	> summary(diag(vcov(obj1)))
      	Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
	-1.029e-20  6.670e+05  3.502e+06  7.781e+06  1.075e+07  4.768e+07 

That negative variance (though essentially 0) doesn't bode well. Indeed,
requiring nearly 1000 iterations for apparent convergence is in itself an
indication of problems.

> options(contrasts=c('contr.treatment','contr.poly'))
> obj2<-multinom(class~.,sample,weights=wts2,maxit=1000)
> 
> predict(obj1,type='probs')[1:5,]
> predict(obj2,type='probs')[1:5,]
> 
> appreciate any suggestions!
> 

I don't know anything about your application, so I hesitate to make
recommendations, but (assuming that the model you've fit makes sense) in the
abstract I'd suggest collecting a lot more data or fitting a much simpler
model.

Regards,
 John

> 
> 
> --- John Fox <jfox at mcmaster.ca> wrote:
> 
> > Dear chip,
> > 
> > > -----Original Message-----
> > > From: r-help-bounces at stat.math.ethz.ch 
> > > [mailto:r-help-bounces at stat.math.ethz.ch] On
> > Behalf Of array chip
> > > Sent: Thursday, April 14, 2005 1:00 PM
> > > To: John Fox
> > > Cc: R-help at stat.math.ethz.ch
> > > Subject: RE: [R] multinom and contrasts
> > > 
> > > Dear John,
> > > 
> > > Thanks for the answer! In my own dataset, The
> > > multinom() did not converge even after I had tried
> > to
> > > increase the maximum number of iteration (from
> > default 100 to
> > > 1000). In this situation, there is some bigger
> > diffrenece in
> > > fitted probabilities under different contrasts
> > (e.g. 
> > > 0.9687817 vs. 0.9920816). My question is whether
> > the analysis
> > > (fitted probabilities) is still valid if it does
> > not
> > > converge? and what else can I try about it?
> > > 
> > 
> > If multinom() doesn't converge to a stable solution after 1000 
> > iterations, it's probably safe to say that the problem is 
> > ill-conditioned in some respect. Have you looked at the covariance 
> > matrix of the estimates?
> > 
> > Regards,
> >  John
> > 
> > > Thank you!
> > > 
> > > 
> > > 
> > > 
> > > --- John Fox <jfox at mcmaster.ca> wrote:
> > > > Dear chip,
> > > > 
> > > > The difference is small and is due to
> > computational error.
> > > > 
> > > > Your example:
> > > > 
> > > > > max(abs(zz[1:10,] - yy[1:10,]))
> > > > [1] 2.207080e-05
> > > > 
> > > > Tightening the convergence tolerance in
> > multinom() eliminates the
> > > > difference:
> > > > 
> > > > >
> > > >
> > options(contrasts=c('contr.treatment','contr.poly'))
> > > > >
> > > >
> > >
> >
> xx<-multinom(Type~Infl+Cont,data=housing[-c(1,10,11,22,25,30),],
> > > > reltol=1.0e-12)
> > > > # weights:  20 (12 variable)
> > > > initial  value 91.495428
> > > > iter  10 value 91.124526
> > > > final  value 91.124523
> > > > converged
> > > > > yy<-predict(xx,type='probs')
> > > > >
> > options(contrasts=c('contr.helmert','contr.poly'))
> > > > >
> > > >
> > >
> >
> xx<-multinom(Type~Infl+Cont,data=housing[-c(1,10,11,22,25,30),],
> > > > reltol=1.0e-12)
> > > > # weights:  20 (12 variable)
> > > > initial  value 91.495428
> > > > iter  10 value 91.125287
> > > > iter  20 value 91.124523
> > > > iter  20 value 91.124523
> > > > iter  20 value 91.124523
> > > > final  value 91.124523
> > > > converged
> > > > > zz<-predict(xx,type='probs')
> > > > > max(abs(zz[1:10,] - yy[1:10,]))
> > > > [1] 1.530021e-08
> > > > 
> > > > I hope this helps,
> > > >  John
> > > > 
> > > > --------------------------------
> > > > John Fox
> > > > Department of Sociology
> > > > McMaster University
> > > > Hamilton, Ontario
> > > > Canada L8S 4M4
> > > > 905-525-9140x23604
> > > > http://socserv.mcmaster.ca/jfox
> > > > --------------------------------
> > > > 
> > > > > -----Original Message-----
> > > > > From: r-help-bounces at stat.math.ethz.ch 
> > > > > [mailto:r-help-bounces at stat.math.ethz.ch] On
> > > > Behalf Of array chip
> > > > > Sent: Wednesday, April 13, 2005 6:26 PM
> > > > > To: R-help at stat.math.ethz.ch
> > > > > Subject: [R] multinom and contrasts
> > > > > 
> > > > > Hi,
> > > > > 
> > > > > I found that using different contrasts (e.g.
> > > > > contr.helmert vs. contr.treatment) will
> > generate
> > > > different
> > > > > fitted probabilities from multinomial logistic
> > > > regression
> > > > > using multinom(); while the fitted
> > probabilities
> > > > from binary
> > > > > logistic regression seem to be the same. Why
> > is
> > > > that? and for
> > > > > multinomial logisitc regression, what contrast
> > > > should be
> > > > > used? I guess it's helmert?
> > > > > 
> > > > > here is an example script:
> > > > > 
> > > > > library(MASS)
> > > > > library(nnet)
> > > > > 
> > > > >       #### multinomial logistic
> > > > >
> > > >
> > options(contrasts=c('contr.treatment','contr.poly'))
> > > > >
> > > >
> > >
> >
> xx<-multinom(Type~Infl+Cont,data=housing[-c(1,10,11,22,25,30),])
> > > > > yy<-predict(xx,type='probs')
> > > > > yy[1:10,]
> > > > > 
> > > > >
> > options(contrasts=c('contr.helmert','contr.poly'))
> > > > >
> > > >
> > >
> >
> xx<-multinom(Type~Infl+Cont,data=housing[-c(1,10,11,22,25,30),])
> > > > > zz<-predict(xx,type='probs')
> > > > > zz[1:10,]
> > > > > 
> > > > > 
> > > > >       ##### binary logistic
> > > > >
> > > >
> > options(contrasts=c('contr.treatment','contr.poly'))
> > > > >
> > > >
> > >
> >
> obj.glm<-glm(Cont~Infl+Type,family='binomial',data=housing[-c(
> > > > 1,10,11,22,25,30),])
> > > > > yy<-predict(xx,type='response')
> > > > > 
> > > > >
> > options(contrasts=c('contr.helmert','contr.poly'))
> > > > >
> > > >
> > >
> >
> obj.glm<-glm(Cont~Infl+Type,family='binomial',data=housing[-c(
> > > > 1,10,11,22,25,30),])
> > > > > zz<-predict(xx,type='response')
> > > > > 
> > > > > Thanks
> > > > > 
> > > > > ______________________________________________
> > > > > R-help at 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
> > > > 
> > > >
> > > 
> > > ______________________________________________
> > > R-help at 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
> > 
> > 
> 
> 
> 
> 		
> __________________________________ 
> 
> 
>




More information about the R-help mailing list