HW4

## Loading required package: CDM

## Loading required package: mvtnorm

## **********************************
## ** CDM 8.2-6 (2022-08-25 15:43:23)       
## ** Cognitive Diagnostic Models  **
## **********************************

## * TAM 4.2-21 (2024-02-19 18:52:08)

## Loading required package: stats4

## Loading required package: lattice

## Loading required package: MASS

## Loading required package: msm

## Loading required package: polycor

## 
## Attaching package: 'ltm'

## The following object is masked from 'package:mirt':
## 
##     Science

# 載入測驗資料
data <- read_xlsx("IQ.xlsx")

# 查看資料的結構與前幾行
head(data)

## # A tibble: 6 × 22
##      ID    IQ Item1 Item2 Item3 Item4 Item5 Item6 Item7 Item8 Item9 Item10
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1   112     0     1     1     0     1     1     1     0     1      0
## 2     2    90     1     1     0     0     1     1     1     0     0      0
## 3     3   142     1     1     1     1     1     1     1     1     1      1
## 4     4    86     1     1     0     1     1     0     0     1     1      0
## 5     5    89     1     1     0     1     0     1     0     0     0      1
## 6     6    86     1     1     1     0     0     1     0     0     1      0
## # ℹ 10 more variables: Item11 <dbl>, Item12 <dbl>, Item13 <dbl>, Item14 <dbl>,
## #   Item15 <dbl>, Item16 <dbl>, Item17 <dbl>, Item18 <dbl>, Item19 <dbl>,
## #   Item20 <dbl>

# 提取測驗作答部分 (假設題目欄位為第 3 到第 22 列)
data_irt <- data[, 3:22]

# 使用 Rasch 模型進行能力估計
data_Mirt <- mirt(data_irt,
                  1,
                  itemtype = 'Rasch')

## Iteration: 1, Log-Lik: -3295.881, Max-Change: 0.03525Iteration: 2, Log-Lik: -3295.669, Max-Change: 0.01993Iteration: 3, Log-Lik: -3295.614, Max-Change: 0.01067Iteration: 4, Log-Lik: -3295.597, Max-Change: 0.00592Iteration: 5, Log-Lik: -3295.592, Max-Change: 0.00319Iteration: 6, Log-Lik: -3295.590, Max-Change: 0.00178Iteration: 7, Log-Lik: -3295.589, Max-Change: 0.00111Iteration: 8, Log-Lik: -3295.589, Max-Change: 0.00057Iteration: 9, Log-Lik: -3295.589, Max-Change: 0.00035Iteration: 10, Log-Lik: -3295.588, Max-Change: 0.00029Iteration: 11, Log-Lik: -3295.588, Max-Change: 0.00011Iteration: 12, Log-Lik: -3295.588, Max-Change: 0.00005

# 模型摘要
summary(data_Mirt)

##           F1    h2
## Item1  0.491 0.241
## Item2  0.491 0.241
## Item3  0.491 0.241
## Item4  0.491 0.241
## Item5  0.491 0.241
## Item6  0.491 0.241
## Item7  0.491 0.241
## Item8  0.491 0.241
## Item9  0.491 0.241
## Item10 0.491 0.241
## Item11 0.491 0.241
## Item12 0.491 0.241
## Item13 0.491 0.241
## Item14 0.491 0.241
## Item15 0.491 0.241
## Item16 0.491 0.241
## Item17 0.491 0.241
## Item18 0.491 0.241
## Item19 0.491 0.241
## Item20 0.491 0.241
## 
## SS loadings:  4.824 
## Proportion Var:  0.241 
## 
## Factor correlations: 
## 
##    F1
## F1  1

# 進行每個題目的擬合檢查與繪圖
for (i in 1:length(data_irt)) {
  # 針對 Rasch 模型計算題目擬合指標，並繪製題目特徵曲線
  itemplot <- itemfit(data_Mirt,
                      group.bins = 15,        
                      # 分組的數量（將數據分成15組）
                      empirical.plot = i,     
                      # 繪製第 i 題的經驗曲線
                      empirical.CI = 0.95,    
                      # 設置95%置信區間
                      method = 'ML')         
  #使用最大似然法估計
  print(itemplot) # 輸出結果
}

# 提取潛在能力估計分數 (latent trait scores)
latent_IQ <- as.vector(fscores(data_Mirt)) # 提取每個受試者的能力分數

# 重新計算總分 (測驗作答的總得分)
sum_score <- apply(data_irt, 1, sum)

# 將結果整合成新的資料框
IQ_measures <- data.frame(ID = data[, 1], # 假設第一欄為受試者 ID
                          IQ = data[, 2], # 假設第二欄為 IQ 分數
                          sum_score = sum_score,
                          latent_IQ = latent_IQ)

# 查看整合後的資料
head(IQ_measures)

##   ID  IQ sum_score  latent_IQ
## 1  1 112        14  0.1663636
## 2  2  90        10 -0.6508563
## 3  3 142        19  1.4905519
## 4  4  86         9 -0.8450976
## 5  5  89         9 -0.8450976
## 6  6  86        10 -0.6508563

# IQ 與潛在能力分數的散佈圖
plot(IQ_measures$IQ, IQ_measures$latent_IQ,
     xlab = "IQ",
     ylab = "Latent Trait (Ability)",
     main = "Scatterplot: IQ vs Latent Ability",
     col = "blue", pch = 19)

# IQ 與總得分的散佈圖
plot(IQ_measures$IQ, IQ_measures$sum_score,
     xlab = "IQ",
     ylab = "Sum Score",
     main = "Scatterplot: IQ vs Sum Score",
     col = "red", pch = 19)

# 計算相關係數
cor_IQ_Latent <- cor(IQ_measures$IQ, IQ_measures$latent_IQ)
cat("Correlation between IQ and Latent Ability:", cor_IQ_Latent, "\n")

## Correlation between IQ and Latent Ability: 0.8726615

cor_IQ_SumScore <- cor(IQ_measures$IQ, IQ_measures$sum_score)
cat("Correlation between IQ and Sum Score:", cor_IQ_SumScore, "\n")

## Correlation between IQ and Sum Score: 0.8661066

# 使用 IQ 預測潛在能力分數
lm_IQ_Latent <- lm(latent_IQ ~ IQ, data = IQ_measures)
summary(lm_IQ_Latent)

## 
## Call:
## lm(formula = latent_IQ ~ IQ, data = IQ_measures)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.46208 -0.27398  0.01064  0.25525  1.16320 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.849034   0.158932  -30.51   <2e-16 ***
## IQ           0.048687   0.001578   30.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4073 on 298 degrees of freedom
## Multiple R-squared:  0.7615, Adjusted R-squared:  0.7607 
## F-statistic: 951.7 on 1 and 298 DF,  p-value: < 2.2e-16

# 可視化回歸線
plot(IQ_measures$IQ, IQ_measures$latent_IQ,
     xlab = "IQ",
     ylab = "Latent Trait (Ability)",
     main = "Regression: Latent Ability vs IQ",
     col = "blue", pch = 19)
abline(lm_IQ_Latent, col = "red", lwd = 2)

# 使用 IQ 預測總得分
lm_IQ_sum <- lm(sum_score ~ IQ, data = IQ_measures)
summary(lm_IQ_sum)

## 
## Call:
## lm(formula = sum_score ~ IQ, data = IQ_measures)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1535 -1.2134  0.1599  1.2274  5.2843 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.860013   0.736941  -12.02   <2e-16 ***
## IQ           0.218890   0.007318   29.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.889 on 298 degrees of freedom
## Multiple R-squared:  0.7501, Adjusted R-squared:  0.7493 
## F-statistic: 894.7 on 1 and 298 DF,  p-value: < 2.2e-16

# 可視化回歸線
plot(IQ_measures$IQ, IQ_measures$sum_score,
     xlab = "IQ",
     ylab = "sum_score",
     main = "Regression: sum_score vs IQ",
     col = "red", pch = 19)
abline(lm_IQ_sum, col = "blue", lwd = 2)

# 潛在能力分數的分佈
hist(IQ_measures$latent_IQ,
     breaks = 20,
     col = "skyblue",
     main = "Distribution of Latent Trait (Ability)",
     xlab = "Latent Trait (Ability)")

# IQ 分數的分佈
hist(IQ_measures$IQ,
     breaks = 20,
     col = "pink",
     main = "Distribution of IQ Scores",
     xlab = "IQ")

# 提取 Rasch 模型的題目參數
item_params_mirt <- coef(data_Mirt, IRTpars = TRUE)

# 檢查 item_params 的結構（確保參數名稱與位置正確）
str(item_params_mirt)

## List of 21
##  $ Item1    : num [1, 1:4] 1 -1.97 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item2    : num [1, 1:4] 1 -2.21 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item3    : num [1, 1:4] 1 -1.68 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item4    : num [1, 1:4] 1 -1.45 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item5    : num [1, 1:4] 1 -1.13 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item6    : num [1, 1:4] 1 -1.09 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item7    : num [1, 1:4] 1 -0.592 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item8    : num [1, 1:4] 1 0.152 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item9    : num [1, 1:4] 1 -0.492 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item10   : num [1, 1:4] 1 -0.248 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item11   : num [1, 1:4] 1 -0.0396 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item12   : num [1, 1:4] 1 0.495 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item13   : num [1, 1:4] 1 0.314 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item14   : num [1, 1:4] 1 1.32 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item15   : num [1, 1:4] 1 -1.97 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item16   : num [1, 1:4] 1 -1.52 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item17   : num [1, 1:4] 1 -1.59 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item18   : num [1, 1:4] 1 -1.3 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item19   : num [1, 1:4] 1 -0.78 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ Item20   : num [1, 1:4] 1 -0.329 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:4] "a" "b" "g" "u"
##  $ GroupPars: num [1, 1:2] 0 0.921
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "par"
##   .. ..$ : chr [1:2] "MEAN_1" "COV_11"
##  - attr(*, "class")= chr [1:2] "mirt_list" "list"

# 提取難度參數 b（矩陣形式的第 2 列）
difficulty_params_mirt <- sapply(item_params_mirt[names(item_params_mirt) != "GroupPars"], function(x) x[1, "b"])

# 查看提取結果
print(difficulty_params_mirt)

##       Item1       Item2       Item3       Item4       Item5       Item6 
## -1.96741931 -2.20839960 -1.68150305 -1.45254683 -1.12640448 -1.08821186 
##       Item7       Item8       Item9      Item10      Item11      Item12 
## -0.59196800  0.15249816 -0.49211113 -0.24788802 -0.03958903  0.49502092 
##      Item13      Item14      Item15      Item16      Item17      Item18 
##  0.31405468  1.32038029 -1.96741931 -1.51876765 -1.58700354 -1.30483281 
##      Item19      Item20 
## -0.77988298 -0.32863191

# 1. 執行 Rasch 模型
model_TAM <- tam(data_irt)

## ....................................................
## Processing Data      2024-12-01 21:06:32.363067 
##     * Response Data: 300 Persons and  20 Items 
##     * Numerical integration with 21 nodes
##     * Created Design Matrices   ( 2024-12-01 21:06:32.364285 )
##     * Calculated Sufficient Statistics   ( 2024-12-01 21:06:32.36706 )
## ....................................................
## Iteration 1     2024-12-01 21:06:32.368153
## E Step
## M Step Intercepts   |---
##   Deviance = 6611.5511
##   Maximum item intercept parameter change: 0.275712
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.093087
## ....................................................
## Iteration 2     2024-12-01 21:06:32.37326
## E Step
## M Step Intercepts   |---
##   Deviance = 6592.7838 | Absolute change: 18.7673 | Relative change: 0.00284665
##   Maximum item intercept parameter change: 0.020698
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000403
## ....................................................
## Iteration 3     2024-12-01 21:06:32.373719
## E Step
## M Step Intercepts   |--
##   Deviance = 6592.1091 | Absolute change: 0.6747 | Relative change: 0.00010234
##   Maximum item intercept parameter change: 0.014873
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.001499
## ....................................................
## Iteration 4     2024-12-01 21:06:32.374083
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.7179 | Absolute change: 0.3912 | Relative change: 5.935e-05
##   Maximum item intercept parameter change: 0.011448
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.001963
## ....................................................
## Iteration 5     2024-12-01 21:06:32.374438
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.4908 | Absolute change: 0.2271 | Relative change: 3.445e-05
##   Maximum item intercept parameter change: 0.008797
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.001952
## ....................................................
## Iteration 6     2024-12-01 21:06:32.374795
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.3587 | Absolute change: 0.1321 | Relative change: 2.004e-05
##   Maximum item intercept parameter change: 0.006751
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.001741
## ....................................................
## Iteration 7     2024-12-01 21:06:32.375146
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.2817 | Absolute change: 0.077 | Relative change: 1.168e-05
##   Maximum item intercept parameter change: 0.005177
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.001468
## ....................................................
## Iteration 8     2024-12-01 21:06:32.375483
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.2367 | Absolute change: 0.045 | Relative change: 6.82e-06
##   Maximum item intercept parameter change: 0.003969
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.001197
## ....................................................
## Iteration 9     2024-12-01 21:06:32.375826
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.2105 | Absolute change: 0.0263 | Relative change: 3.99e-06
##   Maximum item intercept parameter change: 0.003042
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000957
## ....................................................
## Iteration 10     2024-12-01 21:06:32.376163
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1951 | Absolute change: 0.0154 | Relative change: 2.34e-06
##   Maximum item intercept parameter change: 0.002332
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000754
## ....................................................
## Iteration 11     2024-12-01 21:06:32.376519
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.186 | Absolute change: 0.009 | Relative change: 1.37e-06
##   Maximum item intercept parameter change: 0.001788
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000589
## ....................................................
## Iteration 12     2024-12-01 21:06:32.376867
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1808 | Absolute change: 0.0053 | Relative change: 8e-07
##   Maximum item intercept parameter change: 0.00137
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000458
## ....................................................
## Iteration 13     2024-12-01 21:06:32.377228
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1777 | Absolute change: 0.0031 | Relative change: 4.7e-07
##   Maximum item intercept parameter change: 0.00105
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000354
## ....................................................
## Iteration 14     2024-12-01 21:06:32.377584
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1758 | Absolute change: 0.0018 | Relative change: 2.8e-07
##   Maximum item intercept parameter change: 0.000805
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000273
## ....................................................
## Iteration 15     2024-12-01 21:06:32.37794
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1748 | Absolute change: 0.0011 | Relative change: 1.6e-07
##   Maximum item intercept parameter change: 0.000617
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.00021
## ....................................................
## Iteration 16     2024-12-01 21:06:32.378275
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1741 | Absolute change: 6e-04 | Relative change: 1e-07
##   Maximum item intercept parameter change: 0.000473
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000162
## ....................................................
## Iteration 17     2024-12-01 21:06:32.378601
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1738 | Absolute change: 4e-04 | Relative change: 6e-08
##   Maximum item intercept parameter change: 0.000363
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 0.000124
## ....................................................
## Iteration 18     2024-12-01 21:06:32.378929
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1735 | Absolute change: 2e-04 | Relative change: 3e-08
##   Maximum item intercept parameter change: 0.000278
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 9.5e-05
## ....................................................
## Iteration 19     2024-12-01 21:06:32.379254
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1734 | Absolute change: 1e-04 | Relative change: 2e-08
##   Maximum item intercept parameter change: 0.000213
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 7.3e-05
## ....................................................
## Iteration 20     2024-12-01 21:06:32.37959
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1733 | Absolute change: 1e-04 | Relative change: 1e-08
##   Maximum item intercept parameter change: 0.000164
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 5.6e-05
## ....................................................
## Iteration 21     2024-12-01 21:06:32.379922
## E Step
## M Step Intercepts   |--
##   Deviance = 6591.1733 | Absolute change: 0 | Relative change: 1e-08
##   Maximum item intercept parameter change: 0.000126
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 4.3e-05
## ....................................................
## Iteration 22     2024-12-01 21:06:32.380261
## E Step
## M Step Intercepts   |-
##   Deviance = 6591.1733 | Absolute change: 0 | Relative change: 0
##   Maximum item intercept parameter change: 9.6e-05
##   Maximum item slope parameter change: 0
##   Maximum regression parameter change: 0
##   Maximum variance parameter change: 3.3e-05
## ....................................................
## Item Parameters
##    xsi.index xsi.label     est
## 1          1     Item1 -1.9669
## 2          2     Item2 -2.2079
## 3          3     Item3 -1.6810
## 4          4     Item4 -1.4521
## 5          5     Item5 -1.1259
## 6          6     Item6 -1.0877
## 7          7     Item7 -0.5915
## 8          8     Item8  0.1529
## 9          9     Item9 -0.4917
## 10        10    Item10 -0.2474
## 11        11    Item11 -0.0392
## 12        12    Item12  0.4954
## 13        13    Item13  0.3145
## 14        14    Item14  1.3208
## 15        15    Item15 -1.9669
## 16        16    Item16 -1.5183
## 17        17    Item17 -1.5865
## 18        18    Item18 -1.3043
## 19        19    Item19 -0.7794
## 20        20    Item20 -0.3282
## ...................................
## Regression Coefficients
##      [,1]
## [1,]    0
## 
## Variance:
##        [,1]
## [1,] 0.9205
## 
## 
## EAP Reliability:
## [1] 0.751
## 
## -----------------------------
## Start:  2024-12-01 21:06:32.362714
## End:  2024-12-01 21:06:32.383236 
## Time difference of 0.02052188 secs

# 1-2.檢驗試題難度 
xsi_params <- model_TAM$xsi[1:1]
# 1-2-1.檢驗試題難度資料匡
print(head(xsi_params))

##             xsi
## Item1 -1.966916
## Item2 -2.207892
## Item3 -1.681005
## Item4 -1.452055
## Item5 -1.125924
## Item6 -1.087732

# 1-3.模型適配度item_fit
msq.itemfit(model_TAM)

## $itemfit
##      item fitgroup    Outfit   Outfit_t     Outfit_p     Infit    Infit_t
## 1   Item1        1 1.1204386  0.7047519 0.4809646480 1.0289322  0.3117560
## 2   Item2        2 1.2031823  0.9934260 0.3205023993 1.0688199  0.6047865
## 3   Item3        3 1.1716886  1.1226887 0.2615697181 1.0544049  0.6410243
## 4   Item4        4 1.0593392  0.4852398 0.6275062691 1.0548604  0.7301201
## 5   Item5        5 0.8787064 -1.1444457 0.2524388474 0.9473591 -0.8188566
## 6   Item6        6 0.9094273 -0.8583668 0.3906899787 0.9498326 -0.7954308
## 7   Item7        7 0.9366212 -0.8214879 0.4113684158 0.9686934 -0.6159221
## 8   Item8        8 1.2332663  3.4907245 0.0004817126 1.1743443  3.5875787
## 9   Item9        9 0.8905202 -1.5490004 0.1213816290 0.9270945 -1.5270609
## 10 Item10       10 0.9078802 -1.4455877 0.1482928891 0.9231328 -1.7013715
## 11 Item11       11 0.8802878 -1.9931768 0.0462420908 0.9093185 -2.0437908
## 12 Item12       12 0.8892420 -1.5991161 0.1097948127 0.9041320 -1.9496410
## 13 Item13       13 0.9022448 -1.5234963 0.1276345787 0.9295124 -1.4961079
## 14 Item14       14 1.0660329  0.5852387 0.5583872516 1.0318691  0.4538056
## 15 Item15       15 1.3530232  1.8297036 0.0672942795 1.1211669  1.1625283
## 16 Item16       16 0.9130229 -0.6064056 0.5442454452 0.9553171 -0.5445193
## 17 Item17       17 0.9646924 -0.2015790 0.8402458775 0.9832814 -0.1758510
## 18 Item18       18 0.9794630 -0.1319635 0.8950131606 0.9908041 -0.1077117
## 19 Item19       19 1.1034706  1.1954391 0.2319155024 1.0473035  0.8829108
## 20 Item20       20 1.0870406  1.2886026 0.1975362900 1.0707601  1.5004652
##        Infit_p
## 1  0.755225997
## 2  0.545320890
## 3  0.521506906
## 4  0.465316774
## 5  0.412868232
## 6  0.426362923
## 7  0.537945912
## 8  0.000333763
## 9  0.126745876
## 10 0.088873244
## 11 0.040974220
## 12 0.051218925
## 13 0.134625533
## 14 0.649968788
## 15 0.245020947
## 16 0.586084158
## 17 0.860410978
## 18 0.914224412
## 19 0.377284485
## 20 0.133493947
## 
## $summary_itemfit
##           fit        M         SD
## Outfit Outfit 1.022480 0.14014153
## Infit   Infit 1.002047 0.07505481
## 
## $time
## [1] "2024-12-01 21:06:32 CST" "2024-12-01 21:06:32 CST"
## 
## $CALL
## msq.itemfit(object = model_TAM)
## 
## attr(,"class")
## [1] "msq.itemfit"

# 1-4.適配度圖示化(以18題為例)
plot(model_TAM,items = 18)

## Iteration in WLE/MLE estimation  1   | Maximal change  0.9488 
## Iteration in WLE/MLE estimation  2   | Maximal change  0.1044 
## Iteration in WLE/MLE estimation  3   | Maximal change  0.0031 
## Iteration in WLE/MLE estimation  4   | Maximal change  2e-04 
## Iteration in WLE/MLE estimation  5   | Maximal change  0 
## ----
##  WLE Reliability= 0.713

## ....................................................
##  Plots exported in png format into folder:
##  /Users/Shared/作業4/Plots

# 2. 提取 EAP 能力估計和標準誤
eap_abilities <- model_TAM$person$EAP        # 提取 EAP 能力分數
eap_se <- model_TAM$person$SD.EAP           # 提取 EAP 標準誤

# 3. 整合能力估計和標準誤到資料框
person_abilities <- data.frame(
  ID = data[, 1], # 第一列是受試者 ID
  IQ = data[, 2], # 第二列是受試者 IQ
  EAP = eap_abilities,
  SE = eap_se
)

# 4. 查看結果
print(head(person_abilities))

##   ID  IQ        EAP        SE
## 1  1 112  0.1667368 0.4682997
## 2  2  90 -0.6502470 0.4413045
## 3  3 142  1.4906969 0.5791548
## 4  4  86 -0.8445068 0.4405079
## 5  5  89 -0.8445068 0.4405079
## 6  6  86 -0.6502470 0.4413045

# 計算相關性
cor_IQ_EAP <- cor(person_abilities$IQ, person_abilities$EAP)
cat("Correlation between IQ and EAP Ability:", cor_IQ_EAP, "\n")

## Correlation between IQ and EAP Ability: 0.8726578

# 能力估計分佈圖
hist(eap_abilities, breaks = 20, col = "skyblue",
     main = "Distribution of EAP Ability Estimates",
     xlab = "EAP Ability")

# 標準誤分佈圖
hist(person_abilities$SE, breaks = 20, col = "pink",
     main = "Distribution of Standard Error",
     xlab = "Standard Error")

# 能力估計與其標準誤關係圖
plot(eap_abilities, eap_se,
     xlab = "EAP Ability",
     ylab = "Standard Error",
     main = "EAP Ability vs. Standard Error",
     col = "red", pch = 19)

# 安裝並載入 ltm 套件
library(ltm)

# 假設 data_irt 是你的作答數據，進行 Rasch 模型分析
model_lTM_rasch <- rasch(data_irt)

# 設定能力值的上下限
# 假設上下限分別設定為 -4 和 4
ability_limits <- seq(-4, 4, by = 1)

# 使用 factor.scores 提取每個人的能力值，並應用上下限
rasch_scores <- factor.scores(
  model_lTM_rasch,
  resp.patterns = data_irt, # 使用原始數據
  qrs = ability_limits      # 設定能力值上下限
)

# 查看得分結果
head(rasch_scores$scores)

## NULL

# 1. 提取 z1 能力估計和標準誤
z1_abilities <- rasch_scores$score.dat$z1      # 提取 EAP 能力分數
z1_se <- rasch_scores$score.dat$se.z1           # 提取 EAP 標準誤

# 2. 整合能力估計和標準誤到資料框
person_abilities_ltm <- data.frame(
  ID = data[, 1], # 第一列是受試者 ID
  IQ = data[, 2], # 第二列是受試者 IQ
  z1 = z1_abilities,
  SE = z1_se
)

# 3. 計算相關性
cor_IQ_z1 <- cor(person_abilities_ltm$IQ, person_abilities_ltm$z1)
cat("Correlation between IQ and z1 Ability:", cor_IQ_z1, "\n")

## Correlation between IQ and z1 Ability: 0.8725787

# TAM 套件的能力估計值
eap_abilities_tam <- model_TAM$person$EAP

# mirt 套件的能力估計值
latent_abilities_mirt <- as.vector(fscores(data_Mirt))

# ltm 套件的能力估計值
latent_abilities_ltm <- rasch_scores$score.dat$z1  

# 整合能力估計值
abilities_df <- data.frame(
  TAM = eap_abilities_tam,
  mirt = latent_abilities_mirt,
  ltm = latent_abilities_ltm
)

# 檢查整合結果
head(abilities_df)

##          TAM       mirt        ltm
## 1  0.1667368  0.1663636  0.1479396
## 2 -0.6502470 -0.6508563 -0.6837206
## 3  1.4906969  1.4905519  1.4872948
## 4 -0.8445068 -0.8450976 -0.8816281
## 5 -0.8445068 -0.8450976 -0.8816281
## 6 -0.6502470 -0.6508563 -0.6837206

library(corrplot)

## corrplot 0.92 loaded

# 計算相關矩陣
correlation_matrix <- cor(abilities_df, use = "complete.obs")

# 查看相關矩陣
print(correlation_matrix)

##            TAM      mirt       ltm
## TAM  1.0000000 1.0000000 0.9999984
## mirt 1.0000000 1.0000000 0.9999983
## ltm  0.9999984 0.9999983 1.0000000

# 繪製相關矩陣圖
corrplot(correlation_matrix, type = "lower", order = "hclust", tl.col = "black", tl.cex = 0.7, tl.srt = 45)

#TAM難度
difficulty_tam <- model_TAM$xsi

# 提取完整參數表
params_table <- mod2values(data_Mirt)

# 查看參數表
head(params_table)

##   group  item class name parnum    value lbound ubound   est prior.type prior_1
## 1   all Item1  dich   a1      1 1.000000   -Inf    Inf FALSE       none     NaN
## 2   all Item1  dich    d      2 1.967419   -Inf    Inf  TRUE       none     NaN
## 3   all Item1  dich    g      3 0.000000      0      1 FALSE       none     NaN
## 4   all Item1  dich    u      4 1.000000      0      1 FALSE       none     NaN
## 5   all Item2  dich   a1      5 1.000000   -Inf    Inf FALSE       none     NaN
## 6   all Item2  dich    d      6 2.208400   -Inf    Inf  TRUE       none     NaN
##   prior_2
## 1     NaN
## 2     NaN
## 3     NaN
## 4     NaN
## 5     NaN
## 6     NaN

# 篩選出難度參數
difficulty_mirt <- params_table[params_table$name == "d", "value"]

# 查看難度參數
print(difficulty_mirt)

##  [1]  1.96741931  2.20839960  1.68150305  1.45254683  1.12640448  1.08821186
##  [7]  0.59196800 -0.15249816  0.49211113  0.24788802  0.03958903 -0.49502092
## [13] -0.31405468 -1.32038029  1.96741931  1.51876765  1.58700354  1.30483281
## [19]  0.77988298  0.32863191

# 提取 ltm 的難度參數
difficulty_ltm <- coef(model_lTM_rasch)

# 整合難度參數
difficulty_df <- data.frame(
  TAM = difficulty_tam$xsi,          # TAM 的難度參數
  mirt = difficulty_mirt,            # mirt 的難度參數
  ltm = difficulty_ltm[, "Dffclt"]   # ltm 的難度參數（通常標為 "Dffclt"）
)

# 查看整合結果
print(difficulty_df)

##                TAM        mirt         ltm
## Item1  -1.96691557  1.96741931 -2.05033928
## Item2  -2.20789226  2.20839960 -2.30144965
## Item3  -1.68100512  1.68150305 -1.75259962
## Item4  -1.45205511  1.45254683 -1.51362832
## Item5  -1.12592363  1.12640448 -1.17389365
## Item6  -1.08773240  1.08821186 -1.13396820
## Item7  -0.59150824  0.59196800 -0.61673793
## Item8   0.15292428 -0.15249816  0.15904464
## Item9  -0.49165566  0.49211113 -0.51281420
## Item10 -0.24744341  0.24788802 -0.25816758
## Item11 -0.03915397  0.03958903 -0.04117492
## Item12  0.49543112 -0.49502092  0.51606104
## Item13  0.31447326 -0.31405468  0.32733970
## Item14  1.32075438 -1.32038029  1.37627790
## Item15 -1.96691557  1.96741931 -2.05017521
## Item16 -1.51827400  1.51876765 -1.58286903
## Item17 -1.58650802  1.58700354 -1.65353771
## Item18 -1.30434576  1.30483281 -1.35982508
## Item19 -0.77941543  0.77988298 -0.81271033
## Item20 -0.32818366  0.32863191 -0.34245868

# 計算相關矩陣
correlation_matrix_diff <- cor(difficulty_df, use = "complete.obs")

# 查看相關矩陣
print(correlation_matrix_diff)

##      TAM mirt ltm
## TAM    1   -1   1
## mirt  -1    1  -1
## ltm    1   -1   1

# 繪製相關矩陣圖
corrplot(correlation_matrix_diff, type = "lower", order = "hclust", tl.col = "black", tl.cex = 0.7, tl.srt = 45)