[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