[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