[R-sig-ME] ordinal regression with MCMCglmm

Jarrod Hadfield j.hadfield at ed.ac.uk
Wed Apr 14 16:41:45 CEST 2010


Dear Thierry,

I think (but you had better check) that it would actually be something  
like

pnorm(-Xb, 0, sqrt(1+v)),
pnorm(cp[1] - Xb,0, sqrt(1+v)) - pnorm(Xb,0, sqrt(1+v))
pnorm(cp[2] - Xb,0, sqrt(1+v)) - pnorm(cp[1] - Xb,0, sqrt(1+v))
1- pnorm(cp[2] - Xb,0, sqrt(1+v))

where v is  the sum of the variance components for the residual and  
random effects.

Alternatively, if the residual variance is set to one and there are no  
random effects

pnorm(-Xb, 0, sqrt(2)),
pnorm(cp[1] - Xb,0, sqrt(2)) - pnorm(Xb,0, sqrt(2))
pnorm(cp[2] - Xb,0, sqrt(2)) - pnorm(cp[1] - Xb,0, sqrt(2))
1- pnorm(cp[2] - Xb,0, sqrt(2))

or if the prediction includes random effects:

pnorm(-(Xb+Zu), 0, sqrt(2)),
pnorm(cp[1] - (Xb+Zu),0, sqrt(2)) - pnorm((Xb+Zu),0, sqrt(2))
pnorm(cp[2] - (Xb+Zu),0, sqrt(2)) - pnorm(cp[1] - (Xb+Zu),0, sqrt(2))
1- pnorm(cp[2] - (Xb+Zu),0, sqrt(2))

Please, do not take this as gospel. I have not got time to check these  
results, they are a hunch.

Cheers,

Jarrod


On 14 Apr 2010, at 15:25, ONKELINX, Thierry wrote:

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


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.




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