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(-rho*sqrt((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)

}