 Email discussion lists for the UK Education and Research communities  ## BUGS@JISCMAIL.AC.UK

#### View:

 Message: [ First | Previous | Next | Last ] By Topic: [ First | Previous | Next | Last ] By Author: [ First | Previous | Next | Last ] Font: Proportional Font  LISTSERV Archives  BUGS Home  BUGS 2006

#### Options  Subscribe or Unsubscribe   Log In   Get Password Subject: problem simulating from censored gamma in rounded aggregate data

From:  Ramon I. Garcia

Date: Fri, 7 Apr 2006 21:18:54 +0100

Content-Type: text/plain

Parts/Attachments:  text/plain (125 lines)
 ```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 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 ```