Dear All, I am trying to reproduce the example that I found online here bit.ly/11VG4ha However, when I run my script (pasted at the end of the email), I notice that there is a factor 2 between the values for the coefficients for the categorical variable female calculated by my script and in the online example. Any idea about where this difference comes from? Besides, how can I extract the p-values for the calculated coefficients of my model and for the relative risk ratios? Best Regards Lorenzo ########################################################################### library(ares) library(foreign) ## See the Stata example at bit.ly/11VG4ha mydata <- read.dta("ats.ucla.edu/stat/data/hsb2.dta") ## IMPORTANT: redefine the base line!!! mydata$ses2 <- relevel(mydata$ses, ref = "middle") mymodel <- multinom(ses2 ~ science+ socst+ female, data=mydata) print(summary(mymodel)) print("The relative risk ratio (RRR) is, ") print(exp(coef(mymodel)))
On May 1, 2013, at 6:48 AM, Lorenzo Isella wrote:> Dear All, > I am trying to reproduce the example that I found online here > > > bit.ly/11VG4ha > > However, when I run my script (pasted at the end of the email), I notice that there is a factor 2 between the values for the coefficients for the categorical variable female calculated by my script and in the online example.You cannot interpret parameters in a regression model unless you understand the coding that is being used. In R the coding is 0/1 for factors being reference to the control level. In some other coding systems the coding is -1/1 and the coefficients for the effect of that difference are twice what they would be for a 0/1 comparison.> Any idea about where this difference comes from? > Besides, how can I extract the p-values for the calculated coefficients of my model and for the relative risk ratios?Generally the summary function will provide those. Read the help page (or the code) if you want to learn how the results are delivered and or return()-ed in the result.> Best Regards > > Lorenzo > > ########################################################################### > library(ares) > library(foreign) > > ## See the Stata example at bit.ly/11VG4ha > > mydata <- read.dta("ats.ucla.edu/stat/data/hsb2.dta") > > > ## IMPORTANT: redefine the base line!!! > > mydata$ses2 <- relevel(mydata$ses, ref = "middle") > > > > mymodel <- multinom(ses2 ~ science+ socst+ female, data=mydata) > > print(summary(mymodel)) > > print("The relative risk ratio (RRR) is, ") > > print(exp(coef(mymodel))) > > ______________________________________________ > R-help at r-project.org mailing list > stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.David Winsemius Alameda, CA, USA
On May 1, 2013, at 15:48 , Lorenzo Isella wrote:> Dear All, > I am trying to reproduce the example that I found online here > > > bit.ly/11VG4ha > > > However, when I run my script (pasted at the end of the email), I notice that there is a factor 2 between the values for the coefficients for the categorical variable female calculated by my script and in the online example. > Any idea about where this difference comes from? > Besides, how can I extract the p-values for the calculated coefficients of my model and for the relative risk ratios? > Best Regards >(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? (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> Lorenzo > > ########################################################################### > library(ares) > library(foreign) > > ## See the Stata example at bit.ly/11VG4ha > > mydata <- read.dta("ats.ucla.edu/stat/data/hsb2.dta") > > > ## IMPORTANT: redefine the base line!!! > > mydata$ses2 <- relevel(mydata$ses, ref = "middle") > > > > mymodel <- multinom(ses2 ~ science+ socst+ female, data=mydata) > > print(summary(mymodel)) > > print("The relative risk ratio (RRR) is, ") > > print(exp(coef(mymodel))) > > ______________________________________________ > R-help at r-project.org mailing list > stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.-- 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