[R-sig-ME] MCMCglmm and left-censored data

[Jarrod Hadfield]

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 Hadfield… 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.