[R-sig-ME] ordinal regression with MCMCglmm

[Malcolm Fairbrother]

Dear list,

I'm picking up on a thread from 2010 (https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003678.html), about ordinal mixed models, fitted with clmm (from the ordinal package) and MCMCglmm.

I would like to fit an ordinal mixed model, with random slopes. Since clmm doesn't do random slopes, I'm trying with MCMCglmm, but I'm not sure I understand how MCMCglmm handles ordinal data.

To illustrate, I'll use the "wine" data from the ordinal package, and a model with only random intercepts.

library(ordinal)
library(MCMCglmm)

data(wine)
str(wine)
summary(fm2 <- clmm2(rating ~ temp + contact, random=judge, data=wine, Hess=TRUE, nAGQ=10))

# to get predicted probabilities for each possible of the five possible responses, for two combinations of covariates:
diff(plogis(c(-Inf,fm2$Theta, Inf) - fm2$beta %*% c(0,0))) # for tempcold and contactno
diff(plogis(c(-Inf,fm2$Theta, Inf) - fm2$beta %*% c(1,1))) # for tempwarm and contactyes

# compare with the raw data (seems to make sense...):
table(wine$rating[wine$temp=="cold" & wine$contact=="no"])
table(wine$rating[wine$temp=="warm" & wine$contact=="yes"])

# now fit the same model, using MCMCglmm:
prior1 <- list(R = list(V = 1, nu = 0.002, fix=1), G = list(G1 = list(V = 1, nu = 0.002)))
MC1 <- MCMCglmm(rating ~ temp + contact, random=~judge, data=wine, family="ordinal", prior=prior1, nitt=130000, thin=100, verbose=F)

# and get predicted probabilities, using the results from MCMCglmm, and the "hunch" approach Jarrod mentioned in the thread above:
diff(pnorm(c(-Inf, 0, colMeans(MC1$CP), Inf)-colMeans(MC1$Sol) %*% c(1,0,0), 0, sqrt(1+sum(colMeans(MC1$VCV))))) # for tempcold and contactno
diff(pnorm(c(-Inf, 0, colMeans(MC1$CP), Inf)-colMeans(MC1$Sol) %*% c(1,1,1), 0, sqrt(1+sum(colMeans(MC1$VCV))))) # for tempwarm and contactyes

The predicted probabilities are similar, but not similar enough that I'm sure this is right. Are the differences due to MCMC error, inappropriate priors, or the way I'm predicting probabilities based on one (or both?) model(s)?

If anyone can offer any clarifications/suggestions, I'd be grateful. (Maybe Thierry figured this out?)

Cheers,
Malcolm



> Date: Wed, 14 Apr 2010 15:41:45 +0100
> From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
> To: "ONKELINX, Thierry" <Thierry.ONKELINX at inbo.be>
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] ordinal regression with MCMCglmm
> Message-ID: <2FE797FF-C1EC-4CA2-8275-C3B0AA9243BE at ed.ac.uk>
> Content-Type: text/plain; charset=US-ASCII; format=flowed; delsp=yes
> 
> 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