5.11 Quanitities in GDINA R package

Code

library(GDINA)
Y <- sim10GDINA$simdat
head(Y)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    0    1    1    0    1    1    1    0     1
## [2,]    1    1    1    1    1    1    1    1    0     1
## [3,]    1    1    1    0    1    1    1    0    1     0
## [4,]    1    0    1    1    1    1    1    1    1     0
## [5,]    0    0    1    1    0    0    1    0    0     0
## [6,]    1    0    0    0    0    0    0    1    0     1

Code

Q <- sim10GDINA$simQ
Q

##       [,1] [,2] [,3]
##  [1,]    1    0    0
##  [2,]    0    1    0
##  [3,]    0    0    1
##  [4,]    1    0    1
##  [5,]    0    1    1
##  [6,]    1    1    0
##  [7,]    1    0    1
##  [8,]    1    1    0
##  [9,]    0    1    1
## [10,]    1    1    1

Let us fit DINA model to the data

Code

est <- GDINA(Y, Q, model = "DINA", verbose = 0)

The logarithm of \(L( {\mathbf{Y}_i}|{\mathbf{\alpha} _c})\) is a matrix of \(N\times 2^K\):

Code

head(indlogLik(est))

##        000   100   010 001   110  101   011   111
## [1,] -13.0  -9.9 -14.2 -10  -9.0 -5.1 -14.1  -5.8
## [2,] -17.7 -14.6 -15.6 -15 -10.5 -9.8 -12.2  -3.9
## [3,] -15.0 -11.8 -12.9 -12 -10.3 -9.7  -7.1  -6.4
## [4,] -15.7 -12.6 -16.9 -13 -11.8 -7.8 -11.1  -5.2
## [5,]  -7.5  -8.8  -8.7  -5 -12.8 -4.0  -8.6 -11.9
## [6,]  -8.4  -5.3  -9.6 -11  -6.7 -9.6 -14.3 -12.6

The logarithm of \(P({\mathbf{\alpha} _c}| {\mathbf{Y}_i})\) is also a matrix of \(N\times 2^K\):

Code

head(indlogPost(est))

##        000    100   010   001  110   101   011     111
## [1,]  -9.3  -5.73 -10.6  -7.1 -4.3 -0.50 -10.4 -0.9696
## [2,] -14.9 -11.38 -13.0 -12.7 -6.7 -6.15  -9.4 -0.0035
## [3,]  -9.9  -6.36  -8.0  -7.7 -4.3 -3.76  -2.0 -0.1953
## [4,] -11.7  -8.15 -13.0  -9.5 -6.8 -2.92  -7.0 -0.0580
## [5,]  -4.5  -5.37  -5.8  -2.3 -8.7 -0.14  -5.5 -7.7685
## [6,]  -4.0  -0.41  -5.3  -6.4 -1.2 -4.31  -9.7 -6.9765

The \(N_j\) can be obtained using the code below

Code

GDINA::extract(est, what = "expectedTotal")

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##  [1,]  274  726    0    0    0    0    0    0
##  [2,]  419  581    0    0    0    0    0    0
##  [3,]  458  542    0    0    0    0    0    0
##  [4,]  139  318  135  408    0    0    0    0
##  [5,]  188  270  232  311    0    0    0    0
##  [6,]  128  291  146  435    0    0    0    0
##  [7,]  139  318  135  408    0    0    0    0
##  [8,]  128  291  146  435    0    0    0    0
##  [9,]  188  270  232  311    0    0    0    0
## [10,]   74  114   65   54  205  177   81  230

The \(R_j\) can be obtained by

Code

GDINA::extract(est, what = "expectedCorrect")

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##  [1,]  8.9  541    0    0    0    0    0    0
##  [2,] 38.5  422    0    0    0    0    0    0
##  [3,] 32.8  484    0    0    0    0    0    0
##  [4,] 22.9   70   38  325    0    0    0    0
##  [5,] 18.3   22   16  222    0    0    0    0
##  [6,] 43.7  185  102  403    0    0    0    0
##  [7,] 15.7  100   36  276    0    0    0    0
##  [8,] 14.6   49   20  301    0    0    0    0
##  [9,] 34.8   78   72  248    0    0    0    0
## [10,] 15.6   26   16   19   74   69   38  195

The score function is the key component for estimating standard errors.

The code below gives the score function for response vector of each individual

\[\mathbf{u}(\tau)=\frac{\partial \ell(\mathbf{y}_i)}{\partial \tau}\]

Code

u <- score(est)
u <- do.call(cbind, u)

The information can be estimated using observed information: \[ FI(\hat{\tau})=N\times J_{\mathbf{y}}(\hat{\tau}) \]

Code

FI <- solve(crossprod(u))
SE <- sqrt(diag(FI))
SE

##                                                                   
## 0.067 0.074 0.031 0.042 0.027 0.035 0.022 0.038 0.013 0.038 0.027 
##                                                         100   010 
## 0.034 0.022 0.038 0.021 0.041 0.020 0.035 0.019 0.038 0.023 0.018 
##   001   110   101   011   111 
## 0.015 0.022 0.020 0.015 0.019