[R-br] Ajuda modelo misto ordinal - MCMCglmm.

[Ivan Bezerra Allaman]

#Boa noite senhores! Após vários dias estudando sobre dados categóricos e pesquisando (vários fóruns com respostas do dono da library - todas incipientes) como interpretar o 'output' da função MCMCglmm venho recorrer a ajuda dos senhores, caso alguém já tenha trabalhado com função MCMCglmm no caso de variáveis dependentes ordinais. Já li e reli todos os pdf's do pacote, coursenotes do Jarrod, enfim, estou esgotado. Para esclarecer a base de dados, os tratamentos (denominado FASES) consistem de 3 níveis (1-pre, 2-própolis, 3-vincris), cujo as variável resposta ordinal tem 3 categorias (1-normal, 2-agudo, 3-cronico). Vide tabela!

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

#Um modelo misto seria:
set.seed(1)
mod1 <- MCMCglmm(ALT.RENAIS ~-1+FASES, random= ~ ANIMAIS,
    family='ordinal',pl=TRUE,data=du)
summary(mod1)
#Aí começa o sofrimento, pois a documentação é insuficiente neste caso. Segundo o próprio Jarrod (fóruns), as médias a posteriori dos coeficientes das covariáveis estão na escala probit. Segundo o meu estudo, estes coeficientes são os scores da distribuição normal padrão. Mais os escores não deveriam corresponder aos cutpoints? Neste caso, teríamos j(categorias da variável resposta)-1 cutpoints, ou seja, 2 cutpoints. O output só me apresenta um cutpoint. Como é possível então calcular as probabilidades com apenas um cutpoint? Segundo a documentação(Vignettes, pag 22), se P(y=k)=F(yk|l(vlatente),1) - F(yk-1|l,1), esse 'um' provavelmente seria a categoria 'um' da variável dependente? Enfim senhores, como posso extrair as probabilidades para as fases referentes a cada categoria da variável dependente? Desde já agradeço a atenção de todos.

#embora um pouco lógico, os resultados são absurdos!
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))

#isso teria algum resultado lógico até certo ponto, mais seria apenas referente a categoria 1 da variável dependente. E as outras?
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])#probabilidade negativa nãooo!



 

\begin{signature}
<<>>=
Prof. Dr. Ivan Bezerra Allaman
Universidade Estadual de Santa Cruz
Departamento de Ciências Exatas e Tecnológicas
Ilhéus/BA - Brasil
Fone: +55 73 3680-5076
E-mail: ivanalaman em yahoo.com.br/ivanalaman em gmail.com
@
\end{signature}
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120320/5728f178/attachment.html>