help with ordinal mixed model (reposted from r-help)

[forwarded to r-sig-mixed]

Ivan Allaman <ivanalaman <at> yahoo.com.br> writes:

Good afternoon, gentlemen!

(Note that we are not all gentlemen here ...)

After several days studying and researching on categorical data
(various forums with answers from the owner of the library - all
incipient) how to interpret the output the function MCMCglmm, come
to enlist the help of you, if someone has already worked with
MCMCglmm function in the case of variables ordinal dependent. I've
read and reread all the pdf's of the package, the coursenotes
Jarrod, finally, I'm exhausted. To clarify the database, the
treatment (called fases) consist of three levels (1-pre, 2-propolis
and 3-vincris) and the ordinal variable response has three
categories (1-normal, 2-agudo, 3 - cronico). See table!

 Thank you for the reproducible example.

 I'm forwarding this to r-sig-mixed-models at r-project.org,
which is really more appropriate.

 Your biggest problem is that your chain is mixing really, really
badly, probably because your data set is too small/the model is
overfitted ...  you may need to use some informative priors to get
things back on track, or you may need to try something simpler ...

 Ben Bolker

##  a few tweaks for prettier code: set factor labels right away
du <- transform(
read.table('http://dl.dropbox.com/u/33619290/Dados/Dtest.txt',
                          header=TRUE),
               FASES=factor(FASES,
                labels=c('Normal','Aguda','Cr?nica')),
               ANIMAIS=factor(ANIMAIS),
               ALT.RENAIS=ordered(ALT.RENAIS,
               labels=c('Pre','Propolis','Vincr')))
summary(du)
du <- na.omit(du)

(tabela <- with(du,table(FASES,ALT.RENAIS)))

## this shows that there really isn't very much information
## in the data set ...

library(ggplot2)
ggplot(du,aes(FASES,ALT.RENAIS,group=ANIMAIS))+
   geom_point()+geom_line()+facet_wrap(~ANIMAIS)


library('MCMCglmm')

#the mixed model:
set.seed(1)
mod1 <- MCMCglmm(ALT.RENAIS ~-1+FASES, random= ~ ANIMAIS,
    family='ordinal',pl=TRUE,data=du)
summary(mod1)

xyplot(mod1$Sol)  ## NOT mixing well

## run for longer ...
mod2 <- MCMCglmm(ALT.RENAIS ~-1+FASES, random= ~ ANIMAIS,
    family='ordinal',pl=TRUE,data=du,nitt=100000)
xyplot(mod2$Sol)
xyplot(mod2$CP)  ## !!!

Then the pain starts, since the documentation is insufficient in
this case.  According to him Jarrod (forums), the a posteriori means
of the coefficients of the covariates are the probit
scale. According to my study, these coefficients are the scores of
standard normal distribution. More scores should not correspond to
cutpoints? In this case, we would have j (response variable
categories) -1 cutpoints, ie, two cutpoints. The output shows me
only one cutpoint. How can then calculate the probabilities with
only one cutpoint? According to the documentation (Vignettes, page
22), if P (y = k) = F (yk | l (vlatente), 1) - F (yk-1 | l, 1), this
'1' would probably be the category '1' of the dependent variable?
Anyway gentlemen, how can I extract the probabilities for the stages
for each category of the dependent variable? I thank everyone's
attention.

From
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003671.html :

[for a four-category model]
If 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)

For an unknown individual in category 1 (i.e. set u=0) the prediction
would be

l = -3.605

## category 1: pnorm(-l) = pnorm(3.605) = 0.999
summary(pnorm(-mod2$Sol[,1]))

## category 2: pnorm(cp[1]-l)-pnorm(-l)
summary(pnorm(mod2$CP-mod2$Sol[,1])-pnorm(-mod2$Sol[,1]))

## category 3: pnorm(cp[1]-l) [OR pnorm(cp[1]-l,lower.tail=FALSE)]
summary(pnorm(mod2$CP,lower.tail=FALSE))

library(ordinal) ## model MUST have an intercept according to ?clm
mod3 <- clmm(ALT.RENAIS ~FASES+(1|ANIMAIS),data=du,link="probit")
summary(mod3)

This is not quite working either -- maybe the model is just plain
overfitted?