I am trying to estimate a repeated measures model with individual-level random intercepts and correlated random slopes partitioned into genetic effects G and environmental effects E. The genetic effects are scaled with fixed row-wise relatedness matrix A, which is relatively sparse. The environmental effects are i.i.d. and meant to simply soak up any unobserved individual heterogeneity not captured by the fixed covariance genetic effects. The response model is specified by the sum G + E, which is making it difficult for identification of the distinct standard deviations and correlations of the G and E effects. The model does a somewhat good job of returning simulated SDs despite sampling quite poorly (with low bulk and tail ESS), but the correlations seem to jump around haphazardly.

I’ve experimented with implementing basic prior regularization and some constraints on the parameters, but I’m hoping I’ve missed more effective strategies for enhancing identification. My code is below. Note that I take advantage of a matrix normal parameterization of the genetic effects to avoid direct computation of the Kronecker product G %x% A.

```
data {
int<lower=1> N; //number of observations
int<lower=0> I; //number of individuals
int<lower=0> M; //number of genetic effects
int<lower=0> E; //number of environment effects
int<lower=1> id_j[N]; //index of individual observations
matrix[I,I] Amat; //genetic relatedness matrix
vector[N] X1; //covariate
real Y[N]; //normal response
}
transformed data{
matrix[I,I] LA = cholesky_decompose(Amat); //lower-triangle A matrix
real conditional_mean0 = 0; //fixed global intercept
}
parameters {
real<lower=-1,upper=1> conditional_mean1; //expectation for subject slopes
vector<lower=0>[M] sdG_U; //genetic effect SDs
vector<lower=0>[E] sdE_U; //environment effect SDs
real<lower=0> sd_R; //residual SD
cholesky_factor_corr[M] LG_u; //genetic effect correlations
cholesky_factor_corr[E] LE_u; //environment effect correlations
matrix[I,M] std_devG; //genetic deviations
matrix[I,E] std_devE; //environment deviations
}
transformed parameters {
matrix[I, M] G_re; //genetic random effects (G_intercept, G_slope)
matrix[I, E] E_re; //environment random effects (E_intercept, E_slope)
vector[I] P0; //phenotypic random intercepts (G + E)
vector[I] P1; //phenotypic random slopes (G + E)
//matrix normal sampling of genetic effects
G_re = LA * std_devG * diag_pre_multiply(sdG_U, LG_u) ;
//environment effects
E_re = std_devE * diag_pre_multiply(sdE_U, LE_u) ;
//phenotypic effects (centered on global intercept and slope)
P0 = rep_vector(conditional_mean0, I) + col(G_re,1) + col(E_re,1);
P1 = rep_vector(conditional_mean1, I) + col(G_re,2) + col(E_re,2);
}
model{
conditional_mean1 ~ std_normal();
to_vector(std_devG)~std_normal();
to_vector(std_devE)~std_normal();
LG_u ~ lkj_corr_cholesky(3);
LE_u ~ lkj_corr_cholesky(3);
to_vector(sdG_U) ~ exponential(3);
to_vector(sdE_U) ~ exponential(3);
sd_R ~ exponential(3); //residual sd
//response Y ~ G_int + E_int + (G_slope + E_slope)*X1
Y ~ normal( P0[id_j] + P1[id_j].*X1, sd_R);
}
```

Here are posterior plots indicative of the identification issues for the G and E SDs and correlations.

I simulated G and E correlations of large magnitude and opposite sign (-0.5 and 0.5 respectively).

Thanks for your help!