How to fit models with multivariate outcome with INLA

[yetanotherd...@gmail.com]

Hello,

I am trying to fit a model with three outcomes. Since my data is multilevel and bayesian regression modeling gives more flexibility, I am trying to fit a multilevel model using INLA. 

I have three outcomes. I was wondering if I could fit a single multivariate model. I went through help pages with examples of multiple likelihood and also referred to another paper that uses bivariate outcome and applied INLA: "A latent process model for forecasting multiple time series in environmental public health surveillance", Morrison et al 2016.

My code snippet for implement is as follows:

n <- nrow(data_in) 

Y <- matrix(NA, n*3, 3)

Y[1:n, 1] <- data_in$Outcome1 

Y[(1+n):(3*n), 1] <- NA # faked observations

Y[(1+n):(2*n), 2] <- data_in$Outcome2   # actual observations

Y[c(1:n,((2*n+1):(3*n)) ), 2] <- NA # faked observations

Y[(1+2*n):(3*n), 3] <- data_in$Outcome3   # actual observations

Y[c(1:(2*n) ), 3] <- NA # faked observations

x <- data_in[,c("P_ID", "var1","var2","var3")]

X <- rbind(x,x,x)

Y <- data.frame(Y)

names(Y) <- c("Outcome1","Outcome2","Outcome3")

formula_tv = list(Outcome1, Outcome2, Outcome3) ~ var1 + var2 + var3 +  f(P_ID, model="iid")

system.time( 

result_tv <- inla(formula_tv, family = c("Gaussian","Gaussian","Gaussian"), data = cbind(Y,X), control.compute=list(dic=TRUE) )

)

summary(result_tv)

I guess the above model assumes that outcome is a trivariate normal distribution and returns single set of coefficients for the latent field.

I was wondering if there is a way to get three sets of coefficients, one for each outcome; similar to a multivariate regression model, that also accounts for covariance between outcomes.

Thank you,

Karthik S

[Haakon Bakka]

In INLA, there is one formula for everything, which means that you have to hack together a formula, e.g.

~ var1.1 + var1.2 + var1.3

where you use many NA values to separate var1 into these three var's

To get covariance between outcomes, I suggest trying to share the iid model between the three familys, e.g. P_ID = rep(1:n, 3).

Does anyone have a simple code example for this?

Kind regards,

Haakon

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.