The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

HMDCM

library(hmcdm)

Load the spatial rotation data

N = dim(Design_array)[1]
J = nrow(Q_matrix)
K = ncol(Q_matrix)
L = dim(Design_array)[3]

(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 
#>  29  29 103 143  46
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.07685 0.17243
#>  0.18488 0.07224
#>  0.11893 0.09148
#>  0.11457 0.25188
#>  0.23254 0.14521
#>    ... 45 more items
#> 
#> Transition Parameters:
#>    lambdas_EAP
#> λ0     -1.3192
#> λ1      1.8989
#> λ2      0.2050
#> λ3      0.1283
#> 
#> Class Probabilities:
#>      pis_EAP
#> 0000  0.1459
#> 0001  0.1697
#> 0010  0.1465
#> 0011  0.2256
#> 0100  0.2020
#>    ... 11 more classes
#> 
#> Deviance Information Criterion (DIC): 19338.99 
#> 
#> Posterior Predictive P-value (PPP):
#> M1: 0.4997
#> M2:  0.49
#> total scores:  0.6243
a <- summary(output_HMDCM)
a$ss_EAP
#>             [,1]
#>  [1,] 0.07685255
#>  [2,] 0.18488425
#>  [3,] 0.11892786
#>  [4,] 0.11456749
#>  [5,] 0.23254269
#>  [6,] 0.18006215
#>  [7,] 0.12160810
#>  [8,] 0.21652203
#>  [9,] 0.19235236
#> [10,] 0.14879741
#> [11,] 0.20244524
#> [12,] 0.15301554
#> [13,] 0.14864318
#> [14,] 0.16782918
#> [15,] 0.08299005
#> [16,] 0.12658188
#> [17,] 0.17201443
#> [18,] 0.11441438
#> [19,] 0.21313259
#> [20,] 0.15480563
#> [21,] 0.23201792
#> [22,] 0.22040852
#> [23,] 0.18641166
#> [24,] 0.15644341
#> [25,] 0.15068218
#> [26,] 0.14069133
#> [27,] 0.15763156
#> [28,] 0.26286838
#> [29,] 0.20866321
#> [30,] 0.16014609
#> [31,] 0.20695261
#> [32,] 0.18302578
#> [33,] 0.16749500
#> [34,] 0.19857339
#> [35,] 0.17062041
#> [36,] 0.11672392
#> [37,] 0.19205306
#> [38,] 0.11277445
#> [39,] 0.16949512
#> [40,] 0.15137924
#> [41,] 0.22840721
#> [42,] 0.17458218
#> [43,] 0.11333775
#> [44,] 0.16533007
#> [45,] 0.15727484
#> [46,] 0.22936398
#> [47,] 0.13081081
#> [48,] 0.18149035
#> [49,] 0.22000283
#> [50,] 0.19266879
a$lambdas_EAP
#>          [,1]
#> λ0 -1.3191652
#> λ1  1.8988731
#> λ2  0.2049675
#> λ3  0.1282824
mean(a$PPP_total_scores)
#> [1] 0.6232245
mean(upper.tri(a$PPP_item_ORs))
#> [1] 0.49
mean(a$PPP_item_means)
#> [1] 0.4968571

(3) Evaluate the accuracy of estimated parameters

Attribute-wise agreement rate between true and estimated alphas

AAR_vec <- numeric(L)
for(t in 1:L){
  AAR_vec[t] <- mean(Alphas[,,t]==a$Alphas_est[,,t])
}
AAR_vec
#> [1] 0.9092857 0.9200000 0.9521429 0.9685714 0.9735714

Pattern-wise agreement rate between true and estimated alphas

PAR_vec <- numeric(L)
for(t in 1:L){
  PAR_vec[t] <- mean(rowSums((Alphas[,,t]-a$Alphas_est[,,t])^2)==0)
}
PAR_vec
#> [1] 0.7028571 0.7314286 0.8371429 0.8942857 0.9200000

(4) Evaluate the fit of the model to the observed response

a$DIC
#>              Transition Response_Time Response    Joint    Total
#> D_bar          2133.712            NA 15036.72 1323.997 18494.42
#> D(theta_bar)   1891.812            NA 14468.75 1289.299 17649.86
#> DIC            2375.612            NA 15604.68 1358.694 19338.99

head(a$PPP_total_scores)
#>           [,1]      [,2]       [,3]      [,4]      [,5]
#> [1,] 0.3857143 1.0000000 0.67142857 0.4714286 0.8000000
#> [2,] 0.5000000 0.8714286 0.02857143 1.0000000 0.2285714
#> [3,] 0.5571429 0.5000000 1.00000000 1.0000000 0.7714286
#> [4,] 0.3857143 0.8000000 0.81428571 0.9714286 0.5142857
#> [5,] 0.8857143 0.5428571 0.88571429 0.5428571 0.7714286
#> [6,] 0.9285714 0.6142857 0.92857143 0.1142857 1.0000000
head(a$PPP_item_means)
#> [1] 0.3714286 0.5142857 0.4714286 0.5285714 0.3571429 0.4000000
head(a$PPP_item_ORs)
#>      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
#> [1,]   NA 0.6857143 0.9857143 0.5428571 0.3714286 0.6857143 0.3857143 0.5285714
#> [2,]   NA        NA 0.3285714 0.8000000 0.2571429 0.7571429 0.6428571 0.6285714
#> [3,]   NA        NA        NA 0.4000000 0.7285714 0.1285714 0.8142857 0.9142857
#> [4,]   NA        NA        NA        NA 0.9142857 0.6571429 0.9285714 0.4142857
#> [5,]   NA        NA        NA        NA        NA 0.5714286 0.5285714 0.9142857
#> [6,]   NA        NA        NA        NA        NA        NA 0.8571429 0.5285714
#>           [,9]     [,10]     [,11]      [,12]      [,13]     [,14]     [,15]
#> [1,] 0.4285714 0.1285714 0.6285714 0.74285714 0.57142857 0.4857143 0.6428571
#> [2,] 0.5571429 0.5428571 0.8000000 0.80000000 0.17142857 0.3857143 0.4428571
#> [3,] 0.8285714 0.2857143 0.8285714 0.18571429 0.37142857 0.7000000 0.9142857
#> [4,] 0.6714286 0.2428571 0.4000000 0.80000000 0.37142857 0.5142857 0.4571429
#> [5,] 0.5857143 0.5571429 0.5428571 0.01428571 0.21428571 0.6142857 0.2000000
#> [6,] 0.8000000 0.6857143 0.8571429 0.60000000 0.01428571 0.3428571 0.2428571
#>          [,16]      [,17]     [,18]     [,19]     [,20]     [,21]      [,22]
#> [1,] 0.7571429 0.75714286 0.9285714 0.6000000 0.5571429 0.2285714 0.25714286
#> [2,] 0.4571429 0.24285714 0.8857143 0.1857143 0.7857143 0.5142857 0.34285714
#> [3,] 0.6714286 0.25714286 0.8571429 0.7571429 0.7571429 0.6571429 0.42857143
#> [4,] 0.6142857 0.62857143 0.7571429 0.6857143 0.7285714 0.6571429 0.50000000
#> [5,] 0.5285714 0.02857143 0.9571429 0.2285714 0.7428571 0.6428571 0.42857143
#> [6,] 0.4428571 0.62857143 0.8000000 0.5285714 0.6000000 0.2142857 0.08571429
#>          [,23]      [,24]      [,25]      [,26]     [,27]      [,28]     [,29]
#> [1,] 0.5000000 0.01428571 0.08571429 0.04285714 0.8714286 0.08571429 0.8000000
#> [2,] 0.2857143 0.31428571 0.52857143 0.18571429 0.1285714 0.42857143 0.1857143
#> [3,] 0.8142857 0.82857143 0.32857143 0.41428571 0.9285714 0.20000000 0.8571429
#> [4,] 0.2000000 0.38571429 0.15714286 0.00000000 1.0000000 0.88571429 0.8857143
#> [5,] 0.5142857 0.20000000 0.54285714 0.17142857 0.9857143 0.54285714 0.6857143
#> [6,] 0.2000000 0.28571429 0.01428571 0.10000000 0.5142857 0.07142857 0.1857143
#>          [,30]      [,31]     [,32]     [,33]     [,34]      [,35]      [,36]
#> [1,] 0.1142857 0.18571429 0.4285714 0.8714286 0.2714286 0.45714286 0.47142857
#> [2,] 0.2000000 0.45714286 0.3714286 0.7714286 0.2571429 0.52857143 0.47142857
#> [3,] 0.4428571 0.04285714 0.5285714 0.6000000 0.7000000 0.08571429 0.07142857
#> [4,] 0.7428571 0.64285714 0.2142857 0.1571429 0.1285714 0.44285714 0.27142857
#> [5,] 0.8714286 0.67142857 0.7142857 0.7571429 0.8285714 0.60000000 0.92857143
#> [6,] 0.4714286 0.18571429 0.3142857 0.5857143 0.1857143 0.54285714 0.47142857
#>          [,37]      [,38]     [,39]     [,40]      [,41]     [,42]     [,43]
#> [1,] 0.3142857 0.15714286 0.7142857 0.3285714 0.45714286 0.9857143 0.1000000
#> [2,] 0.3142857 0.64285714 0.4714286 0.6571429 0.07142857 0.5285714 0.5857143
#> [3,] 0.3285714 0.04285714 0.2142857 0.6285714 0.51428571 0.7285714 0.5000000
#> [4,] 0.2285714 0.72857143 0.1142857 0.8000000 0.78571429 0.5857143 0.2000000
#> [5,] 0.9857143 0.45714286 0.4571429 0.4571429 0.62857143 0.2428571 0.2857143
#> [6,] 0.1571429 0.88571429 0.2000000 0.6285714 0.14285714 0.7000000 0.8857143
#>          [,44]     [,45]     [,46]      [,47]     [,48]      [,49]     [,50]
#> [1,] 0.4142857 0.9428571 0.5428571 0.57142857 0.2571429 0.17142857 0.6571429
#> [2,] 0.7000000 0.3571429 0.5857143 0.28571429 0.2571429 0.40000000 0.3428571
#> [3,] 0.4857143 0.8142857 0.7142857 0.91428571 0.4428571 0.17142857 0.8571429
#> [4,] 0.7000000 0.5571429 0.3142857 0.04285714 0.1285714 0.78571429 0.6714286
#> [5,] 0.4000000 0.8714286 0.7142857 0.85714286 0.5428571 0.54285714 0.4428571
#> [6,] 0.8428571 0.5714286 0.4571429 0.04285714 0.2000000 0.08571429 0.4285714
library(bayesplot)
pp_check(output_HMDCM)

pp_check(output_HMDCM, plotfun="dens_overlay", type="item_mean")

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 with `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.

pp_check(output_HMDCM, plotfun="scatter_avg", type="total_score")

pp_check(output_HMDCM, plotfun="error_scatter_avg", type="total_score")

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)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.