[R] multinom and contrasts
array chip
arrayprofile at yahoo.com
Thu Apr 14 22:35:27 CEST 2005
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".
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)
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!
--- 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
>
>
__________________________________
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: sample.txt
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20050414/bc7de7fc/sample.txt
More information about the R-help
mailing list