Usted puede caber multinivel GLMM con una distribución de Poisson (con la sobre-dispersión) el uso de R en múltiples formas. Pocos R paquetes: lme4, MCMCglmm, arm, etc. Una buena referencia para ver es Gelman y Hill (2007)
Voy a dar un ejemplo de esto rjags paquete en R. Es una interfaz entre R y JAGS (como OpenBUGS o WinBUGS).
$$n_{ij} \sim \mathrm{Poisson}(\theta_{ij})$$
$$\log \theta_{ij} = \beta_0 + \beta_1 \mbox{ } \mathtt{Tratamiento}_{i} + \delta_{ij}$$
$$\delta_{ij} \sim N(0, \sigma^2_{\epsilon})$$
$$i=1 \ldots I, \quad j = 1\ldots J$$
$\mathtt{Tratamiento}_i = 0 \mbox{ o } 1, \ldots, J-1 \mbox{ si } i^{th} \mbox{ observación pertenece al grupo de tratamiento } 1 \mbox {o, } 2, \ldots, J$
El $\delta_{ij}$ parte en el código de los modelos anteriores de la sobredispersión. Pero no hay nadie que te detenga de modelado de correlación entre los individuos (que no creo que las personas son realmente independientes), y dentro de los individuos (medidas repetidas). Asimismo, la tasa parámetro, pueden ser escaladas por alguna otra constante como en rate models. Por favor, consulte Gelman y Hill (2007) para más referencia. Aquí es el JAGS código para el modelo simple:
data{
for (i in 1:I){
ncount[i,1] <- obsTrt1[i]
ncount[i,2] <- obsTrt2[i]
## notice I have only 2 treatments and I individuals
}
}
model{
for (i in 1:I){
nCount[i, 1] ~ dpois( means[i, 1] )
nCount[i, 2] ~ dpois( means[i, 2] )
log( means[i, 1] ) <- mu + b * trt1[i] + disp[i, 1]
log( means[i, 2] ) <- mu + b * trt2[i] + disp[i, 2]
disp[i, 1] ~ dnorm( 0, tau)
disp[i, 2] ~ dnorm( 0, tau)
}
mu ~ dnorm( 0, 0.001)
b ~ dnorm(0, 0.001)
tau ~ dgamma( 0.001, 0.001)
}
Aquí es el R código para implementar el uso (digamos que es el nombre: overdisp.bug)
dataFixedEffect <- list("I" = 10,
"obsTrt1" = obsTrt1 , #vector of n_i1
"obsTrt2" = obsTrt2, #vector of n_i2
"trt1" = trt1, #vector of 0
"trt2" = trt2, #vector of 1
)
initFixedEffect <- list(mu = 0.0 , b = 0.0, tau = 0.01)
simFixedEffect <- jags.model(file = "overdisp.bug",
data = dataFixedEffect,
inits = initFixedEffect,
n.chains = 4,
n.adapt = 1000)
sampleFixedEffect <- coda.samples(model = simFixedEffect,
variable.names = c("mu", "b", "means"),
n.iter = 1000)
meansTrt1 <- as.matrix(sampleFixedEffect[ , 2:11])
meansTrt2 <- as.matrix(sampleFixedEffect[ , 12:21])
Usted puede jugar con sus parámetros posteriores y se pueden introducir parámetros más para hacer el modelado más preciso (nos gusta pensar que este). Básicamente, usted consigue la idea.
Para más detalles sobre el uso de rjags y JAGS, véase Juan Myles Blanco de la página
