Re: [R] problem with polr ?

From: John Fox <jfox_at_mcmaster.ca>
Date: Sat 11 Jun 2005 - 07:40:01 EST


Dear Marc,

> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Marc Girondot
> Sent: Friday, June 10, 2005 3:44 PM
> To: r-help@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@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@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



R-help@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 Received on Sat Jun 11 07:51:57 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:31 EST