Hello Bugs users,
I am analyzing an aggregated dataset with severe rounding and I am
interesting in determining how the rounding of the sample variances may
affect parameters estimates. My data consists of sample means and variances
of reponses from various dose groups. My model for the expected value of the
sample means is a nonlinear function of dose and parameters. My model is
# data : x, s2, ngroups, n, d, c, R
# param : lA,lD,gam,P
# P = c(.5)
# lA = c(1)
# lD = c(2)
# gam = c(1)
# reals2 = s2
# sigma2 = c(1,1,1,1,1)
#
for (i in 1:ngroups) {
g[i] <- A*(1 - (1-P)*(1-exp(log((1-R-P)/(1-P))*pow(d[i]/D,gam))))
prec[i] <- 1/(sigma2[i]/n[i])
x[i] ~ dnorm(g[i],prec[i])
alpha[i] <- (n[i]-1)/2
beta[i] <- (n[i]-1)/(2*sigma2[i])
lower[i] <- s2[i] - c
upper[i] <- s2[i] + c
reals2[i]~ dgamma(alpha[i],beta[i])I(lower[i],upper[i])
sigma2[i] ~ dgamma(1.0,1.0)
}
P ~ dunif(0.001, 0.89)
gam ~ dunif(0.001, 28)
D <- exp(lD)
A <- exp(lA)
lA ~ dnorm(0,1.0E-6)
lD ~ dnorm(0,1.0E-6)
Unfortunately, I get the following error,
cannot sample from interval censored gamma full conditional reals2[1]
Does anyone know why this occurs or trick around this.
I am aware of the 'classical rounding' and 'berckson rounding'
solutions to accomodate rounding suggested in the manual but they don't seem
appropriate for my problem.
Any suggestions would be great.
Thanks,
Ramon
Phd Student UNC-Biostatistics
National Center for Computational Toxicology
US Environmental Protection Agency
This is the R code I am using to call winbugs from
R through fitbugs
-----------------------------------------------------------------------
source("fitbugs.R")
library(R2WinBUGS)
x <- c(13.8, 14.4, 13.5, 13.3, 5.2)# means
n <- c(6, 6, 6, 6, 6) # sample sizes
d <- c(0, 0.25, 0.54, 1, 1.6) # dose values
s2 <- c(0.5, 0.3, 0.6, 0.5, 2.4) # rounded sample variances
c <- .05 # precision value by which number is
# rounded i.e. .5 if to nearest whole
# number
ngroups <- length(x)
R <- 0.1 # fixed value
roundexpbugs <- function() {
# data : x, s2, ngroups, n, d, c, R
# param : lA,lD,gam,P
# P = c(.5)
# lA = c(1)
# lD = c(2)
# gam = c(1)
# reals2 = s2
# sigma2 = c(1,1,1,1,1)
#
for (i in 1:ngroups) {
g[i] <- A*(1 - (1-P)*(1-exp(log((1-R-P)/(1-P))*pow(d[i]/D,gam))))
prec[i] <- 1/(sigma2[i]/n[i])
x[i] ~ dnorm(g[i],prec[i])
alpha[i] <- (n[i]-1)/2
beta[i] <- (n[i]-1)/(2*sigma2[i])
lower[i] <- s2[i] - c
upper[i] <- s2[i] + c
reals2[i]~ dgamma(alpha[i],beta[i])%_%I(lower[i],upper[i])
sigma2[i] ~ dgamma(1.0,1.0)
}
P ~ dunif(0.001, 0.89)
gam ~ dunif(0.001, 28)
D <- exp(lD)
A <- exp(lA)
lA ~ dnorm(0,1.0E-6)
lD ~ dnorm(0,1.0E-6)
}
n.burn <- 1e4 # number of iterations burned
n.samp <- 1e4 # sample size
n.iter <- n.burn + n.samp
n.thin <- 1 # thinning
n.chain <- 1 # number of chains
bugsout <- fitbugs(roundexpbugs,n.chain=n.chain,n.iter=n.iter,
n.thin=n.thin,n.burnin=n.burn,debug=TRUE)
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
