apply the same nimble model to multiple dataset

[Juan Xie]

Hello Everyone,

I  fit a Poisson Gaussian process model using nimble and I want to apply that model to different data. 

I read the discussion on https://groups.google.com/g/nimble-users/c/k63xGcuNbpc and in my first attempt, I tried using R's assignment syntax to assign new data to the compiled model. To check if my attempt was correct or not, I tried a naïve way, i.e., I fit the model using one dataset, save the results, then clear every thing,  then fit the model using a second dataset. I compared the results from the two attempts, which are different. Could you please help and see what might be wrong? (I thought my second attempt should be correct, but not sure why the first attempt gave different results).

Here are my code:

############ first attempt

## 

test.data <- data.frame (gene1 = gene1, gene2 = gene2, x = x, y = y )

## nimble 

expcov <- nimbleFunction(    
  run = function(dists = double(2), rho = double(0), alpha = double(0)) {
    returnType(double(2))
    n <- dim(dists)[1]
    result <- matrix(nrow = n, ncol = n, init = FALSE)
    alpha2 <- alpha*alpha
    for(i in 1:n){
      for(j in 1:n)
        result[i, j] <- alpha2*exp(-dists[i,j]/rho)
    }
    for(i in 1:n)
      result[i,i] <- result[i,i] + 1e-6
    return(result)
  })
cExpcov <- compileNimble(expcov)


code <- nimbleCode({
  mu0 ~ dnorm(0, sd = 100)
  alpha ~ dnorm(0, sd = 0.1)  
  rho ~ dinvgamma(5, 5)
 
  mu[1:N] <- mu0*ones[1:N]
  cov[1:N, 1:N] <- expcov(dists[1:N, 1:N], rho, alpha)
  s[1:N] ~ dmnorm(mu[1:N], cov = cov[1:N, 1:N])
  # likelihood
  for(i in 1:N) {
    lambda[i] <- exp(s[i])
    y[i] ~ dpois(lambda[i])
  }
})

dists <- as.matrix(dist(cbind(test.data$x, test.data$y)))
n <- nrow(test.data)
constants <- list(N = n, dists = dists, ones = rep(1, n))

data <- list(y = test.data$gene1)
inits <- list(alpha = 1, rho = 0.2)

inits$cov <- cExpcov(dists, inits$rho, inits$alpha)
inits$s <-  t(chol(inits$cov)) %*% rnorm(n)
inits$s <- inits$s[ , 1]  # so can give nimble a vector rather than one-column matrix

model <- nimbleModel(code, constants = constants, data = data, inits = inits)
cModel <- compileNimble(model)
conf <- configureMCMC(model)
MCMC <- buildMCMC(conf)
cMCMC <- compileNimble(MCMC, project = cModel)

gene1.rst  <- runMCMC(cMCMC, niter = 10000, nburnin = 5000,setSeed = 12345)

## provide new data

cModel$y <- test.data$gene2

gene2.rst <- runMCMC(cMCMC,niter = 10000, nburnin = 5000, setSeed = 12345)

### for my second attempt, after runMCMC for gene1, I clear everything and then repeat the the whole process, except that I set  data <- list(y = test.data$gene2)

Any comments or suggestions would be greatly appreciated.

[Chris Paciorek]

Hi Juan,

I think your questions have been addressed in the thread that includes posts from Bill Chen and from you, but let us know if you still have questions.

-chris

--
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/0d644aec-8630-4463-9a92-14a498ad30f1n%40googlegroups.com.

[Juan Xie]

Hi Chris,

Thank you very much for the reply. I saw your suggestion on the other post and will try as suggested.

Best,

Juan