[R] Multinomial model and p-values

[Duncan Mackay]

Hi Luciano

There are a number of types of ordinal regression and you need to be 
specific about that.
There are a number of ordinal packages to do ordinal regression.

MASS::polr
arm::bayespolr
ordinal
VGAM
repolr
geepack
etc

Each of them has specific requirements about coding of the variables 
and these MUST be adhered to.

polr may be better suited to your dataset. -- No data: cannot say.

Regards

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au

At 05:21 3/07/2013, you wrote:
>Hello everyone,
>
>I have a dataset which consists of "Pathology scores" (Absent, Mild, Severe)
>as outcome variable, and two main effects: Age (two factors: twenty / thirty
>days old) and Treatment Group (four factors: infected without ATB; infected
>+ ATB1; infected + ATB2; infected + ATB3).
>
>First I tried to fit an ordinal regression model, which seems more
>appropriate given the characteristics of my dependent variable (ordinal).
>However, the assumption of odds proportionality was severely violated
>(graphically), which prompted me to use a multinomial model instead, using
>the "nnet" package.
>
>
>First I chose the outcome level that I need to use as baseline category:
>
>Data$Path <- relevel(Data$Path, ref = "Absent")
>
>Then, I needed to set baseline categories for the independent variables:
>
>Data$Age <- relevel(Data$Age, ref = "Twenty")
>Data$Treat <- relevel(Data$Treat, ref = "infected without ATB")
>
>The model:
>
>test <- multinom(Path ~ Treat + Age, data = Data)
># weights:  18 (10 variable)
>initial  value 128.537638
>iter  10 value 80.623608
>final  value 80.619911
>converged
>
> > summary1 <- summary(test1)
> > summary1
>
>Call:
>multinom(formula = Jej_fact ~ Treat + Age, data = Data)
>
>Coefficients:
>          (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3
>AgeThirty
>Moderate   -2.238106   -1.1738540      -1.709608       -1.599301
>2.684677
>Severe     -1.544361   -0.8696531      -2.991307       -1.506709
>1.810771
>
>Std. Errors:
>          (Intercept)    infected+ATB1   infected+ATB2   infected+ATB3
>AgeThirty
>Moderate   0.7880046    0.8430368       0.7731359       0.7718480
>0.8150993
>Severe     0.6110903    0.7574311       1.1486203       0.7504781
>0.6607360
>
>Residual Deviance: 161.2398
>AIC: 181.2398
>
>For a while, I could not find a way to get the p-values for the model and
>estimates when using nnet:multinom. Yesterday I came across a post where the
>author put forward a similar issue regarding estimation of p-values for
>coefficients
>(http://stats.stackexchange.com/questions/9715/how-to-set-up-and-estimate-a-
>multinomial-logit-model-in-r).
>
>There, one blogger suggested that getting p-values from the summary() result
>of multinom() is pretty easy, by first getting the t values as follows:
>
>pt(abs(summary1$coefficients / summary1$standard.errors), df=nrow(Data)-10,
>lower=FALSE)
>
>          (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3
>AgeThirty
>Moderate 0.002670340   0.08325396      0.014506395     0.02025858
>0.0006587898
>Severe   0.006433581   0.12665278      0.005216581     0.02352202
>0.0035612114
>
>I AM NOT a statistician, so don't be baffled by a silly question! I this
>procedure correct?
>
>Cheers for now,
>
>Luciano
>
>______________________________________________
>R-help at 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.