[R-sig-ME] Error in phylogenetic ordinal model with MCMCglmm()

[roee maor]

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