Mixed effect modeling using brms

Modeling
1

Please share your Stan program and accompanying data if possible.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

model {
  vector[N] mu = alpha + beta * x;
  y ~ normal(mu, sigma); 
}

instead of

model{
vector[N] mu = alpha+beta*x;
y~normal(mu,sigma);
}

To include mathematical notation in your post put LaTeX syntax between two $ symbols, e.g.,
p(\theta | y) \propto p(\theta) p(y | \theta).

Don’t forget to add relevant tags to your topic (top right of this form) for application area and/or class of models you work with.

I am running mixed effect code using brms with 37 predictor variables as random terms, but i am only able to take 18 only at a time, taking 19th predictor it shows me the following error as :
Chain 2: Exception: model2a1c16b33b75__namespace::write_array: Cor_1 is not positive definite. (in ‘string’, line 142, column 2 to column 66)

Kindly help

2

Welcome to the Stan Discourse. I hope you find it to be a friendly, helpful group.

It’s not immediately obvious to me what the issue is. Could you share some additional details of your problem?

  1. Could you share some data (real or simulated) and code that demonstrates the issue?
  2. Does the issue occur regardless of which 19 predictors you use? Or is it specific to some combination of variables? For example, if you randomly reorder the predictors, does the issue still occur once you add the 19th?
3 
  FCP_R_24 ~Age + Gender+ Carb_24+ Fiber_24+ protein_g+ fat_g+ mono_fat_g+ poly_fat_g+ trans_fat_g+ sat_fat_g+ omega_three_fatty_acids_g+ sodium_mg+ sugars_g+ added_sugars_g+ vitamin_d_iu+ calcium_mg+ potassium_mg+ zinc_mg+ vitamin_a_iu+ vitamin_e_mg+ vitamin_d2d3_ug+ vitamin_c_mg+ thiamin_mg+ riboflavin_mg+ niacin_mg+ vitamin_b6_mg + vitamin_b12_ug+ vitamin_k_ug+ folate_dfe_ug+ cholesterol_mg+ caffeine_mg+ choline_mg+ pantothenic_acid_mg+ phosphorus_mg+ selenium_mcg+(1 +Age + Gender+Carb_24+ Fiber_24+ protein_g+ fat_g+ mono_fat_g+ poly_fat_g+ trans_fat_g+ sat_fat_g+ omega_three_fatty_acids_g+ sodium_mg+sugars_g+ added_sugars_g+ vitamin_d_iu+ calcium_mg + potassium_mg+ zinc_mg + vitamin_a_iu + vitamin_e_mg| Patient_ID),

These are my predictors, on population level, it shows me the estimate of all predictors, but once i try to run their random effects, it shows me the error like
Chain 2: Exception: model2a1c16b33b75__namespace::write_array: Cor_1 is not positive definite. (in ‘string’, line 142, column 2 to column 66)

4

It might be helpful to use brms::make_stancode to take a look at the actual Stan code that is generated by brms for your model, so that you can clearly see what is happening under the hood.

The error message that you’re encountering involves Cor_1. I believe this refers to the correlation matrix of the random effects for the first grouping variable. In your case, you only have one grouping variable, namely Patient_ID. Thus, Cor_1 contains the correlations between each of your random effects (i.e., random intercept and slopes).
This is typically saved in the generated quantities block like so:

corr_matrix[M_1] Cor_1;
Cor_1 = multiply_lower_tri_self_transpose(L_1);

where M_1 is the number of random effects terms (which in your case is 20, with 1 random intercept and 19 random slopes), and L_1 is the estimated Cholesky factor of the correlation matrix.

The error message is telling you that “Cor_1 is not positive definite.”
In principle, if L_1 is a valid Cholesky factor, then multiply_lower_tri_self_transpose should always produce a valid correlation matrix. So the likely root cause of the problem is with the estimation of L_1.

You have a relatively large number of random effects, so presumably it’s difficult to estimate their correlations. It’s possible that some draws of L_1 correspond to degenerate or nearly singular matrices.

Do you have enough data to estimate your random effects structure? Note that the default prior on random effects correlation matrices in brms is lkj_corr_cholesky(1), which means that a priori all correlation matrices are equally likely (see Details in documentation of brms::set_prior).

You could simplify things by modelling one or more of the random effects as uncorrelated, by splitting the random effects terms in your model formula. So in pseudocode, changing (1 + x1 + x2 | group) into this: (1 | group) + (0 + x1 + x2 | group) would remove the correlations between the random intercept and random slopes, but would keep the correlations between the random slopes, and this: (1 | group) + (0 + x1 | group) + (0 + x2 | group) would remove all correlations. But I can’t judge whether that actually makes theoretical sense for your data and model.

Hope that helps.

2 Likes