LQ factorization in Nimble

[Eva Murphy]

Dear all,

I am fairly new at using Nimble, but I'd like to use it to estimate a factor model we're currently working on. In our factor modeling we have a loadings matrix (matrix of weights) that connects the outcomes to the factors. In our case, this loadings matrix is a sparse matrix with some none zero elements for which we have priors. What I'd like to accomplish is to decompose the loadings matrix into L and Q matrices and use L to fit the factor model. The code I wrote (given below the error message) errors out when I run the runMCMC function with the error message included below. It looks like that maybe my nimble Functions are not written correctly. I was hoping that I could get  some help on to fix this issue. 

Thank you very much!

Eva

## The error

running chain 1...
  [Note] Ignoring non-NA values in `inits` for data nodes: y.T.
Error in nimQR.R(x = c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,  :
  unused argument (x = c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0))

### Part of my code

.....

load.mat = matrix(c(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,0, 0,0,0,1, 0,0,0,0), nrow= 6, ncol = 4, byrow = T)

nimQR.Q = nimbleFunction(
  run = function(A = double(2)){
    returnType(double(2))
    A1 = t(A)
    m <- dim(A1)[1]
    n <- dim(A1)[2]
    Q <- matrix(0, m, n)
    R <- matrix(0, n, n)
    for (k in 1:n) {
      Q[1:m, k] <- A1[1:m, k]
      if (k > 1) {
        for (i in 1:(k - 1)) {
          R[i, k] <- t(Q[1:m, i]) %*% Q[1:m, k]
          Q[1:m, k] <- Q[1:m, k] - R[i, k] * Q[1:m, i]
        }
      }
      R[k, k] <- sum(abs(Q[, k])^2)^(1/2)
      Q[1:m, k] <- Q[1:m, k]/R[k, k]
    }
    return(Q)
})

nimQR.R = nimbleFunction(
  run = function(A = double(2)){
    returnType(double(2))
    A1 = t(A)
    m <- dim(A1)[1]
    n <- dim(A1)[2]
    Q <- matrix(0, m, n)
    R <- matrix(0, n, n)
    for (k in 1:n) {
      Q[1:m, k] <- A1[1:m, k]
      if (k > 1) {
        for (i in 1:(k - 1)) {
          R[i, k] <- t(Q[1:m, i]) %*% Q[1:m, k]
          Q[1:m, k] <- Q[1:m, k] - R[i, k] * Q[1:m, i]
        }
      }
      R[k, k] <- sum(abs(Q[1:m, k])^2)^(1/2)
      Q[1:m, k] <- Q[1:m, k]/R[k, k]
    }
    return(t(R))
   
  })

cNimQR.Q = nimbleRcall(prototype = function(x = double(2)){},
                       Rfun = 'nimQR.Q',
                       returnType = double(2)) #compileNimble(nimQR.Q)
cNimQR.R = nimbleRcall(prototype = function(x = double(2)){},
                       Rfun = 'nimQR.R',
                       returnType = double(2)) # compileNimble(nimQR.R)
model_code=nimbleCode({
 
  for (t in T0:T0){
    for (i in 1:n){
      cen[n*(t-T0)+i] ~ dinterval(y.T[n*(t-T0)+i],bd[n*(t-T0)+i,1:2])
      y.D[n*(t-T0)+i] ~ dpois(ED[(n*(t-T0)+i)]*lambdaD[n*(t-T0)+i])
      y.E[n*(t-T0)+i] ~ dpois(EE[(n*(t-T0)+i)]*lambdaE[n*(t-T0)+i])
      y.T[n*(t-T0)+i] ~ dpois(ET[(n*(t-T0)+i)]*lambdaT[n*(t-T0)+i])
      y.B[n*(t-T0)+i] ~ dpois(EB[(n*(t-T0)+i)]*lambdaB[n*(t-T0)+i])
      y.C[n*(t-T0)+i] ~ dpois(EC[(n*(t-T0)+i)]*lambdaC[n*(t-T0)+i])
      y.I[n*(t-T0)+i] ~ dpois(EI[(n*(t-T0)+i)]*lambdaI[n*(t-T0)+i])
     
      log(lambdaD[n*(t-T0)+i]) <- inprod(L[1,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VD[(n*(t-T0)+i)]
      log(lambdaE[n*(t-T0)+i]) <- inprod(L[2,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VE[(n*(t-T0)+i)]
      log(lambdaT[n*(t-T0)+i]) <- inprod(L[3,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VT[(n*(t-T0)+i)]
      log(lambdaB[n*(t-T0)+i]) <- inprod(L[4,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VB[(n*(t-T0)+i)]
      log(lambdaC[n*(t-T0)+i]) <- inprod(L[5,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VC[(n*(t-T0)+i)]
      log(lambdaI[n*(t-T0)+i]) <- inprod(L[6,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VI[(n*(t-T0)+i)]
     
     
      VD[n*(t-T0)+i] ~ dnorm(0,tau.VD)
      VE[n*(t-T0)+i] ~ dnorm(0,tau.VE)
      VT[n*(t-T0)+i] ~ dnorm(0,tau.VT)
      VB[n*(t-T0)+i] ~ dnorm(0,tau.VB)
      VC[n*(t-T0)+i] ~ dnorm(0,tau.VC)
      VI[n*(t-T0)+i] ~ dnorm(0,tau.VI)
      mu[n*(t-T0)+i,1:F] <- mu0[t-T0+1,1:F]
    }
   
   
   
    # ICAR prior for the spatial factors
    U[(n*(t-T0)+1):(n*(t-T0)+n),1] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
    U[(n*(t-T0)+1):(n*(t-T0)+n),2] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
    U[(n*(t-T0)+1):(n*(t-T0)+n),3] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
    U[(n*(t-T0)+1):(n*(t-T0)+n),4] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
   
   
    #tau.U[t-(T0-1),1] ~ dgamma(0.5,0.5)
    #tau.U[t-(T0-1),2] ~ dgamma(0.5,0.5)
    #tau.U[t-(T0-1),3] ~ dgamma(0.5,0.5)
    #tau.U[t-(T0-1),4] ~ dgamma(0.5,0.5)
   
  }  
 
  for (t in (T0+1):(T0+TT-1)){
    for (i in 1:n){
      cen[n*(t-T0)+i] ~ dinterval(y.T[n*(t-T0)+i],bd[n*(t-T0)+i,1:2])
      y.D[n*(t-T0)+i] ~ dpois(ED[(n*(t-T0)+i)]*lambdaD[n*(t-T0)+i])
      y.E[n*(t-T0)+i] ~ dpois(EE[(n*(t-T0)+i)]*lambdaE[n*(t-T0)+i])
      y.T[n*(t-T0)+i] ~ dpois(ET[(n*(t-T0)+i)]*lambdaT[n*(t-T0)+i])
      y.B[n*(t-T0)+i] ~ dpois(EB[(n*(t-T0)+i)]*lambdaB[n*(t-T0)+i])
      y.C[n*(t-T0)+i] ~ dpois(EC[(n*(t-T0)+i)]*lambdaC[n*(t-T0)+i])
      y.I[n*(t-T0)+i] ~ dpois(EI[(n*(t-T0)+i)]*lambdaI[n*(t-T0)+i])
     

      log(lambdaD[n*(t-T0)+i]) <- inprod(L[1,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VD[(n*(t-T0)+i)]
      log(lambdaE[n*(t-T0)+i]) <- inprod(L[2,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VE[(n*(t-T0)+i)]
      log(lambdaT[n*(t-T0)+i]) <- inprod(L[3,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VT[(n*(t-T0)+i)]
      log(lambdaB[n*(t-T0)+i]) <- inprod(L[4,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VB[(n*(t-T0)+i)]
      log(lambdaC[n*(t-T0)+i]) <- inprod(L[5,1:F],U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VC[(n*(t-T0)+i)]
      log(lambdaI[n*(t-T0)+i]) <- inprod(L[6,1:F], U[n*(t-T0)+i,1:F]+mu[n*(t-T0)+i,1:F]) + VI[(n*(t-T0)+i)]
      VD[n*(t-T0)+i] ~ dnorm(0,tau.VD)
      VE[n*(t-T0)+i] ~ dnorm(0,tau.VE)
      VT[n*(t-T0)+i] ~ dnorm(0,tau.VT)
      VB[n*(t-T0)+i] ~ dnorm(0,tau.VB)
      VC[n*(t-T0)+i] ~ dnorm(0,tau.VC)
      VI[n*(t-T0)+i] ~ dnorm(0,tau.VI)
      mu[n*(t-T0)+i,1:F] <- mu0[t-T0+1,1:F]
    }
   
   
    # ICAR prior for the spatial factors
    U[(n*(t-T0)+1):(n*(t-T0)+n),1] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
    U[(n*(t-T0)+1):(n*(t-T0)+n),2] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
    U[(n*(t-T0)+1):(n*(t-T0)+n),3] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
    U[(n*(t-T0)+1):(n*(t-T0)+n),4] ~ dcar_normal(adj[], weights[], num[], 1, zero_mean=1)
   
    #tau.U[t-(T0-1),1] ~ dgamma(0.5,0.5)
    #tau.U[t-(T0-1),2] ~ dgamma(0.5,0.5)
    #tau.U[t-(T0-1),3] ~ dgamma(0.5,0.5)
    #tau.U[t-(T0-1),4] ~ dgamma(0.5,0.5)
   
  }
 
  #Priors
 
  ### loadings matrix entries
  for(f in 2:K){
    load.mat[f,1] ~ dnorm(0,sd=10)
  }
  load.mat[1,2] ~ dnorm(0,sd=10)
  load.mat[4,3] ~ dnorm(0,sd=10)
  load.mat[6,4] ~ dnorm(0,sd=10)
 
  for(t in T0:(T0+TT-1)){
    mu0[t-T0+1,1] ~ dflat()
    mu0[t-T0+1,2] ~ dflat()
    mu0[t-T0+1,3] ~ dflat()
    mu0[t-T0+1,4] ~ dflat()
  }
 
  tau.VD ~ dgamma(0.5,0.5)
  tau.VE ~ dgamma(0.5,0.5)
  tau.VT ~ dgamma(0.5,0.5)
  tau.VB ~ dgamma(0.5,0.5)
  tau.VC ~ dgamma(0.5,0.5)
  tau.VI ~ dgamma(0.5,0.5)
  eta[1] ~ dunif(0,1)
  eta[2] ~ dunif(0,1)
  eta[3] ~ dunif(0,1)
  eta[4] ~ dunif(0,1)
 
  #Compute L and Q from the LQ factorization
  L[1:6, 1:6]<- cNimQR.R(load.mat[1:6, 1:4])
  Q[1:4, 1:6]<- cNimQR.Q(load.mat[1:6, 1:4])
 
 
})

mod_constants=list(bd=bd,n=n,TT=TT,T0=min(Year),num=num,adj=adj,weights=weights,K=6,F=4)

mod_data=list(cen=cen,y.D=y.D,y.E=y.E,y.T=y.T,y.B =y.B,y.C=y.C,y.I=y.I,ED=ED,EE=EE,ET=ET,EB=EB,EC=EC,EI=EI)


#glm0 = glm(outcomedata$AllDeaths~1,offset=log(outcomedata$TotalPopulation))

mod_inits=list(U=matrix(0,n*TT,4),load.mat=load.mat,mu0=matrix(0,TT,4),VD=rep(0,n*TT),VE=rep(0,n*TT),VT=rep(0,n*TT),VB=rep(0,n*TT),VC=rep(0,n*TT),VI=rep(0,n*TT),y.T=floor(.5*(bd[,1]+bd[,2])))


# Build the model.
nim_model <- nimbleModel(model_code, mod_constants,mod_data,mod_inits)
compiled_model <- compileNimble(nim_model,resetFunctions = TRUE)

# Set up samplers.
mcmc_conf <- configureMCMC(nim_model,monitors=c("mu0","load.mat", "L", "Q","mu",
                                                "U","eta","VD","VE","VT","VB","VC","VI","lambdaD","lambdaE","lambdaT","lambdaB","lambdaC","lambdaI"),useConjugacy = TRUE)

mod_mcmc<-buildMCMC(mcmc_conf)
compiled_mcmc<-compileNimble(mod_mcmc, project = nim_model,resetFunctions = TRUE)

# Run the model
MCS=5000
st<-Sys.time()
samples=runMCMC(compiled_mcmc,inits=mod_inits,
                nchains = 1, nburnin=floor(MCS/2),niter = MCS,samplesAsCodaMCMC = TRUE,thin=50,
                summary = FALSE, WAIC = FALSE,progressBar=TRUE) 

[Eva Murphy]

I believe I fixed my error by changing the argument of the nimbleRcall function to match the argument of the nimbleFunction.

Eva

[Chris Paciorek]

Hi Eva,

I don't think you need `nimbleRcall` at all here. The case you would need `nimbleRcall` is when you want to use some functionality in R that nimble can't compile.

Here you are writing nimbleFunctions to do the QR calculation and nimble will just compile those nimbleFunctions. So if you just remove `cNimQR.Q` and  `cNimQR.R` and in your model code call `nimQR.Q` and `nimQR.R`, it should work.

On an efficiency note, if you start to work with bigger matrices, you will probably want to have a single function that returns both Q and R together rather than doing the same calculation twice. It's a bit of a hassle to do this because a nimble model can work with lists, but you could have a `nimQR` nimbleFunction that returns Q and R as part of a larger matrix (e.g., in your example you could have the result be a 10x6 matrix) and then in your model code extract the two pieces as needed (the first 4 rows and the last 6 rows in your example).

-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/2cf40ad1-e7e0-42e9-a986-154aeaa96892n%40googlegroups.com.

[Murphy, Eva]

Hi Chris,

Thank you for getting back to me so quickly. Indeed, not using `nimbleRcall` fixed my issue. 

Also, thank you for suggesting to use a larger matrix as the output of my nimbleFunction that includes both Q and R matrices. Initially, I wanted to use a list for the output, but I couldn't figure out how to do it. Binding the two matrices into one will work just fine.

Eva 

--

Eva Murphy, PhD

Statistics Department

Wake Forest University

https://evamurphy100.github.io/

[Chris Paciorek]

Great, glad things are working now.

As a note for any future readers, there was a typo in my response - I meant to say "a nimble model cannot work with lists".

-chris