rRUM_indept

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 and response times based on the rRUM model

tau <- numeric(K)
for(k in 1:K){
  tau[k] <- runif(1,.2,.6)
}
R = matrix(0,K,K)
# Initial alphas
p_mastery <- c(.5,.5,.4,.4)
Alphas_0 <- matrix(0,N,K)
for(i in 1:N){
  for(k in 1:K){
    prereqs <- which(R[k,]==1)
    if(length(prereqs)==0){
      Alphas_0[i,k] <- rbinom(1,1,p_mastery[k])
    }
    if(length(prereqs)>0){
      Alphas_0[i,k] <- prod(Alphas_0[i,prereqs])*rbinom(1,1,p_mastery)
    }
  }
}
Alphas <- sim_alphas(model="indept",taus=tau,N=N,L=L,R=R,alpha0=Alphas_0)
table(rowSums(Alphas[,,5]) - rowSums(Alphas[,,1])) # used to see how much transition has taken place
#> 
#>   0   1   2   3   4 
#>  35  89 123  90  13
Smats <- matrix(runif(J*K,.1,.3),c(J,K))
Gmats <- matrix(runif(J*K,.1,.3),c(J,K))
# Simulate rRUM parameters
r_stars <- Gmats / (1-Smats)
pi_stars <- apply((1-Smats)^Q_matrix, 1, prod)

Y_sim <- sim_hmcdm(model="rRUM",Alphas,Q_matrix,Design_array,
                   r_stars=r_stars,pi_stars=pi_stars)

(2) Run the MCMC to sample parameters from the posterior distribution

output_rRUM_indept = hmcdm(Y_sim,Q_matrix,"rRUM_indept",Design_array,
                           100,30,R = R)
#> 0
output_rRUM_indept
#> 
#> Model: rRUM_indept 
#> 
#> Sample Size: 350
#> Number of Items: 
#> Number of Time Points: 
#> 
#> Chain Length: 100, burn-in: 50
summary(output_rRUM_indept)
#> 
#> Model: rRUM_indept 
#> 
#> Item Parameters:
#>  r_stars1_EAP r_stars2_EAP r_stars3_EAP r_stars4_EAP pi_stars_EAP
#>        0.3362       0.5892       0.6147       0.5751       0.7412
#>        0.6162       0.4175       0.5610       0.5986       0.7796
#>        0.6981       0.6128       0.6977       0.1720       0.8454
#>        0.5485       0.6919       0.1878       0.5383       0.8108
#>        0.5194       0.1871       0.6001       0.6343       0.6670
#>    ... 45 more items
#> 
#> Transition Parameters:
#>    taus_EAP
#> τ1   0.5209
#> τ2   0.5267
#> τ3   0.3309
#> τ4   0.3452
#> 
#> Class Probabilities:
#>      pis_EAP
#> 0000 0.06715
#> 0001 0.13740
#> 0010 0.08418
#> 0011 0.02898
#> 0100 0.04788
#>    ... 11 more classes
#> 
#> Deviance Information Criterion (DIC): 22528.85 
#> 
#> Posterior Predictive P-value (PPP):
#> M1: 0.5152
#> M2:  0.49
#> total scores:  0.6134
a <- summary(output_rRUM_indept)
head(a$r_stars_EAP)
#>           [,1]      [,2]      [,3]      [,4]
#> [1,] 0.3361732 0.5891916 0.6146525 0.5750689
#> [2,] 0.6161881 0.4175006 0.5609841 0.5985672
#> [3,] 0.6980576 0.6128256 0.6977122 0.1719537
#> [4,] 0.5484959 0.6918715 0.1877874 0.5382910
#> [5,] 0.5194453 0.1871043 0.6001193 0.6342858
#> [6,] 0.6541725 0.1623346 0.1288330 0.5433937

(3) Check for parameter estimation accuracy

(cor_pistars <- cor(as.vector(pi_stars),as.vector(a$pi_stars_EAP)))
#> [1] 0.9418397
(cor_rstars <- cor(as.vector(r_stars*Q_matrix),as.vector(a$r_stars_EAP*Q_matrix)))
#> [1] 0.9243346

AAR_vec <- numeric(L)
for(t in 1:L){
  AAR_vec[t] <- mean(Alphas[,,t]==a$Alphas_est[,,t])
}
AAR_vec
#> [1] 0.8471429 0.8950000 0.9328571 0.9628571 0.9657143

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.5028571 0.6428571 0.7714286 0.8628571 0.8685714

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

a$DIC
#>              Transition Response_Time Response    Joint    Total
#> D_bar          2098.313            NA 17888.47 1858.738 21845.52
#> D(theta_bar)   2055.570            NA 17249.38 1857.245 21162.19
#> DIC            2141.056            NA 18527.56 1860.232 22528.85
head(a$PPP_total_scores)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.30 0.52 0.98 0.52 0.82
#> [2,] 0.74 0.60 0.76 0.20 0.20
#> [3,] 0.74 0.22 0.54 0.50 0.74
#> [4,] 0.22 0.86 0.76 0.78 0.76
#> [5,] 0.68 0.92 0.94 0.16 0.36
#> [6,] 0.94 0.72 0.16 0.86 0.26
head(a$PPP_item_means)
#> [1] 0.52 0.50 0.48 0.42 0.48 0.60
head(a$PPP_item_ORs)
#>      [,1] [,2] [,3] [,4] [,5]       [,6]     [,7]      [,8]      [,9]     [,10]
#> [1,]   NA 0.54 0.40 0.32 0.64 0.51020408 0.840000 0.1200000 0.5000000 0.3000000
#> [2,]   NA   NA 0.52 0.58 0.56 0.42857143 0.380000 0.4800000 0.4600000 0.6800000
#> [3,]   NA   NA   NA 0.98 0.96 0.06122449 0.220000 0.6400000 1.0000000 0.6200000
#> [4,]   NA   NA   NA   NA 0.10 0.93877551 0.020000 0.7200000 0.7600000 0.8200000
#> [5,]   NA   NA   NA   NA   NA 0.10204082 0.820000 0.6600000 0.2000000 0.1400000
#> [6,]   NA   NA   NA   NA   NA         NA 0.122449 0.3061224 0.8367347 0.9183673
#>          [,11]    [,12]     [,13]      [,14]     [,15]     [,16]     [,17]
#> [1,] 0.0600000 0.680000 0.2400000 0.86000000 0.1800000 0.6200000 0.2000000
#> [2,] 0.8600000 0.960000 0.9600000 0.42000000 0.3600000 0.6600000 0.2600000
#> [3,] 0.1200000 0.060000 0.1600000 0.60000000 0.7200000 0.3600000 0.6800000
#> [4,] 0.7400000 0.320000 0.6400000 0.54000000 0.6800000 0.3000000 0.9800000
#> [5,] 0.3600000 0.520000 0.7000000 0.18000000 0.0600000 0.5400000 0.4600000
#> [6,] 0.2244898 0.755102 0.1020408 0.08163265 0.3877551 0.5102041 0.6734694
#>          [,18]     [,19]     [,20]     [,21]      [,22]     [,23]     [,24]
#> [1,] 0.1000000 0.4800000 0.3000000 0.6800000 0.90000000 0.2200000 0.1200000
#> [2,] 0.3800000 0.9000000 0.5000000 0.6800000 0.46000000 0.1800000 0.3200000
#> [3,] 0.6600000 0.8200000 0.6000000 0.7800000 0.74000000 0.4400000 0.4200000
#> [4,] 0.4600000 0.2800000 0.6000000 1.0000000 0.36000000 0.4000000 0.5200000
#> [5,] 0.4800000 0.2400000 0.6800000 0.6400000 0.66000000 0.9000000 0.1200000
#> [6,] 0.8979592 0.7959184 0.2857143 0.8367347 0.08163265 0.3061224 0.6326531
#>          [,25]     [,26]     [,27]     [,28]     [,29]     [,30]     [,31]
#> [1,] 0.6800000 0.5000000 0.0800000 0.9000000 0.7200000 0.6000000 0.4200000
#> [2,] 0.2400000 0.4800000 0.1000000 0.1400000 0.4400000 0.3000000 0.0200000
#> [3,] 0.2400000 0.2200000 0.1000000 0.9800000 0.7400000 0.7200000 0.3400000
#> [4,] 0.6200000 0.5800000 0.5000000 0.4400000 0.6400000 0.0800000 0.1200000
#> [5,] 0.8000000 1.0000000 0.8400000 0.6200000 0.0600000 0.9400000 0.0000000
#> [6,] 0.5102041 0.9795918 0.7755102 0.5918367 0.7959184 0.1428571 0.1632653
#>           [,32]     [,33]    [,34]     [,35]     [,36]    [,37]     [,38]
#> [1,] 0.74000000 0.5600000 0.660000 0.8400000 0.7800000 0.860000 0.3200000
#> [2,] 0.28000000 0.4200000 0.200000 0.7400000 0.9400000 0.700000 0.9200000
#> [3,] 0.84000000 0.0400000 0.120000 0.9600000 0.4400000 0.960000 0.5800000
#> [4,] 0.80000000 0.4800000 0.440000 0.6000000 0.1200000 0.800000 0.7400000
#> [5,] 0.28000000 0.5800000 0.800000 0.6000000 0.3800000 0.800000 0.6200000
#> [6,] 0.06122449 0.3061224 0.122449 0.4081633 0.4489796 0.122449 0.6734694
#>          [,39]     [,40]     [,41]     [,42] [,43]      [,44]     [,45] [,46]
#> [1,] 0.9400000 0.9200000 0.2600000 0.1400000  0.74 0.78000000 0.0800000  0.88
#> [2,] 0.4400000 0.9400000 0.7600000 0.5400000  0.00 0.88000000 0.7000000  0.02
#> [3,] 0.0800000 0.5800000 0.8600000 0.9000000  0.42 0.34000000 0.8000000  0.54
#> [4,] 0.1600000 0.8800000 0.4200000 0.7600000  0.70 0.08000000 1.0000000  0.02
#> [5,] 0.9800000 0.2000000 0.9000000 0.7000000  0.12 0.36000000 0.6600000  0.80
#> [6,] 0.1632653 0.6530612 0.3469388 0.2857143  0.00 0.04081633 0.4897959  0.00
#>          [,47]     [,48]     [,49]     [,50]
#> [1,] 0.3200000 0.9000000 0.3200000 0.8200000
#> [2,] 0.7800000 0.1800000 0.6000000 0.9400000
#> [3,] 0.5600000 0.5800000 0.4200000 0.8400000
#> [4,] 0.8800000 0.0400000 0.1400000 0.6600000
#> [5,] 0.7000000 0.4400000 0.0400000 0.2000000
#> [6,] 0.9795918 0.2857143 0.1428571 0.1428571