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
#> 31 45 91 149 34
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.1567 0.2107
#> 0.1502 0.1071
#> 0.1256 0.1590
#> 0.1568 0.1881
#> 0.1617 0.2115
#> ... 45 more items
#>
#> Transition Parameters:
#> lambdas_EAP
#> λ0 -1.05621
#> λ1 1.68180
#> λ2 0.26279
#> λ3 0.06381
#>
#> Class Probabilities:
#> pis_EAP
#> 0000 0.1837
#> 0001 0.2195
#> 0010 0.1551
#> 0011 0.2107
#> 0100 0.1724
#> ... 11 more classes
#>
#> Deviance Information Criterion (DIC): 19093.47
#>
#> Posterior Predictive P-value (PPP):
#> M1: 0.5094
#> M2: 0.49
#> total scores: 0.6252
a <- summary(output_HMDCM)
a$ss_EAP
#> [,1]
#> [1,] 0.15674162
#> [2,] 0.15016132
#> [3,] 0.12563010
#> [4,] 0.15682447
#> [5,] 0.16169633
#> [6,] 0.14186899
#> [7,] 0.17978033
#> [8,] 0.08995010
#> [9,] 0.15205233
#> [10,] 0.18347899
#> [11,] 0.10903263
#> [12,] 0.14199766
#> [13,] 0.16327322
#> [14,] 0.12594233
#> [15,] 0.15531927
#> [16,] 0.16357672
#> [17,] 0.16279474
#> [18,] 0.13806592
#> [19,] 0.17988085
#> [20,] 0.14550419
#> [21,] 0.18846648
#> [22,] 0.14130972
#> [23,] 0.16228990
#> [24,] 0.24586299
#> [25,] 0.15600760
#> [26,] 0.25088389
#> [27,] 0.13268170
#> [28,] 0.13904112
#> [29,] 0.22615164
#> [30,] 0.17379717
#> [31,] 0.10068254
#> [32,] 0.15539065
#> [33,] 0.15064743
#> [34,] 0.16170573
#> [35,] 0.14677115
#> [36,] 0.14972082
#> [37,] 0.11654722
#> [38,] 0.14006834
#> [39,] 0.12842057
#> [40,] 0.18173773
#> [41,] 0.13981362
#> [42,] 0.18991598
#> [43,] 0.16928146
#> [44,] 0.20796312
#> [45,] 0.15699462
#> [46,] 0.10547156
#> [47,] 0.17845330
#> [48,] 0.12230480
#> [49,] 0.09015153
#> [50,] 0.15243902
a$lambdas_EAP
#> [,1]
#> λ0 -1.05621488
#> λ1 1.68179534
#> λ2 0.26278613
#> λ3 0.06381302
mean(a$PPP_total_scores)
#> [1] 0.6265551
mean(upper.tri(a$PPP_item_ORs))
#> [1] 0.49
mean(a$PPP_item_means)
#> [1] 0.5125714(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 2320.978 NA 14716.06 1226.063 18263.10
#> D(theta_bar) 2043.147 NA 14204.08 1185.502 17432.73
#> DIC 2598.810 NA 15228.04 1266.623 19093.47
head(a$PPP_total_scores)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.6857143 0.37142857 0.4285714 0.8857143 0.7857143
#> [2,] 0.4428571 0.27142857 1.0000000 0.3714286 0.9000000
#> [3,] 0.5571429 0.05714286 0.7285714 0.7428571 0.5142857
#> [4,] 0.8142857 0.04285714 0.6142857 0.9571429 0.6142857
#> [5,] 0.6000000 0.72857143 0.1000000 0.8142857 0.5142857
#> [6,] 0.3714286 0.58571429 0.2000000 0.6571429 0.5428571
head(a$PPP_item_means)
#> [1] 0.5857143 0.4428571 0.5000000 0.5285714 0.5571429 0.5857143
head(a$PPP_item_ORs)
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] NA 0.6285714 0.4000000 0.7428571 0.4571429 0.9857143 0.5142857 0.2428571
#> [2,] NA NA 0.4285714 0.8000000 0.4857143 0.3857143 0.7857143 0.6714286
#> [3,] NA NA NA 0.5571429 0.9285714 0.7285714 0.6285714 0.7000000
#> [4,] NA NA NA NA 0.8285714 0.9142857 0.9428571 0.9000000
#> [5,] NA NA NA NA NA 0.6285714 0.5857143 0.4857143
#> [6,] NA NA NA NA NA NA 0.7428571 0.5000000
#> [,9] [,10] [,11] [,12] [,13] [,14] [,15]
#> [1,] 0.8285714 0.8571429 0.9571429 0.8857143 0.6285714 0.85714286 0.8714286
#> [2,] 0.5428571 0.5285714 0.9142857 0.6285714 0.8142857 0.05714286 0.1857143
#> [3,] 0.8428571 0.9571429 0.6428571 0.7000000 0.9571429 0.74285714 0.6857143
#> [4,] 0.7571429 0.5142857 0.5714286 0.2857143 0.2428571 0.11428571 0.4714286
#> [5,] 0.9000000 0.7000000 0.8857143 0.1428571 0.8857143 0.57142857 0.3571429
#> [6,] 0.8571429 0.5142857 0.9000000 0.4428571 0.4571429 0.21428571 0.4142857
#> [,16] [,17] [,18] [,19] [,20] [,21] [,22]
#> [1,] 1.0000000 0.80000000 0.8142857 0.9000000 0.81428571 0.58571429 0.5285714
#> [2,] 0.6857143 0.31428571 0.4285714 0.4571429 0.08571429 0.35714286 0.5571429
#> [3,] 0.9142857 0.74285714 0.7857143 0.7428571 0.98571429 0.28571429 0.8571429
#> [4,] 0.2714286 0.04285714 0.6571429 0.6857143 0.17142857 0.40000000 0.7142857
#> [5,] 0.7857143 0.51428571 0.8714286 0.7428571 0.64285714 0.08571429 0.8714286
#> [6,] 0.5714286 0.62857143 0.5714286 0.8857143 0.67142857 0.68571429 0.8714286
#> [,23] [,24] [,25] [,26] [,27] [,28] [,29]
#> [1,] 0.5714286 0.85714286 0.27142857 0.8857143 0.3142857 0.9285714 0.8142857
#> [2,] 0.7000000 0.07142857 0.04285714 0.6285714 0.2142857 0.4000000 0.1857143
#> [3,] 0.5428571 0.58571429 0.47142857 0.5285714 0.2285714 0.9142857 0.4285714
#> [4,] 0.6285714 0.21428571 0.00000000 0.5000000 0.9714286 0.1857143 0.1285714
#> [5,] 0.8285714 0.45714286 0.20000000 0.3428571 0.8285714 0.3857143 0.1428571
#> [6,] 0.8285714 0.88571429 0.75714286 0.5857143 0.8857143 0.6857143 0.4714286
#> [,30] [,31] [,32] [,33] [,34] [,35] [,36]
#> [1,] 0.90000000 0.1000000 1.0000000 0.2000000 0.5714286 0.3428571 0.45714286
#> [2,] 0.04285714 0.6142857 0.4714286 0.6571429 0.7857143 0.8142857 0.51428571
#> [3,] 0.48571429 0.7142857 0.2428571 0.2428571 0.7000000 0.5714286 0.82857143
#> [4,] 0.08571429 0.8000000 0.3000000 0.4857143 0.8000000 0.1142857 0.04285714
#> [5,] 0.38571429 0.4571429 0.0000000 0.1571429 0.9000000 0.5142857 0.44285714
#> [6,] 0.61428571 0.5714286 0.4428571 0.6571429 0.8857143 0.1571429 0.22857143
#> [,37] [,38] [,39] [,40] [,41] [,42] [,43]
#> [1,] 0.4142857 0.6000000 0.05714286 0.7000000 0.6285714 1.0000000 0.62857143
#> [2,] 0.6714286 0.9428571 0.78571429 0.5714286 0.9571429 0.9857143 0.58571429
#> [3,] 0.9428571 0.8142857 0.04285714 0.7428571 0.9285714 0.8142857 0.60000000
#> [4,] 0.9285714 0.4714286 0.80000000 0.6714286 0.4000000 0.9857143 0.78571429
#> [5,] 0.8000000 0.4714286 0.75714286 0.9000000 0.9000000 0.8000000 0.05714286
#> [6,] 0.7000000 0.5428571 0.55714286 0.7857143 0.8000000 0.9428571 0.60000000
#> [,44] [,45] [,46] [,47] [,48] [,49] [,50]
#> [1,] 0.6857143 0.4857143 0.8571429 0.8571429 0.8714286 0.7285714 0.8142857
#> [2,] 0.8857143 0.9000000 0.9571429 0.8285714 0.2000000 0.8857143 0.4428571
#> [3,] 0.4857143 0.1571429 1.0000000 0.6714286 0.8000000 0.7142857 0.2857143
#> [4,] 0.6428571 0.7142857 0.9142857 0.8571429 0.4857143 0.9857143 0.8857143
#> [5,] 0.3428571 0.8428571 0.5428571 0.8571429 0.4142857 0.2000000 0.8142857
#> [6,] 0.5857143 0.9285714 0.2285714 0.8285714 0.1142857 0.3857143 0.4142857
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)