R Under development (unstable) (2024-08-29 r87078 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > set.seed(7485) > library(mstDIF) > > # function to simulate 3PL data > P3pl <- function(theta, a = rep(1, length(b)), b, c = rep(0, length(b)), + generateResponses = FALSE, seed = 1){ + nItem <- length(b) + nPerson <- length(theta) + A <- rep(1, nPerson) %o% a + B <- rep(1, nPerson) %o% b + C <- rep(1, nPerson) %o% c + Theta <- theta %o% rep(1, nItem) + + Theta_B <- Theta - B + P <- C + (1 - C) * (1 + exp(- A * Theta_B))^(-1) + if(generateResponses){ + set.seed(seed) + resp <- (P > matrix(runif(nPerson * nItem), + ncol = nItem, nrow = nPerson)) * 1 + return(resp) + } else (return(P)) + } > > # item parameters > b <- seq(-2, 1.8, length.out = 3) > a <- runif(length(b), 1, 1.5) > c <- rep(0, length(b)) > > # person parameters > nPerson <- 20 > theta <- rnorm(nPerson) > > # generate responses > resp <- P3pl(theta, a, b, generateResponses = TRUE) > > # create orders > metric <- rnorm(nPerson) > factor <- factor(sample(1:3, size = nPerson, replace = TRUE)) > order <- ordered(sample(1:3, size = nPerson, replace = TRUE)) > > # compute scores > scores <- mstDIF:::get_scores(resp, a, b, c, theta, + slope_intercept = FALSE, return_terms = TRUE) > scores $scores [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.18468225 -0.29352519 0.17191034 -0.08138442 0.7041845 0.16688286 [2,] 0.24976935 -0.13752642 0.10360703 -0.15718775 -0.6362792 0.06725971 [3,] 0.25808191 0.24282391 0.05379964 -0.27972165 0.2904038 0.02636932 [4,] -1.75302701 0.03643510 0.13375030 0.92545002 -0.5366830 0.10198830 [5,] 0.26231754 -0.62929507 0.05773734 -0.26498604 -0.8135188 0.02901723 [6,] 0.18569685 0.17124920 0.17098626 -0.08218499 -0.4199869 0.16481453 [7,] 0.08519532 0.22398270 -0.14471676 -0.02501704 -0.1519493 -0.81380608 [8,] 0.11580700 0.24759733 0.19010357 -0.03870471 -0.2272828 0.38960229 [9,] 0.21991317 -0.05828303 0.13709574 -0.11414900 0.5950554 0.10650713 [10,] 0.19486646 0.24290672 0.03447349 -0.37956746 0.1996025 0.01468697 [11,] 0.18122482 0.18219312 0.17500655 -0.07870552 -0.4068408 0.17411170 [12,] 0.22215194 -0.04513158 0.13474084 -0.11668733 0.5876285 0.10330967 [13,] 0.11464363 0.24729999 0.18873158 -0.03813587 -0.2243038 0.39537570 [14,] 0.25852876 -0.68865553 0.05414770 -0.27836213 -0.8292328 0.02659968 [15,] -0.86864569 0.22788034 0.06387389 0.79869143 0.3338023 0.03332827 [16,] 0.26407992 0.18398173 0.08156276 -0.19817396 0.4045920 0.04707138 [17,] 0.23633298 0.03423168 0.11938525 -0.13482229 0.5381899 0.08422634 [18,] -0.40513722 -1.03611565 0.03762819 0.68411508 -0.9059990 0.01644519 [19,] 0.22423371 -0.03306617 0.13253647 -0.11911593 0.5806450 0.10038748 [20,] 0.19966544 0.13066538 0.15769286 -0.09394865 -0.4619674 0.13855327 $terms $terms$P [,1] [,2] [,3] [1,] 0.9219875 0.6280759 0.11999384 [2,] 0.8493248 0.4324902 0.04836178 [3,] 0.7318676 0.2590168 0.01896034 [4,] 0.8871073 0.5213220 0.07333268 [5,] 0.7459927 0.2744067 0.02086427 [6,] 0.9212201 0.6254055 0.11850665 [7,] 0.9760195 0.8644734 0.41484873 [8,] 0.9628989 0.7972820 0.28013587 [9,] 0.8905804 0.5307415 0.07658186 [10,] 0.6361586 0.1780294 0.01056038 [11,] 0.9245554 0.6371307 0.12519160 [12,] 0.8881472 0.5241173 0.07428279 [13,] 0.9634442 0.7999391 0.28428713 [14,] 0.7331708 0.2603911 0.01912598 [15,] 0.7656005 0.2977248 0.02396404 [16,] 0.8100367 0.3608635 0.03384575 [17,] 0.8707636 0.4800221 0.06056130 [18,] 0.6557712 0.1919218 0.01182459 [19,] 0.8858192 0.5178886 0.07218165 [20,] 0.9099438 0.5879622 0.09962401 $terms$a [,1] [,2] [,3] [1,] 2.3673403 0.46734029 -1.4326597 [2,] 1.6576669 -0.24233311 -2.1423331 [3,] 0.9625167 -0.93748326 -2.8374833 [4,] 1.9761161 0.07611608 -1.8238839 [5,] 1.0327166 -0.86728342 -2.7672834 [6,] 2.3571589 0.45715890 -1.4428411 [7,] 3.5526849 1.65268487 -0.2473151 [8,] 3.1213880 1.22138797 -0.6786120 [9,] 2.0098143 0.10981435 -1.7901857 [10,] 0.5355807 -1.36441929 -3.2644193 [11,] 2.4020903 0.50209026 -1.3979097 [12,] 1.9861097 0.08610969 -1.8138903 [13,] 3.1361234 1.23612336 -0.6638766 [14,] 0.9688923 -0.93110767 -2.8311077 [15,] 1.1345939 -0.76540606 -2.6654061 [16,] 1.3901626 -0.50983742 -2.4098374 [17,] 1.8286873 -0.07131272 -1.9713127 [18,] 0.6178027 -1.28219731 -3.1821973 [19,] 1.9638480 0.06384804 -1.8361520 [20,] 2.2171199 0.31711990 -1.5828801 $terms$b [,1] [,2] [,3] [1,] -1.043222 -1.121177 -1.390762 [2,] -1.043222 -1.121177 -1.390762 [3,] -1.043222 -1.121177 -1.390762 [4,] -1.043222 -1.121177 -1.390762 [5,] -1.043222 -1.121177 -1.390762 [6,] -1.043222 -1.121177 -1.390762 [7,] -1.043222 -1.121177 -1.390762 [8,] -1.043222 -1.121177 -1.390762 [9,] -1.043222 -1.121177 -1.390762 [10,] -1.043222 -1.121177 -1.390762 [11,] -1.043222 -1.121177 -1.390762 [12,] -1.043222 -1.121177 -1.390762 [13,] -1.043222 -1.121177 -1.390762 [14,] -1.043222 -1.121177 -1.390762 [15,] -1.043222 -1.121177 -1.390762 [16,] -1.043222 -1.121177 -1.390762 [17,] -1.043222 -1.121177 -1.390762 [18,] -1.043222 -1.121177 -1.390762 [19,] -1.043222 -1.121177 -1.390762 [20,] -1.043222 -1.121177 -1.390762 $terms$c [,1] [,2] [,3] [1,] 1.084613 1.592164 8.333761 [2,] 1.177406 2.312191 20.677487 [3,] 1.366367 3.860753 52.741674 [4,] 1.127259 1.918200 13.636484 [5,] 1.340496 3.644226 47.928833 [6,] 1.085517 1.598963 8.438345 [7,] 1.024570 1.156774 2.410517 [8,] 1.038531 1.254261 3.569696 [9,] 1.122863 1.884156 13.057923 [10,] 1.571935 5.617050 94.693600 [11,] 1.081601 1.569537 7.987757 [12,] 1.125939 1.907970 13.462068 [13,] 1.037943 1.250095 3.517570 [14,] 1.363939 3.840377 52.284912 [15,] 1.306164 3.358807 41.729193 [16,] 1.234512 2.771131 29.545805 [17,] 1.148417 2.083237 16.512196 [18,] 1.524922 5.210455 84.569506 [19,] 1.128899 1.930917 13.853937 [20,] 1.098969 1.700790 10.037741 > > # compute process > process <- apply(mstDIF:::scale_scores(scores$scores, + decorrelate = TRUE)[order(metric),], + 2, cumsum) > process [,1] [,2] [,3] [,4] [,5] [1,] 2.420033e-04 -3.334786e-01 -1.006781e-01 -2.818289e-01 -2.725019e-01 [2,] 7.404527e-02 -2.365103e-01 -2.790570e-01 -3.512988e-01 -9.589193e-02 [3,] -1.210368e-01 -7.134302e-02 -5.975643e-01 1.935734e-01 9.981391e-02 [4,] 1.275834e-02 1.254864e-01 -3.028179e-01 2.360873e-01 -1.017536e-01 [5,] 9.917085e-02 3.384026e-01 -4.148226e-01 3.013782e-01 -1.746471e-01 [6,] 1.558282e-01 2.450313e-01 -2.763645e-01 2.776883e-01 1.168400e-01 [7,] 2.884038e-01 4.476681e-01 2.711118e-02 3.231950e-01 -8.003818e-02 [8,] 3.471061e-01 3.641271e-01 1.529694e-01 2.996919e-01 2.065794e-01 [9,] 4.640467e-01 3.709941e-01 2.077436e-01 2.405129e-01 -5.336269e-02 [10,] 5.661012e-01 6.017797e-01 1.599709e-01 3.079788e-01 -1.717531e-01 [11,] -3.424868e-01 6.527641e-01 3.862716e-01 3.117870e-01 -3.015950e-01 [12,] -2.733318e-01 6.244143e-01 4.354578e-01 2.866215e-01 -4.371310e-02 [13,] -2.539443e-01 7.702933e-01 9.482611e-02 1.029992e-01 7.482810e-02 [14,] -2.736185e-01 3.940135e-01 -1.394281e-02 -2.107656e-01 -1.963198e-01 [15,] -1.721636e-01 6.246434e-01 -8.422372e-02 -1.430374e-01 -3.131174e-01 [16,] -2.850898e-01 7.629785e-01 -5.225636e-01 -5.245715e-01 -2.305197e-01 [17,] -2.012700e-01 2.218021e-01 -7.402078e-01 3.678912e-02 -4.604837e-01 [18,] -6.493485e-02 3.961034e-01 -4.827298e-01 6.781724e-02 -6.761841e-01 [19,] -1.054789e-02 2.920380e-01 -3.309072e-01 4.379308e-02 -3.794727e-01 [20,] 4.059253e-16 -9.540979e-18 -1.249001e-16 2.524023e-16 4.163336e-17 [,6] [1,] -5.781780e-03 [2,] 1.578973e-02 [3,] 6.572980e-02 [4,] 8.266642e-02 [5,] -7.856663e-01 [6,] -7.919143e-01 [7,] -7.680210e-01 [8,] -7.738898e-01 [9,] -7.878639e-01 [10,] -4.507996e-01 [11,] -4.672783e-01 [12,] -4.691381e-01 [13,] -4.241853e-01 [14,] -4.296219e-01 [15,] -8.059560e-02 [16,] -1.974634e-02 [17,] -1.952819e-03 [18,] -1.154427e-03 [19,] -7.702091e-03 [20,] -1.578598e-16 > > # Compare with strucchange > library(strucchange) Loading required package: zoo Attaching package: 'zoo' The following objects are masked from 'package:base': as.Date, as.Date.numeric Loading required package: sandwich > gefp <- gefp(x = scale(scores$scores, scale = FALSE), fit = NULL, + scores = function(x) {x}, order.by = metric) > gefp$process[-1,] -1.3771 2.420033e-04 -3.334786e-01 -1.006781e-01 -2.818289e-01 -2.725019e-01 -0.7075 7.404527e-02 -2.365103e-01 -2.790570e-01 -3.512988e-01 -9.589193e-02 -0.6888 -1.210368e-01 -7.134302e-02 -5.975643e-01 1.935734e-01 9.981391e-02 -0.415 1.275834e-02 1.254864e-01 -3.028179e-01 2.360873e-01 -1.017536e-01 -0.3943 9.917085e-02 3.384026e-01 -4.148226e-01 3.013782e-01 -1.746471e-01 -0.2534 1.558282e-01 2.450313e-01 -2.763645e-01 2.776883e-01 1.168400e-01 -0.1645 2.884038e-01 4.476681e-01 2.711118e-02 3.231950e-01 -8.003818e-02 -0.1123 3.471061e-01 3.641271e-01 1.529694e-01 2.996919e-01 2.065794e-01 -0.1028 4.640467e-01 3.709941e-01 2.077436e-01 2.405129e-01 -5.336269e-02 -0.0593 5.661012e-01 6.017797e-01 1.599709e-01 3.079788e-01 -1.717531e-01 -0.0538 -3.424868e-01 6.527641e-01 3.862716e-01 3.117870e-01 -3.015950e-01 0.3646 -2.733318e-01 6.244143e-01 4.354578e-01 2.866215e-01 -4.371310e-02 0.3877 -2.539443e-01 7.702933e-01 9.482611e-02 1.029992e-01 7.482810e-02 0.5567 -2.736185e-01 3.940135e-01 -1.394281e-02 -2.107656e-01 -1.963198e-01 0.697 -1.721636e-01 6.246434e-01 -8.422372e-02 -1.430374e-01 -3.131174e-01 0.7632 -2.850898e-01 7.629785e-01 -5.225636e-01 -5.245715e-01 -2.305197e-01 0.7685 -2.012700e-01 2.218021e-01 -7.402078e-01 3.678912e-02 -4.604837e-01 0.8811 -6.493485e-02 3.961034e-01 -4.827298e-01 6.781724e-02 -6.761841e-01 1.1 -1.054789e-02 2.920380e-01 -3.309072e-01 4.379308e-02 -3.794727e-01 1.3587 1.489738e-17 -3.595552e-17 -2.750251e-16 1.048080e-16 -7.145443e-17 -1.3771 -5.781780e-03 -0.7075 1.578973e-02 -0.6888 6.572980e-02 -0.415 8.266642e-02 -0.3943 -7.856663e-01 -0.2534 -7.919143e-01 -0.1645 -7.680210e-01 -0.1123 -7.738898e-01 -0.1028 -7.878639e-01 -0.0593 -4.507996e-01 -0.0538 -4.672783e-01 0.3646 -4.691381e-01 0.3877 -4.241853e-01 0.5567 -4.296219e-01 0.697 -8.059560e-02 0.7632 -1.974634e-02 0.7685 -1.952819e-03 0.8811 -1.154427e-03 1.1 -7.702091e-03 1.3587 -2.648082e-17 > > # Is the cumsum score process the same? > max(abs(gefp$process[-1,] - process)) [1] 5.140333e-14 > stopifnot(all(abs(gefp$process[-1,] - process) < 1e-8)) > > # bootstrap vs sctest > test_b <- bootstrap_sctest(resp = resp, a = a, b = b, nSample = 15, + item_selection = 1, + DIF_covariate = list(m1 = metric, + m2 = metric, + m3 = metric, + f = factor, + o1 = order, + o2 = order), + statistic = c("auto", "CvM", "maxLM", "auto", "auto", "WDMo"), + decorrelate = TRUE, theta = theta) > test_b$p it1 m1 1.0000000 m2 1.0000000 m3 1.0000000 f 0.8666667 o1 0.8000000 o2 0.7333333 > test_b$statistic it1 m1 0.5661012 m2 0.1314926 m3 2.2278215 f 2.4590317 o1 1.9481736 o2 1.2594882 > sctest_m1 <- sctest(x = scale(scores$scores, scale = FALSE), + scores = function(x) {x}, parm = c(1, 4), order.by = metric) > sctest_m2 <- sctest(x = scale(scores$scores, scale = FALSE), + scores = function(x) {x}, parm = c(1, 4), order.by = metric, + functional = "CvM") > sctest_m3 <- sctest(x = scale(scores$scores, scale = FALSE), + scores = function(x) {x}, parm = c(1, 4), order.by = metric, + functional = "maxLM") > sctest_f <- sctest(x = scale(scores$scores, scale = FALSE), + scores = function(x) {x}, parm = c(1, 4), order.by = factor, + functional = "LMuo") > sctest_o <- sctest(x = scale(scores$scores, scale = FALSE), + scores = function(x) {x}, parm = c(1, 4), order.by = order, + functional = "maxLMo") > sctest_o2 <- sctest(x = scale(scores$scores, scale = FALSE), + scores = function(x) {x}, parm = c(1, 4), order.by = order, + functional = "WDMo") > > # Are the test statistics the same > abs(test_b$statistic / rbind(sctest_m1$statistic, + sctest_m2$statistic, + sctest_m3$statistic, + sctest_f$statistic, + sctest_o$statistic, + sctest_o2$statistic) - 1) it1 m1 8.881784e-16 m2 1.554312e-15 m3 4.440892e-15 f 5.551115e-16 o1 8.881784e-16 o2 8.881784e-16 > > stopifnot(all(abs(test_b$statistic / rbind(sctest_m1$statistic, + sctest_m2$statistic, + sctest_m3$statistic, + sctest_f$statistic, + sctest_o$statistic, + sctest_o2$statistic) - 1) < 1e-8)) > > # permutation vs bootstrap > test_p <- permutation_sctest(resp = resp, a = a, b = b, nSample = 15, + item_selection = 1, + DIF_covariate = list(m1 = metric, + m2 = metric, + m3 = metric, + f = factor, + o1 = order, + o2 = order), + statistic = c("auto", "CvM", "maxLM", "auto", "auto", "WDMo"), + decorrelate = TRUE, theta = theta) > test_p$p it1 m1 1.0000000 m2 1.0000000 m3 0.9333333 f 1.0000000 o1 0.8666667 o2 0.7333333 > test_p$statistic it1 m1 0.5661012 m2 0.1314926 m3 2.2278215 f 2.4590317 o1 1.9481736 o2 1.2594882 > > test_b$statistic it1 m1 0.5661012 m2 0.1314926 m3 2.2278215 f 2.4590317 o1 1.9481736 o2 1.2594882 > test_p$statistic it1 m1 0.5661012 m2 0.1314926 m3 2.2278215 f 2.4590317 o1 1.9481736 o2 1.2594882 > > # Are the test statistics the same > stopifnot(all(test_b$statistic == test_p$statistic)) > > > > proc.time() user system elapsed 2.65 0.71 3.28