[R-sig-ME] MCMCglmm Bivariate Binary Model

[Jarrod Hadfield]

Hi Chuck,

The problem you encounter is because the prior specification is not  
proper. Something of the form:

prior = list(R = list(V = diag(2), n = 0, fix = 1),
               G = list(G1 = list(V = diag(2), n = 2)))

would be more appropriate. Note I've changed fix=2 to fix=1 so that  
the whole R structure is fixed as an Identity matrix. I've also  
changed n=1 to n=2 in the G structure because a 2x2 matrix requires a  
minimum of n=2 to be proper.  I use identity matrices in the tutorial  
out of convenience, but better priors should be used if they exist.  
Failing that checking sensitivity to different prior specifications  
would be a good idea.

My guess is that it is the R-structure definition that is creating  
problems because the residual variance for ANYHR was being estimated,  
but there is no information in the data for this parameter (because  
its binary).

I don't recommend bivariate binary models in MCMCglmm yet because you  
can't restrict the residual covariance structure to a correlation  
matrix.   This would fix the two residual variances because they  
cannot be estimated, but would allow the residual correlation between  
the outcomes to be estimated. Fixing this correlation to 0 (using  
idh()) will almost certainly bias the random effect correlation.

One possibility would be to treat your bivariate binary model as a  
3-level multinomial 1/0 0/1 1/1. This would also result in a bivariate  
model (where one of the three categories becomes the base line  
category) but it may be easier to set up the model in a way that  
captures the lack of 0/0 outcome.

Cheers,

Jarrod


Quoting Chuck Cleland <ccleland at optonline.net>:

> Hi All,
>
>   I attempted to fit a bivariate binary model similar to the example in
> the Tutorial:
>
> prior = list(R = list(V = diag(2), n = 0, fix = 2),
>              G = list(G1 = list(V = diag(2), n = 1)))
>
> model5 <- MCMCglmm(cbind(ANYHR, ANYPO) ~ trait - 1, random =
> ~us(trait):PROGRAM, rcov = ~idh(trait):units, family=c("categorical",
> "categorical"), data = radars, prior = prior, verbose = FALSE)
>
> Error in MCMCglmm(cbind(ANYHR, ANYPO) ~ trait - 1, random =
> ~us(trait):PROGRAM,  :
>   ill-conditioned G/R structure: use proper priors if you haven't or
> rescale data if you have
>
>   In this dataset, one of the four cells in the 2x2 table is empty:
>
>    Cell Contents
> |-------------------------|
> |                   Count |
> |-------------------------|
>
> Total Observations in Table:  21335
>
>              | ANYPO
>        ANYHR |        0  |        1  | Row Total |
> -------------|-----------|-----------|-----------|
>            0 |        0  |     8999  |     8999  |
> -------------|-----------|-----------|-----------|
>            1 |     6458  |     5878  |    12336  |
> -------------|-----------|-----------|-----------|
> Column Total |     6458  |    14877  |    21335  |
> -------------|-----------|-----------|-----------|
>
>   Is there a way to change the specification for prior to make this work?
>
> thanks,
>
> Chuck
>
> --
> Chuck Cleland, Ph.D.
> NDRI, Inc. (www.ndri.org)
> 71 West 23rd Street, 8th floor
> New York, NY 10010
> tel: (212) 845-4495 (Tu, Th)
> tel: (732) 512-0171 (M, W, F)
> fax: (917) 438-0894
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>



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