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

[Jarrod Hadfield]

Hi,

You do not seem to specify a prior? With ordinal data the residual/ 
units variance is not identified - you need to fix it at something  
(e.g. 1):

prior=list(R=list(V=1, fix=1), G=list(G1=list(....)))

where I will leave ... to your judgement.

If you don't it will generate nonsense.

Jarrod

On 26 Mar 2012, at 15:53, Ben Bolker wrote:

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