Factors and Multinomial Logistic Regression
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