[R-sig-ME] ordinal regression with MCMCglmm

ONKELINX, Thierry Thierry.ONKELINX at inbo.be
Wed Apr 14 16:25:12 CEST 2010


Dear Jarrod,

I'm working on a similar problem. Does it makes sense to calculate that
for the fixed effects only? Something like this:
pnorm(-Xb), 
pnorm(cp[1] - Xb) - pnorm(Xb)
pnorm(cp[2] - Xb) - pnorm(cp[1] - Xb)
1 - pnorm(cp[2] - Xb)

Best regards,

Thierry
------------------------------------------------------------------------
----
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey
  

> -----Oorspronkelijk bericht-----
> Van: r-sig-mixed-models-bounces at r-project.org 
> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens 
> Jarrod Hadfield
> Verzonden: woensdag 14 april 2010 15:35
> Aan: Kari Ruohonen
> CC: Kari Ruohonen; r-sig-mixed-models at r-project.org
> Onderwerp: Re: [R-sig-ME] ordinal regression with MCMCglmm
> 
> Hi,
> 
> Imagine a latent variable (l) that conforms to the  standard 
> linear model.
> 
> l = Xb+Zu+e
> 
> The probabilities of falling into each of the four categories are:
> 
> 
> pnorm(-l)
> 
> pnorm(cp[1]-l)-pnorm(-l)
> 
>   pnorm(cp[2]-l)-pnorm(cp[1]-l)
> 
> 1-pnorm(cp[2]-l)
> 
> where cp is the vector of cut-points with 2 elements. A 3 
> cut-point model would be over-parameterised (unless the
> intercept is zero, which I presume is what polr does (?).
> 
> The factors don't need to be ordered, the order is obtained 
> from levels(resp).  In the future, I may only allowed ordered 
> factors to stop any accidents.
> 
> Cheers,
> 
> Jarrod
> 
> 
> 
> 
> 
> On 14 Apr 2010, at 08:07, Kari Ruohonen wrote:
> 
> > Hi and thanks for the answer. I tried exactly that model 
> syntax before 
> > posting but the output of the "fixed" part had an unexpected 
> > parameterisation and I thought I misspecified the model 
> somehow. The 
> > parameters I got with the above model are
> > - two cutpoints
> > - intercept
> > - effect of group B
> >
> > I would have expected that instead of the intercept and two 
> cutpoints 
> > I would have had three cutpoints as given by polr (MASS 
> package), for 
> > example. Can you explain me the parameterisation in
> MCMCglmm and how 
> > it connects to the one in polr that uses J-1 ordered cutpoints 
> > (J=number of score classes) without an intercept?
> >
> > Also, I am uncertain do I need to convert the "resp" before 
> MCMCglmm 
> > to an ordered factor (with "ordered")?
> >
> > Many thanks,
> >
> > Kari
> >
> > On Tue, 2010-04-13 at 17:41 +0100, Jarrod Hadfield wrote:
> >> Hi Kari,
> >>
> >> The simplest model is
> >>
> >>
> >> m1<-MCMCglmm(resp~treat, random=~group, family="ordinal", 
> >> data=your.data, prior=prior)
> >>
> >> as with multinomial data with a single realisation, the residual 
> >> variance cannot be estimated from the data. The best 
> option is to fix 
> >> it at some value. most programs fix it at zero but 
> MCMCglmm will fail 
> >> to mix if this is done, so I usually fix it at 1:
> >>
> >>
> >> prior=list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=0)))
> >>
> >> I have left the default prior for the fixed effects (not 
> explicitly 
> >> specified above), and the default prior random effect variance 
> >> structure (G) which has a zero degree of belief parameter. 
> Often this 
> >> requires some/more thought, especially if there are few groups or 
> >> replication within groups is low. Sections 1.2, 1.5 & 8.2 in the 
> >> CourseNotes cover priors for variances.
> >>
> >>
> >> Currently there is no option for specifying priors on the 
> cut- points 
> >> - the prior is flat and improper. The posterior in virtually all 
> >> cases will be proper though.
> >>
> >> Cheers,
> >>
> >> Jarrod
> >>
> >> Quoting Kari Ruohonen <kari.ruohonen at utu.fi>:
> >>
> >>> Hi,
> >>> I am trying to figure out how to fit an ordinal regression model 
> >>> with MCMCglmm. The "MCMCglmm Course notes" has a section on 
> >>> multinomial models but no example of ordinal models. 
> Suppose I have 
> >>> the following data
> >>>
> >>>> data
> >>>   resp treat group
> >>> 1     4     A    1
> >>> 2     4     A    1
> >>> 3     3     A    2
> >>> 4     4     A    2
> >>> 5     2     A    3
> >>> 6     4     A    3
> >>> 7     2     A    4
> >>> 8     2     A    4
> >>> 9     3     A    5
> >>> 10    2     A    5
> >>> 11    1     B    6
> >>> 12    1     B    6
> >>> 13    1     B    7
> >>> 14    2     B    7
> >>> 15    2     B    8
> >>> 16    3     B    8
> >>> 17    2     B    9
> >>> 18    1     B    9
> >>> 19    2     B   10
> >>> 20    2     B   10
> >>>
> >>> and the "resp" is an ordinal response, "treat" is a treatment and 
> >>> "group" is membership to a group. Assume I would like to fit an 
> >>> ordinal model between "resp" and "treat" by having 
> "group" effects 
> >>> as random effects. How would I specify such a model in 
> MCMCglmm? And 
> >>> how would I specify the prior distributions?
> >>>
> >>> All help is greatly appreciated.
> >>>
> >>> regards, Kari
> >>>
> >>> _______________________________________________
> >>> R-sig-mixed-models at r-project.org mailing list 
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>>
> >>>
> >>
> >>
> >>
> >
> >
> >
> 
> 
> --
> The University of Edinburgh is a charitable body, registered 
> in Scotland, with registration number SC005336.
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 

Druk dit bericht a.u.b. niet onnodig af.
Please do not print this message unnecessarily.

Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer 
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in  this message 
and any annex are purely those of the writer and may not be regarded as stating 
an official position of INBO, as long as the message is not confirmed by a duly 
signed document.




More information about the R-sig-mixed-models mailing list