[R-sig-ME] MCMCglmmRAM error

Sara Calhim @@c@|h|m @end|ng |rom gm@||@com
Sat May 23 23:28:32 CEST 2020 

Hi Ben

Thanks for the info and suggestions.

My Bayesian background is using JAGS and therefore the notation of MCMCglmm still confuses me. My second question was about the random argument in MCMCglmm, where both random = ~us(trait):animal and random = ~cors(trait):animal are options in the code provided in Hadfield (2015) Meth Ecol Evol – where the the reduced phylogeny option for MCMCglmm is presented.

I’ll try my luck in the r-sig-phylo list. And also force dichotomous branching to my phylogenetic tree to check if the error disappears.

Cheers,
Sara



> On 23 May 2020, at 00:17, Ben Bolker <bbolker using gmail.com> wrote:
> 
>    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
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

More information about the R-sig-mixed-models mailing list