Extending regression model to multiple subjects with shared sigma

[Nathan Constantine-Cooke]

Hello all, 

My apologies if this has been discussed before, but I was unable to find documentation on the website or in this Google group for this topic and I am still new to Nimble.

I am interested in a model which fits subject-specific means with a population-level sigma. To make things more difficult, the number of observations  are not the same for each subject, so I cannot simply pass a matrix of observations. Whilst it is easy to fit this model by explicitly defining the relationships for each parameter, it is not feasible to do so for a large number of subjects. So my question is ultimately, is there a way to programmatically define each parameter in this situation? 

For simplicity, I have adapted the linear regression example from the website for two subjects. Here I have had to explicitly define the relationships. Is there an approach where I only need to define the relationship for one subject, and can then programmatically generate the parameters and their relationships for the other subjects?

Any and all guidance on this would be greatly appreciated. 

Kind regards, Nathan

library(nimble)

set.seed(1)

p <- 15    # number of explanatory variables

n1 <- 10   # number of observations for subject 1

n2 <- 12 # number of observations for subject 2

X1 <- matrix(rnorm(p * n1), nrow = n1, ncol = p) # explanatory variables

X2 <- matrix(rnorm(p * n2), nrow = n2, ncol = p)

true_betas1 <- c(c(0.1, 0.2, 0.3, 0.4, 0.5), rep(0, p - 5)) # coefficients

true_betas2 <- c(c(0.5, 0.4, 0.3, 0.2, 0.1), rep(0, p - 3))

sigma <- 1

y1 <- rnorm(n1, X1 %*% true_betas1, sigma)

y2 <- rnorm(n2, X2 %*% true_betas2, sigma)

code <- nimbleCode({

  beta10 ~ dnorm(0, sd = 100) # Beta0 for subject 1 etc

  beta11 ~ dnorm(0, sd = 100)

  beta12 ~ dnorm(0, sd = 100)

  beta20 ~ dnorm(0, sd = 100)

  beta21 ~ dnorm(0, sd = 100)

  beta22 ~ dnorm(0, sd = 100)

  sigma ~ dunif(0, 100)

  for (i in 1:n1) {

    y1[i] ~ dnorm(beta10 + beta11 * x11[i] + beta12 * x12[i], sd = sigma)

  }

  for (i in 1:n2) {

    y2[i] ~ dnorm(beta20 + beta21*x21[i] + beta22*x22[i], sd = sigma)

  }

})

## extract data for two predictors and center for better MCMC performance

x11 <- X1[, 1] - mean(X1[, 1])

x12 <- X1[, 2] - mean(X1[, 2])

x21 <- X2[, 1] - mean(X2[, 1])

x22 <- X2[, 2] - mean(X2[, 2])

constants <- list(n1 = n1,

                  n2 = n2,

                  x11 = x11,

                  x12 = x12,

                  x21 = x21,

                  x22 = x22)

data <- list(y1 = y1, y2 = y2)

inits <- list(beta10 = mean(y1),

              beta20 = mean(y2),

              beta11 = 0,

              beta12 = 0,

              beta21 = 0,

              beta22 = 0,

              sigma = 1)

nburnin <- 1000

thin <- 1

niter <- 10000

nimbleMCMC_samples_sim <- nimbleMCMC(code = code,

                                     constants = constants,

                                     data = data,

                                     inits = inits,

                                     nburnin = nburnin,

                                     thin = thin,

                                     niter = niter)

[Daniel Turek]

I'm sorry this is so brief, but I think this is what you're looking for (code below).

Make sure to provide "nObservationsOfEachSubject" in the "constants" list, and also here "y" is a ragged array, where only the values of y[i, ] up to the nObservationsOfEachSubject[i]-th value will be used, and the remaining values in that row of y[i, ] won't be used.

Also, the inprod() function is the dot-product of two vectors, the sum of the element-wise product of the vectors.

Hope this helps,

Daniel

code <- nimbleCode({
    sigma ~ dunif(0, 100)
    for(i in 1:nSubjects) {
        for(j in 1:(nPredictors+1)) {   ## number of predictors, plus 1 for intercept
            beta[i,j] ~ dnorm(0, sd = 100)
        }
        for(k in 1:nObservationsOfEachSubject[i]) {
            y[i,k] ~ dnorm(beta[i,nPredictors+1] + inprod(beta[i,1:nPredictors], X[i,1:nPredictors]), sd = sigma)
        }
    }
})

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/a77f528c-d2c9-4749-b44d-f74bcec3e86fn%40googlegroups.com.

[Nathan Constantine-Cooke]

Thank you very much! This was the approach I needed and the model now works perfectly. 

Kind regards, Nathan