Dear list members
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 MCMglmm
(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-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
