Below is the specification of an SEM originally specified in WinBUGS. The most notable difference being
for(j in 1:P){
z[i,j] ~ ordered_logistic(mu[i,j], thd[j]);
}
replacing
for(j in 1:P){
y[i,j]~dnorm(mu[i,j],psi[j])I(thd[j,z[i,j]],thd[j,z[i,j]+1])
ephat[i,j]<y[i,j]mu[i,j]
}
and the thresholds being estimated as unknown parameters, with the threshold data as priors.
The I()
function is according to Bard similar to step()
. However, I was not able to make this code run in WinBUGS as is, nor was I able to utilize step()
correctly in the Stan model (see commented lines in the measurement equation model). Anyway, I considered I might not really need the intermediate normal distribution y
, if I can estimate mu
directly from z
using the ordered_logistic()
function. This, however, leaves no way of estimating ephat
.

I was therefore wondering if this is really necessary?

Also, the model fit will yield the following warning message:
"Warning: 2 of 2 chains had an EBFMI less than 0.2."
Is this caused by the way my model is specified? Iāve been unable to get this model to run on WinBUGS, and there are indications this may have to do with with the code snippet above with the indicator function, but it makes it difficult to know whether the warning is caused by my adaptations of the original model, or problems with the data or original model. 
Finally, I am trying to run the model using OpenCL on the GPU. This has turned out to be horrendously slow. I am therefore trying to understand what is causing the slowdown: Is there something inherent in this model that makes CPU computations better suited for it, or is it due to model specification not being optimized for GPU?
cmdstanr code
Stan model
data{
int<lower=0> N; // number of observations
int<lower=1> P; // number of variables
array[N, P] int<lower=1, upper=5> z; // observed data
}
transformed data {
cov_matrix[4] R = diag_matrix(rep_vector(8.0, 4)); // prior covariance matrix for xi
vector[4] u = rep_vector(0,4); // prior mean for xi
}
parameters {
array[21] real lam; // pattern coefficients
vector[N] eta; // latent variables
array[N] vector[4] xi; // latent variables
row_vector[4] gam; // coefficients
vector<lower=0>[P] psi; // precisions of y
real<lower=0> psd; // precision of eta
cov_matrix[4] phi;
array[P] ordered[4] thd; // Thresholds
}
transformed parameters {
vector[P] sgm = 1/psi;
real sgd = 1/psd;
matrix[4, 4] phx = inverse(phi);
}
model {
// Local variables
array[N] vector[P] mu;
array[N] vector[P] ephat;
vector[N] nu;
vector[N] dthat;
array[21] real var_lam;
real var_gam;
//priors on precisions
psi ~ gamma(10,8);
psd ~ gamma(10,8);
phi ~ wishart(30, R);
//priors on loadings and coefficients
var_lam[1] = 4.0 * psi[2];
var_lam[2] = 4.0 * psi[4];
var_lam[3] = 4.0 * psi[5];
var_lam[4] = 4.0 * psi[6];
var_lam[5] = 4.0 * psi[7];
var_lam[6] = 4.0 * psi[8];
var_lam[7] = 4.0 * psi[9];
var_lam[8] = 4.0 * psi[11];
var_lam[9] = 4.0 * psi[12];
var_lam[10] = 4.0 * psi[13];
var_lam[11] = 4.0 * psi[14];
var_lam[12] = 4.0 * psi[15];
var_lam[13] = 4.0 * psi[17];
var_lam[14] = 4.0 * psi[18];
var_lam[15] = 4.0 * psi[20];
var_lam[16] = 4.0 * psi[21];
var_lam[17] = 4.0 * psi[22];
var_lam[18] = 4.0 * psi[23];
var_lam[19] = 4.0 * psi[24];
var_lam[20] = 4.0 * psi[25];
var_lam[21] = 4.0 * psi[26];
lam ~ normal(0.8,var_lam);
var_gam = 4.0 * psd;
gam ~ normal(0.6,var_gam);
for(i in 1:N){
//measurement equation model
mu[i,1] = eta[i];
mu[i,2] = lam[1] * eta[i];
mu[i,3] = xi[i,1];
mu[i,4] = lam[2] * xi[i,1];
mu[i,5] = lam[3] * xi[i,1];
mu[i,6] = lam[4] * xi[i,1];
mu[i,7] = lam[5] * xi[i,1];
mu[i,8] = lam[6] * xi[i,1];
mu[i,9] = lam[7] * xi[i,1];
mu[i,10] = xi[i,2];
mu[i,11] = lam[8] * xi[i,2];
mu[i,12] = lam[9] * xi[i,2];
mu[i,13] = lam[10] * xi[i,2];
mu[i,14] = lam[11] * xi[i,2];
mu[i,15] = lam[12] * xi[i,2];
mu[i,16] = xi[i,3];
mu[i,17] = lam[13] * xi[i,3];
mu[i,18] = lam[14] * xi[i,3];
mu[i,19] = xi[i,4];
mu[i,20] = lam[15] * xi[i,4];
mu[i,21] = lam[16] * xi[i,4];
mu[i,22] = lam[17] * xi[i,4];
mu[i,23] = lam[18] * xi[i,4];
mu[i,24] = lam[19] * xi[i,4];
mu[i,25] = lam[20] * xi[i,4];
mu[i,26] = lam[21] * xi[i,4];
for(j in 1:P){
// y[i,j] ~ normal(mu[i,j], psi[j]) * step(thd[j, z[i,j]]  y[i,j]);
z[i,j] ~ ordered_logistic(mu[i,j], thd[j]);
}
// ephat[i] = y[i]  mu[i];
// structural equation model
nu[i] = gam * xi[i];
} // end of i
dthat = eta  nu;
xi ~ multi_normal(u, phi);
eta ~ normal(nu, psd);
} //end of model
R code
library(cmdstanr)
library(dplyr)
library(glue)
output_dir < glue::glue("{getwd()}/Chapter6/stanoutput")
if( !dir.exists(output_dir) ) dir.create(output_dir)
model2 < paste0(getwd(), "/Chapter6/ch6Stanmodelwithunknownthresholds.stan")
model2 < cmdstan_model(model2,
force_recompile = TRUE
)
data < list(
N = 338,
P = 26,
z = read.csv(
"./Chapter6/ch6WinBUGSdata.dat",
header = FALSE,
nrows = 338,
skip = 2,
strip.white = TRUE
) %>%
.[,1:(ncol(.)  1)] %>%
as.matrix() %>%
unname()
)
Ļ_priors < read.csv(
"./Chapter6/ch6WinBUGSdata.dat",
header = FALSE,
# nrows = 338,
skip = 342,
strip.white = TRUE
) %>%
.[,1:(ncol(.)  1)] %>%
as.matrix() %>%
unname()
thd < matrix(
c(
2.517,1.245,0.444, 0.848,
1.447,0.420, 0.119, 1.245,
1.671,0.869,0.194, 0.679,
1.642,0.869,0.293, 0.332,
1.671,0.827, 0.052, 0.756,
1.769,1.098,0.469, 0.255,
1.490,0.670,0.082, 0.880,
1.933,0.880,0.317, 1.008,
1.587,0.624, 0.000, 1.008,
1.983,1.348,0.348, 1.045,
1.983,1.229,0.247, 0.869,
2.262,1.426, 0.037, 1.330,
2.371,1.295,0.224, 0.651,
2.039,1.112,0.149, 1.169,
2.262,1.198,0.309, 1.198,
2.176,1.537,0.717, 0.597,
1.447,0.786, 0.119, 1.008,
2.039,1.769,0.661, 0.642,
2.262,1.468, 0.015, 1.214,
2.039,1.406, 0.000, 1.140,
1.702,1.058, 0.149, 0.902,
2.262,1.426,0.309, 0.971,
1.702,0.615, 0.179, 1.229,
2.262,1.671,1.033, 0.420,
2.262,1.468,0.689, 1.045,
2.176,1.537,0.880, 0.661
),
nrow = 26,
ncol = 4,
byrow = TRUE
)
inits < list(
list(
lam = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
psi = c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
psd = 1.0,
gam = c(1.0, 1.0, 1.0, 1.0),
phi = structure(
.Data = c(1.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0),
.Dim = c(4,4)),
xi = Ļ_priors,
thd = thd
),
list(
lam = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5),
psi = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5),
psd = 0.6,
gam = c(0.0, 0.0, 0.0, 0.0),
phi = structure(
.Data = c(0.5, 0.0, 0.0, 0.0,
0.0, 0.5, 0.0, 0.0,
0.0, 0.0, 0.5, 0.0,
0.0, 0.0, 0.0, 0.5),
.Dim = c(4,4)),
xi = Ļ_priors,
thd = thd
)
)
fit2 < model2$sample(
data = data,
seed = 69,
refresh = 500,
init = inits,
output_dir = output_dir,
sig_figs = 5,
chains = n.chains,
parallel_chains = getOption("mc.cores", 1),
opencl_ids = NULL,
iter_warmup = n.burnin,
iter_sampling = n.iter,
save_warmup = TRUE
)
data.rdata (4.5 KB)
WinBUGS model
This is a WinBUGS program for the real example in Chapter 6, Section 6.6.2.
Model: Structural Equation Model with Ordered Categorical Variables
Data Set Names: YO.dat, and XI.dat, where XI.dat are input initial values for xi.
Sample Size: N=338
model{
for(i in 1:N){
#measurement equation model
for(j in 1:P){
y[i,j]~dnorm(mu[i,j],psi[j])I(thd[j,z[i,j]],thd[j,z[i,j]+1])
ephat[i,j]<y[i,j]mu[i,j]
}
mu[i,1]<eta[i]
mu[i,2]<lam[1]*eta[i]
mu[i,3]<xi[i,1]
mu[i,4]<lam[2]*xi[i,1]
mu[i,5]<lam[3]*xi[i,1]
mu[i,6]<lam[4]*xi[i,1]
mu[i,7]<lam[5]*xi[i,1]
mu[i,8]<lam[6]*xi[i,1]
mu[i,9]<lam[7]*xi[i,1]
mu[i,10]<xi[i,2]
mu[i,11]<lam[8]*xi[i,2]
mu[i,12]<lam[9]*xi[i,2]
mu[i,13]<lam[10]*xi[i,2]
mu[i,14]<lam[11]*xi[i,2]
mu[i,15]<lam[12]*xi[i,2]
mu[i,16]<xi[i,3]
mu[i,17]<lam[13]*xi[i,3]
mu[i,18]<lam[14]*xi[i,3]
mu[i,19]<xi[i,4]
mu[i,20]<lam[15]*xi[i,4]
mu[i,21]<lam[16]*xi[i,4]
mu[i,22]<lam[17]*xi[i,4]
mu[i,23]<lam[18]*xi[i,4]
mu[i,24]<lam[19]*xi[i,4]
mu[i,25]<lam[20]*xi[i,4]
mu[i,26]<lam[21]*xi[i,4]
#structural equation model
xi[i,1:4]~dmnorm(u[1:4],phi[1:4,1:4])
eta[i]~dnorm(nu[i],psd)
nu[i]<gam[1]*xi[i,1]+gam[2]*xi[i,2]+gam[3]*xi[i,3]+gam[4]*xi[i,4]
dthat[i]<eta[i]nu[i]
}# end of i
for(i in 1:4){u[i]<0.0}
#priors on loadings and coefficients
var.lam[1]<4.0*psi[2] var.lam[2]<4.0*psi[4] var.lam[3]<4.0*psi[5]
var.lam[4]<4.0*psi[6] var.lam[5]<4.0*psi[7] var.lam[6]<4.0*psi[8]
var.lam[7]<4.0*psi[9] var.lam[8]<4.0*psi[11] var.lam[9]<4.0*psi[12]
var.lam[10]<4.0*psi[13] var.lam[11]<4.0*psi[14] var.lam[12]<4.0*psi[15]
var.lam[13]<4.0*psi[17] var.lam[14]<4.0*psi[18] var.lam[15]<4.0*psi[20]
var.lam[16]<4.0*psi[21] var.lam[17]<4.0*psi[22] var.lam[18]<4.0*psi[23]
var.lam[19]<4.0*psi[24] var.lam[20]<4.0*psi[25] var.lam[21]<4.0*psi[26]
for(i in 1:21){lam[i]~dnorm(0.8,var.lam[i])}
var.gam<4.0*psd
gam[1]~dnorm(0.6,var.gam) gam[2]~dnorm(0.6,var.gam)
gam[3]~dnorm(0.4,var.gam) gam[4]~dnorm(0.4,var.gam)
#priors on precisions
for(j in 1:P){
psi[j]~dgamma(10,8)
sgm[j]<1/psi[j]
}
psd~dgamma(10,8)
sgd<1/psd
phi[1:4,1:4]~dwish(R[1:4,1:4], 30)
phx[1:4,1:4]<inverse(phi[1:4,1:4])
} #end of model
Data Set
list(N=338, P=26,
R=structure(
.Data=c(8.0, 0.0, 0.0, 0.0,
0.0, 8.0, 0.0, 0.0,
0.0, 0.0, 8.0, 0.0,
0.0, 0.0, 0.0, 8.0),
.Dim=c(4,4)),
thd=structure(
.Data=c(200.000,2.517,1.245,0.444, 0.848,200.000,
200.000,1.447,0.420, 0.119, 1.245,200.000,
200.000,1.671,0.869,0.194, 0.679,200.000,
200.000,1.642,0.869,0.293, 0.332,200.000,
200.000,1.671,0.827, 0.052, 0.756,200.000,
200.000,1.769,1.098,0.469, 0.255,200.000,
200.000,1.490,0.670,0.082, 0.880,200.000,
200.000,1.933,0.880,0.317, 1.008,200.000,
200.000,1.587,0.624, 0.000, 1.008,200.000,
200.000,1.983,1.348,0.348, 1.045,200.000,
200.000,1.983,1.229,0.247, 0.869,200.000,
200.000,2.262,1.426, 0.037, 1.330,200.000,
200.000,2.371,1.295,0.224, 0.651,200.000,
200.000,2.039,1.112,0.149, 1.169,200.000,
200.000,2.262,1.198,0.309, 1.198,200.000,
200.000,2.176,1.537,0.717, 0.597,200.000,
200.000,1.447,0.786, 0.119, 1.008,200.000,
200.000,2.039,1.769,0.661, 0.642,200.000,
200.000,2.262,1.468, 0.015, 1.214,200.000,
200.000,2.039,1.406, 0.000, 1.140,200.000,
200.000,1.702,1.058, 0.149, 0.902,200.000,
200.000,2.262,1.426,0.309, 0.971,200.000,
200.000,1.702,0.615, 0.179, 1.229,200.000,
200.000,2.262,1.671,1.033, 0.420,200.000,
200.000,2.262,1.468,0.689, 1.045,200.000,
200.000,2.176,1.537,0.880, 0.661,200.000),
.Dim=c(26,6)),
z=structure(
.Data=c(paste YO.dat here),
.Dim=c(338,26)))
Two different Initial Values
list(
lam=c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
psi=c(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
psd=1.0,
gam=c(1.0, 1.0, 1.0, 1.0),
phi=structure(
.Data=c(1.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 1.0),
.Dim=c(4,4)),
xi=structure(
.Data=c(paste XI.dat here),
.Dim=c(338,4)))
list(
lam=c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5),
psi=c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5),
psd=0.6,
gam=c(0.0, 0.0, 0.0, 0.0),
phi=structure(
.Data=c(0.5, 0.0, 0.0, 0.0,
0.0, 0.5, 0.0, 0.0,
0.0, 0.0, 0.5, 0.0,
0.0, 0.0, 0.0, 0.5),
.Dim=c(4,4)),
xi=structure(
.Data=c(paste XI.dat here),
.Dim=c(338,4)))