R help es - Aug 2010 - gibbs

Hola a todos:

 

Tengo una duda despecto como usar gibbs y r. Me explicaré en un ejemplo
donde todos los datos se encuentran en “library(MCMCglmm)”. Al ejecutar el
código R que escribo en este correo, obtengo las soluciones (Solucion.m).
Pues bien, entiendo que utilizando HPDinterval (….) puedo obtener por
ejemplo la región donde se ubica el 95% de la probabilidad para …. Pues, me
gustaría conocer aparte del valor de solución, el valor de “confiabilidad”,
error, error estándar, …, para cada solución. Obteniendo un resultado final
como “Diert2,  3,64…, s.e. de ... ” (agregado a la lalida de Solucion.m del
código r que escribo a continuación). Por las dudas, HPDinterval de la
librería coda.

 

¿Alguien me puede ayudar al respecto?

 

Desde ya muchas gracias y el código R.

 

 

library(MCMCglmm)

library(coda)

data(ChickWeight)

head(ChickWeight)

xyplot(weight ~ Time | Chick, data = ChickWeight)

prior.m <- list(R = list(V=1, n=0.0002), G=list(G1=list(V= diag(3),n=3)))

m <- MCMCglmm(weight ~ Diet + poly(Time, 2, raw=TRUE),

                      random = ~us(1+ poly (Time, 2, raw=TRUE)):Chick,

                      data = ChickWeight, verbose = TRUE, pr=TRUE, pl=TRUE,

                      prior = prior.m, saveX = TRUE, saveZ = TRUE)


pollos<-posterior.mode(m$Sol[,c(7:56)])

W <- cBind(m$X, m$Z)

prediccion <- W %*% posterior.mode(m$Sol)

xyplot(weight+prediccion ~Time |Chick, data=ChickWeight)

Solucion<-posterior.mode(m$Sol)

Solucion

Solucion.m <-as.matrix(Solucion)

Solucion.m


	[[alternative HTML version deleted]]

Javier Marcuzzi vas escriure el dia dc, 11 ago 2010:
> Tengo una duda despecto como usar gibbs y r. Me explicar? en un ejemplo
> donde todos los datos se encuentran en ?library(MCMCglmm)?. Al ejecutar el
> c?digo R que escribo en este correo, obtengo las soluciones (Solucion.m).
> Pues bien, entiendo que utilizando HPDinterval (?.) puedo obtener por
> ejemplo la regi?n donde se ubica el 95% de la probabilidad para ?. Pues, me
> gustar?a conocer aparte del valor de soluci?n, el valor de ?confiabilidad?,
> error, error est?ndar, ?, para cada soluci?n. Obteniendo un resultado final
> como ?Diert2,  3,64?, s.e. de ... ? (agregado a la lalida de Solucion.m del
> c?digo r que escribo a continuaci?n). Por las dudas, HPDinterval de la
> librer?a coda.
> 
m es un objeto de tipo "mcmc". Se trata simplemente de una gran matrix
con
las simulaciones en las filas y los parámetros estimados en las columnas.
De esta manera, usando quantile() puedes obtener todo lo que quieras:

quantile(m$Sol, c(0.025, 0.05, 0.5, 0.95, 0.975) te daría la mediana de tu
estimador, y los intervalos de credibilidad al 90 y al 95 %.

Un simple sd(m$Sol) devuelve la desviación.

Un apply() aplicado a quantile() siempre es muy útil:
apply(m, 2, quantile, c(0.025, 0.5, 0.975))

Y el error para cada solución:
apply(m, 2, sd)

Un saludo,

-- 
-  Xavier  -