#Load library
library(MBESS)
library(lavaan)
## This is lavaan 0.5-23.1097
## lavaan is BETA software! Please report any bugs.
##
## Attaching package: 'lavaan'
## The following object is masked from 'package:MBESS':
##
## cor2cov
library(semPlot)
library(readr)
library(psych)
library(CTT)
##
## Attaching package: 'CTT'
## The following object is masked from 'package:psych':
##
## polyserial
library(GPArotation)
library(psychometric)
## Loading required package: multilevel
## Loading required package: nlme
## Loading required package: MASS
##
## Attaching package: 'psychometric'
## The following object is masked from 'package:psych':
##
## alpha
library(cocron)
library(knitr)
library(ltm)
## Loading required package: msm
## Loading required package: polycor
##
## Attaching package: 'polycor'
## The following object is masked from 'package:CTT':
##
## polyserial
## The following object is masked from 'package:psych':
##
## polyserial
##
## Attaching package: 'ltm'
## The following object is masked from 'package:cocron':
##
## cronbach.alpha
## The following object is masked from 'package:psych':
##
## factor.scores
library(mirt)
## Loading required package: stats4
## Loading required package: lattice
##
## Attaching package: 'mirt'
## The following object is masked from 'package:ltm':
##
## Science
## The following object is masked from 'package:nlme':
##
## fixef
#Import data
SCORES <- read_csv("C:/Goran/01 FAKS/4TT/Project/Data/SCORES.csv")
## Parsed with column specification:
## cols(
## .default = col_integer()
## )
## See spec(...) for full column specifications.
data <- SCORES
data
## # A tibble: 105 × 25
## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 1 1 1 0 1 1 1 0 1 1 0 1
## 2 1 1 1 1 1 1 1 1 1 1 1 1
## 3 0 1 0 1 1 1 1 1 1 0 1 1
## 4 1 1 1 1 1 1 1 0 1 0 1 1
## 5 1 1 1 1 0 0 1 0 1 0 0 1
## 6 1 1 1 1 1 1 1 1 1 0 1 0
## 7 1 1 1 1 1 1 1 1 1 1 1 1
## 8 1 1 1 1 1 1 1 0 0 1 0 1
## 9 1 1 0 0 0 1 1 0 0 0 0 1
## 10 0 1 1 0 1 1 1 0 0 0 0 0
## # ... with 95 more rows, and 13 more variables: X13 <int>, X14 <int>,
## # X15 <int>, X16 <int>, X17 <int>, X18 <int>, X19 <int>, X20 <int>,
## # X21 <int>, X22 <int>, X23 <int>, X24 <int>, X25 <int>
#Convert to cathegorical data
data[,c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12", "X13", "X14", "X15", "X16", "X17", "X18", "X19", "X20", "X21", "X22", "X23", "X24", "X25")] <-
lapply(data[,c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10", "X11", "X12", "X13", "X14", "X15", "X16", "X17", "X18", "X19", "X20", "X21", "X22", "X23", "X24", "X25")],
ordered)
#CFA
one.model = "science =~ X13 + X3 + X8 + X11 + X4 + X15 + X17 + X22 + X25"
#Run the model
one.fit = cfa(one.model,
data = data,
estimator = "wlsmv")
semPaths(one.fit, whatLabels = "std", layout = "tree", intercepts = TRUE, residuals = TRUE)
#Modindicies, cutoff = >6.645 (rule of thumb)
modindices(one.fit, standardized = TRUE, power = FALSE,
delta = 0.1,
alpha = 0.05,
high.power = 0.75,
sort. = TRUE,
minimum.value = 6.635)
## [1] lhs op rhs mi mi.scaled epc sepc.lv
## [8] sepc.all sepc.nox
## <0 rows> (or 0-length row.names)
#Summary
summary(one.fit, standardized = TRUE, fit.measures = T, rsquare = TRUE)
## lavaan (0.5-23.1097) converged normally after 19 iterations
##
## Number of observations 105
##
## Estimator DWLS Robust
## Minimum Function Test Statistic 16.112 23.446
## Degrees of freedom 27 27
## P-value (Chi-square) 0.951 0.661
## Scaling correction factor 0.882
## Shift parameter 5.169
## for simple second-order correction (Mplus variant)
##
## Model test baseline model:
##
## Minimum Function Test Statistic 193.536 155.229
## Degrees of freedom 36 36
## P-value 0.000 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000 1.000
## Tucker-Lewis Index (TLI) 1.092 1.040
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.000
## 90 Percent Confidence Interval 0.000 0.000 0.000 0.063
## P-value RMSEA <= 0.05 0.992 0.889
##
## Robust RMSEA NA
## 90 Percent Confidence Interval 0.000 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.093 0.093
##
## Weighted Root Mean Square Residual:
##
## WRMR 0.598 0.598
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Robust.sem
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## science =~
## X13 1.000 0.732 0.732
## X3 0.793 0.229 3.461 0.001 0.581 0.581
## X8 0.697 0.223 3.120 0.002 0.510 0.510
## X11 0.745 0.207 3.600 0.000 0.546 0.546
## X4 0.899 0.212 4.232 0.000 0.658 0.658
## X15 0.760 0.234 3.247 0.001 0.556 0.556
## X17 0.898 0.217 4.140 0.000 0.658 0.658
## X22 0.686 0.190 3.616 0.000 0.502 0.502
## X25 0.713 0.206 3.464 0.001 0.522 0.522
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X13 0.000 0.000 0.000
## .X3 0.000 0.000 0.000
## .X8 0.000 0.000 0.000
## .X11 0.000 0.000 0.000
## .X4 0.000 0.000 0.000
## .X15 0.000 0.000 0.000
## .X17 0.000 0.000 0.000
## .X22 0.000 0.000 0.000
## .X25 0.000 0.000 0.000
## science 0.000 0.000 0.000
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X13|t1 -0.876 0.142 -6.184 0.000 -0.876 -0.876
## X3|t1 -1.026 0.150 -6.861 0.000 -1.026 -1.026
## X8|t1 -0.108 0.123 -0.874 0.382 -0.108 -0.108
## X11|t1 -0.084 0.123 -0.680 0.497 -0.084 -0.084
## X4|t1 -0.253 0.124 -2.038 0.042 -0.253 -0.253
## X15|t1 0.623 0.132 4.720 0.000 0.623 0.623
## X17|t1 0.060 0.123 0.486 0.627 0.060 0.060
## X22|t1 0.594 0.131 4.532 0.000 0.594 0.594
## X25|t1 -0.405 0.127 -3.196 0.001 -0.405 -0.405
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X13 0.464 0.464 0.464
## .X3 0.663 0.663 0.663
## .X8 0.740 0.740 0.740
## .X11 0.702 0.702 0.702
## .X4 0.567 0.567 0.567
## .X15 0.690 0.690 0.690
## .X17 0.567 0.567 0.567
## .X22 0.748 0.748 0.748
## .X25 0.727 0.727 0.727
## science 0.536 0.185 2.897 0.004 1.000 1.000
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## X13 1.000 1.000 1.000
## X3 1.000 1.000 1.000
## X8 1.000 1.000 1.000
## X11 1.000 1.000 1.000
## X4 1.000 1.000 1.000
## X15 1.000 1.000 1.000
## X17 1.000 1.000 1.000
## X22 1.000 1.000 1.000
## X25 1.000 1.000 1.000
##
## R-Square:
## Estimate
## X13 0.536
## X3 0.337
## X8 0.260
## X11 0.298
## X4 0.433
## X15 0.310
## X17 0.433
## X22 0.252
## X25 0.273
fitmeasures(one.fit)
## npar fmin
## 18.000 0.077
## chisq df
## 16.112 27.000
## pvalue chisq.scaled
## 0.951 23.446
## df.scaled pvalue.scaled
## 27.000 0.661
## chisq.scaling.factor baseline.chisq
## 0.882 193.536
## baseline.df baseline.pvalue
## 36.000 0.000
## baseline.chisq.scaled baseline.df.scaled
## 155.229 36.000
## baseline.pvalue.scaled baseline.chisq.scaling.factor
## 0.000 1.324
## cfi tli
## 1.000 1.092
## nnfi rfi
## 1.092 0.889
## nfi pnfi
## 0.917 0.688
## ifi rni
## 1.065 1.069
## cfi.scaled tli.scaled
## 1.000 1.040
## cfi.robust tli.robust
## NA NA
## nnfi.scaled nnfi.robust
## 1.040 NA
## rfi.scaled nfi.scaled
## 0.799 0.849
## ifi.scaled rni.scaled
## 0.849 1.023
## rni.robust rmsea
## NA 0.000
## rmsea.ci.lower rmsea.ci.upper
## 0.000 0.000
## rmsea.pvalue rmsea.scaled
## 0.992 0.000
## rmsea.ci.lower.scaled rmsea.ci.upper.scaled
## 0.000 0.063
## rmsea.pvalue.scaled rmsea.robust
## 0.889 NA
## rmsea.ci.lower.robust rmsea.ci.upper.robust
## 0.000 NA
## rmsea.pvalue.robust rmr
## NA 0.085
## rmr_nomean srmr
## 0.093 0.093
## srmr_bentler srmr_bentler_nomean
## 0.085 0.093
## srmr_bollen srmr_bollen_nomean
## 0.085 0.093
## srmr_mplus srmr_mplus_nomean
## 0.085 0.093
## cn_05 cn_01
## 259.926 304.140
## gfi agfi
## 0.952 0.920
## pgfi mfi
## 0.571 1.054
#Import data 2
ctt.data <- read.csv("C:/Goran/01 FAKS/4TT/Project/Data/SCORES.csv", header = TRUE)
str(ctt.data)
## 'data.frame': 105 obs. of 25 variables:
## $ X1 : int 1 1 0 1 1 1 1 1 1 0 ...
## $ X2 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ X3 : int 1 1 0 1 1 1 1 1 0 1 ...
## $ X4 : int 0 1 1 1 1 1 1 1 0 0 ...
## $ X5 : int 1 1 1 1 0 1 1 1 0 1 ...
## $ X6 : int 1 1 1 1 0 1 1 1 1 1 ...
## $ X7 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ X8 : int 0 1 1 0 0 1 1 0 0 0 ...
## $ X9 : int 1 1 1 1 1 1 1 0 0 0 ...
## $ X10: int 1 1 0 0 0 0 1 1 0 0 ...
## $ X11: int 0 1 1 1 0 1 1 0 0 0 ...
## $ X12: int 1 1 1 1 1 0 1 1 1 0 ...
## $ X13: int 0 1 0 1 1 1 0 1 0 1 ...
## $ X14: int 1 1 1 1 1 1 1 1 1 1 ...
## $ X15: int 0 1 1 1 0 0 0 1 1 0 ...
## $ X16: int 1 1 0 1 0 1 1 1 1 1 ...
## $ X17: int 0 1 1 1 0 1 1 0 0 1 ...
## $ X18: int 0 1 1 1 0 1 1 1 1 1 ...
## $ X19: int 0 1 1 1 1 1 1 1 1 1 ...
## $ X20: int 1 1 1 1 0 1 1 1 1 1 ...
## $ X21: int 0 1 1 1 1 1 1 1 1 1 ...
## $ X22: int 0 0 0 0 0 0 1 0 0 0 ...
## $ X23: int 0 0 1 1 1 1 1 0 1 0 ...
## $ X24: int 0 0 1 0 0 1 1 1 0 0 ...
## $ X25: int 0 1 1 1 1 1 0 1 1 0 ...
#Matrix 2
responses <- as.matrix(ctt.data[,c( -1, -2, -5, -6, -7, -9, -10, -12, -14, -16, -18, -19, -20, -21, -23, -24)])
responses
## X3 X4 X8 X11 X13 X15 X17 X22 X25
## [1,] 1 0 0 0 0 0 0 0 0
## [2,] 1 1 1 1 1 1 1 0 1
## [3,] 0 1 1 1 0 1 1 0 1
## [4,] 1 1 0 1 1 1 1 0 1
## [5,] 1 1 0 0 1 0 0 0 1
## [6,] 1 1 1 1 1 0 1 0 1
## [7,] 1 1 1 1 0 0 1 1 0
## [8,] 1 1 0 0 1 1 0 0 1
## [9,] 0 0 0 0 0 1 0 0 1
## [10,] 1 0 0 0 1 0 1 0 0
## [11,] 1 1 1 0 1 0 1 0 1
## [12,] 1 1 0 0 0 0 1 1 0
## [13,] 1 0 1 1 1 0 1 1 0
## [14,] 1 0 0 1 1 0 1 0 1
## [15,] 1 0 1 0 1 0 1 0 1
## [16,] 1 0 0 0 1 0 0 1 1
## [17,] 1 1 1 1 1 0 0 0 1
## [18,] 1 1 0 0 1 0 0 0 0
## [19,] 1 1 0 1 1 0 0 0 1
## [20,] 0 0 1 1 1 0 1 1 0
## [21,] 1 1 1 1 1 0 1 1 1
## [22,] 1 1 1 0 1 0 1 1 1
## [23,] 1 1 0 1 1 1 1 0 1
## [24,] 1 1 1 1 1 1 0 1 1
## [25,] 1 0 1 1 1 0 0 1 1
## [26,] 1 1 1 0 1 1 1 0 0
## [27,] 1 0 1 1 1 0 1 1 1
## [28,] 1 0 0 1 1 0 1 1 1
## [29,] 1 1 0 0 0 0 0 0 0
## [30,] 1 1 1 1 1 1 1 1 1
## [31,] 1 1 1 1 1 1 0 1 0
## [32,] 1 0 0 0 1 0 0 0 0
## [33,] 1 0 1 1 1 1 1 0 1
## [34,] 1 1 0 1 1 0 0 0 0
## [35,] 0 0 0 0 1 0 0 0 0
## [36,] 1 1 1 0 1 0 1 0 1
## [37,] 1 1 1 0 1 1 0 0 0
## [38,] 1 1 1 1 1 1 1 0 1
## [39,] 1 0 1 1 1 0 0 0 0
## [40,] 1 0 0 0 1 0 1 0 1
## [41,] 1 1 1 1 0 0 0 0 1
## [42,] 1 0 0 0 0 0 0 0 0
## [43,] 1 1 1 0 1 0 0 0 0
## [44,] 1 0 0 1 1 0 1 0 0
## [45,] 1 1 0 0 0 0 0 1 1
## [46,] 1 0 0 0 1 0 0 0 0
## [47,] 1 1 1 1 1 0 1 0 1
## [48,] 1 1 1 1 1 1 1 0 1
## [49,] 0 0 0 0 0 0 0 0 1
## [50,] 1 0 0 0 1 0 1 0 1
## [51,] 1 1 1 1 1 0 1 1 1
## [52,] 1 0 1 0 1 1 0 0 1
## [53,] 1 1 0 1 1 1 1 0 1
## [54,] 0 0 0 0 1 0 0 0 1
## [55,] 0 0 1 0 0 0 0 0 0
## [56,] 1 0 1 0 0 0 0 0 0
## [57,] 1 1 1 1 1 0 0 0 0
## [58,] 1 0 1 1 1 0 0 0 1
## [59,] 1 1 1 1 1 0 0 0 0
## [60,] 0 0 0 0 0 0 0 0 0
## [61,] 1 1 1 0 1 0 0 0 0
## [62,] 1 1 0 0 1 0 0 0 0
## [63,] 1 1 1 1 1 0 0 0 1
## [64,] 0 0 0 1 0 0 1 0 0
## [65,] 1 1 1 1 1 0 1 1 1
## [66,] 1 1 0 0 1 1 1 1 1
## [67,] 1 0 0 0 1 0 0 0 1
## [68,] 1 1 1 0 1 0 1 1 1
## [69,] 1 0 0 1 1 0 0 0 1
## [70,] 0 1 1 1 1 0 1 1 1
## [71,] 1 1 1 1 1 0 0 0 1
## [72,] 1 1 1 1 1 0 1 0 0
## [73,] 1 1 1 1 1 1 0 0 1
## [74,] 1 1 0 1 1 1 1 0 1
## [75,] 1 1 0 1 1 1 1 1 1
## [76,] 1 0 1 1 1 0 0 0 1
## [77,] 1 1 1 0 1 1 0 0 1
## [78,] 1 1 1 0 1 0 1 1 1
## [79,] 1 1 1 1 1 1 1 1 0
## [80,] 1 0 0 0 0 0 0 0 0
## [81,] 1 1 1 0 1 1 1 1 1
## [82,] 1 1 1 1 1 0 1 0 1
## [83,] 0 0 1 0 1 0 0 0 1
## [84,] 1 1 1 1 1 1 1 0 1
## [85,] 0 1 1 1 0 0 0 0 0
## [86,] 1 1 1 0 1 0 1 0 1
## [87,] 0 1 0 0 1 0 0 1 1
## [88,] 1 1 1 1 1 0 0 1 1
## [89,] 0 0 0 1 1 0 1 0 0
## [90,] 1 1 1 0 1 1 1 1 1
## [91,] 1 0 1 1 1 1 1 0 1
## [92,] 0 0 0 0 0 0 0 0 0
## [93,] 0 0 0 0 1 0 1 1 1
## [94,] 1 1 1 1 1 0 0 0 1
## [95,] 1 1 0 1 1 0 1 1 1
## [96,] 1 0 1 0 1 0 0 0 1
## [97,] 1 0 1 1 0 0 0 0 1
## [98,] 1 0 0 0 1 1 0 0 0
## [99,] 1 1 0 0 0 0 0 0 1
## [100,] 1 0 0 0 0 0 0 0 1
## [101,] 1 1 0 1 1 1 1 1 1
## [102,] 1 0 0 1 1 0 0 0 0
## [103,] 1 1 0 0 1 0 0 0 1
## [104,] 1 1 0 1 1 0 0 0 1
## [105,] 1 1 0 1 1 0 1 0 0
dimnames(responses) <- NULL
(N <- dim(responses) [2])
## [1] 9
(K <- dim(responses) [1])
## [1] 105
#Alpha Ver. 1
cronbachs.alpha <-
function(X){
X <- data.matrix(responses)
n <- ncol(X) #number of items
k <- nrow(X) #number of examinees
#formula
alphaC <- (n/(n - 1))*(1 - sum(apply(
X, 2, var))/var(rowSums(X)))
return(list("Cronbach's alpha" = alphaC,
"Number of items" = n,
"Number of examinees" = k))
}
dump("cronbachs.alpha", file = "cronbachs.alpha.R")
#Alpha Ver. 2
cronbachs.alpha(responses)
## $`Cronbach's alpha`
## [1] 0.6774392
##
## $`Number of items`
## [1] 9
##
## $`Number of examinees`
## [1] 105
#KR-20
KR20 <-
function(X){
X <- data.matrix(X)
k <- ncol(X)
# Person total score variances
SX <- var(rowSums(X))
# item means
IM <- colMeans(X)
return(((k/(k - 1))*((SX - sum(IM*(1 - IM)))/SX)))
}
dump("KR20", file = "KR20.R")
KR20(responses)
## [1] 0.6817017
# Kuder-Richardson formula 21
KR21 <-
function(X){
X <- data.matrix(X)
n <- ncol(X)
return((n/(n-1))*((var(rowSums(X)) - n*(sum(colMeans(X))/n) *
(1-(sum(colMeans(X))/n))))/var(rowSums(X)))
}
dump("KR21", file = "KR21.R")
KR21(responses)
## [1] 0.6041985
#Covariances
covar = cov(responses, method = "pearson")
covar
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.130402930 0.05384615 0.02582418 0.02435897 0.04761905 0.02179487
## [2,] 0.053846154 0.24230769 0.06538462 0.05192308 0.03846154 0.05961538
## [3,] 0.025824176 0.06538462 0.25054945 0.07307692 0.03708791 0.02692308
## [4,] 0.024358974 0.05192308 0.07307692 0.25128205 0.04487179 0.02948718
## [5,] 0.047619048 0.03846154 0.03708791 0.04487179 0.15567766 0.03205128
## [6,] 0.021794872 0.05961538 0.02692308 0.02948718 0.03205128 0.19743590
## [7,] 0.015567766 0.04807692 0.03708791 0.07051282 0.05311355 0.05448718
## [8,] 0.004029304 0.04423077 0.03131868 0.02435897 0.02426740 0.01217949
## [9,] 0.024175824 0.04423077 0.03406593 0.03076923 0.04945055 0.04423077
## [,7] [,8] [,9]
## [1,] 0.01556777 0.004029304 0.02417582
## [2,] 0.04807692 0.044230769 0.04423077
## [3,] 0.03708791 0.031318681 0.03406593
## [4,] 0.07051282 0.024358974 0.03076923
## [5,] 0.05311355 0.024267399 0.04945055
## [6,] 0.05448718 0.012179487 0.04423077
## [7,] 0.25183150 0.078754579 0.04945055
## [8,] 0.07875458 0.201831502 0.03791209
## [9,] 0.04945055 0.037912088 0.22747253
#the problem with corr tables is that they dont have a cutoff value. We could risk it and declare 0.1 and lower as being a bad value
#Correlation
cor(x = responses, method = "pearson")
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.00000000 0.3029196 0.1428684 0.1345658 0.3342134 0.13583059
## [2,] 0.30291963 1.0000000 0.2653660 0.2104244 0.1980295 0.27255973
## [3,] 0.14286836 0.2653660 1.0000000 0.2912412 0.1877900 0.12105003
## [4,] 0.13456577 0.2104244 0.2912412 1.0000000 0.2268713 0.13238520
## [5,] 0.33421342 0.1980295 0.1877900 0.2268713 1.0000000 0.18281811
## [6,] 0.13583059 0.2725597 0.1210500 0.1323852 0.1828181 1.00000000
## [7,] 0.08590681 0.1946247 0.1476490 0.2803060 0.2682484 0.24435782
## [8,] 0.02483659 0.2000076 0.1392715 0.1081643 0.1369038 0.06101287
## [9,] 0.14036962 0.1883981 0.1426951 0.1286979 0.2627807 0.20871178
## [,7] [,8] [,9]
## [1,] 0.08590681 0.02483659 0.1403696
## [2,] 0.19462474 0.20000756 0.1883981
## [3,] 0.14764904 0.13927150 0.1426951
## [4,] 0.28030596 0.10816426 0.1286979
## [5,] 0.26824843 0.13690384 0.2627807
## [6,] 0.24435782 0.06101287 0.2087118
## [7,] 1.00000000 0.34932230 0.2066101
## [8,] 0.34932230 1.00000000 0.1769370
## [9,] 0.20661013 0.17693704 1.0000000
#CTT
#Reliability
reliability(responses, itemal=TRUE)
##
## Number of Items
## 9
##
## Number of Examinees
## 105
##
## Coefficient Alpha
## 0.677
str(reliability(responses,itemal=TRUE), vec.len = 15) #vec.len numeric (>= 0) indicating how many 'first few' elements are displayed of each vector.
## List of 9
## $ nItem : int 9
## $ nPerson : int 105
## $ alpha : num 0.677
## $ scaleMean : num 5.01
## $ scaleSD : num 2.19
## $ alphaIfDeleted: num [1:9(1d)] 0.663 0.634 0.655 0.651 0.641 0.659 0.635 0.665 0.655
## $ pBis : num [1:9(1d)] 0.292 0.426 0.335 0.355 0.415 0.314 0.42 0.283 0.332
## $ bis : num [1:9(1d)] 0.419 0.516 0.405 0.432 0.562 0.443 0.515 0.404 0.409
## $ itemMean : num [1:9] 0.848 0.6 0.543 0.533 0.81 0.267 0.476 0.276 0.657
## - attr(*, "class")= chr "reliability"
#Alpha and Omega
ci.reliability(data = responses, type = "1", interval.type = "bca", B = 1000)
## $est
## [1] 0.6774392
##
## $se
## [1] 0.04554559
##
## $ci.lower
## [1] 0.5801351
##
## $ci.upper
## [1] 0.7594718
##
## $conf.level
## [1] 0.95
##
## $type
## [1] "alpha"
##
## $interval.type
## [1] "bca bootstrap"
ci.reliability(data = responses, type = "3", interval.type = "bca", B = 1000) #Omega
## $est
## [1] 0.6806404
##
## $se
## [1] 0.04758102
##
## $ci.lower
## [1] 0.576643
##
## $ci.upper
## [1] 0.7545325
##
## $conf.level
## [1] 0.95
##
## $type
## [1] "omega"
##
## $interval.type
## [1] "bca bootstrap"
ci.reliability(data = responses, type = "5", interval.type = "bca", B = 1000) #Categorical omega!
## $est
## [1] 0.7043666
##
## $se
## [1] 0.04793346
##
## $ci.lower
## [1] 0.484424
##
## $ci.upper
## [1] 0.7620391
##
## $conf.level
## [1] 0.95
##
## $type
## [1] "categorical omega"
##
## $interval.type
## [1] "bca bootstrap"
#Item analysis
knitr::kable(item.exam(responses, y = NULL, discrim = TRUE),
align = "c",
caption = "Item Analysis")
Item Analysis
| 0.3611135 |
0.4394712 |
0.2923590 |
0.8476190 |
0.2857143 |
NA |
0.1579415 |
0.1050708 |
NA |
| 0.4922476 |
0.6010540 |
0.4260097 |
0.6000000 |
0.6285714 |
NA |
0.2944551 |
0.2087013 |
NA |
| 0.5005491 |
0.5301980 |
0.3352222 |
0.5428571 |
0.6000000 |
NA |
0.2641234 |
0.1669942 |
NA |
| 0.5012804 |
0.5470220 |
0.3552830 |
0.5333333 |
0.6000000 |
NA |
0.2729025 |
0.1772463 |
NA |
| 0.3945601 |
0.5584000 |
0.4148869 |
0.8095238 |
0.4571429 |
NA |
0.2192707 |
0.1629164 |
NA |
| 0.4443376 |
0.4913282 |
0.3144119 |
0.2666667 |
0.5142857 |
NA |
0.2172735 |
0.1390382 |
NA |
| 0.5018282 |
0.5994096 |
0.4198753 |
0.4761905 |
0.6857143 |
NA |
0.2993648 |
0.2096995 |
NA |
| 0.4492566 |
0.4663133 |
0.2831953 |
0.2761905 |
0.4571429 |
NA |
0.2084944 |
0.1266201 |
NA |
| 0.4769408 |
0.5185750 |
0.3318983 |
0.6571429 |
0.6285714 |
NA |
0.2461490 |
0.1575403 |
NA |
#Sample.SD = Standard deviation of the item
#Item.total = Correlation of the item with the total test score
#Item.Tot.woi = Correlation of item with total test score (scored without item)
#Difficulty = Mean of the item (p)
#Discrimination = Discrimination of the item (u-l)/n
#Item.Criterion = Correlation of the item with the Criterion (y)
#Item.Reliab = Item reliability index
#Item.Rel.woi = Item reliability index (scored without item)
#Item.Validity = Item validity index
#Spearman.brown
#If n.or.r="n", the function will return a new reliability; if n.or.r="r", the function
#will return the factor by which the test length must change to achieve a desired level of reliability.
#(input) The new test length or a desired level of reliability, depending on n.or.r
cronbach.alpha(responses, standardized = TRUE)
##
## Standardized Cronbach's alpha for the 'responses' data-set
##
## Items: 9
## Sample units: 105
## alpha: 0.678
H <- cronbach.alpha(responses, standardized = TRUE)
spearman.brown(H$alpha, input = 0.95, n.or.r = "r")
## $n.new
## [1] 9.008542
rel95 <- spearman.brown(H$alpha, input = 0.95, n.or.r = "r") #returns number of items neede to attain 0.95 reliability
rel.Factor <- rel95$n.new*11 #11=number of items in the model
round(rel.Factor)
## [1] 99
#Splithalf coeff.
splitHalf(responses, raw = FALSE, brute = FALSE, n.sample = 1000, covar = FALSE, check.keys = FALSE, key = NULL, use = "complete") # at use = complete/pairwise
## Split half reliabilities
## Call: splitHalf(r = responses, raw = FALSE, brute = FALSE, n.sample = 1000,
## covar = FALSE, check.keys = FALSE, key = NULL, use = "complete")
##
## Maximum split half reliability (lambda 4) = 0.76
## Guttman lambda 6 = 0.67
## Average split half reliability = 0.66
## Guttman lambda 3 (alpha) = 0.68
## Minimum split half reliability (beta) = 0.6
#Striatums H-index
H #standardised alpha
##
## Standardized Cronbach's alpha for the 'responses' data-set
##
## Items: 9
## Sample units: 105
## alpha: 0.678
G <- sqrt(H$alpha/(1-H$alpha))
G
## [1] 1.452277
Striatums = (4*G+1)/3
Striatums
## [1] 2.269703
round(Striatums)
## [1] 2
#Split items
SPLIT.ITEMS <-
function(X, seed = NULL){
# optional fixed seed
if (!is.null(seed)) {set.seed(seed)}
X <- as.matrix(responses)
# if n = 2x, then lengths Y1 = Y2
# if n = 2x+1, then lenths Y1 = Y2+1
n <- ncol(responses)
index <- sample(1:n, ceiling(n/2))
Y1 <- X[, index ]
Y2 <- X[, -index]
return(list(Y1, Y2))
}
dump("SPLIT.ITEMS", file = "SPLIT.ITEMS.R")
#Computes a bootstrapped (1,000 replicates) Pearson product-moment correlation coefficient between two random subsets of examinees.
pearson <-
function(X, seed = NULL, n = NULL){
source("SPLIT.ITEMS.R")
# optional fixed seed
if (!is.null(seed)) {set.seed(seed)}
# the number of bootstrap replicates
if (is.null(n)) {n <- 1e3}
X <- as.matrix(X)
r <- rep(NA, n)
for (i in 1:n) {
# split items
Y <- SPLIT.ITEMS(X)
# total scores
S1 <- as.matrix(rowSums(Y[[1]]))
S2 <- as.matrix(rowSums(Y[[2]]))
# residual scores
R1 <- S1 - mean(S1)
R2 <- S2 - mean(S2)
# Pearson product-moment correlation coefficient
r[i] <- (t(R1) %*% R2) / (sqrt((t(R1) %*% R1)) * sqrt((t(R2) %*% R2)))
}
return(mean(r))
}
dump("pearson", file = "pearson.R")
#Compute the Pearson product-moment correlation coefficient
pearson(responses, seed = 456, n = 1)
## [1] 0.5124339
# compare
# split items
set.seed(456)
Y <- SPLIT.ITEMS(responses)
# Total scores
S1 <- as.matrix(rowSums(Y[[1]]))
S2 <- as.matrix(rowSums(Y[[2]]))
cor(S1, S2)
## [,1]
## [1,] 0.5124339
#Standard Error of measurement
SEM <-
function(X){
source("cronbachs.alpha.R")
X <- data.matrix(responses)
return(sd(rowSums(responses)) * sqrt(1 - cronbachs.alpha(responses)[[1]]))
}
SEM(responses)
## [1] 1.244043
# 90% confidence interval for the true score
knitr::kable(head(cbind(lower_bound = round(rowSums(responses)-1.65* sd(rowSums(responses))*
sqrt(1-KR20(responses)), 2), observed = rowSums(responses),
upper_bound = round(rowSums(responses)+1.65* sd(rowSums(responses))*
sqrt(1-KR20(responses)), 2)), 105),
align = "c",
caption = "Item Scores")
Item Scores
| -1.04 |
1 |
3.04 |
| 5.96 |
8 |
10.04 |
| 3.96 |
6 |
8.04 |
| 4.96 |
7 |
9.04 |
| 1.96 |
4 |
6.04 |
| 4.96 |
7 |
9.04 |
| 3.96 |
6 |
8.04 |
| 2.96 |
5 |
7.04 |
| -0.04 |
2 |
4.04 |
| 0.96 |
3 |
5.04 |
| 3.96 |
6 |
8.04 |
| 1.96 |
4 |
6.04 |
| 3.96 |
6 |
8.04 |
| 2.96 |
5 |
7.04 |
| 2.96 |
5 |
7.04 |
| 1.96 |
4 |
6.04 |
| 3.96 |
6 |
8.04 |
| 0.96 |
3 |
5.04 |
| 2.96 |
5 |
7.04 |
| 2.96 |
5 |
7.04 |
| 5.96 |
8 |
10.04 |
| 4.96 |
7 |
9.04 |
| 4.96 |
7 |
9.04 |
| 5.96 |
8 |
10.04 |
| 3.96 |
6 |
8.04 |
| 3.96 |
6 |
8.04 |
| 4.96 |
7 |
9.04 |
| 3.96 |
6 |
8.04 |
| -0.04 |
2 |
4.04 |
| 6.96 |
9 |
11.04 |
| 4.96 |
7 |
9.04 |
| -0.04 |
2 |
4.04 |
| 4.96 |
7 |
9.04 |
| 1.96 |
4 |
6.04 |
| -1.04 |
1 |
3.04 |
| 3.96 |
6 |
8.04 |
| 2.96 |
5 |
7.04 |
| 5.96 |
8 |
10.04 |
| 1.96 |
4 |
6.04 |
| 1.96 |
4 |
6.04 |
| 2.96 |
5 |
7.04 |
| -1.04 |
1 |
3.04 |
| 1.96 |
4 |
6.04 |
| 1.96 |
4 |
6.04 |
| 1.96 |
4 |
6.04 |
| -0.04 |
2 |
4.04 |
| 4.96 |
7 |
9.04 |
| 5.96 |
8 |
10.04 |
| -1.04 |
1 |
3.04 |
| 1.96 |
4 |
6.04 |
| 5.96 |
8 |
10.04 |
| 2.96 |
5 |
7.04 |
| 4.96 |
7 |
9.04 |
| -0.04 |
2 |
4.04 |
| -1.04 |
1 |
3.04 |
| -0.04 |
2 |
4.04 |
| 2.96 |
5 |
7.04 |
| 2.96 |
5 |
7.04 |
| 2.96 |
5 |
7.04 |
| -2.04 |
0 |
2.04 |
| 1.96 |
4 |
6.04 |
| 0.96 |
3 |
5.04 |
| 3.96 |
6 |
8.04 |
| -0.04 |
2 |
4.04 |
| 5.96 |
8 |
10.04 |
| 4.96 |
7 |
9.04 |
| 0.96 |
3 |
5.04 |
| 4.96 |
7 |
9.04 |
| 1.96 |
4 |
6.04 |
| 4.96 |
7 |
9.04 |
| 3.96 |
6 |
8.04 |
| 3.96 |
6 |
8.04 |
| 4.96 |
7 |
9.04 |
| 4.96 |
7 |
9.04 |
| 5.96 |
8 |
10.04 |
| 2.96 |
5 |
7.04 |
| 3.96 |
6 |
8.04 |
| 4.96 |
7 |
9.04 |
| 5.96 |
8 |
10.04 |
| -1.04 |
1 |
3.04 |
| 5.96 |
8 |
10.04 |
| 4.96 |
7 |
9.04 |
| 0.96 |
3 |
5.04 |
| 5.96 |
8 |
10.04 |
| 0.96 |
3 |
5.04 |
| 3.96 |
6 |
8.04 |
| 1.96 |
4 |
6.04 |
| 4.96 |
7 |
9.04 |
| 0.96 |
3 |
5.04 |
| 5.96 |
8 |
10.04 |
| 4.96 |
7 |
9.04 |
| -2.04 |
0 |
2.04 |
| 1.96 |
4 |
6.04 |
| 3.96 |
6 |
8.04 |
| 4.96 |
7 |
9.04 |
| 1.96 |
4 |
6.04 |
| 1.96 |
4 |
6.04 |
| 0.96 |
3 |
5.04 |
| 0.96 |
3 |
5.04 |
| -0.04 |
2 |
4.04 |
| 5.96 |
8 |
10.04 |
| 0.96 |
3 |
5.04 |
| 1.96 |
4 |
6.04 |
| 2.96 |
5 |
7.04 |
| 2.96 |
5 |
7.04 |
#Item analysis
item.analysis <-
function(responses){
# CRITICAL VALUES
cvpb = 0.20
cvdl = 0.15
cvdu = 0.85
require(CTT, warn.conflicts = FALSE, quietly = TRUE)
(ctt.analysis <- CTT::reliability(responses, itemal = TRUE, NA.Delete = TRUE))
# Mark items that are potentially problematic
item.analysis <- data.frame(item = seq(1:ctt.analysis$nItem),
r.pbis = ctt.analysis$pBis,
bis = ctt.analysis$bis,
item.mean = ctt.analysis$itemMean,
alpha.del = ctt.analysis$alphaIfDeleted)
# code provided by Dr. Gordon Brooks
if (TRUE) {
item.analysis$check <-
ifelse(item.analysis$r.pbis < cvpb |
item.analysis$item.mean < cvdl |
item.analysis$item.mean > cvdu, "???", "")
}
return(item.analysis)
}
dump("item.analysis", file = "item.analysis.R")
knitr::kable(item.analysis(responses),
align = "c",
caption = "Item Analysis")
Item Analysis
| 1 |
0.2923590 |
0.4190647 |
0.8476190 |
0.6627315 |
|
| 2 |
0.4260097 |
0.5159341 |
0.6000000 |
0.6341815 |
|
| 3 |
0.3352222 |
0.4047167 |
0.5428571 |
0.6551620 |
|
| 4 |
0.3552830 |
0.4324898 |
0.5333333 |
0.6505746 |
|
| 5 |
0.4148869 |
0.5615773 |
0.8095238 |
0.6405185 |
|
| 6 |
0.3144119 |
0.4426467 |
0.2666667 |
0.6586211 |
|
| 7 |
0.4198753 |
0.5152756 |
0.4761905 |
0.6354503 |
|
| 8 |
0.2831953 |
0.4038506 |
0.2761905 |
0.6649574 |
|
| 9 |
0.3318983 |
0.4091737 |
0.6571429 |
0.6554052 |
|
#Item Difficulty
item.difficulty <-
function(responses){
# CRITICAL VALUES
cvpb = 0.20
cvdl = 0.15
cvdu = 0.85
require(CTT, warn.conflicts = FALSE, quietly = TRUE)
ctt.analysis <- CTT::reliability(responses, itemal = TRUE, NA.Delete = TRUE)
test_difficulty <- data.frame(item = 1:ctt.analysis$nItem ,
difficulty = ctt.analysis$itemMean)
plot(test_difficulty,
main = "Test Item Difficulty",
type = "p",
pch = 1,
cex = 2.8,
col = "purple",
ylab = "Item Mean (Difficulty)",
xlab = "Item Number",
ylim = c(0, 1),
xlim = c(0, ctt.analysis$nItem))
abline(h = cvdl, col = "tomato")
abline(h = cvdu, col = "tomato")
abline(h = .3, col = "dodgerblue")
abline(h = .7, col = "dodgerblue")
text(diff(range(test_difficulty[, 1]))/2, 0.7,
"maximum information range",
col = "dodger blue",
pos = 3)
text(diff(range(test_difficulty[, 1]))/2, cvdu,
"rule of thumb acceptable range",
col = "tomato",
pos = 3)
outlier2 <- data.matrix(subset(cbind(test_difficulty[, 1],
test_difficulty[, 2]),
subset = ((test_difficulty[, 2] > cvdl &
test_difficulty[, 2] < .3) |
(test_difficulty[, 2] < cvdu &
test_difficulty[, 2] > .7))))
text(outlier2, paste("i", outlier2[,1], sep = ""),
col = "dodgerblue",
cex = .7)
return(test_difficulty[order(test_difficulty$difficulty),])
}
dump("item.difficulty", file = "item.difficulty.R")
knitr::kable(item.difficulty(responses),
align = "c",
caption = "Difficulty")
Difficulty
| 6 |
6 |
0.2666667 |
| 8 |
8 |
0.2761905 |
| 7 |
7 |
0.4761905 |
| 4 |
4 |
0.5333333 |
| 3 |
3 |
0.5428571 |
| 2 |
2 |
0.6000000 |
| 9 |
9 |
0.6571429 |
| 5 |
5 |
0.8095238 |
| 1 |
1 |
0.8476190 |
#Item discrimination
item.discrimination <-
function(responses){
# CRITICAL VALUES
cvpb = 0.20
cvdl = 0.15
cvdu = 0.85
require(CTT, warn.conflicts = FALSE, quietly = TRUE)
ctt.analysis <- CTT::reliability(responses, itemal = TRUE, NA.Delete = TRUE)
item.discrimination <- data.frame(item = 1:ctt.analysis$nItem ,
discrimination = ctt.analysis$pBis)
plot(item.discrimination,
type = "p",
pch = 1,
cex = 3,
col = "purple",
ylab = "Item-Total Correlation",
xlab = "Item Number",
ylim = c(0, 1),
main = "Test Item Discriminations")
abline(h = cvpb, col = "red")
return(item.discrimination[order(item.discrimination$discrimination),])
}
dump("item.discrimination", file = "item.discrimination.R")
knitr::kable(item.discrimination(responses),
align = "c",
caption = "Item Discrimination")
Item Discrimination
| 8 |
8 |
0.2831953 |
| 1 |
1 |
0.2923590 |
| 6 |
6 |
0.3144119 |
| 9 |
9 |
0.3318983 |
| 3 |
3 |
0.3352222 |
| 4 |
4 |
0.3552830 |
| 5 |
5 |
0.4148869 |
| 7 |
7 |
0.4198753 |
| 2 |
2 |
0.4260097 |
#Item total correlation
test_item.total <-
function(responses){
# CRITICAL VALUES
cvpb = 0.20
cvdl = 0.15
cvdu = 0.85
require(CTT, warn.conflicts = FALSE, quietly = TRUE)
ctt.analysis <- CTT::reliability(responses, itemal = TRUE, NA.Delete = TRUE)
test_item.total <- data.frame(item = 1:ctt.analysis$nItem ,
biserial = ctt.analysis$bis)
plot(test_item.total,
main = "Test Item-Total Correlation",
type = "p",
pch = 1,
cex = 2.8,
col = "purple",
ylab = "Item-Total Correlation",
xlab = "Item Number",
ylim = c(0, 1),
xlim = c(0, ctt.analysis$nItem))
abline(h = cvpb, col = "red")
return(test_item.total[order(test_item.total$biserial),])
}
dump("test_item.total", file = "test_item.total.R")
knitr::kable(test_item.total(responses),
align = "c",
caption = "Item Total Correlation")
Item Total Correlation
| 8 |
8 |
0.4038506 |
| 3 |
3 |
0.4047167 |
| 9 |
9 |
0.4091737 |
| 1 |
1 |
0.4190647 |
| 4 |
4 |
0.4324898 |
| 6 |
6 |
0.4426467 |
| 7 |
7 |
0.5152756 |
| 2 |
2 |
0.5159341 |
| 5 |
5 |
0.5615773 |
#IRT
#2PL model
IRTmodel = ltm(responses~z1, IRT.param = TRUE)
summary(IRTmodel)
##
## Call:
## ltm(formula = responses ~ z1, IRT.param = TRUE)
##
## Model Summary:
## log.Lik AIC BIC
## -533.5194 1103.039 1150.81
##
## Coefficients:
## value std.err z.vals
## Dffclt.Item 1 -1.7915 0.5404 -3.3154
## Dffclt.Item 2 -0.3948 0.2123 -1.8595
## Dffclt.Item 3 -0.2161 0.2595 -0.8330
## Dffclt.Item 4 -0.1514 0.2340 -0.6472
## Dffclt.Item 5 -1.2168 0.2894 -4.2048
## Dffclt.Item 6 1.0742 0.3409 3.1512
## Dffclt.Item 7 0.1019 0.1911 0.5329
## Dffclt.Item 8 1.2009 0.4476 2.6829
## Dffclt.Item 9 -0.7721 0.3111 -2.4819
## Dscrmn.Item 1 1.2005 0.4781 2.5109
## Dscrmn.Item 2 1.3671 0.4505 3.0349
## Dscrmn.Item 3 0.9182 0.3337 2.7515
## Dscrmn.Item 4 1.0404 0.3590 2.8980
## Dscrmn.Item 5 1.8007 0.6514 2.7645
## Dscrmn.Item 6 1.2070 0.4501 2.6818
## Dscrmn.Item 7 1.4086 0.4716 2.9870
## Dscrmn.Item 8 0.9520 0.3936 2.4187
## Dscrmn.Item 9 1.0155 0.3561 2.8517
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Convergence: 0
## max(|grad|): <1e-06
## quasi-Newton: BFGS
knitr::kable(coef(IRTmodel),
align = "c",
caption = "")
| Item 1 |
-1.7915281 |
1.2004682 |
| Item 2 |
-0.3948019 |
1.3671194 |
| Item 3 |
-0.2161422 |
0.9182071 |
| Item 4 |
-0.1514169 |
1.0403556 |
| Item 5 |
-1.2167542 |
1.8007489 |
| Item 6 |
1.0741981 |
1.2070130 |
| Item 7 |
0.1018545 |
1.4085975 |
| Item 8 |
1.2008945 |
0.9520480 |
| Item 9 |
-0.7720712 |
1.0154738 |
plot(IRTmodel, type = "ICC") #add argument items = ... for plot one item at a time
plot(IRTmodel, type = "IIC", items = 0) #test information function
person.fit(IRTmodel)
##
## Person-Fit Statistics and P-values
##
## Call:
## ltm(formula = responses ~ z1, IRT.param = TRUE)
##
## Alternative: Inconsistent response pattern under the estimated model
##
## It 1 It 2 It 3 It 4 It 5 It 6 It 7 It 8 It 9 L0 Lz Pr(<Lz)
## 1 0 0 0 0 0 0 0 0 0 -1.8070 1.0561 0.8545
## 2 0 0 0 0 0 0 0 0 1 -3.3318 0.5405 0.7056
## 3 0 0 0 0 0 1 0 0 1 -6.7483 -1.4500 0.0735
## 4 0 0 0 0 1 0 0 0 0 -3.3725 0.7918 0.7858
## 5 0 0 0 0 1 0 0 0 1 -4.0918 0.5697 0.7155
## 6 0 0 0 0 1 0 1 1 1 -7.4424 -1.7916 0.0366
## 7 0 0 0 1 0 0 1 0 0 -6.5543 -1.2796 0.1003
## 8 0 0 0 1 1 0 1 0 0 -6.3465 -0.9414 0.1733
## 9 0 0 1 0 0 0 0 0 0 -3.7189 0.2574 0.6016
## 10 0 0 1 0 1 0 0 0 1 -4.9273 0.1070 0.5426
## 11 0 0 1 1 1 0 1 1 0 -8.1495 -2.3949 0.0083
## 12 0 1 0 0 1 0 0 1 1 -6.7396 -1.2347 0.1085
## 13 0 1 1 1 0 0 0 0 0 -6.7817 -1.3419 0.0898
## 14 0 1 1 1 0 1 1 0 1 -9.1464 -3.2272 0.0006
## 15 0 1 1 1 1 0 1 1 1 -6.0315 -1.0749 0.1412
## 16 1 0 0 0 0 0 0 0 0 -2.3305 1.2616 0.8964
## 17 1 0 0 0 0 0 0 0 1 -3.2912 0.9823 0.837
## 18 1 0 0 0 1 0 0 0 0 -2.9576 1.4478 0.9262
## 19 1 0 0 0 1 0 0 0 1 -3.2335 1.4449 0.9258
## 20 1 0 0 0 1 0 0 1 1 -5.0443 0.1051 0.5419
## 21 1 0 0 0 1 0 1 0 0 -4.4747 0.5259 0.7005
## 22 1 0 0 0 1 0 1 0 1 -4.2526 0.7319 0.7679
## 23 1 0 0 0 1 1 0 0 0 -5.4736 -0.2665 0.3949
## 24 1 0 0 1 1 0 0 0 0 -3.8814 0.9486 0.8286
## 25 1 0 0 1 1 0 0 0 1 -3.7894 1.0982 0.8639
## 26 1 0 0 1 1 0 1 0 0 -4.8883 0.2269 0.5898
## 27 1 0 0 1 1 0 1 0 1 -4.2918 0.6230 0.7334
## 28 1 0 0 1 1 0 1 1 1 -5.2625 -0.2854 0.3877
## 29 1 0 1 0 0 0 0 0 0 -3.7277 0.6466 0.741
## 30 1 0 1 0 1 0 0 0 1 -3.6842 1.1775 0.8805
## 31 1 0 1 0 1 0 1 0 1 -4.2482 0.6725 0.7494
## 32 1 0 1 0 1 1 0 0 1 -5.3848 -0.2063 0.4183
## 33 1 0 1 1 0 0 0 0 1 -5.2009 -0.0577 0.477
## 34 1 0 1 1 1 0 0 0 0 -4.3242 0.6733 0.7496
## 35 1 0 1 1 1 0 0 0 1 -3.9059 0.9802 0.8365
## 36 1 0 1 1 1 0 0 1 1 -5.0530 -0.0399 0.4841
## 37 1 0 1 1 1 0 1 1 0 -5.8818 -0.7425 0.2289
## 38 1 0 1 1 1 0 1 1 1 -4.5525 0.0339 0.5135
## 39 1 0 1 1 1 1 1 0 1 -4.4989 0.0012 0.5005
## 40 1 1 0 0 0 0 0 0 0 -4.0444 0.5445 0.7069
## 41 1 1 0 0 0 0 0 0 1 -4.4765 0.4290 0.666
## 42 1 1 0 0 0 0 0 1 1 -6.4306 -0.9955 0.1597
## 43 1 1 0 0 0 0 1 1 0 -7.6011 -1.9126 0.0279
## 44 1 1 0 0 1 0 0 0 0 -3.7609 1.0787 0.8596
## 45 1 1 0 0 1 0 0 0 1 -3.5535 1.2872 0.901
## 46 1 1 0 0 1 1 0 0 1 -5.0598 0.0026 0.501
## 47 1 1 0 0 1 1 1 1 1 -5.1143 -0.5352 0.2962
## 48 1 1 0 1 1 0 0 0 0 -4.1895 0.7822 0.783
## 49 1 1 0 1 1 0 0 0 1 -3.6083 1.1674 0.8785
## 50 1 1 0 1 1 0 1 0 0 -4.5128 0.3986 0.6549
## 51 1 1 0 1 1 0 1 1 1 -3.8269 0.4188 0.6623
## 52 1 1 0 1 1 1 1 0 1 -3.7197 0.4082 0.6584
## 53 1 1 0 1 1 1 1 1 1 -3.5934 0.1468 0.5583
## 54 1 1 1 0 1 0 0 0 0 -4.0994 0.8550 0.8037
## 55 1 1 1 0 1 0 1 0 1 -3.4136 1.0487 0.8528
## 56 1 1 1 0 1 0 1 1 1 -3.8998 0.4066 0.6579
## 57 1 1 1 0 1 1 0 0 0 -5.6481 -0.4488 0.3268
## 58 1 1 1 0 1 1 0 0 1 -4.6564 0.1581 0.5628
## 59 1 1 1 0 1 1 1 0 0 -5.3935 -0.4803 0.3155
## 60 1 1 1 0 1 1 1 1 1 -3.7434 0.0945 0.5376
## 61 1 1 1 1 0 0 0 0 1 -5.4457 -0.2183 0.4136
## 62 1 1 1 1 0 0 1 1 0 -7.6386 -2.0486 0.0202
## 63 1 1 1 1 1 0 0 0 0 -4.1907 0.7199 0.7642
## 64 1 1 1 1 1 0 0 0 1 -3.2640 1.2590 0.896
## 65 1 1 1 1 1 0 0 1 1 -3.9002 0.5199 0.6984
## 66 1 1 1 1 1 0 1 0 0 -4.0279 0.5802 0.7191
## 67 1 1 1 1 1 0 1 0 1 -2.5174 1.3560 0.9125
## 68 1 1 1 1 1 0 1 1 1 -2.5358 0.9591 0.8313
## 69 1 1 1 1 1 1 0 0 1 -3.8513 0.4770 0.6833
## 70 1 1 1 1 1 1 0 1 0 -5.8317 -0.9276 0.1768
## 71 1 1 1 1 1 1 0 1 1 -3.9665 0.0791 0.5315
## 72 1 1 1 1 1 1 1 0 1 -2.3054 1.0028 0.842
## 73 1 1 1 1 1 1 1 1 0 -4.3317 -0.3014 0.3816
## 74 1 1 1 1 1 1 1 1 1 -1.6581 0.9874 0.8383
item.fit(IRTmodel) #we want the chi^2 to be NON-SIGNIFICANT! But depends on sample size
##
## Item-Fit Statistics and P-values
##
## Call:
## ltm(formula = responses ~ z1, IRT.param = TRUE)
##
## Alternative: Items do not fit the model
## Ability Categories: 10
##
## X^2 Pr(>X^2)
## It 1 7.1329 0.5224
## It 2 7.5840 0.4751
## It 3 13.1364 0.1072
## It 4 11.4881 0.1755
## It 5 10.3397 0.242
## It 6 12.8634 0.1166
## It 7 8.1821 0.4159
## It 8 19.3725 0.013
## It 9 8.6163 0.3757
#3PL model
IRTmodel2 = tpm(data = ctt.data[,c( -1, -2, -5, -6, -7, -9, -10, -12, -14, -16, -18, -19, -20, -21, -23, -24)],
type = "latent.trait", IRT.param = TRUE)
summary(IRTmodel2)
##
## Call:
## tpm(data = ctt.data[, c(-1, -2, -5, -6, -7, -9, -10, -12, -14,
## -16, -18, -19, -20, -21, -23, -24)], type = "latent.trait",
## IRT.param = TRUE)
##
## Model Summary:
## log.Lik AIC BIC
## -533.3024 1120.605 1192.262
##
## Coefficients:
## value std.err z.vals
## Gussng.X3 0.0328 0.6878 0.0478
## Gussng.X4 0.0000 0.0055 0.0077
## Gussng.X8 0.0000 0.0043 0.0053
## Gussng.X11 0.0000 0.0008 0.0022
## Gussng.X13 0.0011 0.0453 0.0235
## Gussng.X15 0.0578 0.0760 0.7610
## Gussng.X17 0.0001 0.0122 0.0085
## Gussng.X22 0.0000 0.0000 0.0002
## Gussng.X25 0.2172 0.4160 0.5221
## Dffclt.X3 -1.7041 1.2134 -1.4044
## Dffclt.X4 -0.3898 0.2113 -1.8446
## Dffclt.X8 -0.2176 0.2649 -0.8213
## Dffclt.X11 -0.1480 0.2322 -0.6373
## Dffclt.X13 -1.2154 0.2963 -4.1020
## Dffclt.X15 1.0664 0.3056 3.4890
## Dffclt.X17 0.1045 0.1923 0.5433
## Dffclt.X22 1.2400 0.4742 2.6151
## Dffclt.X25 -0.2457 1.0920 -0.2250
## Dscrmn.X3 1.2544 0.5967 2.1025
## Dscrmn.X4 1.3823 0.4508 3.0660
## Dscrmn.X8 0.8956 0.3290 2.7225
## Dscrmn.X11 1.0525 0.3587 2.9343
## Dscrmn.X13 1.7991 0.6528 2.7560
## Dscrmn.X15 1.8071 1.2621 1.4318
## Dscrmn.X17 1.4115 0.4687 3.0114
## Dscrmn.X22 0.9118 0.3842 2.3733
## Dscrmn.X25 1.2954 1.0192 1.2710
##
## Integration:
## method: Gauss-Hermite
## quadrature points: 21
##
## Optimization:
## Optimizer: optim (BFGS)
## Convergence: 0
## max(|grad|): 0.022
knitr::kable(coef(IRTmodel2),
align = "c",
caption = "")
| X3 |
0.0328442 |
-1.7041465 |
1.2544473 |
| X4 |
0.0000426 |
-0.3898091 |
1.3822953 |
| X8 |
0.0000229 |
-0.2175817 |
0.8956320 |
| X11 |
0.0000017 |
-0.1480015 |
1.0524593 |
| X13 |
0.0010641 |
-1.2154204 |
1.7990790 |
| X15 |
0.0578283 |
1.0663876 |
1.8070746 |
| X17 |
0.0001037 |
0.1044896 |
1.4115488 |
| X22 |
0.0000000 |
1.2400445 |
0.9118309 |
| X25 |
0.2171827 |
-0.2457209 |
1.2954246 |
plot(IRTmodel2, type = "ICC")
plot(IRTmodel2, type = "IIC", items = 0) #test information function
person.fit(IRTmodel2) #check for outliers in the model
##
## Person-Fit Statistics and P-values
##
## Call:
## tpm(data = ctt.data[, c(-1, -2, -5, -6, -7, -9, -10, -12, -14,
## -16, -18, -19, -20, -21, -23, -24)], type = "latent.trait",
## IRT.param = TRUE)
##
## Alternative: Inconsistent response pattern under the estimated model
##
## X3 X4 X8 X11 X13 X15 X17 X22 X25 L0 Lz Pr(<Lz)
## 1 0 0 0 0 0 0 0 0 0 -2.0000 1.1014 0.8646
## 2 0 0 0 0 0 0 0 0 1 -3.0969 0.6476 0.7414
## 3 0 0 0 0 0 1 0 0 1 -5.8964 -1.0452 0.148
## 4 0 0 0 0 1 0 0 0 0 -3.4351 0.8314 0.7971
## 5 0 0 0 0 1 0 0 0 1 -4.0863 0.5716 0.7162
## 6 0 0 0 0 1 0 1 1 1 -7.4368 -1.8092 0.0352
## 7 0 0 0 1 0 0 1 0 0 -6.6053 -1.2748 0.1012
## 8 0 0 0 1 1 0 1 0 0 -6.3078 -0.9236 0.1778
## 9 0 0 1 0 0 0 0 0 0 -3.8000 0.3142 0.6233
## 10 0 0 1 0 1 0 0 0 1 -4.9344 0.0926 0.5369
## 11 0 0 1 1 1 0 1 1 0 -8.0565 -2.3294 0.0099
## 12 0 1 0 0 1 0 0 1 1 -6.7371 -1.2565 0.1045
## 13 0 1 1 1 0 0 0 0 0 -6.7733 -1.3253 0.0925
## 14 0 1 1 1 0 1 1 0 1 -9.3438 -3.3765 0.0004
## 15 0 1 1 1 1 0 1 1 1 -6.0879 -1.0974 0.1362
## 16 1 0 0 0 0 0 0 0 0 -2.4781 1.2883 0.9012
## 17 1 0 0 0 0 0 0 0 1 -3.2814 0.9883 0.8385
## 18 1 0 0 0 1 0 0 0 0 -2.9810 1.4762 0.93
## 19 1 0 0 0 1 0 0 0 1 -3.2772 1.4073 0.9203
## 20 1 0 0 0 1 0 0 1 1 -5.0345 0.0874 0.5348
## 21 1 0 0 0 1 0 1 0 0 -4.4277 0.5509 0.7092
## 22 1 0 0 0 1 0 1 0 1 -4.2656 0.6847 0.7532
## 23 1 0 0 0 1 1 0 0 0 -5.4812 -0.3100 0.3783
## 24 1 0 0 1 1 0 0 0 0 -3.8578 0.9694 0.8338
## 25 1 0 0 1 1 0 0 0 1 -3.8175 1.0469 0.8524
## 26 1 0 0 1 1 0 1 0 0 -4.8091 0.2558 0.6009
## 27 1 0 0 1 1 0 1 0 1 -4.2687 0.5883 0.7218
## 28 1 0 0 1 1 0 1 1 1 -5.2161 -0.2759 0.3913
## 29 1 0 1 0 0 0 0 0 0 -3.7787 0.6822 0.7524
## 30 1 0 1 0 1 0 0 0 1 -3.7004 1.1384 0.8725
## 31 1 0 1 0 1 0 1 0 1 -4.2310 0.6377 0.7382
## 32 1 0 1 0 1 1 0 0 1 -5.6105 -0.4243 0.3357
## 33 1 0 1 1 0 0 0 0 1 -5.2382 -0.0971 0.4613
## 34 1 0 1 1 1 0 0 0 0 -4.2494 0.7061 0.7599
## 35 1 0 1 1 1 0 0 0 1 -3.8982 0.9376 0.8258
## 36 1 0 1 1 1 0 0 1 1 -5.0144 -0.0430 0.4829
## 37 1 0 1 1 1 0 1 1 0 -5.8102 -0.6869 0.2461
## 38 1 0 1 1 1 0 1 1 1 -4.5463 0.0429 0.5171
## 39 1 0 1 1 1 1 1 0 1 -4.4516 -0.0691 0.4725
## 40 1 1 0 0 0 0 0 0 0 -4.1158 0.5527 0.7098
## 41 1 1 0 0 0 0 0 0 1 -4.5365 0.3869 0.6506
## 42 1 1 0 0 0 0 0 1 1 -6.4353 -1.0174 0.1545
## 43 1 1 0 0 0 0 1 1 0 -7.5057 -1.8595 0.0315
## 44 1 1 0 0 1 0 0 0 0 -3.7215 1.1012 0.8646
## 45 1 1 0 0 1 0 0 0 1 -3.5694 1.2357 0.8917
## 46 1 1 0 0 1 1 0 0 1 -5.2504 -0.2088 0.4173
## 47 1 1 0 0 1 1 1 1 1 -5.0229 -0.5624 0.2869
## 48 1 1 0 1 1 0 0 0 0 -4.1130 0.8065 0.79
## 49 1 1 0 1 1 0 0 0 1 -3.5835 1.1252 0.8697
## 50 1 1 0 1 1 0 1 0 0 -4.4189 0.4375 0.6691
## 51 1 1 0 1 1 0 1 1 1 -3.8183 0.4416 0.6706
## 52 1 1 0 1 1 1 1 0 1 -3.5705 0.3861 0.6503
## 53 1 1 0 1 1 1 1 1 1 -3.4024 0.1699 0.5674
## 54 1 1 1 0 1 0 0 0 0 -4.0155 0.8881 0.8128
## 55 1 1 1 0 1 0 1 0 1 -3.3747 1.0396 0.8507
## 56 1 1 1 0 1 0 1 1 1 -3.9128 0.4169 0.6616
## 57 1 1 1 0 1 1 0 0 0 -5.7863 -0.5916 0.277
## 58 1 1 1 0 1 1 0 0 1 -4.7498 -0.0039 0.4984
## 59 1 1 1 0 1 1 1 0 0 -5.4638 -0.5777 0.2817
## 60 1 1 1 0 1 1 1 1 1 -3.6155 0.0903 0.536
## 61 1 1 1 1 0 0 0 0 1 -5.4547 -0.2596 0.3976
## 62 1 1 1 1 0 0 1 1 0 -7.5414 -1.9732 0.0242
## 63 1 1 1 1 1 0 0 0 0 -4.0958 0.7538 0.7745
## 64 1 1 1 1 1 0 0 0 1 -3.2205 1.2380 0.8921
## 65 1 1 1 1 1 0 0 1 1 -3.8803 0.5328 0.7029
## 66 1 1 1 1 1 0 1 0 0 -3.9761 0.6190 0.732
## 67 1 1 1 1 1 0 1 0 1 -2.5070 1.3750 0.9154
## 68 1 1 1 1 1 0 1 1 1 -2.6727 0.9513 0.8293
## 69 1 1 1 1 1 1 0 0 1 -3.7944 0.3996 0.6553
## 70 1 1 1 1 1 1 0 1 0 -5.8847 -0.9969 0.1594
## 71 1 1 1 1 1 1 0 1 1 -3.8436 0.0673 0.5268
## 72 1 1 1 1 1 1 1 0 1 -2.1034 1.0128 0.8444
## 73 1 1 1 1 1 1 1 1 0 -4.4094 -0.3613 0.3589
## 74 1 1 1 1 1 1 1 1 1 -1.5078 0.9776 0.8359
item.fit(IRTmodel2) #we want the chi^2 to be NON-SIGNIFICANT! But depends on sample size
##
## Item-Fit Statistics and P-values
##
## Call:
## tpm(data = ctt.data[, c(-1, -2, -5, -6, -7, -9, -10, -12, -14,
## -16, -18, -19, -20, -21, -23, -24)], type = "latent.trait",
## IRT.param = TRUE)
##
## Alternative: Items do not fit the model
## Ability Categories: 10
##
## X^2 Pr(>X^2)
## X3 6.7900 0.4511
## X4 8.1744 0.3175
## X8 14.2901 0.0463
## X11 12.5491 0.0839
## X13 10.6302 0.1556
## X15 9.5006 0.2187
## X17 9.0508 0.249
## X22 19.5022 0.0068
## X25 5.7731 0.5665
#Compare models
#anova(IRTmodel, IRTmodel2) # p=1, but I get an error
#If the P-value is non-sig, then the models are the same, thus it is better to use the simpler (2PL model)
#Mirt 2PL
mirtmodel = mirt(data = ctt.data[,c( -1, -2, -5, -6, -7, -9, -10, -12, -14, -16, -18, -19, -20, -21, -23, -24)],
model = 1,
itemtype = "gpcm")
##
Iteration: 1, Log-Lik: -538.256, Max-Change: 0.39741
Iteration: 2, Log-Lik: -534.811, Max-Change: 0.23002
Iteration: 3, Log-Lik: -533.902, Max-Change: 0.13082
Iteration: 4, Log-Lik: -533.596, Max-Change: 0.06178
Iteration: 5, Log-Lik: -533.545, Max-Change: 0.03715
Iteration: 6, Log-Lik: -533.528, Max-Change: 0.02335
Iteration: 7, Log-Lik: -533.521, Max-Change: 0.00804
Iteration: 8, Log-Lik: -533.520, Max-Change: 0.00604
Iteration: 9, Log-Lik: -533.520, Max-Change: 0.00340
Iteration: 10, Log-Lik: -533.520, Max-Change: 0.00133
Iteration: 11, Log-Lik: -533.520, Max-Change: 0.00027
Iteration: 12, Log-Lik: -533.520, Max-Change: 0.00022
Iteration: 13, Log-Lik: -533.520, Max-Change: 0.00017
Iteration: 14, Log-Lik: -533.520, Max-Change: 0.00014
Iteration: 15, Log-Lik: -533.520, Max-Change: 0.00012
Iteration: 16, Log-Lik: -533.520, Max-Change: 0.00028
Iteration: 17, Log-Lik: -533.520, Max-Change: 0.00022
Iteration: 18, Log-Lik: -533.520, Max-Change: 0.00014
Iteration: 19, Log-Lik: -533.520, Max-Change: 0.00003
summary(mirtmodel)
## F1 h2
## X3 0.576 0.332
## X4 0.626 0.392
## X8 0.475 0.225
## X11 0.522 0.272
## X13 0.727 0.528
## X15 0.578 0.335
## X17 0.638 0.407
## X22 0.488 0.238
## X25 0.512 0.263
##
## SS loadings: 2.992
## Proportion Var: 0.332
##
## Factor correlations:
##
## F1
## F1 1
coef(mirtmodel, IRTpars = T) #coefficiants
## $X3
## a b1
## par 1.2 -1.792
##
## $X4
## a b1
## par 1.367 -0.395
##
## $X8
## a b1
## par 0.918 -0.216
##
## $X11
## a b1
## par 1.04 -0.151
##
## $X13
## a b1
## par 1.8 -1.217
##
## $X15
## a b1
## par 1.207 1.074
##
## $X17
## a b1
## par 1.409 0.102
##
## $X22
## a b1
## par 0.952 1.201
##
## $X25
## a b1
## par 1.015 -0.772
##
## $GroupPars
## MEAN_1 COV_11
## par 0 1
plot(mirtmodel, type = "trace") #curves for all items
plot(mirtmodel, type = "info") # test information curve
plot(mirtmodel) #expected total score
itemfit(mirtmodel) #item fit statistics
## item S_X2 df.S_X2 p.S_X2
## 1 X3 0.625 3 0.891
## 2 X4 1.145 4 0.887
## 3 X8 8.791 5 0.118
## 4 X11 3.872 5 0.568
## 5 X13 1.156 3 0.764
## 6 X15 1.190 3 0.755
## 7 X17 1.638 3 0.651
## 8 X22 6.973 4 0.137
## 9 X25 2.271 4 0.686