[R] problem with polr ?

[John Fox]

Dear Marc,

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Marc Girondot
> Sent: Friday, June 10, 2005 3:44 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] problem with polr ?
> 
> I want to fit a multinomial model with logit link.
> For example let this matrix to be analyzed:
> male   female   aborted   factor
> 10     12         1         1.2
> 14     14         4         1.3
> 15     12         3         1.4
> 
> (this is an example, not the true data which are far more complex...)
> 
> I suppose the correct function to analyze these data is polr 
> from MASS library.
> 

Actually, polr fits the proportional-odds logistic regression model (or a
similar ordered probit model), which requires an ordinal response. I don't
believe that it makes sense to think of a < f < m. For the multinomial logit
model, see the multinom function in the nnet package.

> The data have been entered in a text file like this:
> 
> output  factor  n
> m       1.2     10
> f       1.2     12
> a       1.2     1
> m       1.3     14
> f       1.3     14
> a       1.3     4
> m       1.4     15
> f       1.4     12
> a       1.4     3
> 
> However, after having performed the analysis, it appears this 
> is not correct as the 3 outputs per experiment are not linked...
> 
> library(MASS)
> dt.plr <- polr(output ~ factor, data=dt, weights=n)
> >  dt.pr1<-predict(dt.plr, , type="probs")
> >  dt.pr1
>             a         f         m
> 1 0.09987167 0.4578184 0.4423099
> 2 0.09987167 0.4578184 0.4423099
> 3 0.09987167 0.4578184 0.4423099
> 4 0.09437078 0.4477902 0.4578390
> 5 0.09437078 0.4477902 0.4578390
> 6 0.09437078 0.4477902 0.4578390
> 7 0.08914287 0.4374067 0.4734505
> 8 0.08914287 0.4374067 0.4734505
> 9 0.08914287 0.4374067 0.4734505
> 

These are the fitted probabilities of response in each of the three
categories. Note that each row sums to 1, as it should. Of course, the model
likely doesn't make sense.

> 
> ---------------------------Another linked problem
> 
> Also, I don't understand what the meaning of the residual 
> deviance that is displayed.
> Let modify the file so that the model can also be analyzed 
> using binomial model:
> 
> output  factor  n
> m       1.2     10
> f       1.2     12
> a       1.2     0
> m       1.3     14
> f       1.3     14
> a       1.3     0
> m       1.4     15
> f       1.4     12
> a       1.4     0
> 
> dt.plr
> Call:
> polr(formula = output ~ factor, data = dt, weights = n)
> 
> Coefficients:
>    factor
> 2.034848
> 
> Intercepts:
>         a|f        f|m
> -16.306511   2.632410
> 
> Residual Deviance: 106.2304
> AIC: 112.2304
> 
> whereas the corresponding scaled deviance for the binomial 
> model (removing a column) is 1.842e-3...
> 
> 

I'm surprised that you were able to get estimates here at all, since there
are no observations at category a; nevetheless, polr is estimating two
thresholds, one between the never-observed a and f. I expect that this is a
numerical artifact. Note that if you simply remove the rows for a rather
than set the counts to 0, polr will complain that there are only two
categories.

I hope this helps,
 John

> -- 
> 
> __________________________________________________________
> Marc Girondot, Pr
> Laboratoire Ecologie, Systématique et Evolution
> Equipe de Conservation des Populations et des Communautés
> CNRS, ENGREF et Université Paris-Sud 11 , UMR 8079
> Bâtiment 362
> 91405 Orsay Cedex, France
> 
> Tel:  33 1 (0)1.69.15.72.30   Fax: 33 1 (0)1 69 
> 15 56 96   e-mail: marc.girondot at ese.u-psud.fr
> Web: http://www.ese.u-psud.fr/epc/conservation/Marc.html
> Skype: girondot
> Fax in US: 1-425-732-6934
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html