Ordinal mixed model with variance-covariance matrix

Dear Friends,
I am planning to fit an ordinal mixed model, in a genome-wide association study.
I have measurements (y) on a sample of dogs (n=116), and a matrix of genotypes as the covariates (n*p), p~ 11 million. My y vector has coded as: Control(0), Low risk (1), Medium risk (2), and High risk (4), depends on the severity of a disease. I am planning to run a model as follow:

mu = is the overall mean
p = is a single covariate (we want to obtain p-values or LRT for this part)
g = is random part g~N(0,var(g)*G)
e = residual

Here, y is a vector of ordinal values (0,1,2,3,4), p is a fixed covariate, I want to write a loop over all p covariates, then they can enter the model iteratively and obtain p-values for each covariate. The g, is a random part and has a normal distribution g~N(0,var(g)G), G is a n*n covariance matrix between all 116 dogs ( we call it “genomic relationship matrix”), always is a positive definite square matrix. I would appreciate it if you could help me to fit this model and include G matrix as a covariance matrix in the random part.

Kind regards

1 Like

Hi sorry for note getting to you earlier,
I think such a model would be most easily expressed using the brms package, which supports both ordinal outcomes (using the cumulative family, see https://journals.sagepub.com/doi/10.1177/2515245918823199 for more background) and fixed correlation matrices (Estimating Phylogenetic Multilevel Models with brms • brms).

I’ll however note that with such a big dataset, fitting with brms (or any MCMC method) might be computationally very demanding. You are however likely to also be able to fit such a model with R-INLA, which doesn’t use MCMC and uses an approximation instead and is thus quite a bit faster, while in high-data regime the approximation tends to be very good.

p-values and likelihood-ratio tests don’t have a direct counterpart in the Bayesian paradigm, so I am not sure you will obtain what you need :-) For my current best thinking on how to answer questions about model/hypothesis choice see Hypothesis testing, model selection, model comparison - some thoughts which lists some ways to express such questions when fitting Bayesian models.

Best of luck with your model!

1 Like

Thanks Martin , for the great navigation! do you know is there any R package in the Frequentist paradigm to help me in this case?

I don’t really follow the frequentist packages, so can’t help much. I think support for both random effects AND custom covariance matrix AND ordinal models would be a bit hard to find.

You might want to give brms a try first and maybe it will be fast enough for your purposes. If not, I think inlabru (which is the more recent interface to INLA) coud have you covered - it supports ordinal models and its generic covariance structures should let you specify the phylogenetic matrix (although inlabru works with precision matrix which could change the model a bit).

The ‘ordinal’ package will let you run ordinal models with 1 or 2 random effects in a frequentist framework. But, convergence can be a challenge if there is complexity to the design/data.

Thank you so much Martin, I will follow the R-INLA, hope that can help! I tried the brms, which was the same as mcmcglmm, slow in computation!

Thanks JLC
Unfortunately seems that this package can not handle the variance-covariance matrix between samples for random effects.


Indeed. I found frequentist approaches don’t handle these types of models well. As Martin suggested, implementing your model in ‘brms’ or directly in Stan is likely worthwhile.

1 Like

Thanks, JLC and Martin for all your helps and advises !!

1 Like