[R] Factors and Multinomial Logistic Regression

[peter dalgaard]

On May 1, 2013, at 22:40 , Lorenzo Isella wrote:

> 
>> 
>> (A) The example doesn't run for me. library(ares) is not available on current R versions, but even where it is available, it doesn't provide a multinom() function?
> 
> 
> Apologies, ares is not needed at all. Please find the correct script at the end of the email.

It still doesn't work!!!!!


> 
> mymodel <- multinom(ses2 ~ science+ socst+ female, data=mydata)
Error: could not find function "multinom"
> 
> 
>

Please come up with a _reproduceable_ example. Also show what output you are getting. And show options("contrasts"); either something has changed them from the default or we are not talking about the same multinom() function. Halved coefficients is entirely compatible with using  "contr.sum" for unordered factors.


> 
>> 
>> (B) If I insert library(nnet), to get a multinom(), I get exactly the same result as Stata does!
>> 
>> Did you by any chance diddle with options(contrasts=...)?
>> 
>> -pd
> 
> No, I did not. The point is that if I use a variable female, which has two levels, then I do not reproduce the results of stata for that variable only.
> If instead I define a variable "sex" which assumes the numerical values 0/1, then I reproduce entirely the results by stata.
> 
> Hope this helps.
> 
> Lorenzo
> 
> 
> ##################################################################
> 
> 
> library(foreign)
> 
> ## See the Stata example at http://bit.ly/11VG4ha
> 
> mydata <- read.dta("http://www.ats.ucla.edu/stat/data/hsb2.dta")
> 
> 
> sex <- rep(0, dim(mydata)[1])
> 
> sel <- which(mydata$female=="male")
> 
> sex[sel] <- 1
> 
> mydata$sex <- sex
> 
> ## IMPORTANT: redefine the base line!!!
> 
> mydata$ses2 <- relevel(mydata$ses, ref = "middle")
> 
> 
> ## NB: for some reason, if I use female (a factor assuming two values)
> ## I do not reproduce the results of the example.
> ## I need to use a variable which is numeric and assumes two values
> ## (that is why I introduced the variable sex))
> 
> ## mymodel <- multinom(ses2 ~ science+ socst+ sex, data=mydata)
> 
> 
> mymodel <- multinom(ses2 ~ science+ socst+ female, data=mydata)
> 
> 
> 
> 
> print(summary(mymodel))
> 
> print("The relative risk ratio (RRR) is, ")
> 
> print(exp(coef(mymodel)))
> 
> 
> 

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com