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
#> 41 41 79 140 49
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.17388 0.09328
#> 0.13336 0.14467
#> 0.13961 0.24452
#> 0.16495 0.12158
#> 0.08505 0.11944
#> ... 45 more items
#>
#> Transition Parameters:
#> lambdas_EAP
#> λ0 -0.95282
#> λ1 1.72919
#> λ2 0.21742
#> λ3 0.05902
#>
#> Class Probabilities:
#> pis_EAP
#> 0000 0.1578
#> 0001 0.1951
#> 0010 0.1946
#> 0011 0.1861
#> 0100 0.1921
#> ... 11 more classes
#>
#> Deviance Information Criterion (DIC): 19337.81
#>
#> Posterior Predictive P-value (PPP):
#> M1: 0.5117
#> M2: 0.49
#> total scores: 0.6248
a <- summary(output_HMDCM)
a$ss_EAP
#> [,1]
#> [1,] 0.17388408
#> [2,] 0.13335914
#> [3,] 0.13960816
#> [4,] 0.16494693
#> [5,] 0.08504727
#> [6,] 0.15229472
#> [7,] 0.18367384
#> [8,] 0.15483915
#> [9,] 0.18908900
#> [10,] 0.18639015
#> [11,] 0.16323812
#> [12,] 0.17291495
#> [13,] 0.24159677
#> [14,] 0.26116288
#> [15,] 0.16743593
#> [16,] 0.21614725
#> [17,] 0.11215030
#> [18,] 0.18573990
#> [19,] 0.16437316
#> [20,] 0.23021428
#> [21,] 0.18735544
#> [22,] 0.14007320
#> [23,] 0.10047528
#> [24,] 0.16207539
#> [25,] 0.18020759
#> [26,] 0.21088515
#> [27,] 0.16874801
#> [28,] 0.14283116
#> [29,] 0.18138966
#> [30,] 0.19005069
#> [31,] 0.14604301
#> [32,] 0.11364015
#> [33,] 0.12522509
#> [34,] 0.22717094
#> [35,] 0.17852813
#> [36,] 0.19671082
#> [37,] 0.10432572
#> [38,] 0.18545782
#> [39,] 0.23402973
#> [40,] 0.16456017
#> [41,] 0.12403590
#> [42,] 0.24371591
#> [43,] 0.20245449
#> [44,] 0.16189683
#> [45,] 0.14463502
#> [46,] 0.09887668
#> [47,] 0.15577063
#> [48,] 0.15040923
#> [49,] 0.17701869
#> [50,] 0.18046105
a$lambdas_EAP
#> [,1]
#> λ0 -0.95281601
#> λ1 1.72918944
#> λ2 0.21741915
#> λ3 0.05901882
mean(a$PPP_total_scores)
#> [1] 0.6243102
mean(upper.tri(a$PPP_item_ORs))
#> [1] 0.49
mean(a$PPP_item_means)
#> [1] 0.5028571(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 2167.145 NA 15079.11 1268.443 18514.70
#> D(theta_bar) 1903.833 NA 14549.61 1238.150 17691.59
#> DIC 2430.458 NA 15608.62 1298.736 19337.81
head(a$PPP_total_scores)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.8285714 0.5285714 0.8285714 0.8428571 0.1571429
#> [2,] 0.6285714 0.9571429 0.8000000 0.5428571 0.2000000
#> [3,] 0.7714286 0.9714286 0.1714286 0.4571429 0.9142857
#> [4,] 0.5000000 0.2714286 0.6857143 0.8857143 0.4714286
#> [5,] 0.3714286 0.6428571 0.5571429 0.6714286 0.7428571
#> [6,] 0.4428571 0.7285714 0.2857143 1.0000000 0.8428571
head(a$PPP_item_means)
#> [1] 0.5000000 0.5571429 0.5857143 0.4571429 0.4857143 0.5142857
head(a$PPP_item_ORs)
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] NA 0.8857143 0.5428571 0.5285714 0.8714286 0.6428571 0.4571429 0.3857143
#> [2,] NA NA 0.7714286 0.6714286 0.6142857 0.7000000 0.5571429 0.8857143
#> [3,] NA NA NA 0.8714286 0.5285714 0.6714286 0.5428571 0.7285714
#> [4,] NA NA NA NA 0.5714286 0.9000000 0.4285714 0.3571429
#> [5,] NA NA NA NA NA 0.1285714 0.5000000 0.7428571
#> [6,] NA NA NA NA NA NA 0.4714286 0.8285714
#> [,9] [,10] [,11] [,12] [,13] [,14] [,15]
#> [1,] 0.3571429 0.6857143 0.3571429 0.8857143 0.6857143 1.00000000 0.7142857
#> [2,] 0.3000000 0.3000000 0.4142857 0.4000000 0.2000000 0.57142857 0.2285714
#> [3,] 0.4857143 0.7142857 0.3857143 0.7571429 0.7285714 0.04285714 0.6000000
#> [4,] 0.5285714 0.4571429 0.9857143 0.9571429 0.7000000 0.90000000 0.9857143
#> [5,] 0.3857143 0.5571429 0.3428571 0.5857143 0.8285714 0.82857143 0.2142857
#> [6,] 0.6000000 0.6428571 0.7000000 0.8714286 0.1571429 0.64285714 0.6142857
#> [,16] [,17] [,18] [,19] [,20] [,21] [,22]
#> [1,] 0.9857143 1.0000000 0.4428571 0.5000000 0.7857143 0.8142857 0.07142857
#> [2,] 0.8571429 0.8857143 0.2000000 0.2000000 0.3714286 0.4285714 0.20000000
#> [3,] 0.2571429 0.9285714 0.4428571 0.5285714 0.9285714 0.6857143 0.84285714
#> [4,] 0.9714286 0.9285714 0.7142857 0.5142857 0.1000000 0.2285714 0.35714286
#> [5,] 0.9428571 0.8857143 0.4428571 0.6714286 0.3142857 0.2857143 0.31428571
#> [6,] 0.9000000 0.7714286 0.7285714 0.7714286 0.2571429 0.3142857 0.64285714
#> [,23] [,24] [,25] [,26] [,27] [,28] [,29]
#> [1,] 0.04285714 0.1428571 0.1571429 0.01428571 0.08571429 0.27142857 0.7857143
#> [2,] 0.45714286 0.5000000 0.8571429 0.08571429 0.22857143 0.07142857 0.1000000
#> [3,] 0.88571429 1.0000000 0.9571429 0.87142857 0.82857143 0.61428571 0.7857143
#> [4,] 0.42857143 0.9000000 0.5714286 0.11428571 0.47142857 0.40000000 0.6571429
#> [5,] 0.24285714 0.2285714 0.7142857 0.02857143 0.14285714 0.20000000 0.6428571
#> [6,] 0.61428571 0.9000000 0.9571429 0.21428571 0.18571429 0.67142857 0.4857143
#> [,30] [,31] [,32] [,33] [,34] [,35] [,36]
#> [1,] 0.12857143 0.04285714 0.5142857 0.01428571 0.3571429 0.20000000 0.3000000
#> [2,] 0.00000000 0.57142857 0.9000000 0.27142857 0.3571429 0.32857143 0.6857143
#> [3,] 0.30000000 0.21428571 0.3285714 0.01428571 0.1142857 0.44285714 0.2857143
#> [4,] 0.98571429 0.07142857 0.2714286 0.30000000 0.6142857 0.05714286 0.4857143
#> [5,] 0.02857143 0.10000000 0.5428571 0.21428571 0.5428571 0.12857143 0.1000000
#> [6,] 0.85714286 0.50000000 0.9142857 0.42857143 0.7857143 0.31428571 1.0000000
#> [,37] [,38] [,39] [,40] [,41] [,42] [,43]
#> [1,] 0.04285714 0.27142857 0.5142857 0.2714286 0.1285714 0.18571429 0.2142857
#> [2,] 0.37142857 0.17142857 0.5714286 0.5857143 0.8428571 0.61428571 0.6857143
#> [3,] 0.85714286 0.18571429 0.6285714 0.2428571 0.3857143 0.34285714 0.3285714
#> [4,] 0.24285714 0.05714286 0.4285714 0.8571429 0.4285714 0.07142857 0.7428571
#> [5,] 0.07142857 0.00000000 0.7285714 0.4142857 0.8571429 0.41428571 0.5142857
#> [6,] 0.87142857 0.18571429 0.7857143 0.5285714 0.5714286 0.27142857 0.9285714
#> [,44] [,45] [,46] [,47] [,48] [,49] [,50]
#> [1,] 0.4857143 0.8571429 0.8571429 0.55714286 0.62857143 0.4571429 0.30000000
#> [2,] 0.6000000 0.2857143 0.6857143 0.12857143 0.85714286 0.2428571 0.24285714
#> [3,] 0.4142857 0.6285714 0.4142857 0.04285714 0.61428571 0.6571429 0.81428571
#> [4,] 0.7000000 0.5285714 0.3571429 0.12857143 0.02857143 0.2285714 0.15714286
#> [5,] 0.4000000 0.8428571 0.7714286 0.82857143 0.42857143 0.3714286 0.38571429
#> [6,] 0.6428571 0.5857143 0.6857143 0.10000000 0.22857143 0.1285714 0.08571429
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)