Clustering Variables Problem

[Matthew Carli]

Hi all,

I’m trying to fit a model where variables are grouped into a number of clusters which are then used as indices in an index regression. At first I was interested in the number of clusters being determined by the model, but decided to first work on a model where the number of clusters is fixed for simplicity’s sake. This means that the number of variables in any one cluster isn’t known a priori. In trying to model this I seemed to be breaking some of the nimble indexing conventions (I could not hardcode vector dimensions of the groups, as I did not know them), so I tried to make these group vectors all be the length of the total number of chemicals (J), and then “zero out” chemicals that did not belong. I’m new to nimble, so I’m not sure if what I’m trying to do is even possible, or if there is some nimble trick and/or functionality I should use instead of my current approach. Here is what my model looks like right now:

This nimble function uses a vector of group membership  (gnum) to replace values of chemicals that do not belong in a group:

nimble_which <- nimbleFunction(

    run = function(vec1 = double(1), gnum = double(1), index_num = double(0), val = double(0)) {

        returnType(double(1))

            vec1[which(gnum!=index_num)] <- val

        return(vec1)

})

 

testCode <- nimbleCode({

    for (i in 1:N)           # Index Regression

    {

        y[i]  ~ dbern(max(1e-16, min(1 - 1e-16, p[i])))        

        logit(p[i]) <- beta0 + beta1 * index1[i] + beta2 * index2[i] + beta3 * index3[i]

        index1[i] <- inprod(w1[1:J], q[i,J])                # Forming the indices in 3-index model

        index2[i] <- inprod(w2[1:J], q[i,J])

        index3[i] <- inprod(w3[1:J], q[i,J])

    }

    w1[1:J] <- nimble_which(weight1[1:J], d[1:J], 1, 0) # Zeroing weights of variables not in group

    w2[1:J] <- nimble_which(weight2[1:J], d[1:J], 2, 0)

    w3[1:J] <- nimble_which(weight3[1:J], d[1:J], 3, 0)

 

    weight1[1:J] ~ ddirch(a1[1:J])

    weight2[1:J] ~ ddirch(a2[1:J])

    weight3[1:J] ~ ddirch(a3[1:J])

 

    a1[1:J] <- nimble_which(alpha1[1:J], d[1:J], 1, .0001) # Small value on vars not in group

    a2[1:J] <- nimble_which(alpha2[1:J], d[1:J], 2, .0001)

    a3[1:J] <- nimble_which(alpha3[1:J], d[1:J], 3, .0001)

   

    for (j in 1:J){

        d[j] ~ dcat(pi[1:K])              # vector of K group membership for J variables

    }

   

### Priors

    pi[1:K] ~ ddirch(alpha[1:K])

    # Intercept

    beta0 ~ dnorm(0, tau0)

    tau0 <- 1/sqrt(sigma0)

    sigma0 ~ dunif(0,100)

    # Index Effects

    beta1 ~ dnorm(0, tau1)

    tau1 <- 1/pow(sigma1,2)

    sigma1 ~ dunif(0,100)

    beta2 ~ dnorm(0, tau2)

    tau2 <- 1/pow(sigma2,2)

    sigma2 ~ dunif(0,100)

    beta3 ~ dnorm(0, tau3)

    tau3 <- 1/pow(sigma3,2)

    sigma3 ~ dunif(0,100)

    # Alphas

    for (j in 1:J){

    alpha1[j] <- exp(delta1[j])

    delta1[j] ~ dnorm(0, taud1)

 

    alpha2[j] <- exp(delta2[j])

    delta2[j] ~ dnorm(0, taud2)

 

    alpha3[j] <- exp(delta3[j])

    delta3[j] ~ dnorm(0, taud3)

    }

    taud1 <- 1/sqrt(sigmad1)

    sigmad1 ~ dunif(0,10)

    taud2 <- 1/sqrt(sigmad2)

    sigmad2 ~ dunif(0,10)

    taud3 <- 1/sqrt(sigmad3)

    sigmad3 ~ dunif(0,10)

})

I am able to create the model object, but when I then go to compile it I get the message, Error: Failed to create the shared library. Run 'printErrors()' to see the compilation errors. The printErrors() output is quite long, so I’ve attached it in a text file below my post. Hopefully I was clear when describing my problem, if you need any more info or the data I’m working with please let me know!

Best,

Matt

[Daniel Turek]

Matthew, thanks for your use of NIMBLE.  This sounds like an interesting application, and I'm glad to see you're doing some programming with NIMBLE.

First let me say, without a reproducible example, it's hard to diagnose for sure.  But that said, I think your use of the inprod() function is giving at least some problems here.  I think the two arguments to the inprod() function will both have to be vectors of the same length, in order to compute the inner-product between the two vectors.  That is, instead of:

        index2[i] <- inprod(w2[1:J], q[i,J])

the proper use of inprod() should look more like:

        index2[i] <- inprod(w2[1:J], q[i,1:J])

I didn't see q[] elsewhere in the model, so I'm not sure if q[] is a 2D-array or 3D-array, and I know my change above might not "jive" with your other code. But it's only to serve as an example.

If indeed you wanted to multiple the scalar q[i,J] by each element of the entire vector w2[1:J], then sum the resulting vector, then rather try:

    index2[i] <- sum(w2[1:J] * q[i,J])

Can you try correcting this - either using the inprod() function, or not - and see if it works?

If you're still having problems, then sending a reproducible example (perhaps off-list, if necessary) would help immensely.

Cheers,

Daniel

--
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/20404f66-f9e4-437a-82c0-e8849e7007b4n%40googlegroups.com.