## this is a simplified model
## mu = beta_0jk + beta_1jk*x
## where the betas are correlated and have some mean.

library(mvtnorm)
n = 100
K = 10
rho = 0.5
Sigma = matrix(c(1, rho, rho, 1), 2, 2)
mean = c(1, 2)

## we have-to/can write this as
## mu = mean.0 + mean.1*x + beta.zero.mean_0jk + beta.zero.mean_1jk*x
## its kind of clean to use the model f(.., model="iid") to
## generate gaussians...
hyper.fixed = list(prec = list(initial = -5, fixed=TRUE))
formula = y ~ -1 + f(mean.0,  model="iid", hyper = hyper.fixed) +
                   f(mean.1,  x,  model="iid", hyper = hyper.fixed) +
                   f(beta.0, model="iid2d", n = 2*n) +
                   f(beta.1, x, copy="beta.0")

## simulate some data from the model and generate the indices for the
## formula

idx = 1
mean.0 = c()
mean.1 = c()
beta.0 = c()
beta.1 = c()
x = c()
y = c()

for(i in 1:n) {
    beta = rmvnorm(1, mean=mean, sigma = Sigma)
    for(k in 1:K) {
        xx = rnorm(1)
        mu = beta[1] + beta[2]*xx
        yy = mu + rnorm(1)

        mean.0 = c(mean.0, 1)
        mean.1 = c(mean.1, 1)
        beta.0 = c(beta.0, i)
        beta.1 = c(beta.1, i + n) ## since beta.1's are the second
                                  ## half of the vector; see the doc
                                  ## for this model
        x = c(x, xx)
        y = c(y, yy)
    }
}

r = inla(formula,
        data = data.frame(mean.0, mean.1, beta.0, beta.1, x, y),
        family = "gaussian")

