On May 2, 2013, at 20:33 , Lorenzo Isella wrote:
> On Wed, 01 May 2013 23:49:07 +0200, peter dalgaard <pda...@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....@cbs.dk Priv: pda...@gmail.com
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
