From 27c4b8bd28c1de0605ce7d25ab3197a80ca5c960 Mon Sep 17 00:00:00 2001 From: AngeDakouri <149696288+AngeDakouri@users.noreply.github.com> Date: Mon, 14 Sep 2026 10:46:10 +0200 Subject: [PATCH] Add files via upload Here the stan file related to the SEM model (I add the prior calculations in the section generated quantities) and a pipeline to apply the prior sensitivity. --- BayesSEMsensitivity/BayesSEM.stan | 2270 +++++++++++++++++++++ BayesSEMsensitivity/BayesSEMSensitivity.R | 210 ++ 2 files changed, 2480 insertions(+) create mode 100644 BayesSEMsensitivity/BayesSEM.stan create mode 100644 BayesSEMsensitivity/BayesSEMSensitivity.R diff --git a/BayesSEMsensitivity/BayesSEM.stan b/BayesSEMsensitivity/BayesSEM.stan new file mode 100644 index 0000000..68c641d --- /dev/null +++ b/BayesSEMsensitivity/BayesSEM.stan @@ -0,0 +1,2270 @@ +/* This file is based on LERSIL.stan by Ben Goodrich. + https://github.com/bgoodri/LERSIL */ +functions { // you can use these in R following `rstan::expose_stan_functions("foo.stan")` + // mimics lav_mvnorm_cluster_implied22l(): + matrix calc_W_tilde(matrix sigma_w, vector mu_w, array[] int var1_idx, int p_tilde) { + matrix[p_tilde, p_tilde + 1] out = rep_matrix(0, p_tilde, p_tilde + 1); // first column is mean vector + vector[p_tilde] mu1 = rep_vector(0, p_tilde); + matrix[p_tilde, p_tilde] sig1 = rep_matrix(0, p_tilde, p_tilde); + mu1[var1_idx] = mu_w; + sig1[var1_idx, var1_idx] = sigma_w; + out = append_col(mu1, sig1); + return out; + } + matrix calc_B_tilde(matrix sigma_b, vector mu_b, array[] int var2_idx, int p_tilde) { + matrix[p_tilde, p_tilde + 1] out = rep_matrix(0, p_tilde, p_tilde + 1); + vector[p_tilde] mu2 = rep_vector(0, p_tilde); + matrix[p_tilde, p_tilde] sig2 = rep_matrix(0, p_tilde, p_tilde); + mu2[var2_idx] = mu_b; + sig2[var2_idx, var2_idx] = sigma_b; + out = append_col(mu2, sig2); + return out; + } + vector twolevel_logdens(array[] vector mean_d, array[] matrix cov_d, matrix S_PW, array[] vector YX, array[] int nclus, array[] int clus_size, array[] int clus_sizes, int nclus_sizes, array[] int clus_size_ns, vector impl_Muw, matrix impl_Sigmaw, vector impl_Mub, matrix impl_Sigmab, array[] int ov_idx1, array[] int ov_idx2, array[] int within_idx, array[] int between_idx, array[] int both_idx, int p_tilde, int N_within, int N_between, int N_both){ + matrix[p_tilde, p_tilde + 1] W_tilde; + matrix[p_tilde, p_tilde] W_tilde_cov; + matrix[p_tilde, p_tilde + 1] B_tilde; + matrix[p_tilde, p_tilde] B_tilde_cov; + vector[p_tilde] Mu_WB_tilde; + vector[N_between] Mu_z; + int N_wo_b = p_tilde - N_between; + vector[N_wo_b] Mu_y; + vector[N_wo_b] Mu_w; + vector[N_wo_b] Mu_b; + vector[N_between + N_wo_b] Mu_full; + matrix[N_between, N_between] Sigma_zz; + matrix[N_between, N_between] Sigma_zz_inv; + real Sigma_zz_ld; + matrix[N_wo_b, N_between] Sigma_yz; + matrix[N_wo_b, N_between] Sigma_yz_zi; + matrix[N_wo_b, N_wo_b] Sigma_b; + matrix[N_wo_b, N_wo_b] Sigma_b_z; + matrix[N_wo_b, N_wo_b] Sigma_w; + matrix[N_wo_b, N_wo_b] Sigma_w_inv; + real Sigma_w_ld; + matrix[N_wo_b, N_wo_b] Sigma_j; + matrix[N_wo_b, N_wo_b] Sigma_j_inv; + matrix[N_wo_b, N_between] Sigma_ji_yz_zi; + matrix[N_between, N_between] Vinv_11; + real Sigma_j_ld; + vector[nclus_sizes] L; + vector[nclus_sizes] B; + array[N_between] int bidx; + array[p_tilde - N_between] int notbidx; + real q_zz; + real q_yz; + real q_yyc; + vector[nclus_sizes] P; + vector[nclus_sizes] q_W; + vector[nclus_sizes] L_W; + vector[nclus_sizes] loglik; + vector[nclus_sizes] nperclus = to_vector(clus_sizes) .* to_vector(clus_size_ns); + int cluswise = 0; + int r1 = 1; // running index of rows of YX corresponding to each cluster + // 1. compute necessary vectors/matrices, like lav_mvnorm_cluster_implied22l() of lav_mvnorm_cluster.R + if (nclus[2] == nclus_sizes) cluswise = 1; + W_tilde = calc_W_tilde(impl_Sigmaw, impl_Muw, ov_idx1, p_tilde); + W_tilde_cov = block(W_tilde, 1, 2, p_tilde, p_tilde); + B_tilde = calc_B_tilde(impl_Sigmab, impl_Mub, ov_idx2, p_tilde); + B_tilde_cov = block(B_tilde, 1, 2, p_tilde, p_tilde); + Mu_WB_tilde = rep_vector(0, p_tilde); + if (N_within > 0) { + for (i in 1:N_within) { + Mu_WB_tilde[within_idx[i]] = W_tilde[within_idx[i], 1]; + B_tilde[within_idx[i], 1] = 0; + } + } + if (N_both > 0) { + for (i in 1:N_both) { + Mu_WB_tilde[both_idx[i]] = B_tilde[both_idx[i], 1] + W_tilde[both_idx[i], 1]; + } + } + // around line 71 of lav_mvnorm_cluster.R + if (N_between > 0) { + bidx = between_idx[1:N_between]; + notbidx = between_idx[(N_between + 1):p_tilde]; + Mu_z = to_vector(B_tilde[bidx, 1]); + Mu_y = Mu_WB_tilde[notbidx]; + Mu_w = to_vector(W_tilde[notbidx, 1]); + Mu_b = to_vector(B_tilde[notbidx, 1]); + Sigma_zz = B_tilde_cov[bidx, bidx]; + Sigma_yz = B_tilde_cov[notbidx, bidx]; + Sigma_b = B_tilde_cov[notbidx, notbidx]; + Sigma_w = W_tilde_cov[notbidx, notbidx]; + } else { + Mu_y = Mu_WB_tilde; + Mu_w = W_tilde[,1]; + Mu_b = B_tilde[,1]; + Sigma_b = B_tilde_cov; + Sigma_w = W_tilde_cov; + } + // 2. compute lpdf, around line 203 of lav_mvnorm_cluster + Sigma_w_inv = inverse_spd(Sigma_w); + Sigma_w_ld = log_determinant(Sigma_w); + if (N_between > 0) { + Sigma_zz_inv = inverse_spd(Sigma_zz); + Sigma_zz_ld = log_determinant(Sigma_zz); + Sigma_yz_zi = Sigma_yz * Sigma_zz_inv; + Sigma_b_z = Sigma_b - Sigma_yz * Sigma_yz_zi'; + } else { + Sigma_zz_ld = 0; + Sigma_b_z = Sigma_b; + } + Mu_full = append_row(Mu_z, Mu_y); + for (clz in 1:nclus_sizes) { + int nj = clus_sizes[clz]; + matrix[N_between + N_wo_b, N_between + N_wo_b] Y2Yc = crossprod(to_matrix(mean_d[clz] - Mu_full)'); + matrix[N_between, N_between] Y2Yc_zz; + matrix[N_wo_b, N_between] Y2Yc_yz; + matrix[N_wo_b, N_wo_b] Y2Yc_yy; + array[N_between] int uord_bidx; + array[N_wo_b] int uord_notbidx; + if (!cluswise) Y2Yc += cov_d[clz]; // variability between clusters of same size, will always equal 0 for clusterwise + if (N_between > 0) { + for (k in 1:N_between) { + uord_bidx[k] = k; + } + if (N_wo_b > 0) { + for (k in 1:N_wo_b) { + uord_notbidx[k] = N_between + k; + } + } + Y2Yc_zz = Y2Yc[uord_bidx, uord_bidx]; + Y2Yc_yz = Y2Yc[uord_notbidx, uord_bidx]; + Y2Yc_yy = Y2Yc[uord_notbidx, uord_notbidx]; + } else { + Y2Yc_yy = Y2Yc; + } + Sigma_j = (nj * Sigma_b_z) + Sigma_w; + Sigma_j_inv = inverse_spd(Sigma_j); // FIXME npd exceptions for within-only + Sigma_j_ld = log_determinant(Sigma_j); + L[clz] = Sigma_zz_ld + Sigma_j_ld; + if (N_between > 0) { + Sigma_ji_yz_zi = Sigma_j_inv * Sigma_yz_zi; + Vinv_11 = Sigma_zz_inv + nj * (Sigma_yz_zi' * Sigma_ji_yz_zi); + q_zz = sum(Vinv_11 .* Y2Yc_zz); + q_yz = -nj * sum(Sigma_ji_yz_zi .* Y2Yc_yz); + } else { + q_zz = 0; + q_yz = 0; + } + q_yyc = -nj * sum(Sigma_j_inv .* Y2Yc_yy); + B[clz] = q_zz + 2 * q_yz - q_yyc; + if (cluswise) { + matrix[N_wo_b, nj] Y_j; + if (N_between > 0) { + for (i in 1:nj) { + Y_j[, i] = YX[r1 - 1 + i, notbidx] - mean_d[clz, uord_notbidx]; // would be nice to have to_matrix() here + } + } else { + for (i in 1:nj) { + Y_j[, i] = YX[r1 - 1 + i] - mean_d[clz]; // would be nice to have to_matrix() here + } + } + r1 += nj; // for next iteration through loop + q_W[clz] = sum(Sigma_w_inv .* tcrossprod(Y_j)); + } + } + if (!cluswise) { + q_W = (nperclus - to_vector(clus_size_ns)) * sum(Sigma_w_inv .* S_PW); + } + L_W = (nperclus - to_vector(clus_size_ns)) * Sigma_w_ld; + loglik = -.5 * ((L .* to_vector(clus_size_ns)) + (B .* to_vector(clus_size_ns)) + q_W + L_W); + // add constant, line 300 lav_mvnorm_cluster + P = nperclus * (N_within + N_both) + to_vector(clus_size_ns) * N_between; + loglik += -.5 * (P * log(2 * pi())); + return loglik; + } + /* + Fills in the elements of a coefficient matrix containing some mix of + totally free, free subject to a sign constraint, and fixed elements + @param free_elements vector of unconstrained elements + @param skeleton matrix of the same dimensions as the output whose elements are + positive_infinity(): if output element is totally free + other: if output element is fixed to that number + @return matrix of coefficients + */ + matrix fill_matrix(vector free_elements, matrix skeleton, array[,] int eq_skeleton, int pos_start, int spos_start) { + int R = rows(skeleton); + int C = cols(skeleton); + matrix[R, C] out; + int pos = spos_start; // position of eq_skeleton + int freepos = pos_start; // position of free_elements + int eqelem = 0; + for (c in 1:C) for (r in 1:R) { + real rc = skeleton[r, c]; + if (is_inf(rc)) { // free + int eq = eq_skeleton[pos, 1]; + int wig = eq_skeleton[pos, 3]; + if (eq == 0 || wig == 1) { + out[r,c] = free_elements[freepos]; + freepos += 1; + } else { + eqelem = eq_skeleton[pos, 2]; + out[r,c] = free_elements[eqelem]; + } + pos += 1; + } else out[r,c] = skeleton[r, c]; // fixed, so do not bump pos + } + return out; + } + vector fill_prior(vector free_elements, array[] real pri_mean, array[,] int eq_skeleton) { + int R = dims(eq_skeleton)[1]; + int eqelem = 0; + int pos = 1; + vector[num_elements(pri_mean)] out; + for (r in 1:R) { + if (pos <= num_elements(pri_mean)) { + int eq = eq_skeleton[r, 1]; + int wig = eq_skeleton[r, 3]; + if (eq == 0) { + out[pos] = pri_mean[pos]; + pos += 1; + } else if (wig == 1) { + eqelem = eq_skeleton[r, 2]; + out[pos] = free_elements[eqelem]; + pos += 1; + } + } + } + return out; + } + /* + * This is a bug-free version of csr_to_dense_matrix and has the same arguments + */ + matrix to_dense_matrix(int m, int n, vector w, array[] int v, array[] int u) { + matrix[m, n] out = rep_matrix(0, m, n); + int pos = 1; + for (i in 1:m) { + int start = u[i]; + int nnz = u[i + 1] - start; + for (j in 1:nnz) { + out[i, v[pos]] = w[pos]; + pos += 1; + } + } + return out; + } + // sign function + int sign(real x) { + if (x > 0) + return 1; + else + return -1; + } + // sign-constrain a vector of loadings + vector sign_constrain_load(vector free_elements, int npar, array[,] int sign_mat) { + vector[npar] out; + for (i in 1:npar) { + if (sign_mat[i,1]) { + int lookupval = sign_mat[i,2]; + if (free_elements[lookupval] < 0) { + out[i] = -free_elements[i]; + } else { + out[i] = free_elements[i]; + } + } else { + out[i] = free_elements[i]; + } + } + return out; + } + // sign-constrain a vector of regressions or covariances + vector sign_constrain_reg(vector free_elements, int npar, array[,] int sign_mat, vector load_par1, vector load_par2) { + vector[npar] out; + for (i in 1:npar) { + if (sign_mat[i,1]) { + int lookupval1 = sign_mat[i,2]; + int lookupval2 = sign_mat[i,3]; + if (sign(load_par1[lookupval1]) * sign(load_par2[lookupval2]) < 0) { + out[i] = -free_elements[i]; + } else { + out[i] = free_elements[i]; + } + } else { + out[i] = free_elements[i]; + } + } + return out; + } + // obtain covariance parameter vector for correlation/sd matrices + vector cor2cov(array[] matrix cormat, array[] matrix sdmat, int num_free_elements, array[] matrix matskel, array[,] int wskel, int ngrp) { + vector[num_free_elements] out; + int R = rows(to_matrix(cormat[1])); + int pos = 1; // position of eq_skeleton + int freepos = 1; // position of free_elements + for (g in 1:ngrp) { + for (c in 1:(R-1)) for (r in (c+1):R) { + if (is_inf(matskel[g,r,c])) { + if (wskel[pos,1] == 0) { + out[freepos] = sdmat[g,r,r] * sdmat[g,c,c] * cormat[g,r,c]; + freepos += 1; + } + pos += 1; + } + } + } + return out; + } + // E step of EM algorithm on latent continuous space + array[] matrix estep(array[] vector YXstar, array[] vector Mu, array[] matrix Sigma, array[] int Nobs, array[,] int Obsvar, array[] int startrow, array[] int endrow, array[] int grpnum, int Np, int Ng) { + int p = dims(YXstar)[2]; + array[Ng] matrix[p, p + 1] out; //mean vec + cov mat + matrix[dims(YXstar)[1], p] YXfull; // columns consistenly ordered + matrix[p, p] T2pat; + array[p] int obsidx; + int r1; + int r2; + int grpidx; + int Nmis; + for (g in 1:Ng) { + out[g] = rep_matrix(0, p, p + 1); + } + for (mm in 1:Np) { + obsidx = Obsvar[mm,]; + r1 = startrow[mm]; + r2 = endrow[mm]; + grpidx = grpnum[mm]; + Nmis = p - Nobs[mm]; + if (Nobs[mm] < p) { + matrix[Nobs[mm], Nobs[mm]] Sig22 = Sigma[grpidx, obsidx[1:Nobs[mm]], obsidx[1:Nobs[mm]]]; + matrix[Nmis, Nmis] Sig11 = Sigma[grpidx, obsidx[(Nobs[mm] + 1):p], obsidx[(Nobs[mm] + 1):p]]; + matrix[Nmis, Nobs[mm]] Sig12 = Sigma[grpidx, obsidx[(Nobs[mm] + 1):p], obsidx[1:Nobs[mm]]]; + matrix[Nobs[mm], Nobs[mm]] S22inv = inverse_spd(Sig22); + matrix[Nmis, Nmis] T2p11 = Sig11 - (Sig12 * S22inv * Sig12'); + // partition into observed/missing, compute Sigmas, add to out + for (jj in r1:r2) { + vector[Nmis] ymis; + ymis = Mu[grpidx, obsidx[(Nobs[mm] + 1):p]] + (Sig12 * S22inv * (YXstar[jj, 1:Nobs[mm]] - Mu[grpidx, obsidx[1:Nobs[mm]]])); + for (kk in 1:Nobs[mm]) { + YXfull[jj, obsidx[kk]] = YXstar[jj, kk]; + } + for (kk in (Nobs[mm] + 1):p) { + YXfull[jj, obsidx[kk]] = ymis[kk - Nobs[mm]]; + } + } + T2pat = crossprod(YXfull[r1:r2,]); + // correction for missing cells/conditional covariances + for (jj in 1:Nmis) { + for (kk in jj:Nmis) { + T2pat[obsidx[Nobs[mm] + jj], obsidx[Nobs[mm] + kk]] = T2pat[obsidx[Nobs[mm] + jj], obsidx[Nobs[mm] + kk]] + (r2 - r1 + 1) * T2p11[jj, kk]; + if (kk > jj) { + T2pat[obsidx[Nobs[mm] + kk], obsidx[Nobs[mm] + jj]] = T2pat[obsidx[Nobs[mm] + jj], obsidx[Nobs[mm] + kk]]; + } + } + } + } else { + // complete data + for (jj in r1:r2) { + for (kk in 1:Nobs[mm]) { + YXfull[jj, obsidx[kk]] = YXstar[jj, kk]; + } + } + T2pat = crossprod(YXfull[r1:r2,]); + } + for (i in 1:p) { + out[grpidx,i,1] += sum(YXfull[r1:r2,i]); + } + out[grpidx,,2:(p+1)] += T2pat; + } + return out; + } + matrix sig_inv_update(matrix Sigmainv, array[] int obsidx, int Nobs, int np, real logdet) { + matrix[Nobs + 1, Nobs + 1] out = rep_matrix(0, Nobs + 1, Nobs + 1); + int nrm = np - Nobs; + matrix[nrm, nrm] H; + matrix[nrm, Nobs] A; + if (nrm == 0) { + out[1:Nobs, 1:Nobs] = Sigmainv; + out[Nobs + 1, Nobs + 1] = logdet; + } else { + H = Sigmainv[obsidx[(Nobs + 1):np], obsidx[(Nobs + 1):np]]; + A = Sigmainv[obsidx[(Nobs + 1):np], obsidx[1:Nobs]]; + out[1:Nobs, 1:Nobs] = Sigmainv[obsidx[1:Nobs], obsidx[1:Nobs]] - A' * mdivide_left_spd(H, A); + out[Nobs + 1, Nobs + 1] = logdet + log_determinant(H); + } + return out; + } + real multi_normal_suff(vector xbar, matrix S, vector Mu, matrix Supdate, int N) { + int Nobs = dims(S)[1]; + real out; + // using elementwise multiplication + sum here for efficiency + out = -.5 * N * ( sum(Supdate[1:Nobs, 1:Nobs] .* (S + (xbar - Mu) * (xbar - Mu)')) + Supdate[Nobs + 1, Nobs + 1] + Nobs * log(2 * pi()) ); + if(is_nan(out) || out == positive_infinity()) out = negative_infinity(); + return out; + } + // compute mean vectors and cov matrices for a single group (two-level models) + array[] vector calc_mean_vecs(array[] vector YXstar, array[] vector mean_d, array[] int nclus, array[] int Xvar, array[] int Xbetvar, int Nx, int Nx_between, int p_tilde) { + vector[Nx] ov_mean = rep_vector(0, Nx); + vector[Nx_between] ov_mean_d = rep_vector(0, Nx_between); + int nr = dims(YXstar)[1]; + array[2] vector[p_tilde] out; + for (i in 1:2) out[i] = rep_vector(0, p_tilde); + if (Nx > 0) { + for (i in 1:nr) { + ov_mean += YXstar[i, Xvar[1:Nx]]; + } + ov_mean *= pow(nclus[1], -1); + out[1, 1:Nx] = ov_mean; + } + if (Nx_between > 0) { + for (cc in 1:nclus[2]) { + ov_mean_d += mean_d[cc, 1:Nx_between]; + } + ov_mean_d *= pow(nclus[2], -1); + out[2, 1:Nx_between] = ov_mean_d; + } + return out; + } + array[] matrix calc_cov_mats(array[] vector YXstar, array[] vector mean_d, array[] vector mean_vecs, array[] int nclus, array[] int Xvar, array[] int Xbetvar, int Nx, int Nx_between, int p_tilde) { + matrix[Nx_between, Nx_between] cov_mean_d = rep_matrix(0, Nx_between, Nx_between); + matrix[Nx, Nx] cov_w = rep_matrix(0, Nx, Nx); + matrix[Nx, Nx] cov_w_inv; + int nr = dims(YXstar)[1]; + array[3] matrix[p_tilde, p_tilde] out; + for (i in 1:3) out[i] = rep_matrix(0, p_tilde, p_tilde); + if (Nx > 0) { + for (i in 1:nr) { + cov_w += tcrossprod(to_matrix(YXstar[i, Xvar[1:Nx]] - mean_vecs[1, 1:Nx])); + } + cov_w *= pow(nclus[1], -1); + cov_w_inv[1:Nx, 1:Nx] = inverse_spd(cov_w); + out[2, 1:Nx, 1:Nx] = cov_w; + out[3, 1:Nx, 1:Nx] = cov_w_inv; + out[3, Nx + 1, Nx + 1] = log_determinant(cov_w); // need log_determinant for multi_normal_suff + } + if (Nx_between > 0) { + for (cc in 1:nclus[2]) { + cov_mean_d += tcrossprod(to_matrix(mean_d[cc, 1:Nx_between] - mean_vecs[2, 1:Nx_between])); + } + cov_mean_d *= pow(nclus[2], -1); + out[1, 1:Nx_between, 1:Nx_between] = cov_mean_d; + } + return out; + } + // compute log_lik of fixed.x variables for a single group (two-level models) + vector calc_log_lik_x(array[] vector mean_d, vector ov_mean_d, matrix cov_mean_d, matrix cov_w, matrix cov_w_inv, array[] int nclus, array[] int cluster_size, array[] int Xvar, array[] int Xbetvar, int Nx, int Nx_between) { + vector[nclus[2]] out = rep_vector(0, nclus[2]); + for (cc in 1:nclus[2]) { + if (Nx > 0) { + out[cc] += multi_normal_suff(mean_d[cc, Xvar[1:Nx]], cov_w[1:Nx, 1:Nx], mean_d[cc, Xvar[1:Nx]], cov_w_inv[1:(Nx + 1), 1:(Nx + 1)], cluster_size[cc]); + } + if (Nx_between > 0) { + out[cc] += multi_normal_lpdf(mean_d[cc, 1:Nx_between] | ov_mean_d[1:Nx_between], cov_mean_d[1:Nx_between, 1:Nx_between]); + } + } + return out; + } + // fill covariance matrix with blocks + array[] matrix fill_cov(array[] matrix covmat, array[,] int blkse, array[] int nblk, + array[] matrix mat_1, array[] matrix mat_2, array[] matrix mat_3, + array[] matrix mat_4, array[] matrix mat_5, array[,] int psiorder, + array[,] int psirevord) { + int Ng = dims(covmat)[1]; + array[Ng] matrix[dims(covmat)[2], dims(covmat)[3]] out; + for (g in 1:Ng) { + out[g] = covmat[g, psiorder[g], psiorder[g]]; + } + for (k in 1:sum(nblk)) { + int blkidx = blkse[k, 6]; + int arrayidx = blkse[k, 5]; + int blkgrp = blkse[k, 4]; + int srow = blkse[k, 1]; + int erow = blkse[k, 2]; + if (arrayidx == 1) { + out[blkgrp, srow:erow, srow:erow] = mat_1[blkidx]; + } else if (arrayidx == 2) { + out[blkgrp, srow:erow, srow:erow] = mat_2[blkidx]; + } else if (arrayidx == 3) { + out[blkgrp, srow:erow, srow:erow] = mat_3[blkidx]; + } else if (arrayidx == 4) { + out[blkgrp, srow:erow, srow:erow] = mat_4[blkidx]; + } else { + out[blkgrp, srow:erow, srow:erow] = mat_5[blkidx]; + } + } + for (g in 1:Ng) { + out[g] = out[g, psirevord[g], psirevord[g]]; + } + return out; + } + real eval_priors(vector Lambda_y_free, vector B_free, vector Nu_free, vector Alpha_free, vector sd0, vector lambda_y_primn, array[] real lambda_y_sd, array[] int lambda_y_pri, array[] int lambda_y_blk, vector b_primn, array[] real b_sd, array[] int b_pri, array[] int b_blk, vector nu_primn, array[] real nu_sd, array[] int nu_pri, array[] int nu_blk, vector alpha_primn, array[] real alpha_sd, array[] int alpha_pri, array[] int alpha_blk, array[] int len_free, vector blkparm1, vector blkparm2, int priblklen, int npriblks) { + real out = 0.0; + if (npriblks == 0) { + out += normal_lpdf(Lambda_y_free | lambda_y_primn, lambda_y_sd); + out += normal_lpdf(B_free | b_primn, b_sd); + out += normal_lpdf(Nu_free | nu_primn, nu_sd); + out += normal_lpdf(Alpha_free | alpha_primn, alpha_sd); + } else { + for (i in 1:npriblks) { + vector[priblklen] parvec; + int npars = 1; + if (len_free[1] > 0) { + for (j in 1:len_free[1]) { + if (lambda_y_pri[j] == 0) { + out += normal_lpdf(Lambda_y_free[j] | lambda_y_primn[j], lambda_y_sd[j]); + } else if (lambda_y_blk[j] == i && lambda_y_pri[j] == 1) { + parvec[npars] = Lambda_y_free[j]; + npars += 1; + } + } + } + if (len_free[4] > 0) { + for (j in 1:len_free[4]) { + if (b_pri[j] == 0) { + out += normal_lpdf(B_free[j] | b_primn[j], b_sd[j]); + } else if (b_blk[j] == i && b_pri[j] == 1) { + parvec[npars] = B_free[j]; + npars += 1; + } + } + } + if (len_free[13] > 0) { + for (j in 1:len_free[13]) { + if (nu_pri[j] == 0) { + out += normal_lpdf(Nu_free[j] | nu_primn[j], nu_sd[j]); + } else if (nu_blk[j] == i && nu_pri[j] == 1) { + parvec[npars] = Nu_free[j]; + npars += 1; + } + } + } + if (len_free[14] > 0) { + for (j in 1:len_free[14]) { + if (alpha_pri[j] == 0) { + out += normal_lpdf(Alpha_free[j] | alpha_primn[j], alpha_sd[j]); + } else if (alpha_blk[j] == i && alpha_pri[j] == 1) { + parvec[npars] = Alpha_free[j]; + npars += 1; + } + } + } + out += normal_lpdf(parvec[1:(npars - 1)] | 0, sd0[i]); + out += student_t_lpdf(sd0[i] | blkparm1[i], 0, blkparm2[i]) - log(.5); // left-truncated at 0 + } + } + return out; + } +} +data { + // see p. 2 https://books.google.com/books?id=9AC-s50RjacC + int p; // number of manifest response variables + int p_c; // number of manifest level 2 variables + int q; // number of manifest predictors + int m; // number of latent endogenous variables + int m_c; // number of latent level 2 variables + int n; // number of latent exogenous variables + int Ng; // number of groups + int missing; // are there missing values? + int save_lvs; // should we save lvs? + int Np; // number of group-by-missing patterns combos + array[Ng] int N; // number of observations per group + array[Np] int Nobs; // number of observed variables in each missing pattern + array[Np] int Nordobs; // number of ordinal observed variables in each missing pattern + array[Np, p + q] int Obsvar; // indexing of observed variables + int Ntot; // number of observations across all groups + array[Np] int startrow; // starting row for each missing pattern + array[Np] int endrow; // ending row for each missing pattern + array[Np] int grpnum; // group number for each row of data + int wigind; // do any parameters have approx equality constraint ('wiggle')? + int has_data; // are the raw data on y and x available? + int ord; // are there any ordinal variables? + int multilev; // is this a multilevel dataset? + int Nord; // how many ordinal variables? + array[Nord] int ordidx; // indexing of ordinal variables + array[Np, Nord] int OrdObsvar; // indexing of observed ordinal variables in YXo + int Noent; // how many observed entries of ordinal variables (for data augmentation) + array[p + q - Nord] int contidx; // indexing of continuous variables + array[Nord] int nlevs; // how many levels does each ordinal variable have + array[Ng, 2] int nclus; // number of level 1 + level 2 observations + int p_tilde; // total number of variables + array[Ntot] vector[multilev ? p_tilde : p + q - Nord] YX; // continuous data + array[Ntot, Nord] int YXo; // ordinal data + array[Np] int Nx; // number of fixed.x variables (within) + array[Np] int Nx_between; // number of fixed.x variables (between) + int use_cov; + int pri_only; + int emiter; // number of em iterations for saturated model in ppp (missing data only) + int use_suff; // should we compute likelihood via mvn sufficient stats? + int do_test; // should we do everything in generated quantities? + array[Np] vector[multilev ? p_tilde : p + q - Nord] YXbar; // sample means of continuous manifest variables + array[Np] matrix[multilev ? (p_tilde + 1) : (p + q - Nord + 1), multilev ? (p_tilde + 1) : (p + q - Nord + 1)] S; // sample covariance matrix among all continuous manifest variables NB!! multiply by (N-1) to use wishart lpdf!! + int npriblks; // how many blocks of parameters for prior distributions? + int priblklen; // max number of parameters in a block + vector[npriblks] blkparm1; // parameters of block priors + vector[npriblks] blkparm2; + array[sum(nclus[,2])] int cluster_size; // number of obs per cluster + array[Ng] int ncluster_sizes; // number of unique cluster sizes + array[sum(ncluster_sizes)] int cluster_sizes; // unique cluster sizes + array[sum(ncluster_sizes)] int cluster_size_ns; // number of clusters of each size + array[Np, multilev ? p_tilde : p + q] int Xvar; // indexing of fixed.x variables (within) + array[Np, multilev ? p_tilde : p + q] int Xdatvar; // indexing of fixed.x in data (differs from Xvar when missing) + array[Np, multilev ? p_tilde : p + q] int Xbetvar; // indexing of fixed.x variables (between) + array[sum(ncluster_sizes)] vector[p_tilde] mean_d; // sample means by unique cluster size + array[sum(ncluster_sizes)] matrix[p_tilde, p_tilde] cov_d; // sample covariances by unique cluster size + array[Ng] matrix[p_tilde, p_tilde] cov_w; // observed "within" covariance matrix + array[sum(nclus[,2])] vector[p_tilde] mean_d_full; // sample means/covs by cluster, for clusterwise log-densities + array[sum(nclus[,2])] matrix[p_tilde, p_tilde] cov_d_full; + array[Ng] vector[p_tilde] xbar_w; // data estimates of within/between means/covs (for saturated logl) + array[Ng] vector[p_tilde] xbar_b; + array[Ng] matrix[p_tilde, p_tilde] cov_b; + array[Ng] real gs; // group size constant, for computation of saturated logl + int N_within; // number of within variables + int N_between; // number of between variables + int N_both; // number of variables at both levels + array[2] int N_lev; // number of observed variables at each level + array[N_within] int within_idx; + array[p_tilde] int between_idx; // between indexing, followed by within/both + array[N_lev[1]] int ov_idx1; + array[N_lev[2]] int ov_idx2; + array[N_both] int both_idx; + vector[multilev ? sum(ncluster_sizes) : Ng] log_lik_x; // ll of fixed x variables by unique cluster size + vector[multilev ? sum(nclus[,2]) : Ng] log_lik_x_full; // ll of fixed x variables by cluster + /* sparse matrix representations of skeletons of coefficient matrices, + which is not that interesting but necessary because you cannot pass + missing values into the data block of a Stan program from R */ + int len_w1; // max number of free elements in Lambda_y per grp + array[Ng] int wg1; // number of free elements in Lambda_y per grp + array[Ng] vector[len_w1] w1; // values of free elements in Lambda_y + array[Ng, len_w1] int v1; // index of free elements in Lambda_y + array[Ng, p + 1] int u1; // index of free elements in Lambda_y + array[sum(wg1), 3] int w1skel; + array[sum(wg1), 2] int lam_y_sign; + int len_lam_y; // number of free elements minus equality constraints + array[len_lam_y] real lambda_y_mn; // prior + array[len_lam_y] real lambda_y_sd; + array[len_lam_y] int lambda_y_pri; + array[len_lam_y] int lambda_y_blk; + // same things but for B + int len_w4; + array[Ng] int wg4; + array[Ng] vector[len_w4] w4; + array[Ng, len_w4] int v4; + array[Ng, m + 1] int u4; + array[sum(wg4), 3] int w4skel; + array[sum(wg4), 3] int b_sign; + int len_b; + array[len_b] real b_mn; + array[len_b] real b_sd; + array[len_b] int b_pri; + array[len_b] int b_blk; + // same things but for diag(Theta) + int len_w5; + array[Ng] int wg5; + array[Ng] vector[len_w5] w5; + array[Ng, len_w5] int v5; + array[Ng, p + 1] int u5; + array[sum(wg5), 3] int w5skel; + int len_thet_sd; + array[len_thet_sd] real theta_sd_shape; + array[len_thet_sd] real theta_sd_rate; + int theta_pow; + // same things but for Theta_r + int len_w7; + array[Ng] int wg7; + array[Ng] vector[len_w7] w7; + array[Ng, len_w7] int v7; + array[Ng, p + 1] int u7; + array[sum(wg7), 3] int w7skel; + int len_thet_r; + array[len_thet_r] real theta_r_alpha; + array[len_thet_r] real theta_r_beta; + // for blocks within Theta_r that receive lkj + array[5] int thetanblk; + array[5] int thetadims; + array[sum(thetanblk), 6] int thetablkse; + vector[sum(thetanblk)] thetablkpri; + int len_w8; + array[Ng] int wg8; + array[Ng] vector[len_w8] w8; + array[Ng, len_w8] int v8; + array[Ng, p + 1] int u8; + array[sum(wg8), 3] int w8skel; + array[Ng, p] int thetaorder; + array[Ng, p] int thetarevord; + // same things but for Psi + int len_w9; + array[Ng] int wg9; + array[Ng] vector[len_w9] w9; + array[Ng, len_w9] int v9; + array[Ng, m + 1] int u9; + array[sum(wg9), 3] int w9skel; + int len_psi_sd; + array[len_psi_sd] real psi_sd_shape; + array[len_psi_sd] real psi_sd_rate; + int psi_pow; + // same things but for Psi_r + int len_w10; + array[Ng] int wg10; + array[Ng] vector[len_w10] w10; + array[Ng, len_w10] int v10; + array[Ng, m + 1] int u10; + array[sum(wg10), 3] int w10skel; + array[sum(wg10), 3] int psi_r_sign; + int len_psi_r; + array[len_psi_r] real psi_r_alpha; + array[len_psi_r] real psi_r_beta; + // for blocks within Psi_r that receive lkj + array[5] int psinblk; + array[5] int psidims; + array[sum(psinblk), 6] int psiblkse; + vector[sum(psinblk)] psiblkpri; + int len_w11; + array[Ng] int wg11; + array[Ng] vector[len_w11] w11; + array[Ng, len_w11] int v11; + array[Ng, m + 1] int u11; + array[sum(wg11), 3] int w11skel; + array[Ng, m] int psiorder; + array[Ng, m] int psirevord; + // same things but for Nu + int len_w13; + array[Ng] int wg13; + array[Ng] vector[len_w13] w13; + array[Ng, len_w13] int v13; + array[Ng, use_cov ? 1 : p + q + 1] int u13; + array[sum(wg13), 3] int w13skel; + int len_nu; + array[len_nu] real nu_mn; + array[len_nu] real nu_sd; + array[len_nu] int nu_pri; + array[len_nu] int nu_blk; + // same things but for Alpha + int len_w14; + array[Ng] int wg14; + array[Ng] vector[len_w14] w14; + array[Ng, len_w14] int v14; + array[Ng, use_cov ? 1 : m + n + 1] int u14; + array[sum(wg14), 3] int w14skel; + array[sum(wg14), 3] int alph_sign; + int len_alph; + array[len_alph] real alpha_mn; + array[len_alph] real alpha_sd; + array[len_alph] int alpha_pri; + array[len_alph] int alpha_blk; + // same things but for Tau + int len_w15; + array[Ng] int wg15; + array[Ng] vector[len_w15] w15; + array[Ng, len_w15] int v15; + array[Ng, sum(nlevs) - Nord + 1] int u15; + array[sum(wg15), 3] int w15skel; + int len_tau; + array[len_tau] real tau_mn; + array[len_tau] real tau_sd; + // Level 2 matrices start here!! + // Lambda + int len_w1_c; + array[Ng] int wg1_c; + array[Ng] vector[len_w1_c] w1_c; + array[Ng, len_w1_c] int v1_c; + array[Ng, p_c + 1] int u1_c; + array[sum(wg1_c), 3] int w1skel_c; + array[sum(wg1_c), 2] int lam_y_sign_c; + int len_lam_y_c; + array[len_lam_y_c] real lambda_y_mn_c; + array[len_lam_y_c] real lambda_y_sd_c; + // same things but for B + int len_w4_c; + array[Ng] int wg4_c; + array[Ng] vector[len_w4_c] w4_c; + array[Ng, len_w4_c] int v4_c; + array[Ng, m_c + 1] int u4_c; + array[sum(wg4_c), 3] int w4skel_c; + array[sum(wg4_c), 3] int b_sign_c; + int len_b_c; + array[len_b_c] real b_mn_c; + array[len_b_c] real b_sd_c; + // same things but for diag(Theta) + int len_w5_c; + array[Ng] int wg5_c; + array[Ng] vector[len_w5_c] w5_c; + array[Ng, len_w5_c] int v5_c; + array[Ng, p_c + 1] int u5_c; + array[sum(wg5_c), 3] int w5skel_c; + int len_thet_sd_c; + array[len_thet_sd_c] real theta_sd_shape_c; + array[len_thet_sd_c] real theta_sd_rate_c; + int theta_pow_c; + // same things but for Theta_r + int len_w7_c; + array[Ng] int wg7_c; + array[Ng] vector[len_w7_c] w7_c; + array[Ng, len_w7_c] int v7_c; + array[Ng, p_c + 1] int u7_c; + array[sum(wg7_c), 3] int w7skel_c; + int len_thet_r_c; + array[len_thet_r_c] real theta_r_alpha_c; + array[len_thet_r_c] real theta_r_beta_c; + // for blocks within Theta_r that receive lkj + array[5] int thetanblk_c; + array[5] int thetadims_c; + array[sum(thetanblk_c), 6] int thetablkse_c; + vector[sum(thetanblk_c)] thetablkpri_c; + int len_w8_c; + array[Ng] int wg8_c; + array[Ng] vector[len_w8_c] w8_c; + array[Ng, len_w8_c] int v8_c; + array[Ng, p_c + 1] int u8_c; + array[sum(wg8_c), 3] int w8skel_c; + array[Ng, p_c] int thetaorder_c; + array[Ng, p_c] int thetarevord_c; + // same things but for Psi + int len_w9_c; + array[Ng] int wg9_c; + array[Ng] vector[len_w9_c] w9_c; + array[Ng, len_w9_c] int v9_c; + array[Ng, m_c + 1] int u9_c; + array[sum(wg9_c), 3] int w9skel_c; + int len_psi_sd_c; + array[len_psi_sd_c] real psi_sd_shape_c; + array[len_psi_sd_c] real psi_sd_rate_c; + int psi_pow_c; + // same things but for Psi_r + int len_w10_c; + array[Ng] int wg10_c; + array[Ng] vector[len_w10_c] w10_c; + array[Ng, len_w10_c] int v10_c; + array[Ng, m_c + 1] int u10_c; + array[sum(wg10_c), 3] int w10skel_c; + array[sum(wg10_c), 3] int psi_r_sign_c; + int len_psi_r_c; + array[len_psi_r_c] real psi_r_alpha_c; + array[len_psi_r_c] real psi_r_beta_c; + // for blocks within Psi_r that receive lkj + array[5] int psinblk_c; + array[5] int psidims_c; + array[sum(psinblk_c), 6] int psiblkse_c; + vector[sum(psinblk_c)] psiblkpri_c; + int len_w11_c; + array[Ng] int wg11_c; + array[Ng] vector[len_w11_c] w11_c; + array[Ng, len_w11_c] int v11_c; + array[Ng, m_c + 1] int u11_c; + array[sum(wg11_c), 3] int w11skel_c; + array[Ng, m_c] int psiorder_c; + array[Ng, m_c] int psirevord_c; + // same things but for Nu + int len_w13_c; + array[Ng] int wg13_c; + array[Ng] vector[len_w13_c] w13_c; + array[Ng, len_w13_c] int v13_c; + array[Ng, p_c + 1] int u13_c; + array[sum(wg13_c), 3] int w13skel_c; + int len_nu_c; + array[len_nu_c] real nu_mn_c; + array[len_nu_c] real nu_sd_c; + // same things but for Alpha + int len_w14_c; + array[Ng] int wg14_c; + array[Ng] vector[len_w14_c] w14_c; + array[Ng, len_w14_c] int v14_c; + array[Ng, m_c + 1] int u14_c; + array[sum(wg14_c), 3] int w14skel_c; + array[sum(wg14_c), 3] int alph_sign_c; + int len_alph_c; + array[len_alph_c] real alpha_mn_c; + array[len_alph_c] real alpha_sd_c; +} +transformed data { // (re)construct skeleton matrices in Stan (not that interesting) + array[Ng] matrix[p, m] Lambda_y_skeleton; + array[Ng] matrix[m, m] B_skeleton; + array[Ng] matrix[p, p] Theta_skeleton; + array[Ng] matrix[p, p] Theta_r_skeleton; + array[Ng] matrix[p, p] Theta_r_skeleton_f; + array[Ng] matrix[m, m] Psi_skeleton; + array[Ng] matrix[m, m] Psi_r_skeleton; + array[Ng] matrix[m, m] Psi_r_skeleton_f; + array[Ng] matrix[p, 1] Nu_skeleton; + array[Ng] matrix[m, 1] Alpha_skeleton; + array[Ng] matrix[sum(nlevs) - Nord, 1] Tau_skeleton; + array[Np] vector[ord ? 0 : (p + q)] YXbarstar; + array[Np] matrix[ord ? 0 : (p + q), ord ? 0 : (p + q)] Sstar; + array[Ng] matrix[p_c, m_c] Lambda_y_skeleton_c; + array[Ng] matrix[m_c, m_c] B_skeleton_c; + array[Ng] matrix[p_c, p_c] Theta_skeleton_c; + array[Ng] matrix[p_c, p_c] Theta_r_skeleton_c; + array[Ng] matrix[p_c, p_c] Theta_r_skeleton_f_c; + array[Ng] matrix[m_c, m_c] Psi_skeleton_c; + array[Ng] matrix[m_c, m_c] Psi_r_skeleton_c; + array[Ng] matrix[m_c, m_c] Psi_r_skeleton_f_c; + array[Ng] matrix[p_c, 1] Nu_skeleton_c; + array[Ng] matrix[m_c, 1] Alpha_skeleton_c; + matrix[m, m] I = diag_matrix(rep_vector(1, m)); + matrix[m_c, m_c] I_c = diag_matrix(rep_vector(1, m_c)); + int Ncont = p + q - Nord; + array[max(nclus[,2]) > 1 ? max(nclus[,2]) : 0] int intone; + array[Ng,2] int g_start1; + array[Ng,2] int g_start4; + array[Ng,2] int g_start5; + array[Ng,2] int g_start7; + array[Ng,2] int g_start9; + array[Ng,2] int g_start10; + array[Ng,2] int g_start13; + array[Ng,2] int g_start14; + array[Ng,2] int g_start15; + array[Ng,2] int g_start1_c; + array[Ng,2] int g_start4_c; + array[Ng,2] int g_start5_c; + array[Ng,2] int g_start7_c; + array[Ng,2] int g_start9_c; + array[Ng,2] int g_start10_c; + array[Ng,2] int g_start13_c; + array[Ng,2] int g_start14_c; + array[15] int len_free; + array[15] int pos; + array[15] int len_free_c; + array[15] int pos_c; + for (i in 1:15) { + len_free[i] = 0; + pos[i] = 1; + len_free_c[i] = 0; + pos_c[i] = 1; + } + for (g in 1:Ng) { + Lambda_y_skeleton[g] = to_dense_matrix(p, m, w1[g], v1[g,], u1[g,]); + B_skeleton[g] = to_dense_matrix(m, m, w4[g], v4[g,], u4[g,]); + Theta_skeleton[g] = to_dense_matrix(p, p, w5[g], v5[g,], u5[g,]); + Theta_r_skeleton[g] = to_dense_matrix(p, p, w7[g], v7[g,], u7[g,]); + Theta_r_skeleton_f[g] = to_dense_matrix(p, p, w8[g], v8[g,], u8[g,]); + Psi_skeleton[g] = to_dense_matrix(m, m, w9[g], v9[g,], u9[g,]); + Psi_r_skeleton[g] = to_dense_matrix(m, m, w10[g], v10[g,], u10[g,]); + Psi_r_skeleton_f[g] = to_dense_matrix(m, m, w11[g], v11[g,], u11[g,]); + if (!use_cov) { + Nu_skeleton[g] = to_dense_matrix((p + q), 1, w13[g], v13[g,], u13[g,]); + Alpha_skeleton[g] = to_dense_matrix((m + n), 1, w14[g], v14[g,], u14[g,]); + } + Tau_skeleton[g] = to_dense_matrix(sum(nlevs) - Nord, 1, w15[g], v15[g,], u15[g,]); + Lambda_y_skeleton_c[g] = to_dense_matrix(p_c, m_c, w1_c[g], v1_c[g,], u1_c[g,]); + B_skeleton_c[g] = to_dense_matrix(m_c, m_c, w4_c[g], v4_c[g,], u4_c[g,]); + Theta_skeleton_c[g] = to_dense_matrix(p_c, p_c, w5_c[g], v5_c[g,], u5_c[g,]); + Theta_r_skeleton_c[g] = to_dense_matrix(p_c, p_c, w7_c[g], v7_c[g,], u7_c[g,]); + Theta_r_skeleton_f_c[g] = to_dense_matrix(p_c, p_c, w8_c[g], v8_c[g,], u8_c[g,]); + Psi_skeleton_c[g] = to_dense_matrix(m_c, m_c, w9_c[g], v9_c[g,], u9_c[g,]); + Psi_r_skeleton_c[g] = to_dense_matrix(m_c, m_c, w10_c[g], v10_c[g,], u10_c[g,]); + Psi_r_skeleton_f_c[g] = to_dense_matrix(m_c, m_c, w11_c[g], v11_c[g,], u11_c[g,]); + Nu_skeleton_c[g] = to_dense_matrix(p_c, 1, w13_c[g], v13_c[g,], u13_c[g,]); + Alpha_skeleton_c[g] = to_dense_matrix(m_c, 1, w14_c[g], v14_c[g,], u14_c[g,]); + // count free elements in Lambda_y_skeleton + g_start1[g,1] = len_free[1] + 1; + g_start1[g,2] = pos[1]; + for (i in 1:p) { + for (j in 1:m) { + if (is_inf(Lambda_y_skeleton[g,i,j])) { + if (w1skel[pos[1],2] == 0 || w1skel[pos[1],3] == 1) len_free[1] += 1; + pos[1] += 1; + } + } + } + // same thing but for B_skeleton + g_start4[g,1] = len_free[4] + 1; + g_start4[g,2] = pos[4]; + for (i in 1:m) { + for (j in 1:m) { + if (is_inf(B_skeleton[g,i,j])) { + if (w4skel[pos[4],2] == 0 || w4skel[pos[4],3] == 1) len_free[4] += 1; + pos[4] += 1; + } + } + } + // same thing but for Theta_skeleton + g_start5[g,1] = len_free[5] + 1; + g_start5[g,2] = pos[5]; + for (i in 1:p) { + if (is_inf(Theta_skeleton[g,i,i])) { + if (w5skel[pos[5],2] == 0 || w5skel[pos[5],3] == 1) len_free[5] += 1; + pos[5] += 1; + } + } + // same thing but for Theta_r_skeleton + g_start7[g,1] = len_free[7] + 1; + g_start7[g,2] = pos[7]; + for (i in 1:(p-1)) { + for (j in (i+1):p) { + if (is_inf(Theta_r_skeleton[g,j,i])) { + if (w7skel[pos[7],2] == 0 || w7skel[pos[7],3] == 1) len_free[7] += 1; + pos[7] += 1; + } + if (is_inf(Theta_r_skeleton_f[g,j,i])) { + if (w8skel[pos[8],2] == 0 || w8skel[pos[8],3] == 1) len_free[8] += 1; + pos[8] += 1; + } + } + } + // same thing but for Psi_skeleton + g_start9[g,1] = len_free[9] + 1; + g_start9[g,2] = pos[9]; + for (i in 1:m) { + if (is_inf(Psi_skeleton[g,i,i])) { + if (w9skel[pos[9],2] == 0 || w9skel[pos[9],3] == 1) len_free[9] += 1; + pos[9] += 1; + } + } + // same thing but for Psi_r_skeleton + g_start10[g,1] = len_free[10] + 1; + g_start10[g,2] = pos[10]; + for (i in 1:(m-1)) { + for (j in (i+1):m) { + if (is_inf(Psi_r_skeleton[g,j,i])) { + if (w10skel[pos[10],2] == 0 || w10skel[pos[10],3] == 1) len_free[10] += 1; + pos[10] += 1; + } + if (is_inf(Psi_r_skeleton_f[g,j,i])) { + if (w11skel[pos[11],2] == 0 || w11skel[pos[11],3] == 1) len_free[11] += 1; + pos[11] += 1; + } + } + } + if (!use_cov) { + // same thing but for Nu_skeleton + // pos = len_free13 + 1; + g_start13[g,1] = len_free[13] + 1; + g_start13[g,2] = pos[13]; + for (i in 1:(p+q)) { + if (is_inf(Nu_skeleton[g,i,1])) { + if (w13skel[pos[13],2] == 0 || w13skel[pos[13],3] == 1) len_free[13] += 1; + pos[13] += 1; + } + } + // same thing but for Alpha_skeleton + g_start14[g,1] = len_free[14] + 1; + g_start14[g,2] = pos[14]; + for (i in 1:(m+n)) { + if (is_inf(Alpha_skeleton[g,i,1])) { + if (w14skel[pos[14],2] == 0 || w14skel[pos[14],3] == 1) len_free[14] += 1; + pos[14] += 1; + } + } + } + // same thing but for Tau_skeleton + g_start15[g,1] = len_free[15] + 1; + g_start15[g,2] = pos[15]; + for (i in 1:(sum(nlevs) - Nord)) { + if (is_inf(Tau_skeleton[g,i,1])) { + if (w15skel[pos[15],2] == 0 || w15skel[pos[15],3] == 1) len_free[15] += 1; + pos[15] += 1; + } + } + // now level 2 + // count free elements in Lambda_y_skeleton + g_start1_c[g,1] = len_free_c[1] + 1; + g_start1_c[g,2] = pos_c[1]; + for (i in 1:p_c) { + for (j in 1:m_c) { + if (is_inf(Lambda_y_skeleton_c[g,i,j])) { + if (w1skel_c[pos_c[1],2] == 0 || w1skel_c[pos_c[1],3] == 1) len_free_c[1] += 1; + pos_c[1] += 1; + } + } + } + // same thing but for B_skeleton + g_start4_c[g,1] = len_free_c[4] + 1; + g_start4_c[g,2] = pos_c[4]; + for (i in 1:m_c) { + for (j in 1:m_c) { + if (is_inf(B_skeleton_c[g,i,j])) { + if (w4skel_c[pos_c[4],2] == 0 || w4skel_c[pos_c[4],3] == 1) len_free_c[4] += 1; + pos_c[4] += 1; + } + } + } + // same thing but for Theta_skeleton + g_start5_c[g,1] = len_free_c[5] + 1; + g_start5_c[g,2] = pos_c[5]; + for (i in 1:p_c) { + if (is_inf(Theta_skeleton_c[g,i,i])) { + if (w5skel_c[pos_c[5],2] == 0 || w5skel_c[pos_c[5],3] == 1) len_free_c[5] += 1; + pos_c[5] += 1; + } + } + // same thing but for Theta_r_skeleton + g_start7_c[g,1] = len_free_c[7] + 1; + g_start7_c[g,2] = pos_c[7]; + for (i in 1:(p_c-1)) { + for (j in (i+1):p_c) { + if (is_inf(Theta_r_skeleton_c[g,j,i])) { + if (w7skel_c[pos_c[7],2] == 0 || w7skel_c[pos_c[7],3] == 1) len_free_c[7] += 1; + pos_c[7] += 1; + } + if (is_inf(Theta_r_skeleton_f_c[g,j,i])) { + if (w8skel_c[pos_c[8],2] == 0 || w8skel_c[pos_c[8],3] == 1) len_free_c[8] += 1; + pos_c[8] += 1; + } + } + } + // same thing but for Psi_skeleton + g_start9_c[g,1] = len_free_c[9] + 1; + g_start9_c[g,2] = pos_c[9]; + for (i in 1:m_c) { + if (is_inf(Psi_skeleton_c[g,i,i])) { + if (w9skel_c[pos_c[9],2] == 0 || w9skel_c[pos_c[9],3] == 1) len_free_c[9] += 1; + pos_c[9] += 1; + } + } + // same thing but for Psi_r_skeleton + g_start10_c[g,1] = len_free_c[10] + 1; + g_start10_c[g,2] = pos_c[10]; + for (i in 1:(m_c-1)) { + for (j in (i+1):m_c) { + if (is_inf(Psi_r_skeleton_c[g,j,i])) { + if (w10skel_c[pos_c[10],2] == 0 || w10skel_c[pos_c[10],3] == 1) len_free_c[10] += 1; + pos_c[10] += 1; + } + if (is_inf(Psi_r_skeleton_f_c[g,j,i])) { + if (w11skel_c[pos_c[11],2] == 0 || w11skel_c[pos_c[11],3] == 1) len_free_c[11] += 1; + pos_c[11] += 1; + } + } + } + // same thing but for Nu_skeleton + // pos = len_free13 + 1; + g_start13_c[g,1] = len_free_c[13] + 1; + g_start13_c[g,2] = pos_c[13]; + for (i in 1:p_c) { + if (is_inf(Nu_skeleton_c[g,i,1])) { + if (w13skel_c[pos_c[13],2] == 0 || w13skel_c[pos_c[13],3] == 1) len_free_c[13] += 1; + pos_c[13] += 1; + } + } + // same thing but for Alpha_skeleton + g_start14_c[g,1] = len_free_c[14] + 1; + g_start14_c[g,2] = pos_c[14]; + for (i in 1:m_c) { + if (is_inf(Alpha_skeleton_c[g,i,1])) { + if (w14skel_c[pos_c[14],2] == 0 || w14skel_c[pos_c[14],3] == 1) len_free_c[14] += 1; + pos_c[14] += 1; + } + } + } + // for clusterwise loglik computations + if (max(nclus[,2]) > 1) for (i in 1:max(nclus[,2])) intone[i] = 1; + if (!ord && (use_suff || use_cov)) { + // sufficient stat matrices by pattern, moved to left for missing + for (patt in 1:Np) { + Sstar[patt] = rep_matrix(0, p + q, p + q); + Sstar[patt, 1:Nobs[patt], 1:Nobs[patt]] = S[patt, Obsvar[patt, 1:Nobs[patt]], Obsvar[patt, 1:Nobs[patt]]]; + for (j in 1:Nobs[patt]) { + YXbarstar[patt,j] = YXbar[patt, Obsvar[patt,j]]; + } + } + } +} +parameters { + // free elements (possibly with inequality constraints) for coefficient matrices + vector[len_free[1]] Lambda_y_free; + vector[len_free[4]] B_free; + vector[len_free[5]] Theta_sd_free; + array[thetanblk[1]] corr_matrix[thetadims[1]] Theta_r_mat_1; + array[thetanblk[2]] corr_matrix[thetadims[2]] Theta_r_mat_2; + array[thetanblk[3]] corr_matrix[thetadims[3]] Theta_r_mat_3; + array[thetanblk[4]] corr_matrix[thetadims[4]] Theta_r_mat_4; + array[thetanblk[5]] corr_matrix[thetadims[5]] Theta_r_mat_5; + vector[len_free[7]] Theta_r_free; // to use beta prior + vector[len_free[9]] Psi_sd_free; + array[psinblk[1]] corr_matrix[psidims[1]] Psi_r_mat_1; + array[psinblk[2]] corr_matrix[psidims[2]] Psi_r_mat_2; + array[psinblk[3]] corr_matrix[psidims[3]] Psi_r_mat_3; + array[psinblk[4]] corr_matrix[psidims[4]] Psi_r_mat_4; + array[psinblk[5]] corr_matrix[psidims[5]] Psi_r_mat_5; + vector[len_free[10]] Psi_r_free; + vector[len_free[13]] Nu_free; + vector[len_free[14]] Alpha_free; + vector[len_free[15]] Tau_ufree; + vector[npriblks] sd0; // shrink_t parameters + vector[Noent] z_aug; //augmented ordinal data + vector[len_free_c[1]] Lambda_y_free_c; + vector[len_free_c[4]] B_free_c; + vector[len_free_c[5]] Theta_sd_free_c; + array[thetanblk_c[1]] corr_matrix[thetadims_c[1]] Theta_r_mat_1_c; + array[thetanblk_c[2]] corr_matrix[thetadims_c[2]] Theta_r_mat_2_c; + array[thetanblk_c[3]] corr_matrix[thetadims_c[3]] Theta_r_mat_3_c; + array[thetanblk_c[4]] corr_matrix[thetadims_c[4]] Theta_r_mat_4_c; + array[thetanblk_c[5]] corr_matrix[thetadims_c[5]] Theta_r_mat_5_c; + vector[len_free_c[7]] Theta_r_free_c; // to use beta prior + vector[len_free_c[9]] Psi_sd_free_c; + array[psinblk_c[1]] corr_matrix[psidims_c[1]] Psi_r_mat_1_c; + array[psinblk_c[2]] corr_matrix[psidims_c[2]] Psi_r_mat_2_c; + array[psinblk_c[3]] corr_matrix[psidims_c[3]] Psi_r_mat_3_c; + array[psinblk_c[4]] corr_matrix[psidims_c[4]] Psi_r_mat_4_c; + array[psinblk_c[5]] corr_matrix[psidims_c[5]] Psi_r_mat_5_c; + vector[len_free_c[10]] Psi_r_free_c; + vector[len_free_c[13]] Nu_free_c; + vector[len_free_c[14]] Alpha_free_c; +} +transformed parameters { + array[Ng] matrix[p, m] Lambda_y; + array[Ng] matrix[m, m] B; + array[Ng] matrix[p, p] Theta_sd; + array[Ng] matrix[p, p] T_r_lower; + array[Ng] matrix[p, p] Theta_r; + array[Ng] matrix[p + q, 1] Nu; + array[Ng] matrix[m + n, 1] Alpha; + array[Ng] matrix[p_c, m_c] Lambda_y_c; + array[Ng] matrix[m_c, m_c] B_c; + array[Ng] matrix[p_c, p_c] Theta_sd_c; + array[Ng] matrix[p_c, p_c] T_r_lower_c; + array[Ng] matrix[p_c, p_c] Theta_r_c; + array[Ng] matrix[p_c, 1] Nu_c; + array[Ng] matrix[m_c, 1] Alpha_c; + array[Ng] matrix[sum(nlevs) - Nord, 1] Tau_un; + array[Ng] matrix[sum(nlevs) - Nord, 1] Tau; + vector[len_free[15]] Tau_free; + real tau_jacobian; + array[Ng] matrix[m, m] Psi; + array[Ng] matrix[m, m] Psi_sd; + array[Ng] matrix[m, m] Psi_r_lower; + array[Ng] matrix[m, m] Psi_r; + array[Ng] matrix[m_c, m_c] Psi_c; + array[Ng] matrix[m_c, m_c] Psi_sd_c; + array[Ng] matrix[m_c, m_c] Psi_r_lower_c; + array[Ng] matrix[m_c, m_c] Psi_r_c; + vector[len_free[1]] lambda_y_primn; + vector[len_free[4]] b_primn; + vector[len_free[13]] nu_primn; + vector[len_free[14]] alpha_primn; + vector[len_free[15]] tau_primn; + vector[len_free_c[1]] lambda_y_primn_c; + vector[len_free_c[4]] b_primn_c; + vector[len_free_c[13]] nu_primn_c; + vector[len_free_c[14]] alpha_primn_c; + array[Ng] matrix[p, m] Lambda_y_A; // = Lambda_y * (I - B)^{-1} + array[Ng] matrix[p_c, m_c] Lambda_y_A_c; + array[Ng] vector[p + q] Mu; + array[Ng] matrix[p + q, p + q] Sigma; // model covariance matrix + array[Ng] matrix[p + q, p + q] Sigmainv_grp; // model covariance matrix + array[Ng] real logdetSigma_grp; + array[Np] matrix[p + q + 1, p + q + 1] Sigmainv; // for updating S^-1 by missing data pattern + array[Ng] vector[p_c] Mu_c; + array[Ng] matrix[p_c, p_c] Sigma_c; // level 2 model covariance matrix + array[Ng] matrix[N_both + N_within, N_both + N_within] S_PW; + array[Ntot] vector[p + q] YXstar; + array[Ntot] vector[Nord] YXostar; // ordinal data + for (g in 1:Ng) { + // model matrices + Lambda_y[g] = fill_matrix(Lambda_y_free, Lambda_y_skeleton[g], w1skel, g_start1[g,1], g_start1[g,2]); + B[g] = fill_matrix(B_free, B_skeleton[g], w4skel, g_start4[g,1], g_start4[g,2]); + Theta_sd[g] = fill_matrix(Theta_sd_free, Theta_skeleton[g], w5skel, g_start5[g,1], g_start5[g,2]); + T_r_lower[g] = fill_matrix(Theta_r_free, Theta_r_skeleton[g], w7skel, g_start7[g,1], g_start7[g,2]); + Theta_r[g] = T_r_lower[g] + transpose(T_r_lower[g]) - diag_matrix(rep_vector(1, p)); + if (!use_cov) { + Nu[g] = fill_matrix(Nu_free, Nu_skeleton[g], w13skel, g_start13[g,1], g_start13[g,2]); + Alpha[g] = fill_matrix(Alpha_free, Alpha_skeleton[g], w14skel, g_start14[g,1], g_start14[g,2]); + } + Psi[g] = diag_matrix(rep_vector(0, m)); + if (m > 0) { + Psi_sd[g] = fill_matrix(Psi_sd_free, Psi_skeleton[g], w9skel, g_start9[g,1], g_start9[g,2]); + Psi_r_lower[g] = fill_matrix(Psi_r_free, Psi_r_skeleton[g], w10skel, g_start10[g,1], g_start10[g,2]); + Psi_r[g] = Psi_r_lower[g] + transpose(Psi_r_lower[g]) - diag_matrix(rep_vector(1, m)); + } + // level 2 matrices + Lambda_y_c[g] = fill_matrix(Lambda_y_free_c, Lambda_y_skeleton_c[g], w1skel_c, g_start1_c[g,1], g_start1_c[g,2]); + B_c[g] = fill_matrix(B_free_c, B_skeleton_c[g], w4skel_c, g_start4_c[g,1], g_start4_c[g,2]); + Theta_sd_c[g] = fill_matrix(Theta_sd_free_c, Theta_skeleton_c[g], w5skel_c, g_start5_c[g,1], g_start5_c[g,2]); + T_r_lower_c[g] = fill_matrix(Theta_r_free_c, Theta_r_skeleton_c[g], w7skel_c, g_start7_c[g,1], g_start7_c[g,2]); + Theta_r_c[g] = T_r_lower_c[g] + transpose(T_r_lower_c[g]) - diag_matrix(rep_vector(1, p_c)); + Nu_c[g] = fill_matrix(Nu_free_c, Nu_skeleton_c[g], w13skel_c, g_start13_c[g,1], g_start13_c[g,2]); + Alpha_c[g] = fill_matrix(Alpha_free_c, Alpha_skeleton_c[g], w14skel_c, g_start14_c[g,1], g_start14_c[g,2]); + Psi_c[g] = diag_matrix(rep_vector(0, m_c)); + if (m_c > 0) { + Psi_sd_c[g] = fill_matrix(Psi_sd_free_c, Psi_skeleton_c[g], w9skel_c, g_start9_c[g,1], g_start9_c[g,2]); + Psi_r_lower_c[g] = fill_matrix(Psi_r_free_c, Psi_r_skeleton_c[g], w10skel_c, g_start10_c[g,1], g_start10_c[g,2]); + Psi_r_c[g] = Psi_r_lower_c[g] + transpose(Psi_r_lower_c[g]) - diag_matrix(rep_vector(1, m_c)); + } + } + if (sum(thetanblk) > 0) { + // we need to define a separate parameter for each dimension of correlation matrix, + // so we need all these Theta_r_mats + Theta_r = fill_cov(Theta_r, thetablkse, thetanblk, Theta_r_mat_1, Theta_r_mat_2, Theta_r_mat_3, Theta_r_mat_4, Theta_r_mat_5, thetaorder, thetarevord); + } + if (sum(thetanblk_c) > 0) { + // we need to define a separate parameter for each dimension of correlation matrix, + // so we need all these Theta_r_mats + Theta_r_c = fill_cov(Theta_r_c, thetablkse_c, thetanblk_c, Theta_r_mat_1_c, Theta_r_mat_2_c, Theta_r_mat_3_c, Theta_r_mat_4_c, Theta_r_mat_5_c, thetaorder_c, thetarevord_c); + } + if (sum(psinblk) > 0) { + Psi_r = fill_cov(Psi_r, psiblkse, psinblk, Psi_r_mat_1, Psi_r_mat_2, Psi_r_mat_3, Psi_r_mat_4, Psi_r_mat_5, psiorder, psirevord); + } + if (sum(psinblk_c) > 0) { + Psi_r_c = fill_cov(Psi_r_c, psiblkse_c, psinblk_c, Psi_r_mat_1_c, Psi_r_mat_2_c, Psi_r_mat_3_c, Psi_r_mat_4_c, Psi_r_mat_5_c, psiorder_c, psirevord_c); + } + // see p. 3 https://books.google.com/books?id=9AC-s50RjacC + for (g in 1:Ng) { + if (m > 0) { + Lambda_y_A[g] = mdivide_right(Lambda_y[g], I - B[g]); // = Lambda_y * (I - B)^{-1} + Psi[g] = quad_form_sym(Psi_r[g], Psi_sd[g]); + } + if (!use_cov) { + Mu[g] = to_vector(Nu[g]); + } else if(has_data) { + Mu[g] = YXbar[g]; // doesn't enter in likelihood, just for lppd + loo + } + if (p > 0) { + Sigma[g, 1:p, 1:p] = quad_form_sym(Theta_r[g], Theta_sd[g]); + if (m > 0) { + Sigma[g, 1:p, 1:p] += quad_form_sym(Psi[g], transpose(Lambda_y_A[g])); + if (!use_cov) Mu[g, 1:p] += to_vector(Lambda_y_A[g] * Alpha[g, 1:m, 1]); + } + } + if (m_c > 0) { + Lambda_y_A_c[g] = mdivide_right(Lambda_y_c[g], I_c - B_c[g]); + Psi_c[g] = quad_form_sym(Psi_r_c[g], Psi_sd_c[g]); + } + Mu_c[g] = to_vector(Nu_c[g]); + if (p_c > 0) { + Sigma_c[g, 1:p_c, 1:p_c] = quad_form_sym(Theta_r_c[g], Theta_sd_c[g]); + if (m_c > 0) { + Sigma_c[g, 1:p_c, 1:p_c] += quad_form_sym(Psi_c[g], transpose(Lambda_y_A_c[g])); + Mu_c[g, 1:p_c] += to_vector(Lambda_y_A_c[g] * Alpha_c[g, 1:m_c, 1]); + } + } + if (nclus[g,2] > 1) { + // remove between variables, for likelihood computations + S_PW[g] = cov_w[g, between_idx[(N_between + 1):p_tilde], between_idx[(N_between + 1):p_tilde]]; + } + } + // obtain ordered thresholds; NB untouched for two-level models + if (ord) { + int opos = 1; + int ofreepos = 1; + tau_jacobian = 0; + for (g in 1:Ng) { + int vecpos = 1; + Tau_un[g] = fill_matrix(Tau_ufree, Tau_skeleton[g], w15skel, g_start15[g,1], g_start15[g,2]); + for (i in 1:Nord) { + for (j in 1:(nlevs[i] - 1)) { + real rc = Tau_skeleton[g, vecpos, 1]; + int eq = w15skel[opos, 1]; + int wig = w15skel[opos, 3]; + if (is_inf(rc)) { + if (eq == 0 || wig == 1) { + if (j == 1) { + Tau[g, vecpos, 1] = Tau_un[g, vecpos, 1]; + } else { + Tau[g, vecpos, 1] = Tau[g, (vecpos - 1), 1] + exp(Tau_un[g, vecpos, 1]); + } + Tau_free[ofreepos] = Tau[g, vecpos, 1]; + // this is used if a prior goes on Tau_free, instead of Tau_ufree: + //if (j > 1) { + // tau_jacobian += Tau_un[g, vecpos, 1]; // see https://mc-stan.org/docs/2_24/reference-manual/ordered-vector.html + // } + ofreepos += 1; + } else if (eq == 1) { + int eqent = w15skel[opos, 2]; + Tau[g, vecpos, 1] = Tau_free[eqent]; + } + opos += 1; + } else { + // fixed value + Tau[g, vecpos, 1] = Tau_un[g, vecpos, 1]; + } + vecpos += 1; + } + } + } + } + // prior vectors + if (wigind) { + lambda_y_primn = fill_prior(Lambda_y_free, lambda_y_mn, w1skel); + b_primn = fill_prior(B_free, b_mn, w4skel); + nu_primn = fill_prior(Nu_free, nu_mn, w13skel); + alpha_primn = fill_prior(Alpha_free, alpha_mn, w14skel); + tau_primn = fill_prior(Tau_ufree, tau_mn, w15skel); + lambda_y_primn_c = fill_prior(Lambda_y_free_c, lambda_y_mn_c, w1skel_c); + b_primn_c = fill_prior(B_free_c, b_mn_c, w4skel_c); + nu_primn_c = fill_prior(Nu_free_c, nu_mn_c, w13skel_c); + alpha_primn_c = fill_prior(Alpha_free_c, alpha_mn_c, w14skel_c); + } else { + lambda_y_primn = to_vector(lambda_y_mn); + b_primn = to_vector(b_mn); + nu_primn = to_vector(nu_mn); + alpha_primn = to_vector(alpha_mn); + tau_primn = to_vector(tau_mn); + lambda_y_primn_c = to_vector(lambda_y_mn_c); + b_primn_c = to_vector(b_mn_c); + nu_primn_c = to_vector(nu_mn_c); + alpha_primn_c = to_vector(alpha_mn_c); + } + // NB nothing below this will be used for two level, because we need other tricks to + // compute the likelihood + // continuous responses underlying ordinal data + if (ord) { + int idxvec = 0; + for (patt in 1:Np) { + for (i in startrow[patt]:endrow[patt]) { + for (j in 1:Nordobs[patt]) { + int obspos = OrdObsvar[patt,j]; + int vecpos = YXo[i,obspos] - 1; + idxvec += 1; + if (obspos > 1) vecpos += sum(nlevs[1:(obspos - 1)]) - (obspos - 1); + if (YXo[i,obspos] == 1) { + YXostar[i,obspos] = -30 + (Tau[grpnum[patt], (vecpos + 1), 1] + 30) .* z_aug[idxvec]; + tau_jacobian += log(abs(Tau[grpnum[patt], (vecpos + 1), 1] + 30)); // must add log(U) to tau_jacobian + } else if (YXo[i,obspos] == nlevs[obspos]) { + YXostar[i,obspos] = Tau[grpnum[patt], vecpos, 1] + (30 - Tau[grpnum[patt], vecpos, 1]) .* z_aug[idxvec]; + tau_jacobian += log(abs(30 - Tau[grpnum[patt], vecpos, 1])); + } else { + YXostar[i,obspos] = Tau[grpnum[patt], vecpos, 1] + (Tau[grpnum[patt], (vecpos + 1), 1] - Tau[grpnum[patt], vecpos, 1]) .* z_aug[idxvec]; + tau_jacobian += Tau_un[grpnum[patt], (vecpos + 1), 1]; // jacobian is log(exp(Tau_un)) + } + YXstar[i, ordidx[obspos]] = YXostar[i, obspos]; + } + } + } + } + if (Ncont > 0) { + for (patt in 1:Np) { + for (i in startrow[patt]:endrow[patt]) { + for (j in 1:Ncont) { + YXstar[i, contidx[j]] = YX[i,j]; + } + } + } + } + // move observations to the left + if (missing) { + for (patt in 1:Np) { + for (i in startrow[patt]:endrow[patt]) { + for (j in 1:Nobs[patt]) { + YXstar[i,j] = YXstar[i, Obsvar[patt,j]]; + } + } + } + } + // for computing mvn with sufficient stats + if (!multilev) { + for (g in 1:Ng) { + Sigmainv_grp[g] = inverse_spd(Sigma[g]); + logdetSigma_grp[g] = log_determinant(Sigma[g]); + } + for (patt in 1:Np) { + Sigmainv[patt, 1:(Nobs[patt] + 1), 1:(Nobs[patt] + 1)] = sig_inv_update(Sigmainv_grp[grpnum[patt]], Obsvar[patt,], Nobs[patt], p + q, logdetSigma_grp[grpnum[patt]]); + } + } +} +model { // N.B.: things declared in the model block do not get saved in the output, which is okay here + /* transformed sd parameters for priors */ + vector[len_free[5]] Theta_pri; + vector[len_free[9]] Psi_pri; + vector[len_free_c[5]] Theta_pri_c; + vector[len_free_c[9]] Psi_pri_c; + /* log-likelihood */ + if (multilev && has_data) { + int grpidx; + int r1 = 1; // index clusters per group + int r2 = 0; + int rr1 = 1; // index units per group + int rr2 = 0; + int r3 = 1; // index unique cluster sizes per group + int r4 = 0; + for (mm in 1:Np) { + grpidx = grpnum[mm]; + if (grpidx > 1) { + r1 += nclus[(grpidx - 1), 2]; + rr1 += nclus[(grpidx - 1), 1]; + r3 += ncluster_sizes[(grpidx - 1)]; + } + r2 += nclus[grpidx, 2]; + rr2 += nclus[grpidx, 1]; + r4 += ncluster_sizes[grpidx]; + target += twolevel_logdens(mean_d[r3:r4], cov_d[r3:r4], S_PW[grpidx], YX[rr1:rr2], + nclus[grpidx,], cluster_size[r1:r2], cluster_sizes[r3:r4], + ncluster_sizes[grpidx], cluster_size_ns[r3:r4], Mu[grpidx], + Sigma[grpidx], Mu_c[grpidx], Sigma_c[grpidx], + ov_idx1, ov_idx2, within_idx, between_idx, both_idx, + p_tilde, N_within, N_between, N_both); + if (Nx[grpidx] + Nx_between[grpidx] > 0) target += -log_lik_x; + } + } else if (use_cov && !pri_only) { + for (g in 1:Ng) { + target += wishart_lpdf((N[g] - 1) * Sstar[g] | N[g] - 1, Sigma[g]); + if (Nx[g] > 0) { + array[Nx[g]] int xvars = Xdatvar[g, 1:Nx[g]]; + target += -wishart_lpdf((N[g] - 1) * Sstar[g, xvars, xvars] | N[g] - 1, Sigma[g, xvars, xvars]); + } + } + } else if (has_data && !pri_only) { + array[p + q] int obsidx; + array[p + q] int xidx; + array[p + q] int xdatidx; + int grpidx; + int r1; + int r2; + for (mm in 1:Np) { + obsidx = Obsvar[mm,]; + xidx = Xvar[mm,]; + xdatidx = Xdatvar[mm,]; + grpidx = grpnum[mm]; + r1 = startrow[mm]; + r2 = endrow[mm]; + if (!use_suff) { + target += multi_normal_lpdf(YXstar[r1:r2,1:Nobs[mm]] | Mu[grpidx, obsidx[1:Nobs[mm]]], Sigma[grpidx, obsidx[1:Nobs[mm]], obsidx[1:Nobs[mm]]]); + if (Nx[mm] > 0) { + target += -1.0 * multi_normal_lpdf(YXstar[r1:r2,xdatidx[1:Nx[mm]]] | Mu[grpidx, xidx[1:Nx[mm]]], Sigma[grpidx, xidx[1:Nx[mm]], xidx[1:Nx[mm]]]); + } + } else { + // sufficient stats + target += multi_normal_suff(YXbarstar[mm, 1:Nobs[mm]], Sstar[mm, 1:Nobs[mm], 1:Nobs[mm]], Mu[grpidx, obsidx[1:Nobs[mm]]], Sigmainv[mm, 1:(Nobs[mm] + 1), 1:(Nobs[mm] + 1)], r2 - r1 + 1); + if (Nx[mm] > 0 && !missing) { + target += -1.0 * multi_normal_suff(YXbarstar[mm, xdatidx[1:Nx[mm]]], Sstar[mm, xdatidx[1:Nx[mm]], xdatidx[1:Nx[mm]]], Mu[grpidx, xidx[1:Nx[mm]]], sig_inv_update(Sigmainv[mm], xidx, Nx[mm], p + q, logdetSigma_grp[grpidx]), r2 - r1 + 1); + } else if (Nx[mm] > 0) { + target += -1.0 * multi_normal_lpdf(YXstar[r1:r2,xdatidx[1:Nx[mm]]] | Mu[grpidx, xidx[1:Nx[mm]]], Sigma[grpidx, xidx[1:Nx[mm]], xidx[1:Nx[mm]]]); + } + } + } + if (ord) { + target += tau_jacobian; + } + } + /* prior densities in log-units, first for unbounded parameters that could have shrinkage priors */ + target += eval_priors(Lambda_y_free, B_free, Nu_free, Alpha_free, sd0, lambda_y_primn, lambda_y_sd, lambda_y_pri, lambda_y_blk, b_primn, b_sd, b_pri, b_blk, nu_primn, nu_sd, nu_pri, nu_blk, alpha_primn, alpha_sd, alpha_pri, alpha_blk, len_free, blkparm1, blkparm2, priblklen, npriblks); + target += normal_lpdf(Tau_ufree | tau_primn, tau_sd); + target += normal_lpdf(Lambda_y_free_c | lambda_y_primn_c, lambda_y_sd_c); + target += normal_lpdf(B_free_c | b_primn_c, b_sd_c); + target += normal_lpdf(Nu_free_c | nu_primn_c, nu_sd_c); + target += normal_lpdf(Alpha_free_c | alpha_primn_c, alpha_sd_c); + /* transform sd parameters to var or prec, depending on + what the user wants. */ + Theta_pri = Theta_sd_free; + if (len_free[5] > 0 && theta_pow != 1) { + for (i in 1:len_free[5]) { + Theta_pri[i] = Theta_sd_free[i]^(theta_pow); + target += log(abs(theta_pow)) + (theta_pow - 1)*log(Theta_sd_free[i]); + } + } + Psi_pri = Psi_sd_free; + if (len_free[9] > 0 && psi_pow != 1) { + for (i in 1:len_free[9]) { + Psi_pri[i] = Psi_sd_free[i]^(psi_pow); + target += log(abs(psi_pow)) + (psi_pow - 1)*log(Psi_sd_free[i]); + } + } + target += gamma_lpdf(Theta_pri | theta_sd_shape, theta_sd_rate); + target += gamma_lpdf(Psi_pri | psi_sd_shape, psi_sd_rate); + target += beta_lpdf(.5 * (1 + Theta_r_free) | theta_r_alpha, theta_r_beta) + log(.5) * len_free[7]; // the latter term is the jacobian moving from (-1,1) to (0,1), because beta_lpdf is defined on (0,1) + if (sum(psinblk) > 0) { + for (k in 1:sum(psinblk)) { + int blkidx = psiblkse[k, 6]; + int arrayidx = psiblkse[k, 5]; + if (arrayidx == 1) { + target += lkj_corr_lpdf(Psi_r_mat_1[blkidx] | psiblkpri[k]); + } else if (arrayidx == 2) { + target += lkj_corr_lpdf(Psi_r_mat_2[blkidx] | psiblkpri[k]); + } else if (arrayidx == 3) { + target += lkj_corr_lpdf(Psi_r_mat_3[blkidx] | psiblkpri[k]); + } else if (arrayidx == 4) { + target += lkj_corr_lpdf(Psi_r_mat_4[blkidx] | psiblkpri[k]); + } else { + target += lkj_corr_lpdf(Psi_r_mat_5[blkidx] | psiblkpri[k]); + } + } + } + if (len_free[10] > 0) { + target += beta_lpdf(.5 * (1 + Psi_r_free) | psi_r_alpha, psi_r_beta) + log(.5) * len_free[10]; + } + if (sum(thetanblk) > 0) { + for (k in 1:sum(thetanblk)) { + int blkidx = thetablkse[k, 6]; + int arrayidx = thetablkse[k, 5]; + if (arrayidx == 1) { + target += lkj_corr_lpdf(Theta_r_mat_1[blkidx] | thetablkpri[k]); + } else if (arrayidx == 2) { + target += lkj_corr_lpdf(Theta_r_mat_2[blkidx] | thetablkpri[k]); + } else if (arrayidx == 3) { + target += lkj_corr_lpdf(Theta_r_mat_3[blkidx] | thetablkpri[k]); + } else if (arrayidx == 4) { + target += lkj_corr_lpdf(Theta_r_mat_4[blkidx] | thetablkpri[k]); + } else { + target += lkj_corr_lpdf(Theta_r_mat_5[blkidx] | thetablkpri[k]); + } + } + } + // and the same for level 2 + Theta_pri_c = Theta_sd_free_c; + if (len_free_c[5] > 0 && theta_pow_c != 1) { + for (i in 1:len_free_c[5]) { + Theta_pri_c[i] = Theta_sd_free_c[i]^(theta_pow_c); + target += log(abs(theta_pow_c)) + (theta_pow_c - 1)*log(Theta_sd_free_c[i]); + } + } + Psi_pri_c = Psi_sd_free_c; + if (len_free_c[9] > 0 && psi_pow_c != 1) { + for (i in 1:len_free_c[9]) { + Psi_pri_c[i] = Psi_sd_free_c[i]^(psi_pow_c); + target += log(abs(psi_pow_c)) + (psi_pow_c - 1)*log(Psi_sd_free_c[i]); + } + } + target += gamma_lpdf(Theta_pri_c | theta_sd_shape_c, theta_sd_rate_c); + target += gamma_lpdf(Psi_pri_c | psi_sd_shape_c, psi_sd_rate_c); + if (sum(psinblk_c) > 0) { + for (k in 1:sum(psinblk_c)) { + int blkidx = psiblkse_c[k, 6]; + int arrayidx = psiblkse_c[k, 5]; + if (arrayidx == 1) { + target += lkj_corr_lpdf(Psi_r_mat_1_c[blkidx] | psiblkpri_c[k]); + } else if (arrayidx == 2) { + target += lkj_corr_lpdf(Psi_r_mat_2_c[blkidx] | psiblkpri_c[k]); + } else if (arrayidx == 3) { + target += lkj_corr_lpdf(Psi_r_mat_3_c[blkidx] | psiblkpri_c[k]); + } else if (arrayidx == 4) { + target += lkj_corr_lpdf(Psi_r_mat_4_c[blkidx] | psiblkpri_c[k]); + } else { + target += lkj_corr_lpdf(Psi_r_mat_5_c[blkidx] | psiblkpri_c[k]); + } + } + } else if (len_free_c[10] > 0) { + target += beta_lpdf(.5 * (1 + Psi_r_free_c) | psi_r_alpha_c, psi_r_beta_c) + log(.5) * len_free_c[10]; + } + if (sum(thetanblk_c) > 0) { + for (k in 1:sum(thetanblk_c)) { + int blkidx = thetablkse_c[k, 6]; + int arrayidx = thetablkse_c[k, 5]; + if (arrayidx == 1) { + target += lkj_corr_lpdf(Theta_r_mat_1_c[blkidx] | thetablkpri_c[k]); + } else if (arrayidx == 2) { + target += lkj_corr_lpdf(Theta_r_mat_2_c[blkidx] | thetablkpri_c[k]); + } else if (arrayidx == 3) { + target += lkj_corr_lpdf(Theta_r_mat_3_c[blkidx] | thetablkpri_c[k]); + } else if (arrayidx == 4) { + target += lkj_corr_lpdf(Theta_r_mat_4_c[blkidx] | thetablkpri_c[k]); + } else { + target += lkj_corr_lpdf(Theta_r_mat_5_c[blkidx] | thetablkpri_c[k]); + } + } + } else if (len_free_c[7] > 0) { + target += beta_lpdf(.5 * (1 + Theta_r_free_c) | theta_r_alpha_c, theta_r_beta_c) + log(.5) * len_free_c[7]; + } +} +generated quantities { // these matrices are saved in the output but do not figure into the likelihood + // see p. 34 https://books.google.com/books?id=9AC-s50RjacC + // try to compute the log priors: + real log_prior; + log_prior = 0; + log_prior += + eval_priors( + Lambda_y_free, + B_free, + Nu_free, + Alpha_free, + sd0, + lambda_y_primn, + lambda_y_sd, + lambda_y_pri, + lambda_y_blk, + b_primn, + b_sd, + b_pri, + b_blk, + nu_primn, + nu_sd, + nu_pri, + nu_blk, + alpha_primn, + alpha_sd, + alpha_pri, + alpha_blk, + len_free, + blkparm1, + blkparm2, + priblklen, + npriblks + ); + + vector[len_free[5]] Theta_pri; + vector[len_free[9]] Psi_pri; + vector[len_free_c[5]] Theta_pri_c; + vector[len_free_c[9]] Psi_pri_c; + + log_prior += normal_lpdf(Tau_ufree | tau_primn, tau_sd); + log_prior += normal_lpdf(Lambda_y_free_c | lambda_y_primn_c, lambda_y_sd_c); + log_prior += normal_lpdf(B_free_c | b_primn_c, b_sd_c); + log_prior += normal_lpdf(Nu_free_c | nu_primn_c, nu_sd_c); + log_prior += normal_lpdf(Alpha_free_c | alpha_primn_c, alpha_sd_c); + + + Theta_pri = Theta_sd_free; + if (len_free[5] > 0 && theta_pow != 1) { + for (i in 1:len_free[5]) { + Theta_pri[i] = Theta_sd_free[i]^(theta_pow); + log_prior += log(abs(theta_pow)) + (theta_pow - 1)*log(Theta_sd_free[i]); + } + } + + Psi_pri = Psi_sd_free; + if (len_free[9] > 0 && psi_pow != 1) { + for (i in 1:len_free[9]) { + Psi_pri[i] = Psi_sd_free[i]^(psi_pow); + log_prior += log(abs(psi_pow)) + (psi_pow - 1)*log(Psi_sd_free[i]); + } + } + log_prior += gamma_lpdf(Theta_pri | theta_sd_shape, theta_sd_rate); + log_prior += gamma_lpdf(Psi_pri | psi_sd_shape, psi_sd_rate); + log_prior += beta_lpdf(.5 * (1 + Theta_r_free) | theta_r_alpha, theta_r_beta) + log(.5) * len_free[7]; // the latter term is the + // jacobian moving from (-1,1) to (0,1), because beta_lpdf is defined on (0,1) + + + if (sum(psinblk) > 0) { + for (k in 1:sum(psinblk)) { + int blkidx = psiblkse[k, 6]; + int arrayidx = psiblkse[k, 5]; + if (arrayidx == 1) { + log_prior += lkj_corr_lpdf(Psi_r_mat_1[blkidx] | psiblkpri[k]); + } else if (arrayidx == 2) { + log_prior += lkj_corr_lpdf(Psi_r_mat_2[blkidx] | psiblkpri[k]); + } else if (arrayidx == 3) { + log_prior += lkj_corr_lpdf(Psi_r_mat_3[blkidx] | psiblkpri[k]); + } else if (arrayidx == 4) { + log_prior += lkj_corr_lpdf(Psi_r_mat_4[blkidx] | psiblkpri[k]); + } else { + log_prior += lkj_corr_lpdf(Psi_r_mat_5[blkidx] | psiblkpri[k]); + } + } + } + if (len_free[10] > 0) { + log_prior += beta_lpdf(.5 * (1 + Psi_r_free) | psi_r_alpha, psi_r_beta) + log(.5) * len_free[10]; + } + if (sum(thetanblk) > 0) { + for (k in 1:sum(thetanblk)) { + int blkidx = thetablkse[k, 6]; + int arrayidx = thetablkse[k, 5]; + if (arrayidx == 1) { + log_prior += lkj_corr_lpdf(Theta_r_mat_1[blkidx] | thetablkpri[k]); + } else if (arrayidx == 2) { + log_prior += lkj_corr_lpdf(Theta_r_mat_2[blkidx] | thetablkpri[k]); + } else if (arrayidx == 3) { + log_prior += lkj_corr_lpdf(Theta_r_mat_3[blkidx] | thetablkpri[k]); + } else if (arrayidx == 4) { + log_prior += lkj_corr_lpdf(Theta_r_mat_4[blkidx] | thetablkpri[k]); + } else { + log_prior += lkj_corr_lpdf(Theta_r_mat_5[blkidx] | thetablkpri[k]); + } + } + } + // and the same for level 2 + Theta_pri_c = Theta_sd_free_c; + if (len_free_c[5] > 0 && theta_pow_c != 1) { + for (i in 1:len_free_c[5]) { + Theta_pri_c[i] = Theta_sd_free_c[i]^(theta_pow_c); + log_prior += log(abs(theta_pow_c)) + (theta_pow_c - 1)*log(Theta_sd_free_c[i]); + } + } + Psi_pri_c = Psi_sd_free_c; + if (len_free_c[9] > 0 && psi_pow_c != 1) { + for (i in 1:len_free_c[9]) { + Psi_pri_c[i] = Psi_sd_free_c[i]^(psi_pow_c); + log_prior += log(abs(psi_pow_c)) + (psi_pow_c - 1)*log(Psi_sd_free_c[i]); + } + } + log_prior += gamma_lpdf(Theta_pri_c | theta_sd_shape_c, theta_sd_rate_c); + log_prior += gamma_lpdf(Psi_pri_c | psi_sd_shape_c, psi_sd_rate_c); + if (sum(psinblk_c) > 0) { + for (k in 1:sum(psinblk_c)) { + int blkidx = psiblkse_c[k, 6]; + int arrayidx = psiblkse_c[k, 5]; + if (arrayidx == 1) { + log_prior += lkj_corr_lpdf(Psi_r_mat_1_c[blkidx] | psiblkpri_c[k]); + } else if (arrayidx == 2) { + log_prior += lkj_corr_lpdf(Psi_r_mat_2_c[blkidx] | psiblkpri_c[k]); + } else if (arrayidx == 3) { + log_prior += lkj_corr_lpdf(Psi_r_mat_3_c[blkidx] | psiblkpri_c[k]); + } else if (arrayidx == 4) { + log_prior += lkj_corr_lpdf(Psi_r_mat_4_c[blkidx] | psiblkpri_c[k]); + } else { + log_prior += lkj_corr_lpdf(Psi_r_mat_5_c[blkidx] | psiblkpri_c[k]); + } + } + } else if (len_free_c[10] > 0) { + log_prior += beta_lpdf(.5 * (1 + Psi_r_free_c) | psi_r_alpha_c, psi_r_beta_c) + log(.5) * len_free_c[10]; + } + if (sum(thetanblk_c) > 0) { + for (k in 1:sum(thetanblk_c)) { + int blkidx = thetablkse_c[k, 6]; + int arrayidx = thetablkse_c[k, 5]; + if (arrayidx == 1) { + log_prior += lkj_corr_lpdf(Theta_r_mat_1_c[blkidx] | thetablkpri_c[k]); + } else if (arrayidx == 2) { + log_prior += lkj_corr_lpdf(Theta_r_mat_2_c[blkidx] | thetablkpri_c[k]); + } else if (arrayidx == 3) { + log_prior += lkj_corr_lpdf(Theta_r_mat_3_c[blkidx] | thetablkpri_c[k]); + } else if (arrayidx == 4) { + log_prior += lkj_corr_lpdf(Theta_r_mat_4_c[blkidx] | thetablkpri_c[k]); + } else { + log_prior += lkj_corr_lpdf(Theta_r_mat_5_c[blkidx] | thetablkpri_c[k]); + } + } + } else if (len_free_c[7] > 0) { + log_prior += beta_lpdf(.5 * (1 + Theta_r_free_c) | theta_r_alpha_c, theta_r_beta_c) + log(.5) * len_free_c[7]; + } + + + + // sign constraints and correlations + vector[len_free[1]] ly_sign; + vector[len_free[4]] bet_sign; + vector[len_free[14]] al_sign; + array[Ng] matrix[p, p] Thetmat; + array[Ng] matrix[p, p] Thet; + array[Ng] matrix[m, m] PSmat; + array[Ng] matrix[m, m] PS; + vector[len_free[8]] Theta_cov; + vector[len_free[5]] Theta_var; + vector[len_free[10]] P_r; + vector[len_free[11]] Psi_cov; + vector[len_free[9]] Psi_var; + // level 2 + vector[len_free_c[1]] ly_sign_c; + vector[len_free_c[4]] bet_sign_c; + vector[len_free_c[14]] al_sign_c; + array[Ng] matrix[p_c, p_c] Thetmat_c; + array[Ng] matrix[p_c, p_c] Thet_c; + array[Ng] matrix[m_c, m_c] PSmat_c; + array[Ng] matrix[m_c, m_c] PS_c; + vector[len_free_c[8]] Theta_cov_c; + vector[len_free_c[5]] Theta_var_c; + vector[len_free_c[10]] P_r_c; + vector[len_free_c[11]] Psi_cov_c; + vector[len_free_c[9]] Psi_var_c; + // loglik + ppp + vector[multilev ? sum(nclus[,2]) : (use_cov ? Ng : Ntot)] log_lik; // for loo, etc + vector[multilev ? sum(nclus[,2]) : (use_cov ? Ng : Ntot)] log_lik_sat; // for ppp + array[Ntot] vector[multilev ? p_tilde : p + q] YXstar_rep; // artificial data + vector[multilev ? sum(nclus[,2]) : (use_cov ? Ng : Ntot)] log_lik_rep; // for loo, etc + vector[multilev ? sum(nclus[,2]) : (use_cov ? Ng : Ntot)] log_lik_rep_sat; // for ppp + array[Ng] matrix[p + q, p + q + 1] satout; + array[Ng] matrix[p + q, p + q + 1] satrep_out; + array[Ng] vector[p + q] Mu_sat; + array[Ng] matrix[p + q, p + q] Sigma_sat; + array[Ng] matrix[p + q, p + q] Sigma_sat_inv_grp; + array[Ng] real logdetS_sat_grp; + array[Np] matrix[p + q + 1, p + q + 1] Sigma_sat_inv; + array[Ng] vector[p + q] Mu_rep_sat; + array[Ng] matrix[p + q, p + q] Sigma_rep_sat; + array[Ng] matrix[p + q, p + q] Sigma_rep_sat_inv_grp; + array[Np] matrix[p + q + 1, p + q + 1] Sigma_rep_sat_inv; + array[Ng] real logdetS_rep_sat_grp; + matrix[p + q, p + q] zmat; + array[sum(nclus[,2])] vector[p_tilde] mean_d_rep; + vector[multilev ? sum(nclus[,2]) : Ng] log_lik_x_rep; + array[Ng] matrix[N_both + N_within, N_both + N_within] S_PW_rep; + array[Ng] matrix[p_tilde, p_tilde] S_PW_rep_full; + array[Ng] vector[p_tilde] ov_mean_rep; + array[Ng] vector[p_tilde] xbar_b_rep; + array[Ng] matrix[N_between, N_between] S2_rep; + array[Ng] matrix[p_tilde, p_tilde] S_B_rep; + array[Ng] matrix[p_tilde, p_tilde] cov_b_rep; + real ppp; + // first deal with sign constraints: + ly_sign = sign_constrain_load(Lambda_y_free, len_free[1], lam_y_sign); + bet_sign = sign_constrain_reg(B_free, len_free[4], b_sign, Lambda_y_free, Lambda_y_free); + al_sign = sign_constrain_reg(Alpha_free, len_free[14], alph_sign, Lambda_y_free, rep_vector(1, len_free[1])); + if (len_free[10] > 0) { + P_r = sign_constrain_reg(Psi_r_free, len_free[10], psi_r_sign, Lambda_y_free, Lambda_y_free); + } + ly_sign_c = sign_constrain_load(Lambda_y_free_c, len_free_c[1], lam_y_sign_c); + bet_sign_c = sign_constrain_reg(B_free_c, len_free_c[4], b_sign_c, Lambda_y_free_c, Lambda_y_free_c); + al_sign_c = sign_constrain_reg(Alpha_free_c, len_free_c[14], alph_sign_c, Lambda_y_free_c, rep_vector(1, len_free_c[1])); + if (len_free_c[10] > 0) { + P_r_c = sign_constrain_reg(Psi_r_free_c, len_free_c[10], psi_r_sign_c, Lambda_y_free_c, Lambda_y_free_c); + } + for (g in 1:Ng) { + matrix[p, p] Thtmp = fill_matrix(Theta_r_free, Theta_r_skeleton[g], w7skel, g_start7[g,1], g_start7[g,2]); + Thetmat[g] = Thtmp + transpose(Thtmp) - diag_matrix(rep_vector(1, p)); + if (m > 0) { + matrix[m, m] Pstmp = fill_matrix(P_r, Psi_r_skeleton[g], w10skel, g_start10[g,1], g_start10[g,2]); + PSmat[g] = Pstmp + transpose(Pstmp) - diag_matrix(rep_vector(1, m)); + } + if (p_c > 0) { + matrix[p_c, p_c] Thtmp_c = fill_matrix(Theta_r_free_c, Theta_r_skeleton_c[g], w7skel_c, g_start7_c[g,1], g_start7_c[g,2]); + Thetmat_c[g] = Thtmp_c + transpose(Thtmp_c) - diag_matrix(rep_vector(1, p_c)); + } + if (m_c > 0) { + matrix[m_c, m_c] Pstmp_c = fill_matrix(P_r_c, Psi_r_skeleton_c[g], w10skel_c, g_start10_c[g,1], g_start10_c[g,2]); + PSmat_c[g] = Pstmp_c + transpose(Pstmp_c) - diag_matrix(rep_vector(1, m_c)); + } + } + if (sum(thetanblk) > 0) { + Thetmat = fill_cov(Thetmat, thetablkse, thetanblk, Theta_r_mat_1, Theta_r_mat_2, Theta_r_mat_3, Theta_r_mat_4, Theta_r_mat_5, thetaorder, thetarevord); + } + if (sum(thetanblk_c) > 0) { + Thetmat_c = fill_cov(Thetmat_c, thetablkse_c, thetanblk_c, Theta_r_mat_1_c, Theta_r_mat_2_c, Theta_r_mat_3_c, Theta_r_mat_4_c, Theta_r_mat_5_c, thetaorder_c, thetarevord_c); + } + if (sum(psinblk) > 0) { + PSmat = fill_cov(PSmat, psiblkse, psinblk, Psi_r_mat_1, Psi_r_mat_2, Psi_r_mat_3, Psi_r_mat_4, Psi_r_mat_5, psiorder, psirevord); + } + if (sum(psinblk_c) > 0) { + PSmat_c = fill_cov(PSmat_c, psiblkse_c, psinblk_c, Psi_r_mat_1_c, Psi_r_mat_2_c, Psi_r_mat_3_c, Psi_r_mat_4_c, Psi_r_mat_5_c, psiorder_c, psirevord_c); + } + for (g in 1:Ng) { + Thet[g] = quad_form_sym(Thetmat[g], Theta_sd[g]); + Thet_c[g] = quad_form_sym(Thetmat_c[g], Theta_sd_c[g]); + PS[g] = quad_form_sym(PSmat[g], Psi_sd[g]); + PS_c[g] = quad_form_sym(PSmat_c[g], Psi_sd_c[g]); + } + // off-diagonal covariance parameter vectors, from cor/sd matrices: + if (p > 0 && len_free[8] > 0) { + /* iden is created so that we can re-use cor2cov, even though + we don't need to multiply to get covariances */ + array[Ng] matrix[p, p] iden; + for (g in 1:Ng) { + iden[g] = diag_matrix(rep_vector(1, p)); + } + Theta_cov = cor2cov(Thet, iden, len_free[8], Theta_r_skeleton_f, w8skel, Ng); + } + Theta_var = Theta_sd_free .* Theta_sd_free; + if (m > 0 && len_free[11] > 0) { + array[Ng] matrix[m, m] iden; + for (g in 1:Ng) { + iden[g] = diag_matrix(rep_vector(1, m)); + } + Psi_cov = cor2cov(PS, iden, len_free[11], Psi_r_skeleton_f, w11skel, Ng); + } else { + Psi_cov = P_r; + } + Psi_var = Psi_sd_free .* Psi_sd_free; + // and for level 2 + if (p_c > 0 && len_free_c[7] > 0) { + array[Ng] matrix[p, p] iden; + for (g in 1:Ng) { + iden[g] = diag_matrix(rep_vector(1, p)); + } + Theta_cov_c = cor2cov(Thet_c, iden, len_free_c[7], Theta_r_skeleton_c, w7skel_c, Ng); + } + Theta_var_c = Theta_sd_free_c .* Theta_sd_free_c; + if (m_c > 0 && len_free_c[11] > 0) { + array[Ng] matrix[m_c, m_c] iden_c; + for (g in 1:Ng) { + iden_c[g] = diag_matrix(rep_vector(1, m_c)); + } + Psi_cov_c = cor2cov(PS_c, iden_c, len_free_c[11], Psi_r_skeleton_f_c, w11skel_c, Ng); + } else { + Psi_cov_c = P_r_c; + } + Psi_var_c = Psi_sd_free_c .* Psi_sd_free_c; + { // log-likelihood + array[p + q] int obsidx; + array[p + q] int xidx; + array[p + q] int xdatidx; + int r1; + int r2; + int r3; + int r4; + int rr1; + int rr2; + int grpidx; + int clusidx; + if (do_test && use_cov) { + for (g in 1:Ng) { + Sigma_rep_sat[g] = wishart_rng(N[g] - 1, Sigma[g]); + } + } else if (do_test && has_data) { + // generate level 2 data, then level 1 + if (multilev) { + array[p_tilde - N_between] int notbidx; + notbidx = between_idx[(N_between + 1):p_tilde]; + r1 = 1; + rr1 = 1; + clusidx = 1; + r2 = 1; + for (gg in 1:Ng) { + matrix[p_c, p_c] Sigma_c_chol = cholesky_decompose(Sigma_c[gg]); + matrix[p + q, p + q] Sigma_chol = cholesky_decompose(Sigma[gg]); + S_PW_rep[gg] = rep_matrix(0, N_both + N_within, N_both + N_within); + S_PW_rep_full[gg] = rep_matrix(0, p_tilde, p_tilde); + S_B_rep[gg] = rep_matrix(0, p_tilde, p_tilde); + ov_mean_rep[gg] = rep_vector(0, p_tilde); + for (cc in 1:nclus[gg, 2]) { + vector[p_c] YXstar_rep_c; + vector[p_tilde] YXstar_rep_tilde; + YXstar_rep_c = multi_normal_cholesky_rng(Mu_c[gg], Sigma_c_chol); + YXstar_rep_tilde = calc_B_tilde(Sigma_c[gg], YXstar_rep_c, ov_idx2, p_tilde)[,1]; + for (ii in r1:(r1 + cluster_size[clusidx] - 1)) { + vector[N_within + N_both] Ywb_rep; + Ywb_rep = multi_normal_cholesky_rng(Mu[gg], Sigma_chol); + YXstar_rep[ii] = YXstar_rep_tilde; + for (ww in 1:(p_tilde - N_between)) { + YXstar_rep[ii, notbidx[ww]] += Ywb_rep[ww]; + } + ov_mean_rep[gg] += YXstar_rep[ii]; + } + for (jj in 1:p_tilde) { + mean_d_rep[clusidx, jj] = mean(YXstar_rep[r1:(r1 + cluster_size[clusidx] - 1), jj]); + } + r1 += cluster_size[clusidx]; + clusidx += 1; + } // cc + ov_mean_rep[gg] *= pow(nclus[gg, 1], -1); + xbar_b_rep[gg] = ov_mean_rep[gg]; + r1 -= nclus[gg, 1]; // reset for S_PW + clusidx -= nclus[gg, 2]; + if (N_between > 0) { + S2_rep[gg] = rep_matrix(0, N_between, N_between); + for (ii in 1:N_between) { + xbar_b_rep[gg, between_idx[ii]] = mean(mean_d_rep[clusidx:(clusidx + nclus[gg, 2] - 1), between_idx[ii]]); + } + } + for (cc in 1:nclus[gg, 2]) { + for (ii in r1:(r1 + cluster_size[clusidx] - 1)) { + S_PW_rep_full[gg] += tcrossprod(to_matrix(YXstar_rep[ii] - mean_d_rep[clusidx])); + } + S_B_rep[gg] += cluster_size[clusidx] * tcrossprod(to_matrix(mean_d_rep[clusidx] - ov_mean_rep[gg])); + if (N_between > 0) { + S2_rep[gg] += tcrossprod(to_matrix(mean_d_rep[clusidx, between_idx[1:N_between]] - xbar_b_rep[gg, between_idx[1:N_between]])); + } + r1 += cluster_size[clusidx]; + clusidx += 1; + } + S_PW_rep_full[gg] *= pow(nclus[gg, 1] - nclus[gg, 2], -1); + S_B_rep[gg] *= pow(nclus[gg, 2] - 1, -1); + S2_rep[gg] *= pow(nclus[gg, 2], -1); + // mods to between-only variables: + if (N_between > 0) { + array[N_between] int betonly = between_idx[1:N_between]; + S_PW_rep_full[gg, betonly, betonly] = rep_matrix(0, N_between, N_between); + // Y2: mean_d_rep; Y2c: mean_d_rep - ov_mean_rep + for (ii in 1:N_between) { + for (jj in 1:(N_both + N_within)) { + S_B_rep[gg, between_idx[ii], between_idx[(N_between + jj)]] *= (gs[gg] * nclus[gg, 2] * pow(nclus[gg, 1], -1)); + S_B_rep[gg, between_idx[(N_between + jj)], between_idx[ii]] = S_B_rep[gg, between_idx[ii], between_idx[(N_between + jj)]]; + } + } + S_B_rep[gg, betonly, betonly] = rep_matrix(0, N_between, N_between); + for (cc in 1:nclus[gg, 2]) { + S_B_rep[gg, betonly, betonly] += tcrossprod(to_matrix(mean_d_rep[cc, betonly] - ov_mean_rep[gg, betonly])); + } + S_B_rep[gg, betonly, betonly] *= gs[gg] * pow(nclus[gg, 2], -1); + } + cov_b_rep[gg] = pow(gs[gg], -1) * (S_B_rep[gg] - S_PW_rep_full[gg]); + if (N_between > 0) { + cov_b_rep[gg, between_idx[1:N_between], between_idx[1:N_between]] = S2_rep[gg]; + } + rr1 = r1 - nclus[gg, 1]; + r2 = clusidx - nclus[gg, 2]; + Mu_rep_sat[gg] = rep_vector(0, N_within + N_both); + if (N_within > 0) { + for (j in 1:N_within) { + xbar_b_rep[gg, within_idx[j]] = 0; + Mu_rep_sat[gg, within_idx[j]] = ov_mean_rep[gg, within_idx[j]]; + } + } + S_PW_rep[gg] = S_PW_rep_full[gg, notbidx, notbidx]; + if (Nx[gg] > 0 || Nx_between[gg] > 0) { + array[2] vector[p_tilde] mnvecs; + array[3] matrix[p_tilde, p_tilde] covmats; + mnvecs = calc_mean_vecs(YXstar_rep[rr1:(r1 - 1)], mean_d_rep[r2:(clusidx - 1)], nclus[gg], Xvar[gg], Xbetvar[gg], Nx[gg], Nx_between[gg], p_tilde); + covmats = calc_cov_mats(YXstar_rep[rr1:(r1 - 1)], mean_d_rep[r2:(clusidx - 1)], mnvecs, nclus[gg], Xvar[gg], Xbetvar[gg], Nx[gg], Nx_between[gg], p_tilde); + log_lik_x_rep[r2:(clusidx - 1)] = calc_log_lik_x(mean_d_rep[r2:(clusidx - 1)], + mnvecs[2], covmats[1], + covmats[2], covmats[3], + nclus[gg], cluster_size[r2:(clusidx - 1)], + Xvar[gg], Xbetvar[gg], Nx[gg], Nx_between[gg]); + } // Nx[gg] > 0 + } // gg + } else { + for (mm in 1:Np) { + matrix[Nobs[mm], Nobs[mm]] Sigma_chol; + obsidx = Obsvar[mm,]; + xidx = Xvar[mm,]; + xdatidx = Xdatvar[mm,]; + grpidx = grpnum[mm]; + r1 = startrow[mm]; + r2 = endrow[mm]; + Sigma_chol = cholesky_decompose(Sigma[ grpidx, obsidx[1:Nobs[mm]], obsidx[1:Nobs[mm]] ]); + for (jj in r1:r2) { + YXstar_rep[jj, 1:Nobs[mm]] = multi_normal_cholesky_rng(Mu[grpidx, obsidx[1:Nobs[mm]]], Sigma_chol); + } + } + if (missing) { + // start values for Mu and Sigma + for (g in 1:Ng) { + Mu_sat[g] = rep_vector(0, p + q); + Mu_rep_sat[g] = Mu_sat[g]; + Sigma_sat[g] = diag_matrix(rep_vector(1, p + q)); + Sigma_rep_sat[g] = Sigma_sat[g]; + } + for (jj in 1:emiter) { + satout = estep(YXstar, Mu_sat, Sigma_sat, Nobs, Obsvar, startrow, endrow, grpnum, Np, Ng); + satrep_out = estep(YXstar_rep, Mu_rep_sat, Sigma_rep_sat, Nobs, Obsvar, startrow, endrow, grpnum, Np, Ng); + // M step + for (g in 1:Ng) { + Mu_sat[g] = satout[g,,1]/N[g]; + Sigma_sat[g] = satout[g,,2:(p + q + 1)]/N[g] - Mu_sat[g] * Mu_sat[g]'; + Mu_rep_sat[g] = satrep_out[g,,1]/N[g]; + Sigma_rep_sat[g] = satrep_out[g,,2:(p + q + 1)]/N[g] - Mu_rep_sat[g] * Mu_rep_sat[g]'; + } + } + } else { + // complete data; Np patterns must only correspond to groups + for (mm in 1:Np) { + array[3] int arr_dims = dims(YXstar); + matrix[endrow[mm] - startrow[mm] + 1, arr_dims[2]] YXsmat; // crossprod needs matrix + matrix[endrow[mm] - startrow[mm] + 1, arr_dims[2]] YXsrepmat; + r1 = startrow[mm]; + r2 = endrow[mm]; + grpidx = grpnum[mm]; + for (jj in 1:(p + q)) { + Mu_sat[grpidx,jj] = mean(YXstar[r1:r2,jj]); + Mu_rep_sat[grpidx,jj] = mean(YXstar_rep[r1:r2,jj]); + } + for (jj in r1:r2) { + YXsmat[jj - r1 + 1] = (YXstar[jj] - Mu_sat[grpidx])'; + YXsrepmat[jj - r1 + 1] = (YXstar_rep[jj] - Mu_rep_sat[grpidx])'; + } + Sigma_sat[grpidx] = crossprod(YXsmat)/N[grpidx]; + Sigma_rep_sat[grpidx] = crossprod(YXsrepmat)/N[grpidx]; + // FIXME? Sigma_sat[grpidx] = tcrossprod(YXsmat); does not throw an error?? + } + } + for (g in 1:Ng) { + Sigma_sat_inv_grp[g] = inverse_spd(Sigma_sat[g]); + logdetS_sat_grp[g] = log_determinant(Sigma_sat[g]); + Sigma_rep_sat_inv_grp[g] = inverse_spd(Sigma_rep_sat[g]); + logdetS_rep_sat_grp[g] = log_determinant(Sigma_rep_sat[g]); + } + for (mm in 1:Np) { + Sigma_sat_inv[mm, 1:(Nobs[mm] + 1), 1:(Nobs[mm] + 1)] = sig_inv_update(Sigma_sat_inv_grp[grpnum[mm]], Obsvar[mm,], Nobs[mm], p + q, logdetS_sat_grp[grpnum[mm]]); + Sigma_rep_sat_inv[mm, 1:(Nobs[mm] + 1), 1:(Nobs[mm] + 1)] = sig_inv_update(Sigma_rep_sat_inv_grp[grpnum[mm]], Obsvar[mm,], Nobs[mm], p + q, logdetS_rep_sat_grp[grpnum[mm]]); + } + } + } + // compute log-likelihoods + if (multilev) { // multilevel + r1 = 1; + r3 = 1; + r2 = 0; + r4 = 0; + for (mm in 1:Np) { + grpidx = grpnum[mm]; + if (grpidx > 1) { + r1 += nclus[(grpidx - 1), 2]; + r3 += nclus[(grpidx - 1), 1]; + } + r2 += nclus[grpidx, 2]; + r4 += nclus[grpidx, 1]; + log_lik[r1:r2] = twolevel_logdens(mean_d_full[r1:r2], cov_d_full[r1:r2], S_PW[grpidx], YX[r3:r4], + nclus[grpidx,], cluster_size[r1:r2], cluster_size[r1:r2], + nclus[grpidx,2], intone[1:nclus[grpidx,2]], Mu[grpidx], + Sigma[grpidx], Mu_c[grpidx], Sigma_c[grpidx], + ov_idx1, ov_idx2, within_idx, between_idx, both_idx, + p_tilde, N_within, N_between, N_both); + if (Nx[grpidx] + Nx_between[grpidx] > 0) log_lik[r1:r2] -= log_lik_x_full[r1:r2]; + } + } + zmat = rep_matrix(0, p + q, p + q); + // reset for 2-level loglik: + rr1 = 1; + r3 = 1; + rr2 = 0; + r4 = 0; + for (mm in 1:Np) { + obsidx = Obsvar[mm,]; + xidx = Xvar[mm, 1:(p + q)]; + xdatidx = Xdatvar[mm, 1:(p + q)]; + grpidx = grpnum[mm]; + r1 = startrow[mm]; + r2 = endrow[mm]; + if (use_cov) { + log_lik[mm] = wishart_lpdf((N[mm] - 1) * Sstar[mm] | N[mm] - 1, Sigma[mm]); + if (do_test) { + log_lik_sat[mm] = -log_lik[mm] + wishart_lpdf((N[mm] - 1) * Sstar[mm] | N[mm] - 1, Sstar[mm]); + log_lik_rep[mm] = wishart_lpdf(Sigma_rep_sat[mm] | N[mm] - 1, Sigma[mm]); + log_lik_rep_sat[mm] = wishart_lpdf(Sigma_rep_sat[mm] | N[mm] - 1, pow(N[mm] - 1, -1) * Sigma_rep_sat[mm]); + } + if (Nx[mm] > 0) { + array[Nx[mm]] int xvars = xdatidx[1:Nx[mm]]; + log_lik[mm] += -wishart_lpdf((N[mm] - 1) * Sstar[mm, xvars, xvars] | N[mm] - 1, Sigma[mm, xvars, xvars]); + if (do_test) { + log_lik_sat[mm] += wishart_lpdf((N[mm] - 1) * Sstar[mm, xvars, xvars] | N[mm] - 1, Sigma[mm, xvars, xvars]); + log_lik_sat[mm] += -wishart_lpdf((N[mm] - 1) * Sstar[mm, xvars, xvars] | N[mm] - 1, Sstar[mm, xvars, xvars]); + log_lik_rep[mm] += -wishart_lpdf(Sigma_rep_sat[mm, xvars, xvars] | N[mm] - 1, Sigma[mm, xvars, xvars]); + log_lik_rep_sat[mm] += -wishart_lpdf(Sigma_rep_sat[mm, xvars, xvars] | N[mm] - 1, pow(N[mm] - 1, -1) * Sigma_rep_sat[mm, xvars, xvars]); + } + } + } else if (has_data && !multilev) { + for (jj in r1:r2) { + log_lik[jj] = multi_normal_suff(YXstar[jj, 1:Nobs[mm]], zmat[1:Nobs[mm], 1:Nobs[mm]], Mu[grpidx, obsidx[1:Nobs[mm]]], Sigmainv[mm], 1); + // TODO efficiency can be improved by getting inverse/chol of Sigma outside loop + if (Nx[mm] > 0 && !missing) { + log_lik[jj] += -1.0 * multi_normal_suff(YXstar[jj, xdatidx[1:Nx[mm]]], zmat[1:Nx[mm], 1:Nx[mm]], Mu[grpidx, xidx[1:Nx[mm]]], sig_inv_update(Sigmainv[grpidx], xidx, Nx[mm], p + q, logdetSigma_grp[grpidx]), 1); + } else if (Nx[mm] > 0) { + log_lik[jj] += -1.0 * multi_normal_lpdf(YXstar[jj, xdatidx[1:Nx[mm]]] | Mu[grpidx, xidx[1:Nx[mm]]], Sigma[grpidx, xidx[1:Nx[mm]], xidx[1:Nx[mm]]]); + } + } + } + // saturated and y_rep likelihoods for ppp + if (do_test) { + if (multilev) { + // compute clusterwise log_lik_rep for grpidx + if (grpidx > 1) { + rr1 += nclus[(grpidx - 1), 2]; + r3 += nclus[(grpidx - 1), 1]; + } + rr2 += nclus[grpidx, 2]; + r4 += nclus[grpidx, 1]; + // NB: cov_d is 0 when we go cluster by cluster. + // otherwise it is covariance of cluster means by each unique cluster size + // because we go cluster by cluster here, we can reuse cov_d_full everywhere + log_lik_rep[rr1:rr2] = twolevel_logdens(mean_d_rep[rr1:rr2], cov_d_full[rr1:rr2], + S_PW_rep[grpidx], YXstar_rep[r3:r4], + nclus[grpidx,], cluster_size[rr1:rr2], + cluster_size[rr1:rr2], nclus[grpidx,2], + intone[1:nclus[grpidx,2]], Mu[grpidx], + Sigma[grpidx], Mu_c[grpidx], Sigma_c[grpidx], + ov_idx1, ov_idx2, within_idx, between_idx, + both_idx, p_tilde, N_within, N_between, N_both); + log_lik_sat[rr1:rr2] = twolevel_logdens(mean_d_full[rr1:rr2], cov_d_full[rr1:rr2], + S_PW[grpidx], YX[r3:r4], + nclus[grpidx,], cluster_size[rr1:rr2], + cluster_size[rr1:rr2], nclus[grpidx,2], + intone[1:nclus[grpidx,2]], xbar_w[grpidx, ov_idx1], + S_PW[grpidx], xbar_b[grpidx, ov_idx2], cov_b[grpidx, ov_idx2, ov_idx2], + ov_idx1, ov_idx2, within_idx, between_idx, + both_idx, p_tilde, N_within, N_between, N_both); + log_lik_rep_sat[rr1:rr2] = twolevel_logdens(mean_d_rep[rr1:rr2], cov_d_full[rr1:rr2], + S_PW_rep[grpidx], YXstar_rep[r3:r4], + nclus[grpidx,], cluster_size[rr1:rr2], + cluster_size[rr1:rr2], nclus[grpidx,2], + intone[1:nclus[grpidx,2]], Mu_rep_sat[grpidx], + S_PW_rep[grpidx], xbar_b_rep[grpidx, ov_idx2], + cov_b_rep[grpidx, ov_idx2, ov_idx2], + ov_idx1, ov_idx2, + within_idx, between_idx, both_idx, p_tilde, + N_within, N_between, N_both); + if (Nx[grpidx] + Nx_between[grpidx] > 0) { + log_lik_rep[rr1:rr2] -= log_lik_x_rep[rr1:rr2]; + log_lik_sat[rr1:rr2] -= log_lik_x_full[rr1:rr2]; + log_lik_rep_sat[rr1:rr2] -= log_lik_x_rep[rr1:rr2]; + } + // we subtract log_lik here so that _sat always varies and does not lead to + // problems with rhat and neff computations + log_lik_sat[rr1:rr2] -= log_lik[rr1:rr2]; + } else if (!use_cov) { + r1 = startrow[mm]; + r2 = endrow[mm]; + for (jj in r1:r2) { + // log_lik_rep, _sat, _rep_sat + log_lik_rep[jj] = multi_normal_suff(YXstar_rep[jj, 1:Nobs[mm]], zmat[1:Nobs[mm], 1:Nobs[mm]], Mu[grpidx, obsidx[1:Nobs[mm]]], Sigmainv[mm], 1); + log_lik_sat[jj] = multi_normal_suff(YXstar[jj, 1:Nobs[mm]], zmat[1:Nobs[mm], 1:Nobs[mm]], Mu_sat[grpidx, obsidx[1:Nobs[mm]]], Sigma_sat_inv[mm], 1); + log_lik_rep_sat[jj] = multi_normal_suff(YXstar_rep[jj, 1:Nobs[mm]], zmat[1:Nobs[mm], 1:Nobs[mm]], Mu_rep_sat[grpidx, obsidx[1:Nobs[mm]]], Sigma_rep_sat_inv[mm], 1); + // TODO efficiency can be improved by getting inverse/chol of Sigma outside loop + if (Nx[mm] > 0 && !missing) { + log_lik_rep[jj] += -1.0 * multi_normal_suff(YXstar_rep[jj, xdatidx[1:Nx[mm]]], zmat[1:Nx[mm], 1:Nx[mm]], Mu[grpidx, xidx[1:Nx[mm]]], sig_inv_update(Sigmainv[grpidx], xidx, Nx[mm], p + q, logdetSigma_grp[grpidx]), 1); + log_lik_sat[jj] += -1.0 * multi_normal_suff(YXstar[jj, xdatidx[1:Nx[mm]]], zmat[1:Nx[mm], 1:Nx[mm]], Mu_sat[grpidx, xidx[1:Nx[mm]]], sig_inv_update(Sigma_sat_inv[grpidx], xidx, Nx[mm], p + q, logdetS_sat_grp[grpidx]), 1); + log_lik_rep_sat[jj] += -1.0 * multi_normal_suff(YXstar_rep[jj, xdatidx[1:Nx[mm]]], zmat[1:Nx[mm], 1:Nx[mm]], Mu_rep_sat[grpidx, xidx[1:Nx[mm]]], sig_inv_update(Sigma_rep_sat_inv[grpidx], xidx, Nx[mm], p + q, logdetS_rep_sat_grp[grpidx]), 1); + } else if (Nx[mm] > 0) { + log_lik_rep[jj] += -1.0 * multi_normal_lpdf(YXstar_rep[jj, xdatidx[1:Nx[mm]]] | Mu[grpidx, xidx[1:Nx[mm]]], Sigma[grpidx, xidx[1:Nx[mm]], xidx[1:Nx[mm]]]); + log_lik_sat[jj] += -1.0 * multi_normal_lpdf(YXstar[jj, xdatidx[1:Nx[mm]]] | Mu_sat[grpidx, xidx[1:Nx[mm]]], Sigma_sat[grpidx, xidx[1:Nx[mm]], xidx[1:Nx[mm]]]); + log_lik_rep_sat[jj] += -1.0 * multi_normal_lpdf(YXstar_rep[jj, xdatidx[1:Nx[mm]]] | Mu_rep_sat[grpidx, xidx[1:Nx[mm]]], Sigma_rep_sat[grpidx, xidx[1:Nx[mm]], xidx[1:Nx[mm]]]); + } + } + // we subtract log_lik here so that _sat always varies and does not lead to + // problems with rhat and neff computations + log_lik_sat[r1:r2] -= log_lik[r1:r2]; + } + } + } + if (do_test) { + ppp = step((-sum(log_lik_rep) + sum(log_lik_rep_sat)) - (sum(log_lik_sat))); + } else { + ppp = 0; + } + } +} // end with a completely blank line (not even whitespace) \ No newline at end of file diff --git a/BayesSEMsensitivity/BayesSEMSensitivity.R b/BayesSEMsensitivity/BayesSEMSensitivity.R new file mode 100644 index 0000000..6c28bb7 --- /dev/null +++ b/BayesSEMsensitivity/BayesSEMSensitivity.R @@ -0,0 +1,210 @@ +## LIBRARIES: +library(blavaan) +library(posterior) +library(shinystan) +library(priorsense) +library(rstan) + +## DEFINE THE BAYES SEM MODEL: +RI_CLPM_ALSPAC_Z <- ' +# ############################################################ +# # # 1. RANDOM INTERCEPTS (STABLE TRAITS) +# ############################################################ +# # General trend over the time +# +# RI_BMI =~ 1*bmi_12.5y_z + 1*bmi_14y_z + 1*bmi_16y_z + 1*bmi_18y_z + 1*bmi_24y_z +# RI_BE =~ 1*Binge_14 + 1*Binge_16 + 1*Binge_18 + 1*Binge_24 +# +# ############################################################ +# # 2. WITHIN-PERSON +# ############################################################ +# # Create latent factors per each wave that contains only the variation +# +# wBMI_12.5 =~ 1*bmi_12.5y_z +# wBMI_14 =~ 1*bmi_14y_z +# wBMI_16 =~ 1*bmi_16y_z +# wBMI_18 =~ 1*bmi_18y_z +# wBMI_24 =~ 1*bmi_24y_z +# +# wBE_14 =~ 1*Binge_14 +# wBE_16 =~ 1*Binge_16 +# wBE_18 =~ 1*Binge_18 +# wBE_24 =~ 1*Binge_24 +# +# # Force the residual variance of the observed variable to 0 +# # and force the variance in within factors +# bmi_12.5y_z ~~ 0*bmi_12.5y_z +# bmi_14y_z ~~ 0*bmi_14y_z +# bmi_16y_z ~~ 0*bmi_16y_z +# bmi_18y_z ~~ 0*bmi_18y_z +# bmi_24y_z ~~ 0*bmi_24y_z +# +# Binge_14 ~~ 0*Binge_14 +# Binge_16 ~~ 0*Binge_16 +# Binge_18 ~~ 0*Binge_18 +# Binge_24 ~~ 0*Binge_24 +# +# ############################################################ +# # 3. CROSS-LAG AND AUTOREGRESSIVE DYNAMICS +# ############################################################ +# +# BMI PREDICTION (Weight-driven effects) +# wBMI_14 ~ a1*wBMI_12.5 +# wBMI_16 ~ a2*wBMI_14 + d2*wBE_14 +# wBMI_18 ~ a3*wBMI_16 + d3*wBE_16 +# wBMI_24 ~ a4*wBMI_18 + d4*wBE_18 +# +# Binge Eating PREDICTION (Eating-driven effects) +# wBE_14 ~ c1*wBMI_12.5 +# wBE_16 ~ b1*wBE_14 + c2*wBMI_14 +# wBE_18 ~ b2*wBE_16 + c3*wBMI_16 +# wBE_24 ~ b3*wBE_18 + c4*wBMI_18 +# +# ############################################################ +# # 4. VARIANCE AND COVARIANCE +# ############################################################ +# +# # STABLE TRAITS +# RI_BMI ~~ RI_BE +# +# # cross-sectional correlation between deviations (residuals) +# wBMI_14 ~~ wBE_14 +# wBMI_16 ~~ wBE_16 +# wBMI_18 ~~ wBE_18 +# wBMI_24 ~~ wBE_24 +# +# # Define the correlations between RI and first measurements: +# RI_BE ~~ prior("beta(3,3)")*wBE_14 +# RI_BMI ~~ prior("beta(3,3)")*wBMI_12.5 +# ' + + + +## FIT THE MODEL: +fit_RI_CLPM_Z <- bsem( + RI_CLPM_ALSPAC_Z, + data = df, + ordered = c("Binge_14", "Binge_16", "Binge_18", "Binge_24"), + burnin = 6000, + sample = 13000, + n.chains = 4, + bcontrol = list( + control = list( + adapt_delta = 0.99, + max_treedepth = 15)), + dp = dpriors( + beta = "normal(0,0.3)", + psi = "gamma(2,1)", + rho = "beta(2,2)" + ) + ) + + summary(fit_RI_CLPM_Z, standardized = TRUE, fit.measures = TRUE) + lavaan::fitMeasures(fit_RI_CLPM_Z) + + +## STAN-LIKE DATASET CONFIGURATION +## AND SAVE THE PARAMETER IN WHICH I AM INTERESTED IN: + stan_data_list <- fit_RI_CLPM_Z@external$mcmcdata + + pars_save <- c( + "B_free", + "Alpha_free", + "Psi_sd_free", + "Psi_r_mat_1", + "Psi_r_mat_2", + "Nu_free", + "Tau_ufree", + "log_prior", + "log_lik") + + SEM_stan_model <- stan( + file = "BayesSEM.stan", + data = stan_data_list, + chains = 4, + iter = 16000, + warmup = 4000, + cores = 1, + pars = pars_save, + include = TRUE, + sample_file = "Stan_model/chain") + + print(SEM_stan_model, pars = "log_prior") + + psd <- priorsense:::create_priorsense_data.stanfit( + SEM_stan_model, + log_prior_name = "log_prior" + ) + + # If I dind't put the parameter to save I remove the ones that are not useful and contains Nas values + # any_nan <- sapply(psd$draws, function(z) any(is.nan(z))) + # bad_vars <- names(any_nan)[any_nan] + # bad_prefixes <- unique(sub("\\[.*", "", bad_vars)) + # + # draws_clean <- posterior::subset_draws( + # psd$draws, + # variable = bad_prefixes, + # exclude = TRUE + # ) + # + # sum(sapply(draws_clean, function(z) any(is.nan(z)))) + # + # psd2 <- psd + # psd2$draws <- draws_clean + # + # ps <- priorsense:::powerscale_sensitivity.priorsense_data( + # psd2, + # log_prior_name = "log_prior" + # ) + # + + + # Check the sensitivity of the parameter I am interested in: + ps_seq <- powerscale_sequence( + SEM_stan_model, + component = "prior", + variable = c("B_free[2]","B_free[4]","B_free[6]", + "B_free[8]", "B_free[10]","B_free[11]", + "B_free[14]"), + log_prior_name = "log_prior", + lower_alpha = 0.7,upper_alpha = 1.5,length = 50) + + powerscale_plot_dens( + ps_seq, + variable = c(#"B_free[2]","B_free[4]", + #"B_free[6]","B_free[8]", + #"B_free[10]","B_free[11]", + "B_free[14]"# + )) + + powerscale_plot_quantities( + ps_seq, + variable = c(#"B_free[2]","B_free[4]", + #"B_free[6]","B_free[8]", + #"B_free[10]","B_free[11]", + "B_free[14]"# + ), + quantity = c("mean", "sd", "quantile")) + +# Try to tweak the priors +# (In this case the normal distribution related to the coefficients) + stan_data_sens_2 <- stan_data_list +# Now I have N(0, 0.7). +# NOTE: I notice that in the stan model I used the N(0, 0.3^2). +# SO IMPORTANT TO FIX B_SD ACCORDINGLY: + stan_data_sens_2$b_sd <- rep(sqrt(0.7), 14) + SEM_stan_model@stanmodel + + SEM_stan_model_NEW <- sampling( + SEM_stan_model@stanmodel, + data = stan_data_sens, + chains = 4, + iter = 16000, + warmup = 4000, + thin = 1, + cores=1, + seed = 1234, + sample_file = "Stan_sens_03/chain") + +cat("Sampling Completed. Cleaning the memory...\n") +gc()