8 Models for Cross-Classified Multiple Membership Data
Models of section 8.2.1
Cross-Classified-CTCM-1-Model
R Code for jags
#### Loading of data set
load("chapter_8.RData")#
#### Loading packages
library(R2jags)
library(coda)
library(mcmcplots)
#### prepare data for jags
dat<-chapter_8[order(chapter_8$Target,chapter_8$Rater),]
Target <- as.numeric(as.factor(dat$Target))
Rater <- as.numeric(as.factor(dat$Rater))
Y111 <- dat$Y111
Y211 <- dat$Y211
Y311 <- dat$Y311
Y121 <- dat$Y121
Y221 <- dat$Y221
Y321 <- dat$Y321
Y112 <- dat$Y112
Y212 <- dat$Y212
Y312 <- dat$Y312
Y122 <- dat$Y122
Y222 <- dat$Y222
Y322 <- dat$Y322
# create indices for loops in c4model.bug
id.c <- table(dat$Target)
index <- cumsum(as.vector(id.c))
# length(unique(Target))
# length(unique(Rater))
N_prof <- 170 # Number of students
N_stud <- 150 # Number of teacher
n<-length(dat$Y111) # Number of total observations
M <- diag(2) # identity matrix
# data set for jags
dat<-list("Target", "Rater",
"Y111", "Y211","Y311","Y121", "Y221", "Y321",
"Y112", "Y212", "Y312","Y122", "Y222", "Y322",
"n", "N_prof","N_stud", "index", "M")
# parameters that should be sampled by jags
parameters<-c("lamT1","lamT2","mu",
"epsL11", "epsL12", "epsL21", "epsL22",
"var_Trait1","var_Trait2","var_R1", "var_R2",
"var_UR1", "var_UR2", "var_Inter1", "var_Inter2",
"cor.T", "cor.R", "cor.UR", "cor.I",
"con11", "rsp11", "con21", "rsp21", "ind1", "icc1",
"rel11", "rel21",
"con12", "rsp12", "con22", "rsp22", "ind2", "icc2",
"rel12", "rel22")
# create BUGS model
# C4 model for structurally different and non-independent interchangeable raters
# the set of interchangeable raters servers as reference rater group
# 3 indicators, 2 trait, 2 methods
# n = number of total observations i (rater times targets)
# N_stud = number of targets t (here students)
# N_prof = number of raters r (here professors)
# M = identity matrix
# index = index running over unique professors' self-ratings
# Y111 - Y211 - Y311 student ratings of teaching quality
# Y121 - Y221 - Y321 student ratings of engagement
# Y112 - Y212 - Y312 self-reported teaching quality
# Y122 - Y222 - Y322 self-reported engagement
c4model <- "
model{
for(i in 1:n)
{
# Measurement model for non-independent interchangeable raters
# level-1: the interaction-level
Y111[i]~dnorm(tau_Y111[i], theta.Y11)
Y211[i]~dnorm(tau_Y211[i], theta.Y11)
Y311[i]~dnorm(tau_Y311[i], theta.Y11)
Y121[i]~dnorm(tau_Y121[i], theta.Y21)
Y221[i]~dnorm(tau_Y221[i], theta.Y21)
Y321[i]~dnorm(tau_Y321[i], theta.Y21)
tau_Y111[i]<-Trait1[Target[i],1]+URater1[Rater[i],1]+Inter[i,1]
tau_Y211[i]<-Trait2[Target[i],1]+URater2[Rater[i],1]+Inter[i,1]
tau_Y311[i]<-Trait3[Target[i],1]+URater3[Rater[i],1]+Inter[i,1]
tau_Y121[i]<-Trait4[Target[i],2]+URater4[Rater[i],2]+Inter[i,2]
tau_Y221[i]<-Trait5[Target[i],2]+URater5[Rater[i],2]+Inter[i,2]
tau_Y321[i]<-Trait6[Target[i],2]+URater6[Rater[i],2]+Inter[i,2]
Inter[i, 1:2] ~ dmnorm(mu_Inter, Omega_Inter)
}
for(t in 1:N_prof)
{
# Measurement model for structurally different raters
# level-2a: target level
Y112[index[t]]~dnorm(tau_Y112[t], theta.Y12)
Y212[index[t]]~dnorm(tau_Y212[t], theta.Y12)
Y312[index[t]]~dnorm(tau_Y312[t], theta.Y12)
Y122[index[t]]~dnorm(tau_Y122[t], theta.Y22)
Y222[index[t]]~dnorm(tau_Y222[t], theta.Y22)
Y322[index[t]]~dnorm(tau_Y322[t], theta.Y22)
tau_Y112[t]<-lamT1*Trait[t,1] + CRater[t,1] + mu
tau_Y212[t]<-lamT1*Trait[t,1] + CRater[t,1] + mu
tau_Y312[t]<-lamT1*Trait[t,1] + CRater[t,1] + mu
tau_Y122[t]<-lamT2*Trait[t,2] + CRater[t,2] + mu
tau_Y222[t]<-lamT2*Trait[t,2] + CRater[t,2] + mu
tau_Y322[t]<-lamT2*Trait[t,2] + CRater[t,2] + mu
# Latent regression analysis
Trait1[t,1]<-Trait[t,1] + mu
Trait2[t,1]<-Trait[t,1] + mu
Trait3[t,1]<-Trait[t,1] + mu
Trait4[t,2]<-Trait[t,2] + mu
Trait5[t,2]<-Trait[t,2] + mu
Trait6[t,2]<-Trait[t,2] + mu
Trait[t, 1:2] ~ dmnorm(mu_Trait, Omega_Trait)
CRater[t, 1:2] ~ dmnorm(mu_CRater, Omega_CRater)
}
for(r in 1:N_stud)
{
# level-2b: rater-level
URater1[r,1]<-URater[r,1]
URater2[r,1]<-URater[r,1]
URater3[r,1]<-URater[r,1]
URater4[r,2]<-URater[r,2]
URater5[r,2]<-URater[r,2]
URater6[r,2]<-URater[r,2]
URater[r, 1:2] ~ dmnorm(mu_URater, Omega_URater)
}
# Hyper-Priors
mu_Inter[1] <- 0
mu_Inter[2] <- 0
Omega_Inter ~ dwish(M,2)
Sigma_Inter <- inverse(Omega_Inter)
mu_Trait[1] <- 0
mu_Trait[2] <- 0
Omega_Trait ~ dwish(M,2)
Sigma_Trait <- inverse(Omega_Trait)
mu_CRater[1] <- 0
mu_CRater[2] <- 0
Omega_CRater ~ dwish(M,2)
Sigma_CRater <- inverse(Omega_CRater)
mu_URater[1] <- 0
mu_URater[2] <- 0
Omega_URater ~ dwish(M,2)
Sigma_URater <- inverse(Omega_URater)
# Loadings
lamT1~dnorm(1,10)
lamT2~dnorm(1,10)
# Intercepts
# because the manifest variables were simulated with mean of zero (centered),
# we used a strong prior around zero. Change this accordingly in applications!
mu~dnorm(0,0.001)
# Error variables
# Due to equality restrictions,
# we used the same labels for some error variances
# change this to freely estimate all variances individually
theta.Y11 ~ dgamma(.001,.001)
theta.Y21 ~ dgamma(.001,.001)
theta.Y12 ~ dgamma(.001,.001)
theta.Y22 ~ dgamma(.001,.001)
# Error variances
# inverse gamma prior on factor variance
epsL11<-1/theta.Y11
epsL12<-1/theta.Y21
epsL21<-1/theta.Y12
epsL22<-1/theta.Y22
# Factor variances
var_Trait1<-Sigma_Trait[1,1]
var_Trait2<-Sigma_Trait[2,2]
var_R1<-Sigma_CRater[1,1]
var_R2<-Sigma_CRater[2,2]
var_UR1<-Sigma_URater[1,1]
var_UR2<-Sigma_URater[2,2]
var_Inter1<-Sigma_Inter[1,1]
var_Inter2<-Sigma_Inter[2,2]
# Latent correlations
cor.T<-Sigma_Trait[2,1]/(sqrt(var_Trait1)*sqrt(var_Trait2))
cor.R<-Sigma_CRater[2,1]/(sqrt(var_R1)*sqrt(var_R2))
cor.UR<-Sigma_URater[2,1]/(sqrt(var_UR1)*sqrt(var_UR2))
cor.I<-Sigma_Inter[2,1]/(sqrt(var_Inter1)*sqrt(var_Inter2))
# variance coefficients
con11 = var_Trait1/(var_Trait1+var_UR1+var_Inter1)
rsp11 = var_UR1/(var_Trait1+var_UR1+var_Inter1)
con21 = lamT1^2*var_Trait1/(lamT1^2*var_Trait1+var_R1)
rsp21 = var_R1/(lamT1^2*var_Trait1)
ind1 = var_Inter1/(var_Trait1+var_UR1+var_Inter1)
icc1 = 1-ind1
rel11 = (var_Trait1+var_UR1+var_Inter1)/ (var_Trait1+var_UR1+var_Inter1+epsL11)
rel21 = (lamT1^2*var_Trait1+var_R1)/(lamT1^2*var_Trait1+var_R1+epsL21)
con12 = var_Trait2/(var_Trait2+var_UR2+var_Inter2);
rsp12 = var_UR2/(var_Trait2+var_UR2+var_Inter2);
con22 = lamT2^2*var_Trait2/(lamT2^2*var_Trait2+var_R2);
rsp22 = var_R2/(lamT2^2*var_Trait2+var_R2);
ind2 = var_Inter2/(var_Trait2+var_UR2+var_Inter2);
icc2 = 1-ind2;
rel12 = (var_Trait2+var_UR2+var_Inter2)/ (var_Trait2+var_UR2+var_Inter2+epsL12);
rel22 = (lamT2^2*var_Trait2+var_R2)/(lamT2^2*var_Trait2+var_R2+epsL22);
}
"
writeLines(c4model, con = "c4model.bug")
# run jags
set.seed(1234)
c4model<-jags(data=dat,inits=NULL, parameters.to.save = parameters,
n.chains=3, n.iter=7500, n.burnin=2500, n.thin = 5,
model.file="c4model.bug")
# parallel computing
c4model<-jags.parallel(dat, inits=NULL, parameters,n.chains=3,
n.iter=7500, n.burnin=2500, n.thin = 5,
model.file="c4model.bug")
print(c4model)
plot(c4model)
mcmcplot(c4model)
effectiveSize(c4model)
# Traceplots
traplot(c4model, parms=parameters)
# Inspect densities
denplot(c4model, parms=parameters)R Code for rstan
#### Loading packages
library(rstan)
library(tidyverse)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#### prepare data for stan
dat <- chapter_8 %>%
as_tibble %>%
arrange(Target, Rater) %>%
mutate(Rater = as.integer(as.factor(Rater)))
dat_standata <- list(
N_targets = n_distinct(dat$Target),
N_raters = n_distinct(dat$Rater),
N_interactions = n_distinct(sprintf("%i_%i", dat$Target, dat$Rater)),
K_strdiff_1 = 3,
K_strdiff_2 = 3,
K_interch_1 = 3,
K_interch_2 = 3,
targets = dat$Target,
raters = dat$Rater,
Y_interch_1 = dat %>% select(ends_with("11")) %>% as.list,
Y_interch_2 = dat %>% select(ends_with("21")) %>% as.list,
Y_strdiff_1 = dat %>%
group_by(Target) %>%
summarize(across(ends_with("2"), unique)) %>%
select(ends_with("12")) %>%
as.list,
Y_strdiff_2 = dat %>%
group_by(Target) %>%
summarize(across(ends_with("2"), unique)) %>%
select(ends_with("22")) %>%
as.list
)
c4model_stan <- "
data {
int<lower=1> N_targets;
int<lower=1> N_raters;
int<lower=1> N_interactions;
int<lower=1> K_strdiff_1;
int<lower=1> K_strdiff_2;
int<lower=1> K_interch_1;
int<lower=1> K_interch_2;
int<lower=1, upper=N_targets> targets[N_interactions]; // Target indices
int<lower=1, upper=N_raters> raters[N_interactions]; // Rater indices
vector[N_targets] Y_strdiff_1[K_strdiff_1];
vector[N_targets] Y_strdiff_2[K_strdiff_2];
vector[N_interactions] Y_interch_1[K_interch_1];
vector[N_interactions] Y_interch_2[K_interch_2];
}
transformed data {
vector[2] zero_2 = rep_vector(0.0, 2);
}
parameters {
vector<lower=0.0>[2] sd_T;
corr_matrix[2] rho_T;
vector<lower=0.0>[2] sd_R1;
corr_matrix[2] rho_R1;
vector<lower=0.0>[2] sd_R2;
corr_matrix[2] rho_R2;
vector<lower=0.0>[2] sd_I;
corr_matrix[2] rho_I;
vector<lower=0.0>[2] sd_residual_strdiff;
vector<lower=0.0>[2] sd_residual_interch;
array[2] vector[N_targets] Targ;
array[2] vector[N_raters] R1;
array[2] vector[N_targets] R2;
array[2] vector[N_interactions] I;
vector[K_interch_1] intercepts_1;
vector[K_interch_2] intercepts_2;
vector<lower=0.0>[2] loading_T;
}
transformed parameters {
matrix[2, 2] cov_T = quad_form_diag(rho_T, sd_T);
matrix[2, 2] cov_R1 = quad_form_diag(rho_R1, sd_R1);
matrix[2, 2] cov_R2 = quad_form_diag(rho_R2, sd_R2);
matrix[2, 2] cov_I = quad_form_diag(rho_I, sd_I);
}
model {
array[2] vector[N_targets] conditional_means_strdiff;
array[2] vector[N_interactions] conditional_means_interch;
for (i in 1:2) {
conditional_means_strdiff[i] = loading_T[i] * Targ[i] + R2[i];
conditional_means_interch[i] = Targ[i][targets] + R1[i][raters] + I[i];
}
for (t in 1:N_targets) [Targ[1][t], Targ[2][t]]' ~ multi_normal(zero_2, cov_T);
for (r in 1:N_raters) [R1[1][r], R1[2][r]]' ~ multi_normal(zero_2, cov_R1);
for (r in 1:N_targets) [R2[1][r], R2[2][r]]' ~ multi_normal(zero_2, cov_R2);
for (i in 1:N_interactions) [I[1][i], I[2][i]]' ~ multi_normal(zero_2, cov_I);
// sd_T ~ lognormal(0, 1);
// rho_T ~ lkj_corr(2);
// sd_R1 ~ lognormal(0, 1);
// rho_R1 ~ lkj_corr(2);
// sd_R2 ~ lognormal(0, 1);
// rho_R2 ~ lkj_corr(2);
// sd_I ~ lognormal(0, 1);
// rho_I ~ lkj_corr(2);
//
// loading_T ~ lognormal(0, 1);
//
// sd_residual_strdiff ~ lognormal(0, 1);
// sd_residual_interch ~ lognormal(0, 1);
//
// intercepts_1 ~ normal(0, 10);
// intercepts_2 ~ normal(0, 10);
for (k in 1:K_strdiff_1) Y_strdiff_1[k] ~ normal(intercepts_1[k] + conditional_means_strdiff[1], sd_residual_strdiff[1]);
for (k in 1:K_strdiff_2) Y_strdiff_2[k] ~ normal(intercepts_2[k] + conditional_means_strdiff[2], sd_residual_strdiff[2]);
for (k in 1:K_interch_1) Y_interch_1[k] ~ normal(conditional_means_interch[1], sd_residual_interch[1]);
for (k in 1:K_interch_2) Y_interch_2[k] ~ normal(conditional_means_interch[2], sd_residual_interch[2]);
}
generated quantities {
vector<lower=0.0>[2] sq_loading_T = square(loading_T);
vector<lower=0.0>[2] var_residual_strdiff = square(sd_residual_strdiff);
vector<lower=0.0>[2] var_residual_interch = square(sd_residual_interch);
array[2] real<lower=0.0, upper=1.0> con1;
array[2] real<lower=0.0, upper=1.0> con2;
array[2] real<lower=0.0, upper=1.0> rsp1;
array[2] real<lower=0.0, upper=1.0> rsp2;
array[2] real<lower=0.0, upper=1.0> ind;
array[2] real<lower=0.0, upper=1.0> icc;
array[2] real<lower=0.0, upper=1.0> gen;
array[2] real<lower=0.0, upper=1.0> rel1;
array[2] real<lower=0.0, upper=1.0> rel2;
for (i in 1:2) con1[i] = cov_T[i, i] / (cov_T[i, i] + cov_R1[i, i] + cov_I[i, i]);
for (i in 1:2) con2[i] = sq_loading_T[i] * cov_T[i, i] / (sq_loading_T[i] * cov_T[i, i] + cov_R2[i, i]);
for (i in 1:2) rsp1[i] = cov_R1[i, i] / (cov_T[i, i] + cov_R1[i, i] + cov_I[i, i]);
for (i in 1:2) rsp2[i] = cov_R2[i, i] / (sq_loading_T[i] * cov_T[i, i] + cov_R2[i, i]);
for (i in 1:2) ind[i] = cov_I[i, i] / (cov_T[i, i] + cov_R1[i, i] + cov_I[i, i]);
for (i in 1:2) icc[i] = 1 - ind[i];
for (i in 1:2) gen[i] = cov_I[i, i] / (cov_R1[i, i] + cov_I[i, i]);
for (i in 1:2) rel1[i] = (cov_T[i, i] + cov_R1[i, i] + cov_I[i, i]) / (cov_T[i, i] + cov_R1[i, i] + cov_I[i, i] + var_residual_interch[i]);
for (i in 1:2) rel2[i] = (sq_loading_T[i] * cov_T[i, i] + cov_R2[i, i]) / (sq_loading_T[i] * cov_T[i, i] + cov_R2[i, i] + var_residual_strdiff[i]);
// Compute correlations based on covariances
real corr_T = cov_T[1, 2] / sqrt(cov_T[1, 1] * cov_T[2, 2]);
real corr_R1 = cov_R1[1, 2] / sqrt(cov_R1[1, 1] * cov_R1[2, 2]);
real corr_R2 = cov_R2[1, 2] / sqrt(cov_R2[1, 1] * cov_R2[2, 2]);
real corr_I = cov_I[1, 2] / sqrt(cov_I[1, 1] * cov_I[2, 2]);
}
"
writeLines(c4model_stan, con = "c4model.stan")
#### fit c4model with stan
c4model_stan <- stan(
file = "c4model.stan",
data = dat_standata,
pars = c("sd_T", "sd_R1", "sd_R2", "sd_I",
"intercepts_1", "intercepts_2", "loading_T",
"var_residual_strdiff", "var_residual_interch",
"corr_T", "corr_R1", "corr_R2", "corr_I",
"con1", "con2", "rsp1", "rsp2",
"ind", "icc", "gen", "rel1", "rel2"),
include = TRUE,
chains = 3,
cores = 3,
iter = 4000,
seed = 1234)
#### print results
print(c4model_stan, digits = 3)
#### traceplots
traceplot(c4model_stan, pars=c("sd_T", "sd_R1", "sd_R2", "sd_I",
"intercepts_1", "intercepts_2", "loading_T",
"var_residual_strdiff", "var_residual_interch"))
traceplot(c4model_stan, pars=c("con1", "con2", "rsp1", "rsp2",
"ind", "icc", "gen", "rel1", "rel2"))MPlus Code