Hi Tim, You specified G to be a diagonal matrix (no correlation between the random intercept and random slope). This is fine, but make sure that it is what you intended. Also, you drew a single sample for the random effects (rather than 200 samples). This means that everyone in the data set has the same intercept and slope, hence the zero variances. Try the code below. I replicated what you were trying to do, but without a compound symmetric structure for R. I left that off because you should generally avoid specifying a compound symmetric structure for R while also including random effects in the model (the random intercept and the residual within-cluster covariance end up being redundant). An AR(1) structure might be preferable if you want both R-side random effects and G-side random effects in the same model. Also, I switched the sampling of the random effects and the error to be performed in blocks. This decreases the memory requirements for the simulation, which will be helpful if you want to test the simulation for a large number of subjects. Best, Steve proc iml; /*************************************************************************** Specify model parameters ***************************************************************************/ n_l2 = 200; *Number of people; n_l1 = 4; *Number of time points per person; n = n_l1 * n_l2; seed = 314159; *Random value seed; call randseed(seed); *Specifies a seed for all random number generation; *Specify rixed intercept and slope; fixed = {16.7611 0.6602}; *Specify block of G matrix; *Treating G as block-diagonal reduces memory requirements; g_block = {1.8430 0, 0 0.02228}; *Specify block of R matrix; *Treating R as block-diagonal reduces memory requirements; *R is currently set to have a block-diagonal structure; r_block = 1.8787 # i(n_l1); /*************************************************************************** Set up matrices and draw random values ***************************************************************************/ *Draw intercepts and slopes from a multivariate normal distribution; l2_random = randnormal(n_l2, fixed, g_block); l1_random = j(n,2,0); l2_id = j(n,1,0); do i = 1 to (n_l2); *Generate cluster ID variable; l2_id[n_l1*(i-1)+1:n_l1*(i-1)+n_l1,1]=i; *Copy random values from level-2 matrix to level-1 matrix; do k = 1 to n_l1; l1_random[n_l1*(i-1)+k,]=l2_random[i,]; end; end; *Generate level-1 ID variable; l1_id = (1:n)`; *Generate age variable; age = repeat(do(8,14,2)`,n_l2,1); *Draw level-1 errors from a multivariate normal distribution; errors = colvec(randnormal(n_l2, j(n_l1,1,0), r_block)); *Generate dependant variable; dv = l1_random[,1] + l1_random[,2] # age + errors; output = l2_id || l1_id || age || dv; names = {l2_id l1_id age dv}; create multilevelData from output [colname = names]; append from output; close multilevelData; quit; proc mixed data=multilevelData; model dv = age / s; random int age / subject=l2_id type=un; run;
... View more