Hi:

Thanks for the helpful comments from my last post. I am writing today with a new question.

I'd like to specify a bivariate Poisson model with two random effects. The responses are counts of clinical episodes of two different malaria species during an 8-week interval. There are a handful of time-dependent covariates that I would like to include on the RHS of the equation as well as a random effect for child (repeated measures) and one for village (clustering). There's also the complication of needing to add an offset for time-at-risk which subtracts off the time after treatment when a child is not at risk of becoming ill (I will worry about that later). There are 1867 observations from 264 children nested within 11 villages. 

My question is, can I fit this model in MCMCglmm or GLMER? I don't believe I can fit it in GLMER, but I should be able to fit it in MCMCglmm. I'd like to fit a saturated model, so I'd like an estimate of the effects of the covariates on "mixed clinical episodes", i.e. the interaction between counts of episodes between the two species. The two counts are negatively correlated, which I hope can be estimated using this package as well. Does anyone know of a good manuscript that has done something similar? I think what I am truly stuck on is writing down the model. How many parameters should there be on the RHS? An intercept and a slope for each fixed effect for each species? For example, if there is simply age on the RHS, do I end up with 4 parameters: 2 intercepts and 2 slopes for the two response categories? What about Y1*Y2? 

My first attempt at analyzing this is below (which I obviously should not do until I figure out how to write down the model, but perhaps this will help to understand what I am trying to do):

prior.bivar <- list(R = list(V = diag(2), nu = 0, fix=1), 
		 G = list(G1 = list(V = diag(2), n=2),
		 G2 = list(V = diag(2), n=2)))

model.bivar <- MCMCglmm(cbind(PFepisodes, PVepisodes) ~ -1 + trait + trait:Age,
		random = ~ us(trait):child + us(trait):village, 
		rcov = ~ us(trait):units, prior = prior.bivar, 
		family = c("poisson", "poisson"), 
		nitt = 60000, burnin = 10000,
		thin=25, data = MixedDat)

summary(model.bivar)

Thanks in advance!

Katie Colborn
benton@stat.berkeley.edu


	[[alternative HTML version deleted]]

