[R-sig-ME] MCMCglmmRAM error

[Sara Calhim]

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