MCMCglmm and left-censored data

Mon, Jun 3, 2013 7:52 AM

Hi,

If you add:

data$ymax<-ifelse(data$y <= log(0.4), log(0.4), data$y)

and then fit cbind(ymin, ymax) as the response you should get the  
answer you expect. In your simulation ymax is the actual value of y,  
whereas in my simulation actual values less than log(0.4) are missing  
but known to be less than log(0.4). Values greater than log(0.4) are  
uncensored. I think this is the usual definition of left-censoring (?).

Cheers,

Jarrod



Quoting Anders Bjerring Strathe <strathe at sund.ku.dk> on Mon, 3 Jun  
2013 08:25:55 +0000:

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]]

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

MCMCglmm and left-censored data

Thread (2 messages)