Dear all,

I am trying to fit a linear mixed model for left censored Gaussian data that originates from a pedigreed population. I have used MCMCglmm in the past and I want to use it for this problem as well. I have simulated some data in order to understand how MCMCglmm syntax should be specified.

The tutorial states that for left censored data use -Inf in the 1st column and for right censored data use Inf in the second column. If a particular data point is not censored have the same value in both columns.

The R-code for my simple simulation is given below. I have implemented the model in OpenBUGS and it works well by returning estimates that correspond well with the true. I could fit the animal model with OpenBUGS, but I have >10.000 records and it is easy (to compared to BUGS) to fit an animal model with MCMCglmm, thanks to Jarrod Hadﬁeld… I get biased results with MCMCglmm and hence I think that the syntax is not correct and I hope that some of you could show me how it is done…

Thanks a lot in advance,
Anders Strathe

# R-code
sim <- function(N=2, M=100, cv.beta=0.75, cv.err=0.15){
  beta <- rep(rnorm(M, mean=0, sd=cv.beta), each=N)
  err  <- rnorm(M*N, mean=0, sd=cv.err)
  y <- log(0.8) + beta + err
  data.frame(grp=rep(1:M, each=N), y)
}


data <- sim(M=1000)
data$ymin <- ifelse(data$y <= log(0.4), -Inf, data$y)

p.var <- var(data$y, na.rm=TRUE)
prior <- list(G=list(G1=list(V=matrix(p.var/2),n=1)), R=list(V=matrix(p.var/2),n=1))

m1 <- MCMCglmm(cbind(ymin, y) ~ 1, random = ~ grp, family="cengaussian", data=data,
               prior=prior, nitt=20000, thin=1, burnin=10000,
               start= list(G=runif(1, 0, 100), R=runif(1, 0, 100)), verbose=FALSE)

summary(sqrt(m1$VCV))[[1]]

# TRUE values: cv.beta=0.75 (grp);  cv.err=0.15 (units)

# OpenBUGS model:
model
{
    for (i in 1:N) {
        y[i] ~ dnorm(beta[id[i]], tau)C( , LQD[i])
    }
    for (j in 1:M) {
        beta[j] ~ dnorm(alpha, tau.id)
    }
    tau ~ dgamma(0.001, 0.001)
    sigma <- 1/sqrt(tau)
    tau.id ~ dgamma(0.001, 0.001)
    sigma.id <- 1/sqrt(tau.id)
    alpha ~ dnorm(0.00000E+00, 1.00000E-06)
    mu <- exp(alpha)
}






	[[alternative HTML version deleted]]

