12 views

Skip to first unread message

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)

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)

}

}

})

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.

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

Search

Clear search

Close search

Google apps

Main menu