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@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.
