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

Discussion:

roee maor

2018-08-02 13:24:55 UTC

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".

packageVersion("MCMCglmm")

[1] ‘2.25’

R.version.string

[1] "R version 3.4.3 (2017-11-30)"

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

Dexter Locke

2018-08-02 13:39:18 UTC

Permalink

Maybe Response needs to be an ordered factor, not a factor with three levels. Try something like

valid$Response <- as.factor(valid$Response, ordered=T)

See also th clmm package.

HTH, Dexter

Post by 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".

packageVersion("MCMCglmm")

[1] ‘2.25’

R.version.string

[1] "R version 3.4.3 (2017-11-30)"

## 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,

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)
contrasts can be applied only to factors with 2 or more levels
## the shape of the data

str(valid)

$ 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
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Paul Buerkner

2018-08-02 13:45:24 UTC

Permalink

If it doesn't end up working in MCMCglmm, you may also try the brms
package.

See
https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html
for an introduction of
phylogenetic models in brms.

Post by Dexter Locke
Maybe Response needs to be an ordered factor, not a factor with three
levels. Try something like
valid$Response <- as.factor(valid$Response, ordered=T)
See also th clmm package.
HTH, Dexter

packageVersion("MCMCglmm")

[1] ‘2.25’

R.version.string

[1] "R version 3.4.3 (2017-11-30)"

INphylo <- inverseA(mammaltree, nodes="ALL", scale=TRUE) ## phylogeny

with 1399 tips, setting nodes="TIPS" is extremely slow

Post by roee maor

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,

Post by roee maor
+ random = ~ us(trait):Binomial,
+ rcov = ~ us(trait):units,
+ prior = list(R = list(fix=1, V=(1/k) * (I + J), n =

k-1),

Post by roee maor
+ G = list(G1 = list(V = diag(k-1), n =

k-1))),

Post by roee maor
+ 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,

Post by roee maor
+ 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))),

Post by roee maor
+ ginverse = list(Binomial=INphylo$Ainv),
+ burnin = 300000,
+ nitt = 3000000,
+ thin = 2000,
+ family = "ordinal",
+ data = valid,
+ pl = TRUE)
contrasts can be applied only to factors with 2 or more levels
## the shape of the data

str(valid)

$ 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 ...

Post by roee maor
$ 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 ...

Post by roee maor
$ 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
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

roee maor

2018-08-03 11:43:28 UTC

Permalink

Thanks everyone for your quick and helpful responses.
Removing the 'trait' terms indeed solved this problem and I have taken
Jarrod's advice to use family="threshold" instead of "ordinal".

Would it be possible to get a quick summary of the main differences
between ordinal and threshold models? Both seem to be rarely used and
I can't find much documentation on their inner workings in the
MCMCglmm vignette, tutorial or course notes.

Cheers,

If it doesn't end up working in MCMCglmm, you may also try the brms package.
See https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html for an introduction of
phylogenetic models in brms.

Post by Dexter Locke
Maybe Response needs to be an ordered factor, not a factor with three levels. Try something like
valid$Response <- as.factor(valid$Response, ordered=T)
See also th clmm package.
HTH, Dexter

packageVersion("MCMCglmm")

[1] ‘2.25’

R.version.string

[1] "R version 3.4.3 (2017-11-30)"

## 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,

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)
contrasts can be applied only to factors with 2 or more levels
## the shape of the data

str(valid)

_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

HADFIELD Jarrod

2018-08-03 11:58:48 UTC

Permalink

Hi,

With three levels you have 4 cut-points, one of which needs estimating (call it cp.est)

cp<-c(-Inf, 0, cp.est, Inf)

If eta is the fixed+random effect linear predictor, Xb+Zu, then the probability of being in category i is

pnorm(cp[i+1], eta , sqrt(Vr)) - pnorm(cp[i], eta, sqrt(Vr))

when family=âthresholdâ. Vr is the residual (units) variance. When Vr=1 this is standard profit regression. When family=âordinalâ it is:

pnorm(cp[i+1], eta , sqrt(Vr+1)) - pnorm(cp[i], eta, sqrt(Vr+1))

Its not a big deal because the eta cam be rescaled to be equivalent:

eta_ordinal/sqrt(Vr+1) = eta_threshold/sqrt(Vr)

and Vr is fixed.

However, the MCMC algorithm for family=âthresholdâ is more efficient.

Cheers,

Jarrod

On 3 Aug 2018, at 12:43, roee maor <***@gmail.com<mailto:***@gmail.com>> wrote:

Thanks everyone for your quick and helpful responses.
Removing the 'trait' terms indeed solved this problem and I have taken
Jarrod's advice to use family="threshold" instead of "ordinal".

Would it be possible to get a quick summary of the main differences
between ordinal and threshold models? Both seem to be rarely used and
I can't find much documentation on their inner workings in the
MCMCglmm vignette, tutorial or course notes.

Cheers,
RoiOn Thu, 2 Aug 2018 at 14:45, Paul Buerkner <***@gmail.com<mailto:***@gmail.com>> wrote:

If it doesn't end up working in MCMCglmm, you may also try the brms package.

See https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html for an introduction of
phylogenetic models in brms.

2018-08-02 15:39 GMT+02:00 Dexter Locke <***@gmail.com<mailto:***@gmail.com>>:

Maybe Response needs to be an ordered factor, not a factor with three levels. Try something like

valid$Response <- as.factor(valid$Response, ordered=T)

See also th clmm package.

HTH, Dexter

On Aug 2, 2018, at 9:24 AM, roee maor <***@gmail.com<mailto:***@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-***@r-project.org<mailto:R-sig-mixed-***@r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-***@r-project.org<mailto:R-sig-mixed-***@r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-***@r-project.org<mailto:R-sig-mixed-***@r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

roee maor

2018-08-03 12:08:22 UTC

Permalink

That's great Jarrod, thanks very much!
Roi

Post by HADFIELD Jarrod
Hi,
With three levels you have 4 cut-points, one of which needs estimating (call it cp.est)
cp<-c(-Inf, 0, cp.est, Inf)
If eta is the fixed+random effect linear predictor, Xb+Zu, then the probability of being in category i is
pnorm(cp[i+1], eta , sqrt(Vr)) - pnorm(cp[i], eta, sqrt(Vr))
pnorm(cp[i+1], eta , sqrt(Vr+1)) - pnorm(cp[i], eta, sqrt(Vr+1))
eta_ordinal/sqrt(Vr+1) = eta_threshold/sqrt(Vr)
and Vr is fixed.
However, the MCMC algorithm for family=“threshold” is more efficient.
Cheers,
Jarrod
Thanks everyone for your quick and helpful responses.
Removing the 'trait' terms indeed solved this problem and I have taken
Jarrod's advice to use family="threshold" instead of "ordinal".
Would it be possible to get a quick summary of the main differences
between ordinal and threshold models? Both seem to be rarely used and
I can't find much documentation on their inner workings in the
MCMCglmm vignette, tutorial or course notes.
Cheers,
If it doesn't end up working in MCMCglmm, you may also try the brms package.
See https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html for an introduction of
phylogenetic models in brms.
Maybe Response needs to be an ordered factor, not a factor with three levels. Try something like
valid$Response <- as.factor(valid$Response, ordered=T)
See also th clmm package.
HTH, Dexter
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".
packageVersion("MCMCglmm")
[1] ‘2.25’
R.version.string
[1] "R version 3.4.3 (2017-11-30)"
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)
contrasts can be applied only to factors with 2 or more levels
## the shape of the data
str(valid)
$ 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
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________
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.

Pierre de Villemereuil

2018-08-02 13:46:10 UTC

Permalink

Hi,

Contrary to "categorical", the "ordinal" family is not multivariate on the latent scale, so you need to drop everything related to "trait" in the ordinal model and modify the prior to be univariate as well.

I think MCMCglmm is complaining here because "trait" does not have more than one level due to this.

Hope this helps,
Pierre.

packageVersion("MCMCglmm")

[1] ‘2.25’

R.version.string

[1] "R version 3.4.3 (2017-11-30)"

## 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,

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)
contrasts can be applied only to factors with 2 or more levels
## the shape of the data

str(valid)

HADFIELD Jarrod

2018-08-02 13:44:33 UTC

Permalink

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

packageVersion("MCMCglmm")

[1] ‘2.25’

R.version.string

[1] "R version 3.4.3 (2017-11-30)"

## 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,

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)
contrasts can be applied only to factors with 2 or more levels
## the shape of the data

str(valid)

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.