Lecture 08

1 Review Previous Class

1.1 In Previous Class…

1. We introduced a self-reported survey data, called Conspiracy Theories

• The scale has 10 items with 5-point Likert scale response
• Items have varied item difficult and Positive skewed item response distribution

2. We talked about using if factor analysis is a proper method when normality assumption is violated
3. We checked the item characteristics curves

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks rstan::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(here)

here() starts at /Users/jihong/Documents/website-jihong

library(blavaan)

Loading required package: Rcpp
This is blavaan 0.5-8
On multicore systems, we suggest use of future::plan("multicore") or
  future::plan("multisession") for faster post-MCMC computations.

self_color <- c("#DB7093", "#AFEEEE", "#3CB371", "#9370DB", "#FFD700")
root_dir <- "teaching/2024-01-12-syllabus-adv-multivariate-esrm-6553/Lecture07/Code"
dat <- read.csv(here(root_dir, 'conspiracies.csv'))
itemResp <- dat[,1:10]
colnames(itemResp) <- paste0('item', 1:10)
conspiracyItems = itemResp
itemResp |>
  rownames_to_column("ID") |>
  pivot_longer(-ID, names_to = "Item", values_to = "Response") |>
  mutate(Item = factor(Item, levels = paste0('item', 1:10)),
         Response = factor(Response, levels = 1:5)) |>
  ggplot() +
  geom_bar(aes(x = Response, fill = Response, group = Response),
           position = position_stack()) +
  facet_wrap(~ Item, nrow = 1, ncol = 10) +
  theme_bw() +
  scale_fill_manual(values = self_color)

1.2 Today’s Lecture Objectives

1. Dive deep into factor scoring
2. Show how different initial values affect Bayesian model estimation
3. Show how parameterization differs for standardized latent variables vs. marker item scale identification

Lecture 07 Files: Lecture07.R | conspiracies.csv

1.3 Posterior Distribution for Item Parameters

Before moving onto the latent variable, let’s note the posterior distribution of item parameters (for a single item):

f(\mu_i,\lambda_i,\psi_i\mid\boldsymbol{Y})\propto f(\boldsymbol{Y}\mid\mu_i,\lambda_i,\psi_i)f(\mu_i, \lambda_i,\psi_i)

where f(\boldsymbol{Y}\mid\mu_i,\lambda_i,\psi_i) is the (joint) posterior distribution of the parameters for item i conditional on the adta

f(\boldsymbol{Y}\mid\mu_i,\lambda_i,\psi_i) is the distribution we defined for our observed data:

f(\boldsymbol{Y}\mid\mu_i,\lambda_i,\psi_i)\sim N(\mu_i+\lambda_i\theta_p, \psi_i)

f(\mu_i, \lambda_i, \psi_i) is the (joint) prior distribution for each of the parameters, which, are independent:

f(\mu_i, \lambda_i, \psi_i) = f(\mu_i)f(\lambda_i)f(\psi_i)

1.4 Investigating the Latent Variables

The estimated latent variables are then:

save_dir <- "~/Library/CloudStorage/OneDrive-Personal/2024_Spring/ESRM6553 - Advanced Multivariate Modeling/Lecture07"
modelCFA_samples <- readRDS(here(save_dir, "model01.RDS"))
modelCFA_samplesFail <- readRDS(here(save_dir, "model01fail.RDS"))
modelCFA_samplesFix <- readRDS(here(save_dir, "model03fix.RDS"))
model01_blv <- readRDS(here(save_dir, "model01_blv.RDS"))

modelCFA_samples$summary('theta')

# A tibble: 177 × 10
   variable     mean  median    sd   mad     q5     q95  rhat ess_bulk ess_tail
   <chr>       <dbl>   <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 theta[1]   0.0221  0.0222 0.246 0.250 -0.386  0.419   1.00    5122.    5185.
 2 theta[2]   1.53    1.53   0.257 0.254  1.11   1.96    1.00    3754.    4698.
 3 theta[3]   1.71    1.71   0.264 0.263  1.27   2.14    1.00    3809.    4720.
 4 theta[4]  -0.927  -0.924  0.250 0.246 -1.34  -0.517   1.00    4936.    4752.
 5 theta[5]   0.0443  0.0485 0.248 0.244 -0.368  0.449   1.00    5126.    4699.
 6 theta[6]  -0.976  -0.972  0.252 0.249 -1.40  -0.566   1.00    3786.    4650.
 7 theta[7]  -0.336  -0.338  0.254 0.249 -0.754  0.0846  1.00    4421.    4882.
 8 theta[8]  -0.0520 -0.0548 0.248 0.247 -0.460  0.349   1.00    5037.    5503.
 9 theta[9]  -0.785  -0.784  0.244 0.248 -1.19  -0.388   1.00    4633.    5479.
10 theta[10]  0.0477  0.0460 0.248 0.250 -0.357  0.455   1.00    5000.    4896.
# ℹ 167 more rows

1.5 EAP Estimates of Latent Variables

1.6 Density of EAP Estimates

1.7 Density of 500 Posterior draws of θ

1.8 Comparing Posterior Distribution for Individuals

1.9 Comparing EAP Estimate with Posterior SD

1.10 Comparing EAP Estimate with Sum Score

1.11 Comparing EAP Estimate with Factor Score by lavaan

1.12 Posterior Distribution for Person Parameter θ

The posterior distribution of the person parameters (the latent variable for a single person):

f(\theta_p\mid\boldsymbol{Y}) \propto f(\boldsymbol{Y} \mid \theta_p)f(\theta_p)

Here:

• f(\theta_p\mid\boldsymbol{Y}) is the posterior distribution of the latent variable conditional on the observed data;

• f(\boldsymbol{Y} \mid \theta_p) is the model (data) likelihood:

  f(\boldsymbol{Y} \mid \theta_p) = \prod_{i=1}^{I} f(Y_i\mid\theta_p)

  • f(Y_i \mid \theta_p) is one individual item’s data likelihood: f(Y_i \mid \theta_p) \sim N(\mu_i+\lambda_i\theta_p, \psi_i);

• f(\theta_p) \sim N(0,1) is the prior distribution for the latent variable – \theta_p

2 Measurement Model Estimation Fails

2.1 Recall: Stan’s parameters {} Block

parameters {
  vector[nObs] theta;                // the latent variables (one for each person)
  vector[nItems] mu;                 // the item intercepts (one for each item)
  vector[nItems] lambda;             // the factor loadings/item discriminations (one for each item)
  vector<lower=0>[nItems] psi;       // the unique standard deviations (one for each item)
}

Here, the parameterization of \lambda (factor loadings / discrimination parameters) can lead to problems in estimation

• The issue is \lambda_i\theta_p=(-\lambda_i)(-\theta_p)

  • Depending on the random starting values of each of these parameters (per Markov chain), a given chain may converge to different region with others

• To demonstrate, we will start with different random number seend

  • Currently use 09102022: works fine

  • Change to 25102022: big problem

2.2 New Samples Syntax

modelCFA_samplesFail = modelCFA_stan$sample(
  data = modelCFA_data,
  seed = 25102022,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 1000,
  iter_sampling = 2000
)

Convergence fail with maximum of \hat R as:

max(modelCFA_samplesFail$summary(.cores =4)$'rhat')

[1] 1.73562

2.3 Why Convergence Failed

• The issue of exchangeable likelihood: \lambda_i\theta_p=(-\lambda_i)(-\theta_p)

# A tibble: 10 × 10
   variable     mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
   <chr>       <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 lambda[1]  -0.001  0.01  0.743 1.09  -0.843 0.841  1.73     6.08     106.
 2 lambda[2]  -0.002 -0.004 0.873 1.29  -0.972 0.969  1.73     6.10     121.
 3 lambda[3]  -0.002  0.007 0.805 1.19  -0.901 0.899  1.73     6.08     113.
 4 lambda[4]  -0.001 -0.012 0.846 1.25  -0.945 0.941  1.73     6.09     113.
 5 lambda[5]  -0.001 -0.002 0.999 1.47  -1.09  1.09   1.73     6.13     122.
 6 lambda[6]  -0.001  0.008 0.9   1.33  -0.983 0.98   1.73     6.14     103.
 7 lambda[7]  -0.002 -0.002 0.767 1.13  -0.855 0.85   1.73     6.08     118.
 8 lambda[8]  -0.001 -0.002 0.855 1.26  -0.931 0.932  1.73     6.09     108.
 9 lambda[9]  -0.002  0.029 0.863 1.27  -0.963 0.962  1.73     6.12     101.
10 lambda[10] -0.001 -0.003 0.677 0.993 -0.775 0.774  1.73     6.08     105.

2.4 Posterior Trace Plots of λ

Unfortunately, we are unable to extract initial values that were generated automatically by cmdstanr. See here for more details.

2.5 Posterior Density Plots of λ

2.6 Examining Latent Variable

# A tibble: 177 × 10
   variable    mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
   <chr>      <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 theta[1]  -0.002  0.001 0.247 0.245 -0.414 0.405  1.00  3167.     4819. 
 2 theta[2]  -0.001  0.004 1.56  2.27  -1.87  1.87   1.73     6.08    104. 
 3 theta[3]   0      0.029 1.73  2.53  -2.06  2.05   1.73     6.07    105. 
 4 theta[4]  -0.008  0.018 0.962 1.37  -1.26  1.24   1.73     6.07    130. 
 5 theta[5]   0.001  0     0.247 0.249 -0.403 0.410  1.02   177.     4690. 
 6 theta[6]  -0.004  0.05  1.01  1.45  -1.30  1.30   1.73     6.07    106. 
 7 theta[7]  -0.005 -0.002 0.422 0.503 -0.672 0.656  1.60     6.65    104. 
 8 theta[8]  -0.003 -0.001 0.254 0.248 -0.427 0.417  1.03   112.     3799. 
 9 theta[9]  -0.009  0.015 0.825 1.16  -1.11  1.09   1.73     6.06    113. 
10 theta[10] -0.004 -0.004 0.252 0.251 -0.416 0.418  1.02   168.     3352. 
11 theta[11] -0.006 -0.055 1.01  1.45  -1.30  1.29   1.73     6.07     99.4
12 theta[12] -0.004 -0.006 0.356 0.395 -0.571 0.578  1.46     7.62    105. 
13 theta[13] -0.007 -0.025 0.787 1.11  -1.07  1.06   1.73     6.06    106. 
14 theta[14] -0.006 -0.101 1.01  1.45  -1.32  1.29   1.73     6.06     99.0
15 theta[15] -0.005  0.037 0.875 1.24  -1.16  1.14   1.73     6.06     97.9
16 theta[16] -0.004 -0.005 0.259 0.259 -0.434 0.424  1.03    82.2    3490. 
17 theta[17] -0.004 -0.006 0.263 0.263 -0.438 0.430  1.03    88.1    4232. 
18 theta[18] -0.006 -0.017 0.710 0.985 -0.984 0.984  1.73     6.08    101. 
19 theta[19]  0.001 -0.005 1.89  2.77  -2.22  2.23   1.73     6.07    116. 
20 theta[20] -0.007  0.072 0.962 1.37  -1.26  1.25   1.74     6.05    110. 
21 theta[21] -0.008  0     0.350 0.390 -0.571 0.548  1.45     7.81    104. 
22 theta[22] -0.005 -0.051 1.01  1.45  -1.30  1.30   1.73     6.08    109. 
23 theta[23] -0.003  0.032 1.01  1.45  -1.30  1.29   1.73     6.06    114. 
24 theta[24] -0.008 -0.027 1.01  1.43  -1.30  1.30   1.73     6.06    108. 
25 theta[25] -0.005 -0.008 0.602 0.807 -0.883 0.865  1.73     6.07    106. 
26 theta[26] -0.003 -0.009 0.779 1.09  -1.07  1.06   1.73     6.07    107. 
27 theta[27] -0.002 -0.007 0.276 0.275 -0.452 0.451  1.10    24.5     375. 
28 theta[28] -0.002  0.091 1.24  1.80  -1.55  1.55   1.73     6.07     99.7
29 theta[29] -0.004 -0.032 0.928 1.32  -1.23  1.21   1.73     6.06    107. 
30 theta[30] -0.001 -0.057 1.70  2.49  -2.03  2.04   1.73     6.07    112. 
# ℹ 147 more rows

2.7 Posterior Trace Plots of θ

2.8 Posterior Density Plots of θ

2.9 Fixing Convergence

Stan allows starting values to be set via cmdstanr

• Documentation is very lacking, I can show you one method but it may change in the future.

Alternatively:

• In Stan file, restrict \lambda to be positive: vector<lower=0>[nItems] lambda;

• Can also choose prior that has strictly positive range (like log-normal distribution)

• Note: the restriction on the space of \lambda will not permit truely negative values

  • Not ideals as negative \lambda values are informative as a problem with data

  parameters {
    vector[nObs] theta;                // the latent variables (one for each person)
    vector[nItems] mu;                 // the item intercepts (one for each item)
    vector<lower=0>[nItems] lambda;    // the factor loadings/item discriminations (one for each item)
    vector<lower=0>[nItems] psi;       // the unique standard deviations (one for each item)
  }

2.10 Setting Starting (Inital) Values in Stan

Starting values (initial values) are the first values used when an MCMC chain starts

• In Stan, by default, parameters are randomly started between -2 and 2

  • Bounded parameters are transformed so they are unbounded in the algorithm

• What we need:

  • Randomly start all \lambda parameters so that they converge to the \lambda_i\theta_p mode

  • As opposed to the (-\lambda_i)(-\theta_p) mode

2.11 cmdstanr Syntax for Initial Values

Add the init option to the $sample() function of the cmdstanr object:

# set starting values for some of the parameters
modelCFA_samples2fixed = modelCFA_stan$sample(
  data = modelCFA_data,
  seed = 25102022,
  chains = 4,
  parallel_chains = 4,
  iter_warmup = 2000,
  iter_sampling = 2000,
  init = function() list(lambda=rnorm(nItems, mean=10, sd=2))
)

The init option can be specified as a function, here, randomly starting each \lambda following a normal distribution

2.12 Initialization Process

See the lecture R syntax for information on how to confirm starting values are set.

You should find the initial values using $init() function

modelCFA_samplesFix$init()

[[1]]
[[1]]$lambda
 [1]  7.505537  8.652382 13.870335 11.813123 11.473054 10.293251 10.566242
 [8] 11.178553 10.431699 10.310106


[[2]]
[[2]]$lambda
 [1] 11.907653  9.668295 10.510007  8.205055 13.167940 12.785517  6.375484
 [8]  9.567361 11.099712 10.965472


[[3]]
[[3]]$lambda
 [1] 11.521697  9.108684  7.589333 10.602933  6.921710 11.270741 11.405904
 [8]  6.188234 11.877843  9.551016


[[4]]
[[4]]$lambda
 [1]  8.652366 10.891575 12.561234 13.130260  7.598908  9.124610 10.292763
 [8] 10.132038  9.109928  5.315775

2.13 Final Results: Parameters

# Check convergence
max(modelCFA_samplesFix$summary()$'rhat')

[1] 1.006357

print(modelCFA_samplesFix$summary(c("mu", "lambda", "psi"),.cores = 4) |> mutate(across(c(mean, median), \(x) round(x, 3))), n = Inf)

# A tibble: 30 × 10
   variable    mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
   <chr>      <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 mu[1]      2.37   2.37  0.0866 0.0851 2.22  2.51   1.00    1317.    3146.
 2 mu[2]      1.95   1.95  0.0857 0.0846 1.81  2.09   1.01     908.    2329.
 3 mu[3]      1.87   1.87  0.0839 0.0848 1.74  2.01   1.00    1111.    3011.
 4 mu[4]      2.01   2.01  0.0847 0.0844 1.87  2.15   1.00    1040.    2582.
 5 mu[5]      1.98   1.98  0.0859 0.0868 1.84  2.12   1.01     790.    2080.
 6 mu[6]      1.89   1.89  0.0766 0.0768 1.77  2.02   1.01     750.    2168.
 7 mu[7]      1.72   1.72  0.0767 0.0740 1.59  1.85   1.01    1000.    2708.
 8 mu[8]      1.84   1.84  0.0723 0.0719 1.72  1.96   1.01     775.    2012.
 9 mu[9]      1.80   1.81  0.0863 0.0854 1.66  1.95   1.00    1010.    2552.
10 mu[10]     1.52   1.52  0.0807 0.0805 1.39  1.65   1.00    1428.    3074.
11 lambda[1]  0.738  0.737 0.0823 0.0817 0.606 0.877  1.00    3138.    5274.
12 lambda[2]  0.87   0.867 0.0743 0.0743 0.752 0.993  1.00    1903.    3675.
13 lambda[3]  0.802  0.8   0.0777 0.0765 0.678 0.934  1.00    2452.    4289.
14 lambda[4]  0.843  0.839 0.0765 0.0750 0.723 0.971  1.00    2267.    4160.
15 lambda[5]  0.997  0.995 0.0700 0.0687 0.886 1.12   1.00    1393.    3095.
16 lambda[6]  0.899  0.896 0.0637 0.0646 0.800 1.01   1.00    1439.    3092.
17 lambda[7]  0.765  0.764 0.0689 0.0701 0.654 0.879  1.00    2213.    4109.
18 lambda[8]  0.853  0.851 0.0611 0.0612 0.755 0.957  1.00    1443.    2821.
19 lambda[9]  0.861  0.858 0.0784 0.0785 0.736 0.996  1.00    2267.    4043.
20 lambda[10] 0.673  0.672 0.0773 0.0764 0.549 0.800  1.00    3313.    4230.
21 psi[1]     0.891  0.889 0.0494 0.0491 0.814 0.976  1.00   15504.    6185.
22 psi[2]     0.735  0.733 0.0433 0.0436 0.667 0.809  1.00   11583.    5925.
23 psi[3]     0.782  0.78  0.0443 0.0442 0.712 0.858  1.00   13289.    6547.
24 psi[4]     0.757  0.755 0.0445 0.0443 0.687 0.833  1.00   12802.    6653.
25 psi[5]     0.545  0.544 0.0366 0.0369 0.486 0.607  1.00    7541.    6505.
26 psi[6]     0.505  0.503 0.0336 0.0334 0.452 0.563  1.00    9362.    6500.
27 psi[7]     0.686  0.685 0.0404 0.0401 0.623 0.755  1.00   14624.    6491.
28 psi[8]     0.48   0.479 0.0318 0.0315 0.430 0.534  1.00    9793.    7087.
29 psi[9]     0.781  0.779 0.0458 0.0455 0.708 0.858  1.00   11809.    7030.
30 psi[10]    0.839  0.837 0.0467 0.0450 0.766 0.920  1.00   15303.    5988.

2.14 Comparing Results with different inits

Correlation of all parameters across both algorithm runs:

cor(modelCFA_samplesFix$summary(variables = c("mu", "lambda", "psi", "theta"), .cores = 4)$mean,
    modelCFA_samples$summary(variables = c("mu", "lambda", "psi", "theta"), .cores = 4)$mean)

[1] 0.9999963

2.15 CFA model with blavaan

blavaan.model <- ' theta  =~ item1 + item2 + item3 + item4 + item5 + item6 + item7 + item8 + item9 + item10 '
model01_blv <- bcfa(blavaan.model, data=conspiracyItems,
            n.chains = 4, burnin = 1000, sample = 2000,
            target = 'stan', seed = 09102022,
            save.lvs = TRUE, # save sampled latent variable
            std.lv = TRUE,
            bcontrol = list(cores = 4),
            mcmcfile=TRUE # save Stan file
            ) # standardized latent variable

---

2.15.1 blavaan Results

summary(model01_blv)

blavaan 0.5.3 ended normally after 2000 iterations

  Estimator                                      BAYES
  Optimization method                             MCMC
  Number of model parameters                        20

  Number of observations                           177

  Statistic                                 MargLogLik         PPP
  Value                                      -2153.697       0.000

Parameter Estimates:


Latent Variables:
                   Estimate  Post.SD pi.lower pi.upper     Rhat    Prior       
  theta =~                                                                     
    item1             0.738    0.082    0.581    0.902    1.000    normal(0,10)
    item2             0.869    0.077    0.726    1.027    1.001    normal(0,10)
    item3             0.802    0.076    0.659    0.955    1.001    normal(0,10)
    item4             0.842    0.078    0.696    1.002    1.001    normal(0,10)
    item5             0.997    0.071    0.867    1.144    1.001    normal(0,10)
    item6             0.899    0.064    0.780    1.033    1.001    normal(0,10)
    item7             0.764    0.069    0.634    0.911    1.001    normal(0,10)
    item8             0.852    0.061    0.739    0.976    1.001    normal(0,10)
    item9             0.859    0.078    0.709    1.018    1.001    normal(0,10)
    item10            0.673    0.076    0.528    0.826    1.000    normal(0,10)

Variances:
                   Estimate  Post.SD pi.lower pi.upper     Rhat    Prior       
   .item1             0.791    0.088    0.636    0.976    1.000 gamma(1,.5)[sd]
   .item2             0.538    0.064    0.426    0.677    1.000 gamma(1,.5)[sd]
   .item3             0.608    0.071    0.481    0.760    1.000 gamma(1,.5)[sd]
   .item4             0.571    0.068    0.454    0.716    1.000 gamma(1,.5)[sd]
   .item5             0.296    0.040    0.224    0.385    1.000 gamma(1,.5)[sd]
   .item6             0.254    0.034    0.193    0.327    1.000 gamma(1,.5)[sd]
   .item7             0.469    0.055    0.374    0.587    1.000 gamma(1,.5)[sd]
   .item8             0.230    0.030    0.177    0.293    1.000 gamma(1,.5)[sd]
   .item9             0.609    0.073    0.484    0.767    1.000 gamma(1,.5)[sd]
   .item10            0.701    0.079    0.562    0.871    1.000 gamma(1,.5)[sd]
    theta             1.000                                                    

---

2.15.2 Comparing factors score

---

2.15.3 Comparing Parameters

---

2.15.4 Limitation of blavaan 🧛🏼

You can find more details in the Edgar Merkle’s paper.

• Very unreadable Stan file (1891 lines)

• limits in setting up prior distributions (cannot change distribution family for lambda and psi)

• Use rstan rather than cmdstanr

/* 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[dims(covmat)[1]] matrix[dims(covmat)[2], dims(covmat)[3]] out = covmat;
    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];
      }
    }
    return out;
  }
}
data {
  // see https://books.google.com/books?id=9AC-s50RjacC&lpg=PP1&dq=LISREL&pg=PA2#v=onepage&q=LISREL&f=false
  int<lower=0> p; // number of manifest response variables
  int<lower=0> p_c; // number of manifest level 2 variables
  int<lower=0> q; // number of manifest predictors
  int<lower=0> m; // number of latent endogenous variables
  int<lower=0> m_c; // number of latent level 2 variables
  int<lower=0> n; // number of latent exogenous variables
  int<lower=1> Ng; // number of groups
  int<lower=0, upper=1> missing; // are there missing values?
  int<lower=0, upper=1> save_lvs; // should we save lvs?
  int<lower=1> Np; // number of group-by-missing patterns combos
  array[Ng] int<lower=1> N; // number of observations per group
  array[Np] int<lower=1> Nobs; // number of observed variables in each missing pattern
  array[Np] int<lower=0> Nordobs; // number of ordinal observed variables in each missing pattern
  array[Np, p + q] int<lower=0> Obsvar; // indexing of observed variables
  int<lower=1> Ntot; // number of observations across all groups
  array[Np] int<lower=1> startrow; // starting row for each missing pattern
  array[Np] int<lower=1,upper=Ntot> endrow; // ending row for each missing pattern
  array[Np] int<lower=1,upper=Ng> grpnum; // group number for each row of data
  int<lower=0,upper=1> wigind; // do any parameters have approx equality constraint ('wiggle')?
  int<lower=0, upper=1> has_data; // are the raw data on y and x available?
  int<lower=0, upper=1> ord; // are there any ordinal variables?
  int<lower=0, upper=1> multilev; // is this a multilevel dataset?
  int<lower=0> Nord; // how many ordinal variables?
  array[Nord] int<lower=0> ordidx; // indexing of ordinal variables
  array[Np, Nord] int<lower=0> OrdObsvar; // indexing of observed ordinal variables in YXo
  int<lower=0> Noent; // how many observed entries of ordinal variables (for data augmentation)
  array[p + q - Nord] int<lower=0> contidx; // indexing of continuous variables
  array[Nord] int<lower=1> nlevs; // how many levels does each ordinal variable have
  array[Ng, 2] int<lower=1> nclus; // number of level 1 + level 2 observations
  int<lower=0> 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<lower=0> Nx; // number of fixed.x variables (within)
  array[Np] int<lower=0> Nx_between; // number of fixed.x variables (between)
  int<lower=0, upper=1> use_cov;
  int<lower=0, upper=1> pri_only;
  int<lower=0> emiter; // number of em iterations for saturated model in ppp (missing data only)
  int<lower=0, upper=1> use_suff; // should we compute likelihood via mvn sufficient stats?
  int<lower=0, upper=1> 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!!
  array[sum(nclus[,2])] int<lower=1> cluster_size; // number of obs per cluster
  array[Ng] int<lower=1> ncluster_sizes; // number of unique cluster sizes
  array[sum(ncluster_sizes)] int<lower=1> cluster_sizes; // unique cluster sizes
  array[sum(ncluster_sizes)] int<lower=1> cluster_size_ns; // number of clusters of each size
  array[Np, multilev ? p_tilde : p + q] int<lower=0> Xvar; // indexing of fixed.x variables (within)
  array[Np, multilev ? p_tilde : p + q] int<lower=0> Xdatvar; // indexing of fixed.x in data (differs from Xvar when missing)
  array[Np, multilev ? p_tilde : p + q] int<lower=0> 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<lower=0> len_w1; // max number of free elements in Lambda_y per grp
  array[Ng] int<lower=0> 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<lower=1> v1; // index  of free elements in Lambda_y
  array[Ng, p + 1] int<lower=1> u1; // index  of free elements in Lambda_y
  array[sum(wg1), 3] int<lower=0> w1skel;
  array[sum(wg1), 2] int<lower=0> lam_y_sign;
  int<lower=0> len_lam_y; // number of free elements minus equality constraints
  array[len_lam_y] real lambda_y_mn; // prior
  array[len_lam_y] real<lower=0> lambda_y_sd;
  // same things but for B
  int<lower=0> len_w4;
  array[Ng] int<lower=0> wg4;
  array[Ng] vector[len_w4] w4;
  array[Ng, len_w4] int<lower=1> v4;
  array[Ng, m + 1] int<lower=1> u4;
  array[sum(wg4), 3] int<lower=0> w4skel;
  array[sum(wg4), 3] int<lower=0> b_sign;
  int<lower=0> len_b;
  array[len_b] real b_mn;
  array[len_b] real<lower=0> b_sd;
  // same things but for diag(Theta)
  int<lower=0> len_w5;
  array[Ng] int<lower=0> wg5;
  array[Ng] vector[len_w5] w5;
  array[Ng, len_w5] int<lower=1> v5;
  array[Ng, p + 1] int<lower=1> u5;
  array[sum(wg5), 3] int<lower=0> w5skel;
  int<lower=0> len_thet_sd;
  array[len_thet_sd] real<lower=0> theta_sd_shape;
  array[len_thet_sd] real<lower=0> theta_sd_rate;
  int<lower=-2, upper=2> theta_pow;
  // same things but for Theta_r
  int<lower=0> len_w7;
  array[Ng] int<lower=0> wg7;
  array[Ng] vector[len_w7] w7;
  array[Ng, len_w7] int<lower=1> v7;
  array[Ng, p + 1] int<lower=1> u7;
  array[sum(wg7), 3] int<lower=0> w7skel;
  int<lower=0> len_thet_r;
  array[len_thet_r] real<lower=0> theta_r_alpha;
  array[len_thet_r] real<lower=0> theta_r_beta;
  // same things but for Psi
  int<lower=0> len_w9;
  array[Ng] int<lower=0> wg9;
  array[Ng] vector[len_w9] w9;
  array[Ng, len_w9] int<lower=1> v9;
  array[Ng, m + 1] int<lower=1> u9;
  array[sum(wg9), 3] int<lower=0> w9skel;
  int<lower=0> len_psi_sd;
  array[len_psi_sd] real<lower=0> psi_sd_shape;
  array[len_psi_sd] real<lower=0> psi_sd_rate;
  int<lower=-2,upper=2> psi_pow;
  // same things but for Psi_r
  int<lower=0> len_w10;
  array[Ng] int<lower=0> wg10;
  array[Ng] vector[len_w10] w10;
  array[Ng, len_w10] int<lower=1> v10;
  array[Ng, m + 1] int<lower=1> u10;
  array[sum(wg10), 3] int<lower=0> w10skel;
  array[sum(wg10), 3] int<lower=0> psi_r_sign;
  int<lower=0> len_psi_r;
  array[len_psi_r] real<lower=0> psi_r_alpha;
  array[len_psi_r] real<lower=0> psi_r_beta;
  // for blocks within Psi_r that receive lkj
  array[5] int<lower=0> nblk;
  array[5] int<lower=3> psidims;
  array[sum(nblk), 7] int<lower=0> blkse;
  int<lower=0> len_w11;
  array[Ng] int<lower=0> wg11;
  array[Ng] vector[len_w11] w11;
  array[Ng, len_w11] int<lower=1> v11;
  array[Ng, m + 1] int<lower=1> u11;
  array[sum(wg11), 3] int<lower=0> w11skel;
  // same things but for Nu
  int<lower=0> len_w13;
  array[Ng] int<lower=0> wg13;
  array[Ng] vector[len_w13] w13;
  array[Ng, len_w13] int<lower=1> v13;
  array[Ng, use_cov ? 1 : p + q + 1] int<lower=1> u13;
  array[sum(wg13), 3] int<lower=0> w13skel;
  int<lower=0> len_nu;
  array[len_nu] real nu_mn;
  array[len_nu] real<lower=0> nu_sd;
  // same things but for Alpha
  int<lower=0> len_w14;
  array[Ng] int<lower=0> wg14;
  array[Ng] vector[len_w14] w14;
  array[Ng, len_w14] int<lower=0> v14;
  array[Ng, use_cov ? 1 : m + n + 1] int<lower=1> u14;
  array[sum(wg14), 3] int<lower=0> w14skel;
  int<lower=0> len_alph;
  array[len_alph] real alpha_mn;
  array[len_alph] real<lower=0> alpha_sd;
  // same things but for Tau
  int<lower=0> len_w15;
  array[Ng] int<lower=0> wg15;
  array[Ng] vector[len_w15] w15;
  array[Ng, len_w15] int<lower=0> v15;
  array[Ng, sum(nlevs) - Nord + 1] int<lower=1> u15;
  array[sum(wg15), 3] int<lower=0> w15skel;
  int<lower=0> len_tau;
  array[len_tau] real tau_mn;
  array[len_tau] real<lower=0> tau_sd;
  // Level 2 matrices start here!!
  // Lambda
  int<lower=0> len_w1_c;
  array[Ng] int<lower=0> wg1_c;
  array[Ng] vector[len_w1_c] w1_c;
  array[Ng, len_w1_c] int<lower=1> v1_c;
  array[Ng, p_c + 1] int<lower=1> u1_c;
  array[sum(wg1_c), 3] int<lower=0> w1skel_c;
  array[sum(wg1_c), 2] int<lower=0> lam_y_sign_c;
  int<lower=0> len_lam_y_c;
  array[len_lam_y_c] real lambda_y_mn_c;
  array[len_lam_y_c] real<lower=0> lambda_y_sd_c;
  // same things but for B
  int<lower=0> len_w4_c;
  array[Ng] int<lower=0> wg4_c;
  array[Ng] vector[len_w4_c] w4_c;
  array[Ng, len_w4_c] int<lower=1> v4_c;
  array[Ng, m_c + 1] int<lower=1> u4_c;
  array[sum(wg4_c), 3] int<lower=0> w4skel_c;
  array[sum(wg4_c), 3] int<lower=0> b_sign_c;
  int<lower=0> len_b_c;
  array[len_b_c] real b_mn_c;
  array[len_b_c] real<lower=0> b_sd_c;
  // same things but for diag(Theta)
  int<lower=0> len_w5_c;
  array[Ng] int<lower=0> wg5_c;
  array[Ng] vector[len_w5_c] w5_c;
  array[Ng, len_w5_c] int<lower=1> v5_c;
  array[Ng, p_c + 1] int<lower=1> u5_c;
  array[sum(wg5_c), 3] int<lower=0> w5skel_c;
  int<lower=0> len_thet_sd_c;
  array[len_thet_sd_c] real<lower=0> theta_sd_shape_c;
  array[len_thet_sd_c] real<lower=0> theta_sd_rate_c;
  int<lower=-2, upper=2> theta_pow_c;
  // same things but for Theta_r
  int<lower=0> len_w7_c;
  array[Ng] int<lower=0> wg7_c;
  array[Ng] vector[len_w7_c] w7_c;
  array[Ng, len_w7_c] int<lower=1> v7_c;
  array[Ng, p_c + 1] int<lower=1> u7_c;
  array[sum(wg7_c), 3] int<lower=0> w7skel_c;
  int<lower=0> len_thet_r_c;
  array[len_thet_r_c] real<lower=0> theta_r_alpha_c;
  array[len_thet_r_c] real<lower=0> theta_r_beta_c;
  // same things but for Psi
  int<lower=0> len_w9_c;
  array[Ng] int<lower=0> wg9_c;
  array[Ng] vector[len_w9_c] w9_c;
  array[Ng, len_w9_c] int<lower=1> v9_c;
  array[Ng, m_c + 1] int<lower=1> u9_c;
  array[sum(wg9_c), 3] int<lower=0> w9skel_c;
  int<lower=0> len_psi_sd_c;
  array[len_psi_sd_c] real<lower=0> psi_sd_shape_c;
  array[len_psi_sd_c] real<lower=0> psi_sd_rate_c;
  int<lower=-2,upper=2> psi_pow_c;
  // same things but for Psi_r
  int<lower=0> len_w10_c;
  array[Ng] int<lower=0> wg10_c;
  array[Ng] vector[len_w10_c] w10_c;
  array[Ng, len_w10_c] int<lower=1> v10_c;
  array[Ng, m_c + 1] int<lower=1> u10_c;
  array[sum(wg10_c), 3] int<lower=0> w10skel_c;
  array[sum(wg10_c), 3] int<lower=0> psi_r_sign_c;
  int<lower=0> len_psi_r_c;
  array[len_psi_r_c] real<lower=0> psi_r_alpha_c;
  array[len_psi_r_c] real<lower=0> psi_r_beta_c;
  // for blocks within Psi_r that receive lkj
  array[5] int<lower=0> nblk_c;
  array[5] int<lower=3> psidims_c;
  array[sum(nblk_c), 7] int<lower=0> blkse_c;
  int<lower=0> len_w11_c;
  array[Ng] int<lower=0> wg11_c;
  array[Ng] vector[len_w11_c] w11_c;
  array[Ng, len_w11_c] int<lower=1> v11_c;
  array[Ng, m_c + 1] int<lower=1> u11_c;
  array[sum(wg11_c), 3] int<lower=0> w11skel_c;
  // same things but for Nu
  int<lower=0> len_w13_c;
  array[Ng] int<lower=0> wg13_c;
  array[Ng] vector[len_w13_c] w13_c;
  array[Ng, len_w13_c] int<lower=1> v13_c;
  array[Ng, p_c + 1] int<lower=1> u13_c;
  array[sum(wg13_c), 3] int<lower=0> w13skel_c;
  int<lower=0> len_nu_c;
  array[len_nu_c] real nu_mn_c;
  array[len_nu_c] real<lower=0> nu_sd_c;
  // same things but for Alpha
  int<lower=0> len_w14_c;
  array[Ng] int<lower=0> wg14_c;
  array[Ng] vector[len_w14_c] w14_c;
  array[Ng, len_w14_c] int<lower=0> v14_c;
  array[Ng, m_c + 1] int<lower=1> u14_c;
  array[sum(wg14_c), 3] int<lower=0> w14skel_c;
  int<lower=0> len_alph_c;
  array[len_alph_c] real alpha_mn_c;
  array[len_alph_c] real<lower=0> 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[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[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<lower = 0> 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,]);
    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,]);
    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;
 }
      }
    }
    // 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;
 }
      }
    }
    // 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<lower=0>[len_free[5]] Theta_sd_free;
  vector<lower=-1,upper=1>[len_free[7]] Theta_r_free; // to use beta prior
  vector<lower=0>[len_free[9]] Psi_sd_free;
  array[nblk[1]] corr_matrix[psidims[1]] Psi_r_mat_1;
  array[nblk[2]] corr_matrix[psidims[2]] Psi_r_mat_2;
  array[nblk[3]] corr_matrix[psidims[3]] Psi_r_mat_3;
  array[nblk[4]] corr_matrix[psidims[4]] Psi_r_mat_4;
  array[nblk[5]] corr_matrix[psidims[5]] Psi_r_mat_5;
  vector<lower=-1,upper=1>[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<lower=0,upper=1>[Noent] z_aug; //augmented ordinal data
  vector[len_free_c[1]] Lambda_y_free_c;
  vector[len_free_c[4]] B_free_c;
  vector<lower=0>[len_free_c[5]] Theta_sd_free_c;
  vector<lower=-1,upper=1>[len_free_c[7]] Theta_r_free_c; // to use beta prior
  vector<lower=0>[len_free_c[9]] Psi_sd_free_c;
  array[nblk_c[1]] corr_matrix[psidims_c[1]] Psi_r_mat_1_c;
  array[nblk_c[2]] corr_matrix[psidims_c[2]] Psi_r_mat_2_c;
  array[nblk_c[3]] corr_matrix[psidims_c[3]] Psi_r_mat_3_c;
  array[nblk_c[4]] corr_matrix[psidims_c[4]] Psi_r_mat_4_c;
  array[nblk_c[5]] corr_matrix[psidims_c[5]] Psi_r_mat_5_c;
  vector<lower=-1,upper=1>[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(nblk) > 0) {
    // we need to define a separate parameter for each dimension of correlation matrix,
    // so we need all these Psi_r_mats
    Psi_r = fill_cov(Psi_r, blkse, nblk, Psi_r_mat_1, Psi_r_mat_2, Psi_r_mat_3, Psi_r_mat_4, Psi_r_mat_5);
  }
  if (sum(nblk_c) > 0) {
    Psi_r_c = fill_cov(Psi_r_c, blkse_c, nblk_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);
  }
  // see https://books.google.com/books?id=9AC-s50RjacC&lpg=PP1&dq=LISREL&pg=PA3#v=onepage&q=LISREL&f=false
  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] = -10 + (Tau[grpnum[patt], (vecpos + 1), 1] + 10) .* z_aug[idxvec];
     tau_jacobian += log(abs(Tau[grpnum[patt], (vecpos + 1), 1] + 10)); // must add log(U) to tau_jacobian
   } else if (YXo[i,obspos] == nlevs[obspos]) {
     YXostar[i,obspos] = Tau[grpnum[patt], vecpos, 1] + (10 - Tau[grpnum[patt], vecpos, 1]) .* z_aug[idxvec];
     tau_jacobian += log(abs(10 - 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 += -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) {
   target += -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[grpidx], xidx, Nx[mm], p + q, logdetSigma_grp[grpidx]), r2 - r1 + 1);
 }
      }
    }
    if (ord) {
      target += tau_jacobian;
    }
  }
  /* prior densities in log-units */
  target += normal_lpdf(Lambda_y_free | lambda_y_primn, lambda_y_sd);
  target += normal_lpdf(B_free | b_primn, b_sd);
  target += normal_lpdf(Nu_free | nu_primn, nu_sd);
  target += normal_lpdf(Alpha_free | alpha_primn, alpha_sd);
  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(nblk) > 0) {
    for (k in 1:sum(nblk)) {
      int blkidx = blkse[k, 6];
      int arrayidx = blkse[k, 5];
      if (arrayidx == 1) {
 target += lkj_corr_lpdf(Psi_r_mat_1[blkidx] | blkse[k,7]);
      } else if (arrayidx == 2) {
 target += lkj_corr_lpdf(Psi_r_mat_2[blkidx] | blkse[k,7]);
      } else if (arrayidx == 3) {
 target += lkj_corr_lpdf(Psi_r_mat_3[blkidx] | blkse[k,7]);
      } else if (arrayidx == 4) {
 target += lkj_corr_lpdf(Psi_r_mat_4[blkidx] | blkse[k,7]);
      } else {
 target += lkj_corr_lpdf(Psi_r_mat_5[blkidx] | blkse[k,7]);
      }
    }
  }
  if (len_free[10] > 0) {
    target += beta_lpdf(.5 * (1 + Psi_r_free) | psi_r_alpha, psi_r_beta) + log(.5) * len_free[10];
  }
  // 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);
  target += beta_lpdf(.5 * (1 + Theta_r_free_c) | theta_r_alpha_c, theta_r_beta_c) + log(.5) * len_free_c[7];
  if (sum(nblk_c) > 0) {
    for (k in 1:sum(nblk_c)) {
      int blkidx = blkse_c[k, 6];
      int arrayidx = blkse_c[k, 5];
      if (arrayidx == 1) {
 target += lkj_corr_lpdf(Psi_r_mat_1_c[blkidx] | blkse_c[k,7]);
      } else if (arrayidx == 2) {
 target += lkj_corr_lpdf(Psi_r_mat_2_c[blkidx] | blkse_c[k,7]);
      } else if (arrayidx == 3) {
 target += lkj_corr_lpdf(Psi_r_mat_3_c[blkidx] | blkse_c[k,7]);
      } else if (arrayidx == 4) {
 target += lkj_corr_lpdf(Psi_r_mat_4_c[blkidx] | blkse_c[k,7]);
      } else {
 target += lkj_corr_lpdf(Psi_r_mat_5_c[blkidx] | blkse_c[k,7]);
      }
    }
  } 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];
  }
}
generated quantities { // these matrices are saved in the output but do not figure into the likelihood
  // see https://books.google.com/books?id=9AC-s50RjacC&lpg=PP1&dq=LISREL&pg=PA34#v=onepage&q=LISREL&f=false
  // sign constraints and correlations
  vector[len_free[1]] ly_sign;
  vector[len_free[4]] bet_sign;
  array[Ng] matrix[m, m] PSmat;
  array[Ng] matrix[m, m] PS;
  vector[len_free[7]] 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;
  array[Ng] matrix[m_c, m_c] PSmat_c;
  array[Ng] matrix[m_c, m_c] PS_c;
  vector[len_free_c[7]] 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<lower=0, upper=1> 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);
  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);
  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) {
    if (m > 0) {
      PSmat[g] = fill_matrix(P_r, Psi_r_skeleton[g], w10skel, g_start10[g,1], g_start10[g,2]) + transpose(fill_matrix(P_r, Psi_r_skeleton[g], w10skel, g_start10[g,1], g_start10[g,2])) - diag_matrix(rep_vector(1, m));
    }
    if (m_c > 0) {
      PSmat_c[g] = fill_matrix(P_r_c, Psi_r_skeleton_c[g], w10skel_c, g_start10_c[g,1], g_start10_c[g,2]) + transpose(fill_matrix(P_r_c, Psi_r_skeleton_c[g], w10skel_c, g_start10_c[g,1], g_start10_c[g,2])) - diag_matrix(rep_vector(1, m_c));
    }
  }
  if (sum(nblk) > 0) {
    PSmat = fill_cov(PSmat, blkse, nblk, Psi_r_mat_1, Psi_r_mat_2, Psi_r_mat_3, Psi_r_mat_4, Psi_r_mat_5);
  }
  if (sum(nblk_c) > 0) {
    PSmat_c = fill_cov(PSmat_c, blkse_c, nblk_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);
  }
  for (g in 1:Ng) {
    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:
  Theta_cov = cor2cov(Theta_r, Theta_sd, num_elements(Theta_r_free), Theta_r_skeleton, w7skel, Ng);
  Theta_var = Theta_sd_free .* Theta_sd_free;
  if (m > 0 && len_free[11] > 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[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
  Theta_cov_c = cor2cov(Theta_r_c, Theta_sd_c, num_elements(Theta_r_free_c), 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);
   if (Nx[mm] > 0) {
     log_lik[jj] += -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);
   }
 }
      }
      // 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[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);
     // log_lik_sat, log_lik_sat_rep
     if (Nx[mm] > 0) {
       log_lik_rep[jj] += -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] += -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] += -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);
     }
   }
   // 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)

3 Wrapping up

Today, we showed how to model observed data using a normal distribution

• Assumptions of Confirmatory Factor Analysis

  • Not appropriate for our data

  • May not be appropriate for many data sets

• We will have to keep our loading/discrimination parameters positive to ensure converges to the same posterior mode

  • This will continue through the next types of data

4 Next Class

• Next up, categorical distributions for observed data

  • More appropriate for these data as they are discrete categorical responses

4.1 Resources

• Dr. Templin’s slide