[R-sig-ME] Error in phylogenetic ordinal model with MCMCglmm()
HADFIELD Jarrod
j@h@dfield @ending from ed@@c@uk
Thu Aug 2 15:44:33 CEST 2018
Hi,
There is only one ’trait’ in an ordinal model, so the terms with trait in are redundant and causing the problem. You should get rid of all trait and us(trait): terms. With family=“categorical” there are two traits with 3 categories, one for each contrast. Note that is using family=“threshold” will be more efficient than using “ordinal” although they are almost the same model.
Cheers,
Jarrod
> On 2 Aug 2018, at 14:24, roee maor <roeemaor using gmail.com> wrote:
>
> Dear list,
>
> I'm using MCMCglmm to run a phylogenetic model where the response is a
> 3-level ordinal factor (i.e. level 2 is an intermediate phenotype
> between 1 and 3), and the predictors include one factorial (foraging
> habitat), one ordinal (trophic level), and several continuous
> variables.
>
> As far as I know MCMCglmm is the only package that can handle logistic
> models for phylogenetically structured multi-level discrete data, but
> please correct me if that's not the case.
>
> My problem right now is that I can't get MCMCglmm() to work with the
> 'family' argument set to "ordinal", although it does work with
> "categorical".
> Here's the code I'm using:
>
>> packageVersion("MCMCglmm")
> [1] ‘2.25’
>> R.version.string
> [1] "R version 3.4.3 (2017-11-30)"
>
> ## model specifications:
>> INphylo <- inverseA(mammaltree, nodes="ALL", scale=TRUE) ## phylogeny with 1399 tips, setting nodes="TIPS" is extremely slow
>> k <- length(levels(valid$Response))
>> I <- diag(k-1)
>> J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
>
> ## categorical model (unordered response) - runs to completion
>> m1 <- MCMCglmm(Response ~ -1 + trait + ForagingHab + Troph_Lev + Mass + Mean.Diur.Range + Max.Temp.Warmest.M + Temp.Annual.Range + Precip.Driest.Month + PET,
> + random = ~ us(trait):Binomial,
> + rcov = ~ us(trait):units,
> + prior = list(R = list(fix=1, V=(1/k) * (I + J), n = k-1),
> + G = list(G1 = list(V = diag(k-1), n = k-1))),
> + ginverse = list(Binomial=INphylo$Ainv),
> + burnin = 300000,
> + nitt = 3000000,
> + thin = 2000,
> + family = "categorical",
> + data = valid,
> + pl = TRUE)
>
> ## ordinal model and error message
>> m2 <- MCMCglmm(Response ~ -1 + trait + ForagingHab + Troph_Lev + Mass + Mean.Diur.Range + Max.Temp.Warmest.M + Temp.Annual.Range + Precip.Driest.Month + PET,
> + random = ~ us(trait):Binomial,
> + rcov = ~ us(trait):units,
> + prior = list(R = list(fix=1, V=1, n = k-1),
> + G = list(G1 = list(V = diag(k-1), n = k-1))),
> + ginverse = list(Binomial=INphylo$Ainv),
> + burnin = 300000,
> + nitt = 3000000,
> + thin = 2000,
> + family = "ordinal",
> + data = valid,
> + pl = TRUE)
> Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
> contrasts can be applied only to factors with 2 or more levels
>
> ## the shape of the data
>> str(valid)
> 'data.frame': 1399 obs. of 16 variables:
> $ Binomial : chr "Abrocoma_bennettii" "Abrothrix_andinus"
> "Abrothrix_jelskii" "Abrothrix_longipilis" ...
> $ Response : Factor w/ 3 levels "1","2","3": 1 3 2 2 2 3 2 3 1 3 ...
> $ ForagingHab : Factor w/ 7 levels "1","3","4","5",..: 2 2 2 2
> 2 2 2 2 2 2 ...
> $ Troph_Lev : Factor w/ 3 levels "1","2","3": 1 2 2 2 2 3 2 1 1 1 ...
> $ Mass : num 250.5 24.9 34.5 38.9 24.5 ...
> $ Annual.Mean.Temp : num 12.42 7.26 9.18 9.9 8.62 ...
> $ Mean.Diur.Range : num 10.46 13.82 16.16 9.15 7.78 ...
> $ Max.Temp.Warmest.M : num 22 16.6 19.1 19.8 17.2 ...
> $ Min.Temp.Coldest.M : num 3.77 -3.98 -3.24 2.21 1.46 ...
> $ Temp.Annual.Range : num 18.3 20.6 22.4 17.6 15.7 ...
> $ Mean.Temp.Warm.Q : num 16 9.2 11.1 13.8 12.2 ...
> $ Mean.Temp.Cold.Q : num 8.87 4.49 6.39 5.98 4.87 ...
> $ Annual.Precip : num 166 645 558 903 1665 ...
> $ Precip.Driest.Month: num 1.74 5.99 6.31 31.31 104.7 ...
> $ AET : num 213 482 704 455 361 ...
> $ PET : num 1074 1242 1305 677 638 ...
>
>
> I don't understand what factors the error refers to, because there
> sufficient levels in the response even if one is absorbed in the
> intercept.
>
> The R-constraint in the prior is specified as suggested in the
> MCMCglmm tutorial (fix=1, V=1), but the error message is the same
> whether I use this specification or the categorical model
> specification (fix=1, V=(1/k)*(I + J)) .
>
> On a side note - what parameters affect the acceptance rates? The
> categorical models maintain a rate of around 0.3 so I think the mixing
> could be improved.
>
> Any input would be very much appreciated.
>
> Many thanks,
> Roi Maor
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list