[R-sig-ME] MCMCglmmRAM error
Sara Calhim
@@c@|h|m @end|ng |rom gm@||@com
Fri May 22 13:55:09 CEST 2020
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]]
More information about the R-sig-mixed-models
mailing list