[R] multinomial logistic regression with equality constraints?
Jasjeet Singh Sekhon
sekhon at berkeley.edu
Sun Feb 4 05:42:13 CET 2007
Hi Roger,
Walter's command is correct. To match the exact normalization used by
nnet's multinom(), however, you would need to make the coefficients
zero for the first class (i.e., y1) and not the last (i.e., y3).
mr <- multinomRob(list(y2 ~ x1 + x2, y3 ~ x1 + x2, y1~0),data=d,
print.level=1)
The results are:
MNL Estimates:
y1 y2 y3
NA/(Intercept)/(Intercept) 0 -0.6931462 -0.6931462
NA/x1/x1 0 1.3862936 0.6931474
NA/x2/x2 0 0.6931474 1.3862936
Compare to the output from nnet's multinom:
> summary(m1)
Call:
multinom(formula = y ~ x1 + x2, data = d)
Coefficients:
(Intercept) x1 x2
b -0.6931475 1.3862975 0.6931499
c -0.6931475 0.6931499 1.3862975
Also, the MLE MNL objects are in:
mr$mnl
To constrain the x1 coeffs to be equal do:
emr <- multinomRob(list(y2 ~ x1 + x2,y3 ~ x1 + x2, y1~0),data=d, print.level=1,
equality=list(list(y2~x1+0,y3~x1+0)))
To constrain y2's x1 to be equal to y3's x2:
emr2 <- multinomRob(list(y2 ~ x1 + x2,y3 ~ x1 + x2, y1~0),data=d, print.level=1,
equality=list(list(y2~x1+0,y3~x2+0)))
See the multinomRob help file for more details:
http://sekhon.berkeley.edu/robust/multinomRob.html
Any BFGS warnings can be ignored because you are not interested in the
robust estimates (they are comming from LQD estimation and require
changing 'genoud.parms').
Cheers,
Jas.
=======================================
Jasjeet S. Sekhon
Associate Professor
Travers Department of Political Science
Survey Research Center
UC Berkeley
http://sekhon.berkeley.edu/
V: 510-642-9974 F: 617-507-5524
=======================================
Walter Mebane writes:
> Roger,
>
> summary(multinomRob(list(y1 ~ x1 + x2,y2 ~ x1 + x2, y3 ~ 0),data=d,
> print.level=1))
>
> Walter Mebane
>
> Roger Levy writes:
> > Many thanks for pointing this out to me!
> >
> > I'm still a bit confused, however, as to how to use multinomRob. For
> > example I tried to translate the following example using nnet:
> >
> >
> > x1 <- c(1,1,1,1,0,0,0,0,0,0,0,0)
> > x2 <- c(0,0,0,0,1,1,1,1,0,0,0,0)
> > y <- factor(c("a","b","b","c","a","b","c","c","a","a","b","c"))
> > library(nnet)
> > d <- data.frame(x1,x2,y)
> > summary(multinom(y ~ x1 + x2, data=d))
> >
> >
> > into multinomRob as follows:
> >
> >
> > x1 <- c(1,1,1,1,0,0,0,0,0,0,0,0)
> > x2 <- c(0,0,0,0,1,1,1,1,0,0,0,0)
> > y <- factor(c("a","b","b","c","a","b","c","c","a","a","b","c"))
> > y1 <- ifelse(y=="a",1, 0)
> > y2 <- ifelse(y=="b", 1, 0)
> > y3 <- ifelse(y=="c", 1, 0)
> > d <- data.frame(x1,x2,y,y1,y2,y3)
> > summary(multinomRob(list(y1 ~ x1 + x2,y2 ~ x1 + x2, y3 ~ x1 + x2),data=d))
> >
> > but the last command gives me the error message:
> >
> >
> > [1] "multinomMLE: Hessian is not positive definite"
> > Error in obsformation %*% opg : non-conformable arguments
> >
> >
> > though it's not obvious to me why. I also tried a couple other variants:
> >
> >
> > > summary(multinomRob(list(y1 ~ 0,y2 ~ x1 + x2,y3 ~ x1 + x2),data=d))
> > Error in multinomT(Yp = Yp, Xarray = X, xvec = xvec, jacstack =
> > jacstack, :
> > (multinomT): invalid specification of Xarray (regressors not
> > allowed for last category
> > > summary(multinomRob(list(y1 ~ 0,y2 ~ x1 ,y3 ~ x2),data=d))
> > Error in multinomT(Yp = Yp, Xarray = X, xvec = xvec, jacstack =
> > jacstack, :
> > (multinomT): invalid specification of Xarray (regressors not
> > allowed for last category
> >
> >
> > Any advice would be much appreciated!
> >
> >
> > Many thanks,
> >
> > Roger
> >
> > Walter Mebane wrote:
> > > By default, with print.level=0 or greater, the multinomRob program
> > > prints the maximum likelihood estimates with conventional standard
> > > errors before going on to compute the robust estimates.
> > >
> > > Walter Mebane
> > >
> > > Jasjeet Singh Sekhon writes:
> > > >
> > > > Hi Roger,
> > > >
> > > > Yes, multinomRob can handle equality constraints of this type---see
> > > > the 'equality' option. But the function assumes that the outcomes are
> > > > multinomial counts and it estimates overdispersed multinomial logistic
> > > > models via MLE, a robust redescending-M estimator, and LQD which is
> > > > another high breakdown point estimator. It would be a simple matter
> > > > to edit the 'multinomMLE' function to work without counts and to do
> > > > straight MNL instead, but right now it estimates an overdispersed MNL
> > > > model.
> > > >
> > > > Cheers,
> > > > Jas.
> > > >
> > > > =======================================
> > > > Jasjeet S. Sekhon
> > > >
> > > > Associate Professor
> > > > Travers Department of Political Science
> > > > Survey Research Center
> > > > UC Berkeley
> > > >
> > > > http://sekhon.berkeley.edu/
> > > > V: 510-642-9974 F: 617-507-5524
> > > > =======================================
> > > >
> > > >
> > > >
> > > > Roger Levy writes:
> > > > > I'm interested in doing multinomial logistic regression with equality
> > > > > constraints on some of the parameter values. For example, with
> > > > > categorical outcomes Y_1 (baseline), Y_2, and Y_3, and covariates X_1
> > > > > and X_2, I might want to impose the equality constraint that
> > > > >
> > > > > \beta_{2,1} = \beta_{3,2}
> > > > >
> > > > > that is, that the effect of X_1 on the logit of Y_2 is the same as the
> > > > > effect of X_2 on the logit of Y_3.
> > > > >
> > > > > Is there an existing facility or package in R for doing this? Would
> > > > > multinomRob fit the bill?
> > > > >
> > > > > Many thanks,
> > > > >
> > > > > Roger
> > > > >
> > > > >
> > > > > --
> > > > >
> > > > > Roger Levy Email: rlevy at ucsd.edu
> > > > > Assistant Professor Phone: 858-534-7219
> > > > > Department of Linguistics Fax: 858-534-4789
> > > > > UC San Diego Web: http://ling.ucsd.edu/~rlevy
> > > > >
> > > > >
> > >
>
> --
> * - * - * - * - * - * - * - * - * - * - * - * - * - * - * - * - * - *
> Walter R. Mebane, Jr. email: wrm1 at cornell.edu
> Professor office voice: 607/255-3868
> Department of Government cell: 607/592-0546
> Cornell University fax: 607/255-4530
> 217 White Hall WWW: http://macht.arts.cornell.edu/wrm1/
> Ithaca, NY 14853-7901
> * - * - * - * - * - * - * - * - * - * - * - * - * - * - * - * - * - *
More information about the R-help
mailing list