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:
y=mu+p+g+e,
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.
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.
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.
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.