dimension limits on dmnorm?

Help
2010-09-23
2013-12-13
  • Woodrow Setzer

    Woodrow Setzer - 2010-09-23

    I am modeling data from high throughput in vitro bioassays, which are often
    carried out on microtiter plates (say 8 rows and 12 columns). In the data set
    I'm currently analyzing, in addition to plate effects, there are clearly fixed
    row and column effects, with residuals being autocorrelated, based on a set of
    control plates. I am trying to model this with a multivariate normal
    likelihood with an exponential correlation matrix, and an additional parameter
    to allow the correlation to decay faster along one or the other dimension of
    the plate (so, correlation among cells in the same column may decay more
    slowly than correlations among cells in the same row).

    I am simulating data, so I can get an idea of how well this should work for
    the design of my data set, which is a bit 'wacky', but for now, I am just
    working with simulating sets of control plates.

    All goes well for 2 X 3 plates, but for larger, say, 3 X 3,I get an error
    like:

    Welcome to JAGS 2.1.0 on Thu Sep 23 10:31:09 2010

    JAGS is free software and comes with ABSOLUTELY NO WARRANTY

    Loading module: basemod

    Loading module: bugs

    Reading data file data.R

    Compiling model graph

    Declaring variables

    Resolving undeclared variables

    Allocating nodes

    Graph Size: 1283

    Reading parameter file parameters.R

    Error in node y

    Invalid parent values

    Deleting model

    It looks as if I am exceeding some size limit, but getting an inappropriate
    error.

    Am I right? If so, is there a way to expand the size limit? I will be working
    with at least 8 X 12 plates, so the correlation matrix will need to be 96 X
    96.

    Thanks,

    Woody Setzer

    Below are model.bugs, the bugs model, runmods.bugs, a JAGS script for running
    the model, and makedata.R, an R script for creating all the files (other than
    bugs.model) required by runmods.bugs.

    ------------ model.bugs

    var

    Sigma, A, R, C, mu;

    model

    {

    Construct the covariance matrix

    for (i in 1:Ncells) {

    for (j in 1:Ncells) {

    Sigma <-

    sigma2*(equals(i,j) +

    (1 - equals(i,j))*

    (exp(-rho*sqrt((Ix_-Ix)^2 +

    theta*(Jx_-Jx)^2)))

    )

    }

    }

    Omega <- inverse(Sigma)

    for (i in 1:Np) {

    y ~ dmnorm(mu, Omega)

    }

    for (i in 1:Np) {

    Construct the mean vector, mu

    for (j in 1:Ncells) {

    mu <- A_ + R[Ix] + C[Jx]

    }

    }

    for (i in 2:Nr) {

    R_ <- rr_

    }

    R <- -sum(rr)

    for (j in 2:Nc) {

    C <- cc

    }

    C <- -sum(cc)

    Priors

    sigma2 ~ dunif(0.00001, 4)

    for (i in 1:(Nr-1)) {

    rr_ ~ dnorm(0, 0.001)

    }

    for (j in 1:(Nc-1)) {

    cc ~ dnorm(0, 0.001)

    }

    for (i in 1:Np) {

    A_ ~ dnorm(0, 0.001)

    }

    rho ~ dunif(0, 10)

    theta ~ dunif(0, 10)

    }

    ------------ runmod.bugs

    model in model.bugs

    data in data.R

    compile, nchains(1)

    parameters in parameters.R

    parameters to tmp.R

    initialize

    update 2000

    monitor rr

    monitor cc

    monitor rho

    monitor theta

    monitor sigma2

    update 10000

    coda *, stem(JAGS)

    exit

    ------------------ makedata.R

    This script is one of a series that constructs simulated data sets for

    test MCMC methods for estimating dose-response curves for data with

    correlated errors.

    All the action takes place on Nr X Nc microtiter plates. I assume there

    are substantial row and column effects that are constant across plates,

    and that residuals are correlated row-wise and column-wise. Finally,

    plate means vary among plates.

    This simulation only deals with control plates.

    There is a plate mean, a_

    _, p %in% 1:Np, and row effects (R_) and

    column effects (C_)

    Np <- 25 ## Number of plates

    Nc <- 3 ## number of columns on a plate

    Nr <- 2 ## Number of rows on a plate

    Ncells <- Nr * Nc ## Total number of cells on a plate

    A <- rnorm(Np, 0,1)

    rr <- rnorm(Nr - 1, 0, 1)

    R <- c(-sum(rr), rr)

    cc <- rnorm(Nc - 1, 0, 1)

    C <- c(-sum(cc), cc)

    mu = a___

    ___ + R_ + C

    mu <- outer(outer(R,C,"+"), A, "+")

    So, mu is a stack of Nr X Nc matrices. Considered as a vector

    (as.vector(mu)), it is a vector with the Ncells = Nr X Nc wells

    of plate 1 first,

    then the wells of the next plate, etc. Within plates, the wells

    go down rows first: ((1,1), (2,1), (3,1),...,(Nr,1),(1,2),...(Nr,Nc))

    The residuals need to be simulated plate-wise, since wells on the same

    plate are correlated. It will facilitate indexing, if I construct

    two vectors, Ix and Jx, each elements, that map the scaler 1:Ncells

    into the corresponding two dimensional index [Ix_, Jx_]. This is used

    in computing the correlation between elements (i,j) and (i', j'), for

    instance.

    Ix <- as.vector((1:Nr) %o% rep(1, Nc))

    Jx <- as.vector(rep(1, Nr) %o% (1:Nc))

    The parameters for the error distribution

    sigma2 <- 0.045

    rho <- -log(0.5)

    theta <- 1

    Covariance matrix

    Sigma <- sigma2 *

    outer(1:Ncells, 1:Ncells,

    function (i,j) {

    ifelse (i == j,

    1,

    exp(-rhosqrt((Ix_ - Ix)^2 + theta(Jx_-Jx)^2))

    )

    })

    construct a vector of Np replicates from a multivariate normal distribution

    with mean 0 and Sigma given above, then add it to as.vector(mu) to get

    the simulated dataset.

    library(MASS)

    y <- matrix(as.vector(mu) + as.vector(t(mvrnorm(Np, rep(0, Ncells), Sigma))),

    ncol=Np)

    dump(c("y", "Ix", "Jx", "Np","Nr","Nc","Ncells"), file="data.R")

    Save the simulation parameters

    save(A, R, C, rr,cc, sigma2, rho, theta, Np, Nr, Nc, file="simparms.dmp")

    dump(c("A","rr","cc","sigma2","rho","theta"), file="parameters.R")

    if (FALSE) {

    library(coda)

    cn <- read.coda("JAGSchain1.txt","JAGSindex.txt")

    plot(cn)

    heidel.diag(cn)

    effectiveSize(cn)

    }


     
  • Martyn Plummer

    Martyn Plummer - 2010-09-23

    This is known bug that strikes when you invert a large matrix. Omega is not
    quite symmetric and will fail a symmetry test when used as a precision matrix
    for the multivariate normal. The patched version of JAGS (currently cvs only)
    does not do this as 1) the inverse function now uses a Cholesky decomposition
    to do the inverse, preserving symmetry and 2) the symmetry test has been
    removed as it is computationally expensive and fairly redundant. So this will
    be fixed in the next release.

     
  • Woodrow Setzer

    Woodrow Setzer - 2010-09-23

    Thanks. I was taken in by how abruptly the effect happens as the dimension
    increases; 2 X 3 OK, 3 X 3, not. But, if I replace the code

    Omega <- inverse(Sigma)

    with

    O1 <- inverse(Sigma)

    Omega <- 0.5*(O1 + t(O1))

    the error is not triggered, and, for the small example I tried, I get
    reasonable answers. I look forward to the next release!

    Woody

     
  • Woodrow Setzer

    Woodrow Setzer - 2010-09-23

    I'd like to follow up on this. I see that inverse() already has the ability to
    use the Cholesky decomposition for symmetric pd matrices, as long as the spd
    option is set to true. So, I modified the code for Inverse in version 2.1.0 to
    do that, and commented out the test for symmetry in DMNorm.cc. This seems to
    give reasonable answers for my toy problem, and is a bit faster than the
    symmetrization trick and the original code. Do you see a problem with these
    changes? I see from the CVS repository that you've taken a somewhat different
    direction.

     
  • Martyn Plummer

    Martyn Plummer - 2010-09-24

    That's fine.

    In fact WinBUGS/OpenBUGS inverse() function is only defined for SPD matrices.
    I was trying to be more general in JAGS but I did not realise that there would
    be a consequent loss of precision.

     
  • tds151

    tds151 - 2013-12-12

    This bug still seems to exist in JAGS 3.4. I'm running a similar Gaussian process distribution (under a squared exponential covariance) for a parameter set "bb". The covariance matrix is 15 x 15 and I receive an invalid parent value error that appears to relate to lack of SPD. The "Omega" matrix contains squared euclidean distances (with 0's on the diagonal). I tried adding diagonal jitter, but that didn't fix it. Are there any suggestions for making this model work? Please see model, below:

    model
    {
    ################################
    ## Build likelihood for d[i,j], the QCEW - CES count difference for establishment i
    ## at time j
    ################################
    ## likelihood for employment differences from QCEW - CES
    for( i in 1:N ) ## N establishments at T time points / months
    {
    for(j in 1:T)
    {
    d_imp[i,j] ~ dnorm(bb[i,j],tau_d)
    } ## end loop j over T months

                ## draw cluster indicators
        H_g[i]              ~ dcat( beta )
        for(m in 1:M){SC[i,m]       <- equals(m,H_g[i])}
    
        ## GP(DP) prior on trajectories, bb[i,1:T]
        bb[i,1:T]   ~ dmnorm( zeros.T, inverse((1/ksi_star[H_g[i]])*mexp(-        
                                           (1/rho_star[H_g[i]]) * Omega)) )
    } ## end loop i over N establishments
    
    ################################
        ## DP prior on (ksi[s],rho[s]) to borrow strength across error codes
    #################################
    conc_global         ~ dunif(0.5,10)
    
        ## global weights
    beta            ~ ddirich( (conc_global / M) * ones_M )
    
        ## Compute total number of clusters, Mstar.  CL[m] = 1 if n_{m}-1 > 0
        Mstar <- sum(CL[]) 
        for(m in 1:M) 
        {
        n[m]            <- sum(SC[,m])
            CL[m]           <- step(n[m]-1)
        ## prior for locations
        rho_star[m]     ~ dgamma(1.0,1.0)
        ksi_star[m]     ~ dgamma(1.0,1.0)
    } ## end loop m over clusters
    
    ################################
        ## 3rd level priors on covariance matrices
    #################################
    ## tau_d, precision parameter for likelihood
    tau_d           ~ dgamma(1.5,1.5)
    

    } ## end model diff_counts_ss.bugs

     
  • tds151

    tds151 - 2013-12-13

    I've figured out that the problem is the lower numerical precision of mexp() as compared to exp(). In the above script I generate a squared exponential Gaussian process covariance function (of 2 parameters) directly as a matrix, using mexp(). When i re-cast the script to build the matrix, one element-a-time, using exp(), the resulting matrix is able to be inverted and the script is successfully running (though it is still initializing). Is it possible to improve the numerical precision of mexp() to be equivalent to that of exp()?

     

Log in to post a comment.