[R-sig-ME] help with ordinal mixed model (reposted from r-help)

[Ben Bolker]

[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?