{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":45533,"databundleVersionId":5748852,"sourceType":"competition"}],"dockerImageVersionId":30433,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Psychometric analysis of the Jo Wilder online educational game: some insights from the Item Response Theory\n---\n\nIn this notebook I analyse the data from the competition [Predict Student Performance from Game Play – Trace student learning from Jo Wilder online educational game](https://www.kaggle.com/competitions/predict-student-performance-from-game-play) to evaluate the properties of the game from a psychometric perspective leveraging the models of the Item Response Theory (IRT).\n\nThe game ([here for play](https://pbswisconsineducation.org/jowilder/play-the-game/)) is composed of 18 questions (\"items\"). During the game, for each question a player performs a sequence of actions leading to a dichotomous score assessing the correctness of the response (0: fail, 1: hit). In this notebook, I analyze only the final scores of the matrix _subject sessions_ $\\times$ _question items_, evaluating if the test produces consistent measures of individuals.\n\n![image](https://content.pbswisconsineducation.org/wp-content/uploads/2021/06/25221406/jw-facebook-about.png)","metadata":{}},{"cell_type":"markdown","source":"# 1. Theoretical background\n\nWe can consider the game as a performance test with a dichotomous outcome measuring a latent variable. Despite this latent variable has a continuous nature, we can only observe whether it assumes a sufficient level to lead to a success (question hit), or not (question fail).\n\nLet $\\varphi$ be the latent variabile under measurement. By calculating the log-odds of the probability of success in a question ($\\pi$), we can measure the latent trait on a logit scale:\n\n$$ \\varphi = log( \\frac{\\pi}{1 - \\pi} ) $$\n\nWhen $\\varphi = 0$, the player has 50% of probability having success; when $\\varphi < 0$, the probability of failure is greater than the probability of success, and when $\\varphi > 0$, the probability of a hit a question is greater than the probability of fail. Solving the equation for $\\pi$, we obtain the probability of a subject having success in an item:\n\n$$ \\pi = \\frac{e^\\varphi}{1 + e^\\varphi} $$\n\nBasing on this logistic formulation, during decades psychometrician has developed many models, falling under the frameworks of the **Rasch measurement model** and the **Item Response Theory (IRT)**. Basically, different models account for different concepts of what $\\varphi$ is, facing the many ways in which $\\varphi$ can manifest itself.\n\nConsidering the dichotomous case, both Rasch and IRT frameworks assume the latent trait $\\varphi$ depends on two factors: the ability of an individual $\\theta$ and the difficulty of an item $b$, in the following way:\n\n$$ \\varphi = \\theta - b $$\n\nSo, the model formulation for the dichotomous case will be:\n\n$$ \\pi = \\frac{e^{\\theta - b}}{1 + e^{\\theta - b}} $$\n\nThis model is known as dichotomous Rasch model in the Rasch framework, and **1PL** (one-parameter logistic) model in the IRT framework. Furthermore, IRT provides a more complex **2PL** model, which includes a discrimination parameter $a$:\n\n$$ \\pi = \\frac{e^{a (\\theta - b)}}{1 + e^{a(\\theta - b)}} $$\n\nThe **discrimination** is the capacity of an item to distinguish between individuals. Imagine having two students with very distant ability levels. In this case, you do not need to strong discriminating items to understand that they are different. Conversely, when two individuals have similar ability levels, you need to measure the latent trait by using strongly discriminating items. Be aware that introducing a discrimination parameter has several consequences. For example, a person can have a greater probability of having success in a difficult item than in an easy one.\n\nSince this is a simple analysis notebook, I don't dive further into more theoretical or statistical explanations. Anyway, you can read a short but illuminating paper from [Stemler & Naples (2021)](https://doi.org/10.7275/v2gd-4441) for a good introduction.","metadata":{}},{"cell_type":"markdown","source":"# 2. Software environment\n\nIn this notebook I leverage R and the `tidyverse` package system, using the package `mirt` for what concerns Rasch and IRT models, and `psych` to perform correlation analyses.","metadata":{}},{"cell_type":"code","source":"sessionInfo()","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:43:23.609662Z","iopub.execute_input":"2025-08-25T06:43:23.611523Z","iopub.status.idle":"2025-08-25T06:43:23.834892Z","shell.execute_reply":"2025-08-25T06:43:23.833679Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"remotes::install_version(\"psych\", version=\"2.3.9\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-25T06:43:23.836815Z","iopub.execute_input":"2025-08-25T06:43:23.862699Z","iopub.status.idle":"2025-08-25T06:44:06.089234Z","shell.execute_reply":"2025-08-25T06:44:06.070453Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"remotes::install_version(\"mirt\", version=\"1.40\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-25T06:44:06.090963Z","iopub.execute_input":"2025-08-25T06:44:06.091811Z","iopub.status.idle":"2025-08-25T06:47:08.576013Z","shell.execute_reply":"2025-08-25T06:47:08.574659Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"library(psych)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:08.577844Z","iopub.execute_input":"2025-08-25T06:47:08.578995Z","iopub.status.idle":"2025-08-25T06:47:08.641071Z","shell.execute_reply":"2025-08-25T06:47:08.639810Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"library(mirt)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:08.642634Z","iopub.execute_input":"2025-08-25T06:47:08.643473Z","iopub.status.idle":"2025-08-25T06:47:09.419771Z","shell.execute_reply":"2025-08-25T06:47:09.418447Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"library(tidyverse)","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2025-08-25T06:47:09.421510Z","iopub.execute_input":"2025-08-25T06:47:09.422456Z","iopub.status.idle":"2025-08-25T06:47:10.033609Z","shell.execute_reply":"2025-08-25T06:47:10.032385Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"library(doParallel)\nnum_cores <- detectCores()\nsprintf(\"Attached doParallel v%s (%s cores available)\", packageVersion(\"doParallel\"), num_cores)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-25T06:47:10.035221Z","iopub.execute_input":"2025-08-25T06:47:10.036396Z","iopub.status.idle":"2025-08-25T06:47:10.102968Z","shell.execute_reply":"2025-08-25T06:47:10.101806Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 3. Data exploration\n\nData input and preview:","metadata":{}},{"cell_type":"code","source":"game <- \"/kaggle/input/predict-student-performance-from-game-play/\" %>%\n    paste0(\"train_labels.csv\") %>%\n    readr::read_csv(show_col_types = FALSE) %>%\n    dplyr::mutate(\n        correct = as.integer(correct),\n        session_id = strsplit(session_id, \"_\"),\n        id = sapply(session_id, `[[`, 1),\n        item = sapply(session_id, `[[`, 2)\n    ) %>%\n    dplyr::select(-session_id) %>%\n    tidyr::pivot_wider(names_from=\"item\", values_from=\"correct\") %>%\n    tibble::column_to_rownames(\"id\")\n\nitems <- colnames(game)\n\nstr(game)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:10.104648Z","iopub.execute_input":"2025-08-25T06:47:10.105555Z","iopub.status.idle":"2025-08-25T06:47:12.465996Z","shell.execute_reply":"2025-08-25T06:47:12.464841Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"In the following, I print some plots to take a first look at the data.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=6)\n\ngame %>%\n    apply(MARGIN=2, FUN=sum) %>%\n    cbind() %>%\n    as.data.frame() %>%\n    dplyr::rename(\"freq\"=\".\") %>%\n    tibble::rownames_to_column(\"item\") %>%\n    dplyr::mutate(\n        perc=freq/nrow(game)*100,\n        item=substr(item, 2, nchar(item)),\n        item=factor(item, levels=as.character(1:18)),\n        label=sprintf(\"%.1f%%\", perc)\n    ) %>%\n    ggplot(aes(x=item, y=perc)) +\n    geom_segment(aes(x=item, y=0, xend=item, yend=perc)) +\n    geom_point() +\n    geom_text(aes(label=label), position = position_nudge(y=5)) +\n    ggtitle(\"Percentage of correct response by item\") +\n    xlab(\"Item\") + ylab(\"% of hits\") +\n    theme_grey(base_size=18)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:12.467603Z","iopub.execute_input":"2025-08-25T06:47:12.468451Z","iopub.status.idle":"2025-08-25T06:47:13.096424Z","shell.execute_reply":"2025-08-25T06:47:13.095183Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nAt first glance, the percentage of individuals getting success varies enough between items. The game appears to cover a good spectrum of difficulties (they aren't everything easy, nor everything hard).","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=10, repr.plot.height=8)\n\ngame %>%\n    psych::tetrachoric() %>%\n    magrittr::extract2(\"rho\") %>%\n    data.frame() %>%\n    tibble::rownames_to_column(\"from\") %>%\n    tidyr::pivot_longer(cols=tidyselect::all_of(items), names_to=\"to\", values_to=\"r\") %>%\n    dplyr::mutate(\n        from=factor(from, levels=items),\n        to=factor(to, levels=rev(items)),\n        lab=sprintf(\"%.2f\", r),\n        is_diag=from==to,\n        r=ifelse(is_diag, NA, r),\n        lab=ifelse(is_diag, \"\", lab),\n        lab=substr(lab, 2, nchar(lab))\n        \n    ) %>%\n    ggplot(aes(x=from, y=to)) +\n    geom_tile(aes(fill=r)) +\n    geom_text(aes(label=lab), colour=\"#21759B\") +\n    scale_fill_gradientn(colors=rev(hcl.colors(20,\"OrRd\")), limits=c(0,0.5)) +\n    labs(fill=\"Corr\") +\n    ggtitle(\"Tetrachoric correlation matrix between items\") +\n    xlab(NULL) + ylab(NULL) +\n    theme_grey(base_size=18)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:13.098248Z","iopub.execute_input":"2025-08-25T06:47:13.099262Z","iopub.status.idle":"2025-08-25T06:47:14.105632Z","shell.execute_reply":"2025-08-25T06:47:14.104225Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nThe correlations are always positive and they aren't upper than 0.45. This is good notice, because the questions seem to not be redundant. On the whole, however, the correlations are low, but this is expected because of the nature of the game. Two items (16 and 17) have very low correlations with the whole group of items, as if they are measuring a different latent trait.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=6)\n\ngame %>%\n    apply(MARGIN=1, FUN=sum) %>%\n    cbind() %>%\n    as.data.frame() %>%\n    dplyr::rename(\"count\"=\".\") %>%\n    tibble::rownames_to_column(\"person\") %>%\n    dplyr::mutate(count=factor(count, levels=0:18)) %>%\n    dplyr::group_by(count) %>%\n    dplyr::summarise(freq=dplyr::n(), .groups=\"drop\") %>%\n    tidyr::complete(count, fill=list(freq=0L)) %>%\n    dplyr::mutate(perc=freq/nrow(game)*100) %>%\n    ggplot(aes(x=count, y=perc)) +\n    geom_bar(stat=\"identity\") +\n    xlab(\"Number of hits\") +\n    ylab(\"% of individuals\") +\n    ggtitle(\"Distribution of total scores\") +\n    theme_grey(base_size=18)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:14.107737Z","iopub.execute_input":"2025-08-25T06:47:14.108801Z","iopub.status.idle":"2025-08-25T06:47:14.426016Z","shell.execute_reply":"2025-08-25T06:47:14.424793Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nThe distribution of the total score is not symmetric. The scoring leaves a great span to measure abilities under the mode, and fewer scores to measure high abilities.","metadata":{}},{"cell_type":"markdown","source":"# 4. Fitting the Rasch (1PL) model","metadata":{}},{"cell_type":"code","source":"model1pl <- mirt(game, itemtype=\"Rasch\", SE=TRUE)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:14.427932Z","iopub.execute_input":"2025-08-25T06:47:14.428900Z","iopub.status.idle":"2025-08-25T06:47:15.985193Z","shell.execute_reply":"2025-08-25T06:47:15.983812Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Extracting the item parameters:","metadata":{}},{"cell_type":"code","source":"param1pl <- model1pl %>%\n    coef(IRTpars=TRUE, simplify=TRUE) %>% # printSE=TRUE to extract item std. err.\n    magrittr::extract2(\"items\") %>%\n    data.frame() %>%\n    tibble::rownames_to_column(\"item\")\n\nparam1pl","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:15.987021Z","iopub.execute_input":"2025-08-25T06:47:15.988524Z","iopub.status.idle":"2025-08-25T06:47:16.030258Z","shell.execute_reply":"2025-08-25T06:47:16.026395Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The only estimated parameter is $b$ (difficulty), while the other parameters are bounded ($a$ is for 2PL, $g$ is for 3PL and $u$ is for 4PL).\n\nEstimating the person parameters:","metadata":{}},{"cell_type":"code","source":"person1pl <- model1pl %>%\n    fscores(method=\"WLE\", full.scores.SE=TRUE) %>%\n    data.frame()\n\nstr(person1pl)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:47:16.033581Z","iopub.execute_input":"2025-08-25T06:47:16.035726Z","iopub.status.idle":"2025-08-25T06:48:29.148604Z","shell.execute_reply":"2025-08-25T06:48:29.147435Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"_F1_: the first (and unique) latent factor, measured in logits; these are the $\\theta$ parameters.\n\n_SE_F1_: the estimated standard errors of $\\theta$ parameters.","metadata":{}},{"cell_type":"markdown","source":"## 4.1 Evaluation of test calibration","metadata":{}},{"cell_type":"markdown","source":"The histogram below compares the location of item difficulties against the distribution of the abilities of the individuals, putting both on the continuum of the latent trait that the test should measure.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=8)\n\nbin <- 0.25\n\nparam1pl_shifted <- param1pl %>%\n    dplyr::mutate(\n        group=cut(b, breaks=seq(-5,3,by=bin))\n    ) %>%\n    dplyr::group_by(group) %>%\n    dplyr::mutate(\n        y=seq(1,dplyr::n())*(-150)\n    ) %>%\n    dplyr::rename(\"logit\"=\"b\") %>%\n    dplyr::ungroup()\n\nperson1pl %>%\n    dplyr::rename(\"logit\"=\"F1\") %>%\n    ggplot(aes(x=logit)) +\n    geom_histogram(binwidth=bin, fill=\"#61868d\", colour=\"white\") +\n    geom_label(\n        data=param1pl_shifted,\n        mapping=aes(x=logit, y=y, label=item),\n        angle=90\n    ) +\n    theme_gray(base_size=18) +\n    ggtitle(\"Person-item map\") +\n    xlab(\"Latent trait (logits)\") +\n    ylab(\"Number of person\")","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.150280Z","iopub.execute_input":"2025-08-25T06:48:29.151154Z","iopub.status.idle":"2025-08-25T06:48:29.485105Z","shell.execute_reply":"2025-08-25T06:48:29.483730Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nApproximately 2/3 of the spectrum of abilities is covered by items, especially the center part, which is the most dense and requires fine discrimination between subjects. The upper part is poorly covered, so the test may not be capable of discriminating between excellence. The distribution of total scores had already shown a weak tendency towards the ceiling effect.","metadata":{}},{"cell_type":"markdown","source":"## 4.2 Reliability of the test measures\n\nThe overall reliability of the instrument is evaluated by three indices.\n\n- **Separation reliability _R_.** Is a Rasch version of KR-20 or Cronbach's $\\alpha$ indices, calculated as the proportion of the true variance on the observed variance, where the _true_ variance is the genuine value of whatever is being measured. Roughly, good values are greater than 0.7. Range: [0, 1].\n\n- **Separation ratio _G_.** A high variability in persons’ abilities and in item difficulties is necessary so that the measurement is reliable. This index quantifies the number of strata of abilities that the instrument can separate and identify as different. Low values (_G_ < 2) indicate that the instrument may not be sensitive enough to distinguish between individuals with different ability levels. When _G_ = 2, the test can distinguish between high and low performers. Range: [0, Inf).\n\n- **Discernible strata _H_.** This is a transformation of _G_ suggested in the case of asymmetries in outliers. Good values are at least _H_ = 3. When _H_ = 3, the test can distinguish between very high, middle and very low performers. Range: [1/3, Inf).","metadata":{}},{"cell_type":"code","source":"separation <- function(pars, SE) {\n    n <- length(pars)\n    SD2 <- var(pars)*(n-1)/n\n    MSE <- sum(SE^2)/n\n    SA2 <- SD2-MSE\n    R <- 1-MSE/SD2 # SA2/SD2\n    G <- sqrt(SA2/MSE)\n    H <- (4*G+1)/3\n    return(data.frame(\"R\"=R, \"G\"=G, \"H\"=H))\n}\n\nwith(person1pl, separation(F1, SE_F1))","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.486921Z","iopub.execute_input":"2025-08-25T06:48:29.487913Z","iopub.status.idle":"2025-08-25T06:48:29.515686Z","shell.execute_reply":"2025-08-25T06:48:29.514423Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\n_R_ is borderline, and _G_ and _H_ are low, albeit not dramatically. This can be due to a lack of questions capable of measuring the higher part of the latent trait.","metadata":{}},{"cell_type":"markdown","source":"## 4.3 Goodness of fit of the items","metadata":{}},{"cell_type":"markdown","source":"The item fit is evaluated by using the indices [Infit and Outfit](https://www.rasch.org/rmt/rmt162f.htm). The table below reports these values, in two versions: the basic MSQ indices and their _t_ conversion. Since the large sample size, is difficult to interpret _t_ values, so in the next I'll discuss only the MSQ values.","metadata":{}},{"cell_type":"code","source":"(item_stats1pl <- itemfit(model1pl, fit_stats=\"infit\", Theta=person1pl$F1))","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.517436Z","iopub.execute_input":"2025-08-25T06:48:29.518382Z","iopub.status.idle":"2025-08-25T06:48:29.876733Z","shell.execute_reply":"2025-08-25T06:48:29.875677Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The literature indicates various threshold limits to interpret MSQ values. The commonly accepted range for a good fit is [0.6, 1.4] ([Bond & Fox, 2007](https://www.taylorfrancis.com/books/mono/10.4324/9781410614575/applying-rasch-model-trevor-bond-christine-fox)). Other approaches are similar and anyway symmetric around 1.\n\nSo, our items appear really good. **Is everything okay? Not, at all**.\n\nSimulation studies have shown that appropriate critical values depend on the sample size, the number of items, and their difficulty. In an impressive study, [Müller, 2020](https://jsdajournal.springeropen.com/articles/10.1186/s40488-020-00108-7) shows that the range of good fit can be different between Infit and Outfit, and even not centered around 1 for high sample sizes. This is especially true when considering unconditional estimates such as those returned by almost all software packages, mirt included. And, unfortunately, here we have a very large sample size (according to the standards of the Rasch literature).\n\nFor this reason, in the code below I set a simulation to estimate the range of theoretical good fit for items. Stating by the estimated $\\theta$ and $b$, for 500 iterations I simulated a data matrix conforming to this model, I fitted a Rasch model with `mirt` and I calculated Infits and Outfits. At the end, I calculated the quantiles of resulting Infit and Outfit distributions for each item.","metadata":{}},{"cell_type":"code","source":"irt_simulation <- function(param, theta) {\n    check <- FALSE\n    while(!check) {\n        simulated_data <- theta %>%\n            outer(param$b, function(theta, b) 1/(1+exp(-(theta-b)))) %>%\n            sapply(rbinom, n=1, size=1) %>%\n            matrix(nrow=length(theta), dimnames=list(NULL, param$item))\n        check <- simulated_data %>%\n            apply(2, function(x) length(unique(x))>1) %>%\n            all()\n    }\n    model <- simulated_data %>%\n        data.frame() %>%\n        mirt(itemtype=\"Rasch\", verbose=FALSE)\n    simulated_theta <- fscores(model, method=\"WLE\", full.scores.SE=FALSE)\n    stats <- model %>%\n        itemfit(fit_stats=\"infit\", Theta=simulated_theta) %>%\n        dplyr::select(item, outfit, infit)\n    return(stats)\n}","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.878297Z","iopub.execute_input":"2025-08-25T06:48:29.879097Z","iopub.status.idle":"2025-08-25T06:48:29.886785Z","shell.execute_reply":"2025-08-25T06:48:29.885694Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**Warning:** since the code is very slow, I ran it on my machine and I put the resulting data frame in a designed chunk below.","metadata":{}},{"cell_type":"code","source":"if(FALSE) {\n\nregisterDoParallel(num_cores-1)\nsim1pl <- foreach(i=1:500, .combine=\"list\", .packages=c(\"mirt\",\"magrittr\",\"dplyr\")) %dopar% {\n    irt_simulation(param1pl, person1pl$F1)\n}\nstopImplicitCluster()\n\nsim1pl_quantiles <- sim1pl %>%\n    dplyr::bind_rows(.id=\"iter\") %>%\n    tidyr::pivot_longer(\n        cols=all_of(c(\"outfit\",\"infit\")),\n        names_to=\"index\", values_to=\"fit\"\n    ) %>%\n    dplyr::group_by(item, index) %>%\n    dplyr::summarise(\n        Q_med=median(fit),\n        Q_low=quantile(fit, probs=0.05/2, names=FALSE),\n        Q_upp=quantile(fit, probs=1-0.05/2, names=FALSE),\n        .groups=\"drop\"\n    )\n\n}","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.888343Z","iopub.execute_input":"2025-08-25T06:48:29.889163Z","iopub.status.idle":"2025-08-25T06:48:29.896678Z","shell.execute_reply":"2025-08-25T06:48:29.895577Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Dump of simulation results\nsim1pl_quantiles <-\nstructure(list(item = c(\"q1\", \"q1\", \"q10\", \"q10\", \"q11\", \"q11\", \n\"q12\", \"q12\", \"q13\", \"q13\", \"q14\", \"q14\", \"q15\", \"q15\", \"q16\", \n\"q16\", \"q17\", \"q17\", \"q18\", \"q18\", \"q2\", \"q2\", \"q3\", \"q3\", \"q4\", \n\"q4\", \"q5\", \"q5\", \"q6\", \"q6\", \"q7\", \"q7\", \"q8\", \"q8\", \"q9\", \"q9\"\n), index = c(\"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \n\"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \n\"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \n\"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \n\"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \"infit\", \"outfit\", \n\"infit\", \"outfit\"), Q_med = c(0.94893943713553708, 0.88352916198706644, \n0.9561219877471111, 0.91374134029277287, 0.95431532284420617, \n0.89807329881654385, 0.93115869949446739, 0.84719250633111465, \n0.93621093652413945, 0.90906802141329335, 0.94986053304472662, \n0.88622109969382112, 0.95591226858487932, 0.91524552625571087, \n0.94766008100800536, 0.88078647760908557, 0.95124290195093031, \n0.88996357638496582, 0.90087935711842149, 0.81214458551461943, \n0.87482863910237374, 0.79027468847828708, 0.90947487535558125, \n0.8213186727526155, 0.94120973019444398, 0.86570302778199903, \n0.95685372646296252, 0.9111198108704186, 0.94414406692045838, \n0.87187268963070064, 0.94844239538650965, 0.88147199053191194, \n0.95426474701856467, 0.90168778147202233, 0.94805044051234622, \n0.8809411589807572), Q_low = c(0.9383836536614949, 0.86544116260612336, \n0.94603051734979937, 0.89914381238561059, 0.94353194822711706, \n0.88288746304628662, 0.91967455887947536, 0.82132061127154088, \n0.92540971294333429, 0.88834016232842838, 0.93961683548276564, \n0.86874492239886059, 0.94484576342518556, 0.90139093996769526, \n0.93808733298081881, 0.86395389855763072, 0.94145238169194556, \n0.87466730536141202, 0.88772453925286321, 0.75664938276733518, \n0.85950815014280812, 0.71422197160097523, 0.89652422137526966, \n0.77728514313342822, 0.93024168261999673, 0.84534464062632531, \n0.94646212235611626, 0.89674572521221085, 0.93279521330639259, \n0.85302620290264142, 0.93738212084587769, 0.86401349154243801, \n0.94522554031608064, 0.88719952593327789, 0.93724663457587842, \n0.86207883206812719), Q_upp = c(0.95931332798567848, 0.89880938884501582, \n0.96718445376813911, 0.92858818085611405, 0.96475660561059873, \n0.9140955723853782, 0.94216443685827134, 0.87311062465723133, \n0.94626350351714195, 0.93243419853052267, 0.96114947406203566, \n0.90396131990713424, 0.96632810469595576, 0.9307186652291275, \n0.95946870142848262, 0.89876078583794872, 0.96189163462809457, \n0.90599177034313239, 0.91420055564231395, 0.86300020759042972, \n0.88880184041366861, 0.8708154554481975, 0.92154103736776283, \n0.86246419855223155, 0.95230549528116371, 0.88648467221976834, \n0.96751935117133603, 0.92513160671921968, 0.95515415427837924, \n0.89231257144958942, 0.95772689360954621, 0.89831104081191193, \n0.96492092274517138, 0.91655863851321306, 0.95885070570501318, \n0.89881152104162187)), class = c(\"tbl_df\", \"tbl\", \"data.frame\"\n), row.names = c(NA, -36L))","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.898316Z","iopub.execute_input":"2025-08-25T06:48:29.899152Z","iopub.status.idle":"2025-08-25T06:48:29.907632Z","shell.execute_reply":"2025-08-25T06:48:29.906449Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The plot below shows the results, putting on the x-axis the difficulties of items, and on the y-axis the observed Infit and Outfit (those in the table above). The shaded area is defined by the quantiles of simulated data and indicates the goodness of fit limits. The items exceeding the area on the upper side are underfitting, while items exceeding the area on the lower side are overfitting (overfitting is generally indicating redundancy and it's less problematic).\n\nThe output is really different if compared to what is expected from classic literature. Both Infit and Outfit are clearly influenced by the item difficulties but in different ways. Often, good values are less than 1. Anyway, be generous when evaluating the items, because this method is an experimental approach and may not be completely free from bias.","metadata":{}},{"cell_type":"code","source":"item_stats_long <- item_stats1pl %>%\n    tidyr::pivot_longer(\n        cols=all_of(c(\"outfit\",\"infit\")),\n        names_to=\"index\", values_to=\"fit\"\n    ) %>%\n    dplyr::left_join(param1pl, by=\"item\") %>%\n    dplyr::mutate(item_label=substr(item,2,nchar(item)))\n\nitem_stats_long %>%\n    dplyr::left_join(sim1pl_quantiles, by=c(\"item\", \"index\")) %>%\n    dplyr::mutate(item=factor(item, levels=items)) %>%\n    ggplot(aes(x=b, y=fit)) +\n    geom_ribbon(aes(ymin=Q_low, ymax=Q_upp), alpha=0.5, fill=\"#61868d\") +\n    geom_segment(x=-Inf, xend=Inf, y=1, yend=1, colour=\"orange\") +\n    geom_point(pch=21, bg=\"white\", size=8) +\n    geom_text(aes(x=b, y=fit, label=item_label)) +\n    facet_wrap(~index) +\n    theme_gray(base_size=18) +\n    xlab(\"Item difficulty\") + ylab(\"MSQ\")","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:29.909269Z","iopub.execute_input":"2025-08-25T06:48:29.910091Z","iopub.status.idle":"2025-08-25T06:48:30.267933Z","shell.execute_reply":"2025-08-25T06:48:30.266575Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The plot below shows the absolute distance from the goodness-of-fit area of each item.","metadata":{}},{"cell_type":"code","source":"item_stats_long <- item_stats_long %>%\n    dplyr::left_join(sim1pl_quantiles, by=c(\"item\",\"index\")) %>%\n    dplyr::mutate(item=factor(item, levels=items)) %>%\n    dplyr::mutate(\n        dist_from_med=fit-Q_med,\n        dist_from_lim=ifelse(dist_from_med>0, fit-Q_upp, fit-Q_low),\n        dist_from_lim=abs(dist_from_lim),\n        dist_from_lim=ifelse(index==\"outfit\", -dist_from_lim, dist_from_lim)\n    )\n\naxis_brk <- item_stats_long %>%\n    dplyr::pull(dist_from_lim) %>%\n    abs() %>%\n    max() %>%\n    `*`(10) %>%\n    ceiling() %>%\n    `/`(10) %>%\n    {seq(-., ., by=0.1)}\n\naxis_lab <- ifelse(axis_brk<0, -axis_brk, axis_brk) %>% sprintf(fmt=\"%.1f\")\n\nggplot(item_stats_long, aes(x=item, y=dist_from_lim, fill=index)) +\n    geom_bar(stat=\"identity\", colour=\"black\") +\n    scale_y_continuous(\n        limits=range(axis_brk), breaks=axis_brk, labels=axis_lab\n    ) +\n    scale_fill_manual(values=c(\"#b0bf1a\",\"#21759B\"), name=\"Index\") +\n    theme_gray(base_size=18) +\n    xlab(\"Item\") + ylab(\"Absolute distance from the MSQ limit\")","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:30.269980Z","iopub.execute_input":"2025-08-25T06:48:30.270970Z","iopub.status.idle":"2025-08-25T06:48:30.605613Z","shell.execute_reply":"2025-08-25T06:48:30.604291Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nAccording to both Infit and Outfit, the items `q16` and `q17` present the strongest underfit. This appears clear in the barplot, where I show the absolute distance of each item from the goodness-of-fit area. Moreover, the items `q6`, `q8`, `q13`, and `q18` aren't much better. `q18` has a low Infit distance and a strong Outfit distance, indicating that this item is infested by outliers.\n\nAnyway, `q17` and `q16` are the most sources of concern, because they are strongly underfitting (in the Rasch model, underfit is degrading, while overfit is simply unproductive).","metadata":{}},{"cell_type":"markdown","source":"## 4.4 Item characteristic curves for 1PL\n\nItem characteristic curves show the probability of answering correctly to each item along the latent trait. In the plot below, I overlap the curves upon the observed proportions of correct answers, calculated after slitting the sample into 10 groups according to the quantiles of $\\theta$.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=12)\n\nperson1pl %>%\n    dplyr::bind_cols(game) %>%\n    tidyr::pivot_longer(\n        cols=dplyr::all_of(colnames(game)),\n        names_to=\"item\", values_to=\"score\"\n    ) %>%\n    dplyr::rename(\"theta\"=\"F1\") %>%\n    dplyr::mutate(\n        group = cut(theta,\n            breaks=quantile(theta, probs=seq(0,1,by=0.1), names=FALSE),\n            include.lowest=TRUE, right=TRUE\n        )\n    ) %>%\n    dplyr::group_by(item, group) %>%\n    dplyr::summarise(\n        score=mean(score),\n        theta=mean(theta),\n        .groups=\"drop\"\n    ) %>%\n    dplyr::left_join(\n        dplyr::select(param1pl, item, b), by=\"item\"\n    ) %>%\n    dplyr::mutate(\n        latent=theta-b,\n        prob=exp(latent)/(1+exp(latent)),\n        item=factor(item, levels=colnames(game))\n    ) %>%\n    ggplot(aes(x=theta, y=score)) +\n    geom_point(pch=21, bg=\"#cde3e7\", size=3) +\n    geom_line(aes(y=prob)) +\n    facet_wrap(~item, nrow=3) +\n    theme_gray(base_size=18) +\n    xlab(\"Person ability\") +\n    ylab(\"Probability of success\") ","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:30.607463Z","iopub.execute_input":"2025-08-25T06:48:30.608366Z","iopub.status.idle":"2025-08-25T06:48:31.982292Z","shell.execute_reply":"2025-08-25T06:48:31.981061Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nItems `q16` and `q17` present evidence of discrimination issues, because their curves require a different slope to have a better fit of data points. They don't seem to be the only ones. Then, let's go to extend the model with discrimination parameters by fitting the 2PL model.","metadata":{}},{"cell_type":"markdown","source":"# 5. Fitting the 2PL model","metadata":{}},{"cell_type":"code","source":"model2pl <- mirt(game, itemtype=\"2PL\", SE=TRUE)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:31.983909Z","iopub.execute_input":"2025-08-25T06:48:31.984814Z","iopub.status.idle":"2025-08-25T06:48:33.701075Z","shell.execute_reply":"2025-08-25T06:48:33.698043Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Extracting the item parameters:","metadata":{}},{"cell_type":"code","source":"param2pl <- model2pl %>%\n    coef(IRTpars=TRUE, simplify=TRUE) %>%\n    magrittr::extract2(\"items\") %>%\n    data.frame() %>%\n    tibble::rownames_to_column(\"item\")\n\nparam2pl","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:33.704449Z","iopub.execute_input":"2025-08-25T06:48:33.706733Z","iopub.status.idle":"2025-08-25T06:48:33.835861Z","shell.execute_reply":"2025-08-25T06:48:33.834588Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Estimating the person parameters:","metadata":{}},{"cell_type":"code","source":"person2pl <- model2pl %>%\n    fscores(method=\"WLE\", full.scores.SE=TRUE) %>%\n    data.frame()\n\nstr(person2pl)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:48:33.837567Z","iopub.execute_input":"2025-08-25T06:48:33.838473Z","iopub.status.idle":"2025-08-25T06:49:48.808643Z","shell.execute_reply":"2025-08-25T06:49:48.807419Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Scatterplot between item difficulties and item discriminations:","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\n\nparam2pl %>%\n    dplyr::mutate(item=substr(item, 2, nchar(item))) %>%\n    ggplot(aes(x=b, y=a)) +\n    geom_segment(x=-Inf, xend=Inf, y=1, yend=1, colour=\"orange\") +\n    geom_point(pch=21, bg=\"white\", size=8) +\n    geom_text(aes(label=item)) +\n    theme_gray(base_size=18) +\n    xlab(\"2PL item difficulty\") + ylab(\"2PL item discrimination\")","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:49:48.810412Z","iopub.execute_input":"2025-08-25T06:49:48.812003Z","iopub.status.idle":"2025-08-25T06:49:49.035348Z","shell.execute_reply":"2025-08-25T06:49:49.034178Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nItems with $a > 1$ are better able to discriminate between ability levels near the inflection point of the curve. More discriminating items provide greater information about a respondent than do less discriminating items ([Hays et al, 2000](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1815384/)). This is the case of `q6`, which is the most discriminating item (it was misfitting according to 1PL).\n\nItems `q16` and `q17`, and also `q8` and `q9` show low discriminating power ($a < 1$), so they have less capability to differentiate between individuals with near abilities. They can be useful in discriminating between distant subjects.","metadata":{}},{"cell_type":"markdown","source":"## 5.1 Comparison between 2PL and 1PL","metadata":{}},{"cell_type":"markdown","source":"The graph below shows the relation between item difficulties estimated by 1PL and 2PL.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=10, repr.plot.height=8)\n\nparam1pl %>%\n    dplyr::select(item, b) %>%\n    dplyr::left_join(\n        dplyr::select(param2pl, item, a, b),\n        by=\"item\", suffix=c(\"1pl\",\"2pl\")\n    ) %>%\n    dplyr::mutate(item=substr(item, 2, nchar(item))) %>%\n    ggplot(aes(x=b1pl, y=b2pl)) +\n    geom_abline(intercept=0, slope=1) +\n    geom_point(aes(colour=a), size=8) +\n    geom_point(size=8, pch=21, colour=\"black\", bg=\"transparent\") +\n    geom_text(aes(label=item)) +\n    scale_colour_distiller(\n        type=\"div\", name=\"Discrim\",\n        breaks=seq(0.2,1.8, by=0.4),\n        limits=c(0.2, 1.8)\n    ) +\n    theme_gray(base_size=18) +\n    xlab(\"1PL item difficulty\") + ylab(\"2PL item difficulty\")","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:49:49.037092Z","iopub.execute_input":"2025-08-25T06:49:49.038022Z","iopub.status.idle":"2025-08-25T06:49:49.311939Z","shell.execute_reply":"2025-08-25T06:49:49.310784Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\n2PL and 1PL are clearly related, with the exception of `q16` and `q17`, whose difficulties – now we can understand – are strongly influenced by their low discrimination power.","metadata":{}},{"cell_type":"markdown","source":"## 5.2 Item characteristic curves for 2PL","metadata":{}},{"cell_type":"markdown","source":"Below I plot again the item characteristic curves, but corrected according to the 2PL model.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=12)\n\nperson2pl %>%\n    dplyr::bind_cols(game) %>%\n    tidyr::pivot_longer(\n        cols=dplyr::all_of(colnames(game)),\n        names_to=\"item\", values_to=\"score\"\n    ) %>%\n    dplyr::rename(\"theta\"=\"F1\") %>%\n    dplyr::mutate(\n        group = cut(theta,\n            breaks=quantile(theta, probs=seq(0,1,by=0.1), names=FALSE),\n            include.lowest=TRUE, right=TRUE\n        )\n    ) %>%\n    dplyr::group_by(item, group) %>%\n    dplyr::summarise(\n        score=mean(score),\n        theta=mean(theta),\n        .groups=\"drop\"\n    ) %>%\n    dplyr::left_join(\n        dplyr::select(param2pl, item, a, b), by=\"item\"\n    ) %>%\n    dplyr::mutate(\n        latent=a*(theta-b),\n        prob=exp(latent)/(1+exp(latent)),\n        item=factor(item, levels=colnames(game))\n    ) %>%\n    ggplot(aes(x=theta, y=score)) +\n    geom_point(pch=21, bg=\"#cde3e7\", size=3) +\n    geom_line(aes(y=prob)) +\n    facet_wrap(~item, nrow=3) +\n    theme_gray(base_size=18) +\n    xlab(\"Person ability\") +\n    ylab(\"Probability of success\") ","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:49:49.313572Z","iopub.execute_input":"2025-08-25T06:49:49.314422Z","iopub.status.idle":"2025-08-25T06:49:50.373701Z","shell.execute_reply":"2025-08-25T06:49:50.372437Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nBy including a discrimination parameter, the fit seems to have gained a lot, as confirmed by fit indices (lower RMSEA and SRMSR, higher TLI and CFI).","metadata":{}},{"cell_type":"code","source":"list(\"2PL\"=model2pl, \"1PL\"=model1pl) %>% sapply(M2, simplify=FALSE)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:49:50.375910Z","iopub.execute_input":"2025-08-25T06:49:50.376961Z","iopub.status.idle":"2025-08-25T06:49:51.485225Z","shell.execute_reply":"2025-08-25T06:49:51.483999Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6. Unidimensionality check\n\nIn order to evaluate the strength of the components, I run a principal component analysis (PCA) of standardized residuals ([here](https://www.winsteps.com/winman/principalcomponents.htm) for an overview). The idea is that since by IRT analysis, I extract the first, and – theoretically – the only latent component, residuals should be constituted only by random noise. So, PCA should be ineffective, because residuals should not have any common component.\n\nTo evaluate the result, I perforform a simulation study.\n\nI simulate 500 data matrices starting from the model probabilities. For each data matrix, I calculate the eigenvalues, obtaining a total of 500 estimations for each one. Finally, I obtain the 95th percentile of each eigenvalue distribution, plotting these quantiles against the observed ones (I show just the top nine).","metadata":{}},{"cell_type":"markdown","source":"Calculating the standardized residuals on observed data:","metadata":{}},{"cell_type":"code","source":"# Function to get the probability of hit given parameters\nget_prob <- function(theta, param) {\n    1 / (1 + exp(-param$a*(theta-param$b)))\n}\n\n# Get the probability matrix\nprobs2pl <- matrix(nrow=nrow(game), ncol=ncol(game), dimnames=dimnames(game))\nfor(i in 1:nrow(game)) {\n    for(j in 1:ncol(game)) {\n        probs2pl[i,j] <- get_prob(person2pl$F1[i], param2pl[j,])\n    }\n}\n\n# Get the information matrix\ninfo2pl <- sapply(items,\n    function(j) {\n        iteminfo(\n            x=extract.item(model2pl, item=j),\n            Theta=person2pl$F1\n        )\n    }\n)\n\n# Calculate the standardized residuals as raw residuals divided by the square root of information\nstdres2pl <- (game-probs2pl)/sqrt(info2pl)\n\n# Observed eigenvalues\nstdres2pl_eigen <- svd(stdres2pl)$d/sqrt(nrow(stdres2pl))","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:49:51.486931Z","iopub.execute_input":"2025-08-25T06:49:51.487778Z","iopub.status.idle":"2025-08-25T06:50:10.545865Z","shell.execute_reply":"2025-08-25T06:50:10.544352Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Simulating the eigenvalues:","metadata":{}},{"cell_type":"code","source":"# Function to get the probability of hit given parameters\nget_prob <- function(theta, param) {\n    1 / (1 + exp(-param$a*(theta-param$b)))\n}\n\n# Simulation function\nres_simulation <- function() {\n    # Simulate a dataset\n    check <- FALSE\n    while(!check) {\n        simulated_data <- probs2pl %>%\n            sapply(rbinom, n=1, size=1) %>%\n            matrix(ncol=nrow(param2pl), dimnames=list(NULL, param2pl$item))\n        check <- simulated_data %>%\n            apply(2, function(x) length(unique(x))>1) %>%\n            all()\n    }\n    # Fit the model\n    model <- simulated_data %>%\n        data.frame() %>%\n        mirt(itemtype=\"2PL\", verbose=FALSE)\n    # Get the model parameters\n    param <- model %>%\n        coef(IRTpars=TRUE, simplify=TRUE) %>% \n        magrittr::extract2(\"items\") %>%\n        data.frame()\n    # Estimate person parameters\n    theta <- model2pl %>%\n        fscores(method=\"WLE\", full.scores.SE=FALSE) %>% # WLE is too time-consuming\n        data.frame() %>%\n        dplyr::pull(F1)\n    # Get the information matrix\n    info_sq <- sapply(param2pl$item,\n        function(j) {\n            sqrt(iteminfo(\n                x=extract.item(model, item=j),\n                Theta=theta\n            ))\n        }\n    )\n    # Get standardized residuals\n    std_resid <- (simulated_data-probs2pl)/info_sq\n    # Get eigenvalues\n    eigen_val <- svd(std_resid)$d/sqrt(length(theta))\n    return(eigen_val)\n}","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:50:10.547805Z","iopub.execute_input":"2025-08-25T06:50:10.548791Z","iopub.status.idle":"2025-08-25T06:50:10.558513Z","shell.execute_reply":"2025-08-25T06:50:10.557395Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if(FALSE) {\n\ncl <- makeCluster(num_cores - 2)\nclusterSetRNGStream(cl, 666) # set seed\nsimres2pl <- parSapply(cl=cl, X=1:500, FUN=res_simulation)\nstopCluster(cl)\n\nsimres2pl_Q95 <- apply(simres2pl, 2, quantile, probs=0.95, names=FALSE)\n\n}","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:50:10.560212Z","iopub.execute_input":"2025-08-25T06:50:10.561090Z","iopub.status.idle":"2025-08-25T06:50:10.568634Z","shell.execute_reply":"2025-08-25T06:50:10.567508Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Dump of simulation results\nsimres2pl_Q95 <-\nc(3.2036155046821091, 2.9784239546175639, 1.5298754991436261, \n1.4477019462477843, 1.1668104961767276, 1.12673660077358, 1.0018987421583498, \n0.83554603015284079, 0.81654277810232956, 0.8027563285197109, \n0.78666466728479234, 0.77180300942139768, 0.75924150708249427, \n0.72190528903310436, 0.6534037040091073, 0.63532875446232506, \n0.62435787445525259, 0.57803260385381061)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:50:10.570337Z","iopub.execute_input":"2025-08-25T06:50:10.571843Z","iopub.status.idle":"2025-08-25T06:50:10.582233Z","shell.execute_reply":"2025-08-25T06:50:10.580436Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height=8)\n\ntibble::tibble(\n    Component=factor(1:length(items)),\n    Observed=stdres2pl_eigen,\n    Simulated=simres2pl_Q95\n) %>%\ndplyr::filter(as.integer(Component)<=9) %>%\ntidyr::pivot_longer(\n    cols=c(\"Observed\", \"Simulated\"),\n    names_to=\"Origin\", values_to=\"Eigenvalue\"\n) %>%\nggplot(aes(x=Component, y=Eigenvalue, group=Origin, colour=Origin)) +\n    geom_segment(x=-Inf, xend=Inf, y=1, yend=1, colour=\"orange\") +\n    geom_path() +\n    geom_point(size=6) +\n    scale_color_hue(direction=-1) +\n    scale_y_continuous(limits=c(0, 3.5)) +\n    theme_gray(base_size=18)","metadata":{"execution":{"iopub.status.busy":"2025-08-25T06:50:10.584908Z","iopub.execute_input":"2025-08-25T06:50:10.586394Z","iopub.status.idle":"2025-08-25T06:50:10.900907Z","shell.execute_reply":"2025-08-25T06:50:10.899763Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Comment\n\nThe first observed eigenvalue slightly exceeds the 95° percentile of the simulated ones. The elbow method gets suspicious about two components, especially because their eigenvalues exceed 3. The situation is borderline: further investigations are required.","metadata":{}}]}