HMDCM
Load the spatial rotation data
(1) Simulate responses based on the HMDCM model
class_0 <- sample(1:2^K, N, replace = L)
Alphas_0 <- matrix(0,N,K)
for(i in 1:N){
Alphas_0[i,] <- inv_bijectionvector(K,(class_0[i]-1))
}
thetas_true = rnorm(N)
lambdas_true = c(-1, 1.8, .277, .055)
Alphas <- sim_alphas(model="HO_sep",
lambdas=lambdas_true,
thetas=thetas_true,
Q_matrix=Q_matrix,
Design_array=Design_array)
table(rowSums(Alphas[,,5]) - rowSums(Alphas[,,1])) # used to see how much transition has taken place
#>
#> 0 1 2 3 4
#> 33 53 83 134 47
itempars_true <- matrix(runif(J*2,.1,.2), ncol=2)
Y_sim <- sim_hmcdm(model="DINA",Alphas,Q_matrix,Design_array,
itempars=itempars_true)(2) Run the MCMC to sample parameters from the posterior distribution
output_HMDCM = hmcdm(Y_sim,Q_matrix,"DINA_HO",Test_order = Test_order, Test_versions = Test_versions,
chain_length=100,burn_in=30,
theta_propose = 2,deltas_propose = c(.45,.35,.25,.06))
#> 0
output_HMDCM = hmcdm(Y_sim,Q_matrix,"DINA_HO",Design_array,
chain_length=100,burn_in=30,
theta_propose = 2,deltas_propose = c(.45,.35,.25,.06))
#> 0
output_HMDCM
#>
#> Model: DINA_HO
#>
#> Sample Size: 350
#> Number of Items:
#> Number of Time Points:
#>
#> Chain Length: 100, burn-in: 30
summary(output_HMDCM)
#>
#> Model: DINA_HO
#>
#> Item Parameters:
#> ss_EAP gs_EAP
#> 0.1777 0.15443
#> 0.1381 0.30236
#> 0.2758 0.08192
#> 0.1625 0.22738
#> 0.1877 0.19436
#> ... 45 more items
#>
#> Transition Parameters:
#> lambdas_EAP
#> λ0 -1.5656
#> λ1 2.6412
#> λ2 0.1199
#> λ3 0.1878
#>
#> Class Probabilities:
#> pis_EAP
#> 0000 0.1360
#> 0001 0.2129
#> 0010 0.1943
#> 0011 0.2141
#> 0100 0.1387
#> ... 11 more classes
#>
#> Deviance Information Criterion (DIC): 19480.49
#>
#> Posterior Predictive P-value (PPP):
#> M1: 0.5117
#> M2: 0.49
#> total scores: 0.6223
a <- summary(output_HMDCM)
a$ss_EAP
#> [,1]
#> [1,] 0.17773985
#> [2,] 0.13805726
#> [3,] 0.27578995
#> [4,] 0.16252579
#> [5,] 0.18767900
#> [6,] 0.17444797
#> [7,] 0.28794946
#> [8,] 0.13560591
#> [9,] 0.16250225
#> [10,] 0.15875369
#> [11,] 0.13766880
#> [12,] 0.15650249
#> [13,] 0.17087214
#> [14,] 0.11467018
#> [15,] 0.23863191
#> [16,] 0.19535241
#> [17,] 0.22454393
#> [18,] 0.12890438
#> [19,] 0.22218281
#> [20,] 0.18620662
#> [21,] 0.22448626
#> [22,] 0.24910515
#> [23,] 0.13555271
#> [24,] 0.19767606
#> [25,] 0.15442061
#> [26,] 0.24951609
#> [27,] 0.22025459
#> [28,] 0.12448479
#> [29,] 0.21913241
#> [30,] 0.17833893
#> [31,] 0.15369285
#> [32,] 0.11406414
#> [33,] 0.25181884
#> [34,] 0.15563458
#> [35,] 0.18349669
#> [36,] 0.16940022
#> [37,] 0.12277529
#> [38,] 0.15864644
#> [39,] 0.16923324
#> [40,] 0.20026663
#> [41,] 0.16664965
#> [42,] 0.15281167
#> [43,] 0.09040818
#> [44,] 0.22583350
#> [45,] 0.15386080
#> [46,] 0.14948430
#> [47,] 0.19994224
#> [48,] 0.14762701
#> [49,] 0.16946876
#> [50,] 0.15530585
a$lambdas_EAP
#> [,1]
#> λ0 -1.5655547
#> λ1 2.6412092
#> λ2 0.1198932
#> λ3 0.1877814
mean(a$PPP_total_scores)
#> [1] 0.6226286
mean(upper.tri(a$PPP_item_ORs))
#> [1] 0.49
mean(a$PPP_item_means)
#> [1] 0.5157143(3) Evaluate the accuracy of estimated parameters
Attribute-wise agreement rate between true and estimated alphas
(4) Evaluate the fit of the model to the observed response
a$DIC
#> Transition Response_Time Response Joint Total
#> D_bar 1896.874 NA 15452.01 1303.196 18652.07
#> D(theta_bar) 1587.723 NA 14975.35 1260.581 17823.66
#> DIC 2206.024 NA 15928.66 1345.810 19480.49
head(a$PPP_total_scores)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.5285714 0.8142857 0.1714286 0.90000000 0.6142857
#> [2,] 0.5428571 0.3000000 0.7428571 0.48571429 0.3000000
#> [3,] 0.1714286 0.6142857 0.8857143 0.31428571 0.9000000
#> [4,] 0.4571429 0.9428571 0.6428571 0.02857143 1.0000000
#> [5,] 0.3428571 0.3000000 0.3285714 0.85714286 0.8000000
#> [6,] 0.9428571 0.6285714 0.4571429 0.95714286 0.5285714
head(a$PPP_item_means)
#> [1] 0.5142857 0.5142857 0.4714286 0.6000000 0.5285714 0.5142857
head(a$PPP_item_ORs)
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] NA 0.9285714 0.5714286 0.7000000 0.6714286 0.6571429 0.8714286 0.4857143
#> [2,] NA NA 0.4571429 0.6428571 0.8000000 0.5142857 0.1571429 0.2285714
#> [3,] NA NA NA 0.8571429 0.6714286 0.6142857 0.9142857 0.8714286
#> [4,] NA NA NA NA 0.9857143 0.8142857 0.6857143 0.5142857
#> [5,] NA NA NA NA NA 0.7714286 0.1142857 0.6000000
#> [6,] NA NA NA NA NA NA 0.2571429 0.3857143
#> [,9] [,10] [,11] [,12] [,13] [,14] [,15]
#> [1,] 0.3142857 0.5428571 0.4000000 0.90000000 0.9428571 0.9000000 0.94285714
#> [2,] 0.8428571 0.7428571 0.6857143 0.47142857 0.9285714 0.4571429 0.55714286
#> [3,] 0.4142857 0.6571429 0.4000000 0.48571429 0.4428571 0.4428571 0.08571429
#> [4,] 0.6142857 0.5000000 0.6285714 0.01428571 0.9571429 0.3142857 0.57142857
#> [5,] 0.6000000 0.4714286 0.6285714 0.52857143 0.9428571 0.5571429 0.28571429
#> [6,] 0.7000000 0.3428571 0.2000000 0.84285714 0.4428571 0.1000000 0.65714286
#> [,16] [,17] [,18] [,19] [,20] [,21] [,22]
#> [1,] 0.80000000 0.90000000 0.6857143 0.8142857 0.9428571 0.9428571 0.9714286
#> [2,] 0.51428571 0.40000000 0.4428571 0.8428571 0.1571429 0.7285714 0.7714286
#> [3,] 0.30000000 0.01428571 0.3142857 0.4142857 0.3857143 0.1428571 0.7285714
#> [4,] 0.87142857 0.37142857 0.7714286 0.3285714 0.9142857 0.7714286 0.8285714
#> [5,] 0.15714286 0.28571429 0.3428571 0.3571429 0.5857143 0.7714286 0.8714286
#> [6,] 0.07142857 0.14285714 0.4857143 0.9857143 0.2714286 0.3285714 0.4142857
#> [,23] [,24] [,25] [,26] [,27] [,28] [,29]
#> [1,] 0.8142857 0.9571429 0.22857143 0.3857143 0.6000000 0.3285714 0.8142857
#> [2,] 0.9714286 0.9857143 0.08571429 0.8142857 0.7142857 0.8000000 0.5000000
#> [3,] 0.6285714 0.2857143 0.51428571 0.1714286 0.8142857 0.9571429 0.2571429
#> [4,] 0.8142857 0.7142857 0.40000000 0.2571429 0.8285714 0.4142857 0.5285714
#> [5,] 0.8428571 0.6142857 0.25714286 0.3428571 0.7000000 0.3714286 0.1571429
#> [6,] 1.0000000 0.5714286 0.17142857 0.9142857 0.4285714 0.2285714 0.1857143
#> [,30] [,31] [,32] [,33] [,34] [,35] [,36]
#> [1,] 0.5142857 0.75714286 0.6428571 0.3571429 0.2142857 0.24285714 0.00000000
#> [2,] 0.1000000 0.27142857 0.5285714 0.8142857 0.5428571 0.04285714 0.08571429
#> [3,] 0.5714286 0.77142857 0.8428571 0.6857143 0.9428571 0.90000000 0.44285714
#> [4,] 0.1857143 0.01428571 0.7285714 0.7000000 0.9285714 0.94285714 0.25714286
#> [5,] 0.5571429 0.14285714 0.6714286 0.9571429 0.2285714 0.62857143 0.18571429
#> [6,] 0.4142857 0.57142857 0.5000000 0.2857143 0.8142857 0.08571429 0.07142857
#> [,37] [,38] [,39] [,40] [,41] [,42] [,43]
#> [1,] 0.3428571 0.6571429 0.4714286 0.35714286 0.2285714 0.1142857 0.2571429
#> [2,] 0.6000000 0.5714286 0.8857143 0.24285714 0.4142857 0.5428571 0.4285714
#> [3,] 0.6142857 0.6857143 0.7857143 0.78571429 0.4857143 0.5285714 0.9571429
#> [4,] 0.9285714 0.5285714 0.9142857 0.67142857 0.5000000 0.8714286 0.8428571
#> [5,] 0.6857143 0.9000000 0.7857143 0.30000000 0.8428571 0.3714286 0.4000000
#> [6,] 0.7857143 0.4571429 0.7857143 0.02857143 0.1857143 0.3571429 0.6000000
#> [,44] [,45] [,46] [,47] [,48] [,49] [,50]
#> [1,] 0.5000000 0.5000000 0.25714286 0.2571429 0.1714286 0.71428571 0.12857143
#> [2,] 0.4857143 0.4285714 0.08571429 0.2142857 0.3000000 0.07142857 0.08571429
#> [3,] 0.9000000 0.2285714 0.41428571 0.2714286 0.2000000 0.21428571 0.75714286
#> [4,] 0.1714286 0.9285714 0.68571429 0.8857143 0.5000000 0.41428571 0.78571429
#> [5,] 0.7142857 0.7714286 0.17142857 0.7285714 0.1000000 0.07142857 0.07142857
#> [6,] 0.4571429 0.9428571 0.57142857 0.2000000 0.4857143 0.47142857 0.41428571
library(bayesplot)
pp_check(output_HMDCM)
pp_check(output_HMDCM, plotfun="hist", type="item_OR")
#> Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
pp_check(output_HMDCM, plotfun="stat_2d", type="item_mean")
#> Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.
Convergence checking
Checking convergence of the two independent MCMC chains with
different initial values using coda package.
# output_HMDCM1 = hmcdm(Y_sim, Q_matrix, "DINA_HO", Design_array,
# chain_length=100, burn_in=30,
# theta_propose = 2, deltas_propose = c(.45,.35,.25,.06))
# output_HMDCM2 = hmcdm(Y_sim, Q_matrix, "DINA_HO", Design_array,
# chain_length=100, burn_in=30,
# theta_propose = 2, deltas_propose = c(.45,.35,.25,.06))
#
# library(coda)
#
# x <- mcmc.list(mcmc(t(rbind(output_HMDCM1$ss, output_HMDCM1$gs, output_HMDCM1$lambdas))),
# mcmc(t(rbind(output_HMDCM2$ss, output_HMDCM2$gs, output_HMDCM2$lambdas))))
#
# gelman.diag(x, autoburnin=F)