[R] Factors and Multinomial Logistic Regression

[peter dalgaard]

On May 2, 2013, at 20:33 , Lorenzo Isella wrote:

> On Wed, 01 May 2013 23:49:07 +0200, peter dalgaard <pdalgd at gmail.com> wrote:
> 
>> It still doesn't work!!!!!
>> 
> 
> 
> Apologies; since I had already imported nnet in my workspace, the script worked on my machine even without importing it explicitly (see the script at the end of the email).
> Sorry for the confusion.

You still owe us an answer why you thought that this:

Coefficients:
     (Intercept)     science       socst femalefemale
low     1.912288 -0.02356494 -0.03892428   0.81659717
high   -4.057284  0.02292179  0.04300323  -0.03287211

Std. Errors:
     (Intercept)    science      socst femalefemale
low     1.127255 0.02097468 0.01951649    0.3909804
high    1.222937 0.02087182 0.01988933    0.3500151

Residual Deviance: 388.0697 

is at all different from the Stata output. As far as I can tell it is EXACTLY the same!

Apologies for being insistent, but this will come up in Internet searches as "I couldn't make R do what Stata does".

> 
> I now mainly have a question about a definition: I can easily calculate the relative risk ratio (RRR) and its confidence interval (CI) for a given variable of my multinomial regression by exponentiating the variable and its original CI.
> However, how is the standard error on the RRR defined? This is now the only part of the stata calculation which I cannot reproduce.
> Cheers
> 

They would appear just to be delta-method based. 

s.e.(f(thetahat)) =~ f'(thetahat) s.e.(thetahat)

in casu f() is exp() and, e.g., looking at the coef. for female in the "low" table:

> .3909813 * exp(.8166202)
[1] 0.8847277

(It is a pretty useless quantity. Stata itself doesn't use it for much, either.)

> cc <- summary(mymodel)
> exp(cc$coefficients) * cc$standard.errors
     (Intercept)    science      socst femalefemale
low   7.62989469 0.02048619 0.01877141    0.8847053
high  0.02115184 0.02135577 0.02076329    0.3386964



> Lorenzo
> 
> ##############################################################################################
> 
> 
> 
> library(foreign)
> library(nnet)
> ## 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