Extending regression model to multiple subjects with shared sigma

12 views
Skip to first unread message

Nathan Constantine-Cooke

unread,
May 6, 2021, 8:24:31 AMMay 6
to nimble-users
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

unread,
May 6, 2021, 11:21:57 AMMay 6
to Nathan Constantine-Cooke, nimble-users
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

unread,
May 11, 2021, 8:43:12 AMMay 11
to nimble-users
Thank you very much! This was the approach I needed and the model now works perfectly. 

Kind regards, Nathan

Reply all
Reply to author
Forward
0 new messages