Good afternoon, gentlemen! 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!
du <-
transform(read.table('http://dl.dropbox.com/u/33619290/Dados/Dtest.txt',h=T),FASES=factor(FASES),ALT.RENAIS=ordered(ALT.RENAIS))
summary(du)
library('MCMCglmm')
du <- subset(du, ALT.RENAIS != 'NA')
tabela <- table(du[,c(2,4)])
tabela
colnames(tabela) <- c('Normal','Aguda','Cr?nica')
rownames(tabela) <- c('Pre','Propolis','Vincr')
tabela
#the mixed model:
set.seed(1)
mod1 <- MCMCglmm(ALT.RENAIS ~-1+FASES, random= ~ ANIMAIS,
family='ordinal',pl=TRUE,data=du)
summary(mod1)
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.
#Absurd results!
latentv <- mean(mod1$Liab)
cutpoint <- mean(mod1$CP)
pnorm(-(latentv), 0, sqrt(2))
pnorm(cutpoint - (latentv),0, sqrt(2)) - pnorm((latentv),0, sqrt(2))
1- pnorm(cutpoint - (latentv),0, sqrt(2))
#this would have a logical outcome to some extent, more would just be
referring to category '1' of the dependent variable. And the other?
bre <-
c(mean(mod1$Liab),mean(mod1$Sol[,1]),mean(mod1$Sol[,2]),mean(mod1$Sol[,3]))
pnorm(bre[2])-pnorm(bre[1])
pnorm(bre[3])-pnorm(bre[2])
pnorm(bre[4])-pnorm(bre[3])# negative probability?
--
View this message in context:
http://r.789695.n4.nabble.com/Help-ordinal-mixed-model-tp4501943p4501943.html
Sent from the R help mailing list archive at Nabble.com.
