[R-sig-ME] MCMCglmmRAM error

[Ben Bolker]

    I didn't see if this got answered or not.  I don't know the answer 
(sorry): if no-one comes forward here, you might also try the 
r-sig-phylo using r-project.org list ...  Polytomies seem like a good guess - 
shouldn't be too hard to split them with a little bit of random noise 
and see if the problem persists ... ??


    I don't know how cors() works, but in general us() will estimate 
separate variances (<>1) for each component, as well as all the pairwise 
correlations, while cor*() variants will assume the variances are 1 (or 
something fixed): from ?MCMCglmm.

 > ‘corg’  fixes the variances along the diagonal to one and ‘corgh’ 
fixes the variances along the diagonal to those specified in the prior.  
‘cors’ allows correlation submatrices.


    Do you mean cor(trait) or cors(trait)?  The difference between us() 
and cor() is that cor() restricts the variances to 1 while us() 
estimates separate variances for each component ...

On 5/22/20 7:55 AM, Sara Calhim wrote:
> Dear Jarrod (and list)
>
> I am trying to run generalized mixed effects models where
> (i) the dependent variable is a non-ordinal 4-level category
> (ii) I have to control for phylogeny using a rooted, non ultrametric phylo object with n = 32 species and 26 internal nodes; with branch length info; I remove the node labels  that denote bayesian consensus support because of the earlier error "phylogeny tip/node labels are not unique” while running MCMcglmm()
> (iii) fixed effects structure is a single categorical predictor (with 2, 3 or 4 levels depending on the model)
> (iv) I have complete separation between Y categ and some of the X1 categories; there is also a strong 'clumping' of the Y_categ on the phylogenetic tree
>
>> str(dummydat)
> 'data.frame':	32 obs. of  8 variables:
>   $ animal : num  1 2 3 4 5 6 7 8 9 10 ...
>   $ X1     : Factor w/ 4 levels "Habitat1","Habitat2",..: 4 3 4 4 4 3 4 4 3 4 ...
>   $ X2     : Factor w/ 2 levels "Mode1","Mode2": 1 2 1 1 1 1 1 1 1 1 ...
>   $ Y_categ: Factor w/ 4 levels "A","B","C","D": 4 4 4 4 4 4 4 2 2 4 ...
>   $ Y_A    : num  0 0 0 0 0 0 0 0 0 0 ...
>   $ Y_B    : num  0 0 0 0 0 0 0 1 1 0 ...
>   $ Y_C    : num  0 0 0 0 0 0 0 0 0 0 ...
>   $ Y_D    : num  1 1 1 1 1 1 1 0 0 1 ...
>
>> table(dummydat$Y_categ, dummydat$X1)
>      Habitat1 Habitat2 Habitat3 Habitat4
>    A        0        3        1        2
>    B        0        0        2        1
>    C        5        0        0        6
>    D        1        0        3        8
>
>> table(dummydat$Y_categ, dummydat$X2)
>     
>      Mode1 Mode2
>    A     4     2
>    B     2     1
>    C    10     1
>    D    11     1
>
> I have been trying to apply Hadfield's (2015) reduced phylogenetic model (MCMCglmm_RAM2) to fix phylogenetic heritability at 1 but I get the error message that
> "non-reduced nodes do not appear first” in my models mod3 and mod4 below.
>
> I have three questions:
>
> (1) What does this error mean? Are the priors correctly set? Is it due to having polytomies in my phylogeny?
>
> (2) What is the difference between models with random = ~us(trait):animal and random = ~cors(trait):animal?
>
> (3) Alternative models using family = categorical and Y_categ as the dependent variable (mod1 and mod2) run. But even long runs (>1e7 iterations) with considerable thinning failed to give well mixed  posterior chains for X1$Habitat1 and X1$Habitat2 levels of the predictor; also effective sample sizes are tiny. This is why I was trying to used the reduced = T option, with cbind() object as the multinomial dependent variable and family = threshold.
>
>
> Thanks!
> Sara
>
>
> Code:
>
>
>
> # # # model with Y_categ
>
> # prior_cat
>
> k <- length(levels(dummydat$Y_categ))
> I <- diag(k-1)
> J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
>
> prior_cat <- list(R = list(fix=1, V=(1/k) * (I + J), n = k - 1),
>                 G = list(G1 = list(V = diag(k - 1), n = k - 1, fix=2)))
>
> # random = ~us
>
> mod1 <- MCMCglmm(Y_categ ~ -1 + X1, random = ~us(trait):animal,  rcov = ~us(trait):units, data = dummydat, pedigree = dummyphylo, scale = FALSE, family = "categorical", prior = prior_cat, nitt = 1e5)
>
> # random = ~cors
>
> mod2 <- MCMCglmm(Y_categ ~ -1 + X1, random = ~cors(trait):animal,  rcov = ~us(trait):units, data = dummydat, pedigree = dummyphylo, scale = FALSE, family = "categorical", prior = prior_cat)
>
>
> # # # model with cbind(Y_A, Y_B, Y_C, Y_D)
>
> # prior_cbind
>
> prior_cbind <- list(R=list(V=diag(4)*1e-15, fix=1), G=list(G1=list(V=diag(4), nu=0.002, fix=2)))
>
> # random = ~us
>   
> mod3 <- MCMCglmm(cbind(Y_A, Y_B, Y_C, Y_D) ~ -1 + X1, random = ~us(trait):animal,  rcov = ~us(trait):units, data = dummydat, pedigree = dummyphylo, scale = FALSE, family = c(rep( "threshold", 4)), reduced = TRUE, prior = prior_cbind)
>
> # random = ~cors
>
> mod4 <- MCMCglmm(cbind(Y_A, Y_B, Y_C, Y_D) ~ -1 + X1, random = ~cors(trait):animal,  rcov = ~us(trait):units, data = dummydat, pedigree = dummyphylo, scale = FALSE, family = c(rep( "threshold", 4)), reduced = TRUE, prior = prior_cbind)
>
>
>
>
>
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models