# [R] Factors and Multinomial Logistic Regression

peter dalgaard pdalgd at gmail.com
Thu May 2 22:04:26 CEST 2013

```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
>
>
>
> 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

```