---
title: "Detecting item misfit in Rasch models"
author:
  - name: Magnus Johansson
    affiliations:
      - name: RISE Research Institutes of Sweden, Division Built Environment, System Transition
      - name: Centre for Psychiatry Research, Department of Clinical Neuroscience, Karolinska Institutet, & Stockholm Health Care Services, Region Stockholm
    orcid: 0000-0003-1669-592X
    email: pgmj@pm.me
    corresponding: true
keywords:
  - Rasch
  - Psychometrics
  - Item fit
  - Cutoffs
  - Critical values
  - Model fit
abstract: |
  Psychometrics have long relied on rule-of-thumb critical values for goodness of fit metrics. With  powerful personal computers it is both feasible and desirable to use simulation methods to determine appropriate cutoff values. This paper evaluates the use of an R package for Rasch psychometrics that has implemented functions to simplify the process of determining simulation-based cutoff values. Through six simulation studies, comparisons are made between information-weighted conditional item fit ("infit") and item-restscore correlations using Goodman and Kruskal’s $\gamma$. Results indicate the limitations of small samples (n < 500) in correctly detecting item misfit, especially when a larger proportion of items are misfit and/or when misfit items are off-target. Infit with simulation-based cutoffs outperforms item-restscore with sample sizes below 500. Both methods result in problematic rates of false positives with large samples (n >= 1000). Large datasets should be analyzed using nonparametric bootstrap of subsamples with item-restscore to reduce the risk of type-1 errors. Finally, the importance of an iterative analysis process is emphasized, since a situation where several items show underfit will cause other items to show overfit. Underfit items should be removed one at a time, and a re-analysis conducted for each step to avoid erroneously eliminating items.
key-points:
  - Rule-of-thumb critical values for item fit metrics should be replaced with simulation-based values
  - Sample size affects the detection rate and false positive rate
  - More misfitting items increase the false positive rate
  - Mean squared unweighted "outfit" is not useful
date: 2025-03-11
bibliography: references.bib
shorttitle: Detecting item misfit in Rasch models
notebook-links: global
license: "CC BY"
funding: "The author received no specific funding for this work."
citation:
  type: article-journal
  container-title: "Educational Methods & Psychometrics"
  volume: 3
  issue: 18
  doi: "10.61186/emp.2025.5"
  issn: "2943-873X"
editor_options: 
  chunk_output_type: console
---

# Introduction

This paper presents a series of simulations conducted to evaluate methods to detect item misfit due to multidimensionality in Rasch models. First, conditional item infit and outfit [@muller_item_2020] will be under scrutiny. Second, item infit will be compared to the item-restscore method [@kreiner_note_2011;@christensen_item_2013]. Third, a bootstrap method for item-restscore will be presented and tested. This paper is intended for a target group of those who make practical use of Rasch analysis and wish to better understand the expected performance of methods available. As such, we refer readers interested in mathematical and statistical descriptions of the methods to referenced papers detailing this aspect. Only two simple performance metrics will be presented in the results: correct detection rate and false positive rate, both in percentages.

The evaluation of item fit under the Rasch model has, in the majority of published psychometric papers, been conducted using various more or less arbitrary rule-of-thumb critical values. Regarding mean squared (MSQ) item residuals, which should ideally be centered around 1.0, there are two sources often cited. One is the book by Bond and Fox [-@bond_applying_2015], which has garnered around 12 000 citations according to Google Scholar. It contains a table with rule-of-thumb recommendations for various settings, ranging from 0.8–1.2 to 0.5–1.7. Another frequently seen source, which is not an actual peer-reviewed publication and thus lacks citation counts, is the webpage at <https://rasch.org/rmt/rmt162f.htm>, where Mike Linacre states 0.5-1.5 to be "productive for measurement". Neither of these sources seem to rely on simulation studies to support their recommendations. While it is reasonable to accept a non-perfect fit to the Rasch model and also describe what one defines as acceptable levels of misfit, such recommendations would seem less arbitrary if related to simulations showing the range of item fit values found when simulating data that fit the Rasch model.

Müller [-@muller_item_2020] used simulation to show how the range of critical values for conditional item infit varies with sample size. The expected average conditional item infit range was described by Müller as fairly well captured by Smith's rule-of-thumb formula 1±2/$\sqrt{n}$ [@smith_using_1998], where n denotes the sample size. However, the average range does not apply for all items within a dataset, since item location relative to sample mean location also affects the expected model fit for individual items. This means that some items within a set of items varying in location are likely to have item fit values outside Smith's average value range while still fitting the Rasch model. Although primarily affecte by sample size, each item has its variations in the range of expected item fit.

While evaluation of item fit is an essential part of evaluating unidimensionality, it is recommended to use multiple methods. Standardized residuals are frequently analyzed, commonly with principal component analysis (PCA) and an analysis of residual correlations amongst item pairs, often referred to as Yen's Q3. Chou and Wang [-@chou_checking_2010] showed that the critical value for PCA of residuals to support unidimensionality suggested by Smith [-@smith_detecting_2002], using the largest eigenvalue < 1.5, is not generally applicable since it is affected by both test length and sample size. Christensen and colleagues [-@christensen_critical_2017] used simulation methods to illustrate the expected range of residual correlations under different conditions. Both of these papers provide important information about the dubiousness of using rule-of-thumb critical values when the empirical distribution of a statistic is not known, but they leave practitioners without tools to determine appropriate cutoffs to apply in practical analysis work.

It is here proposed that by using parametric bootstrapping one can establish item fit critical cutoff values that are relevant for a specific sample and item set. The procedure uses the estimated properties of the available data and simulates multiple new response datasets that fit the Rasch model to determine the range of plausible item fit values for each item. The R package `easyRasch` [@easyrasch] includes a function to determine item infit and outfit cutoff values using this method and will be tested in the simulation studies in this paper.

Similar developments, moving from rule-of-thumb towards adaptive critical values, have recently taken place in the related field of confirmatory factor analysis. McNeish and Wolf [-@mcneish_direct_2024] have created an R package called `dynamic` that uses simulation to determine appropriate critical values for commonly used model fit metrics for models using ordinal or interval data. 

It is important to note that the conditional item fit described by Müller [-@muller_item_2020] and implemented in the `iarm` R package [@mueller_iarm_2022] should not be confused with the unconditional item fit implemented in software such as Winsteps and RUMM2030, as well as all R packages except `iarm`. Unconditional item fit can result in unreliable item fit in sample sizes as small as 200 with an increasing probability of problems as sample size increases. Readers are strongly recommended to read Müller's paper to fully understand the issues with unconditional item fit. Additionally, the experienced Rasch analyst will perhaps wonder why the Wilson-Hilferty transformed Z statistic (often abbreviated ZSTD), which is based on unconditional MSQ is not included in this analysis. This is also explained in Müller's paper, where she describes both the notorious problems with sample size and shows that conditional item fit makes ZSTD superfluous. The `easyRasch` package, which is used in this paper, uses the `iarm` implementation of conditional item fit.

Currently, there are no published studies on the performance the item-restscore method, as described by Kreiner and Christensen [@kreiner_note_2011;@christensen_item_2013], in detecting misfitting items. Comparing it with an improved version of the long used item infit/outfit methods seemed like a good setting to evaluate item-restscore. The conditional likelihood ratio test [@andersen_goodness_1973] is included in Study 6, since it is a global test of fit that many are likely to be familiar with. As such, it also serves as a point of reference.

There are six simulation studies included in this paper:

1. Conditional item infit and outfit
2. Item-restscore
3. Comparing infit and item-restscore
4. Bootstrapped item-restscore
5. Varying the number of items
6. Conditional likelihood ratio test

# Methods 

```{r}
library(iarm)
library(eRm)
library(ggdist)
library(tidyverse)
library(easyRasch)
library(arrow)
library(showtext)

showtext_auto()

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
rename <- dplyr::rename

theme_rise <- theme_rise(fontfamily = "sans")
```

This is a general description of the methods used. Each study included in this paper has its own brief introduction and method section. A reproducible manuscript with R code and data is available on GitHub: <https://github.com/pgmj/rasch_itemfit>. First, a note on terminology. Non-parametric bootstrapping is synonymous with sampling with replacement. With this method, multiple random samples are drawn from the original response data. Parametric bootstrapping is synonymous with simulation since new data is generated based on person and item parameter estimates made from the original response data. This new data is then used for sampling with replacement.

The simulation of response data used three steps: First, a vector of theta values (person scores on the latent variable's logit scale) was generated using `rnorm(mean = 0, sd = 1.5)`. Second, a set of item locations ranging from -2 to 2 logits were generated for dichotomous items, using `runif(n = 20, min = -2, max = 2)`. The same set of item locations were used for all studies except Study 5, which adds 20 more item locations. Third, the theta values were used to simulate item responses for participants, using `sim.xdim()` from the `eRm` package [@mair_extended_2007], which allows the simulation of multidimensional response data. The sigma matrix used by `sim.xdim()` was specified to use 0.15 on the off-diagonal, where values are between 0 and 1 with lower values indicating stronger multidimensionality. Multiple datasets with 10 000 respondents each were generated using the same item and person parameters, varying the targeting of the misfitting item(s) and number of the misfitting item(s). More details are described in the separate studies. The parametric bootstrapping procedure was implemented using random samples from the simulated datasets. Sample size variations tested are also described in each study. 

The general procedure for the parametric bootstrap was as follows:

1. Estimation of item locations based on simulated item response data, using Conditional Maximum Likelihood [CML, @mair_extended_2007].
2. Estimation of sample theta values using weighted maximum likelihood [@warm_weighted_1989].
3. Simulation of new response data that fit the Rasch model, using the estimated item locations and theta values.
4. Estimation of the dichotomous Rasch model for the new response data using CML.
5. Based on step 4, calculation of conditional item infit and outfit [@muller_item_2020;@mueller_iarm_2022] and/or item-restscore metrics [@kreiner_note_2011;@mueller_iarm_2022].

Steps three and four were iterated over, using sampling with replacement from the estimated theta values as a basis for simulating the response data in step three. Summary statistics were created with a focus on the percentage of correct detection of misfit and false positives. A complete list of software used for the analyses is listed under @sec-addmat.

```{r}
# read pre-generated item locations for reproducibility
items1 <- read_csv("data/rm_items40.csv") %>%
  slice(1:20) %>% 
  pull(location)
```

```{r}
#| eval: false

# we want two dimensional data, lower off diagonal values -> stronger misfit
sigma <- matrix(c(1, 0.15, 0.15, 1), 2)
# a matrix to specify which dimension each item loads on
wmat0 <- matrix(nrow = 20,
               ncol = 2)
# set all items to load on dimension one
wmat0[1:20,1] <- 1
wmat0[1:20,2] <- 0
# item 9, with good targeting in this item set (closest to sample theta mean of 0), is chosen to belong to a second dimension
wmat0[9,1] <- 0
wmat0[9,2] <- 1

# a matrix to specify which dimension each item loads on
wmat1 <- matrix(nrow = 20,
               ncol = 2)
wmat1[1:20,1] <- 1
wmat1[1:20,2] <- 0
wmat1[18,1] <- 0
wmat1[18,2] <- 1

# a matrix to specify which dimension each item loads on
wmat2 <- matrix(nrow = 20,
               ncol = 2)
wmat2[1:20,1] <- 1
wmat2[1:20,2] <- 0
wmat2[13,1] <- 0
wmat2[13,2] <- 1

# generate dichotomous data
simdata0 <- eRm::sim.xdim(10000, items1, cutpoint = "randomized", Sigma = sigma, weightmat = wmat0)
simdata1 <- eRm::sim.xdim(10000, items1, cutpoint = "randomized", Sigma = sigma, weightmat = wmat1)
simdata2 <- eRm::sim.xdim(10000, items1, cutpoint = "randomized", Sigma = sigma, weightmat = wmat2)

simdata <- list(data0 = simdata0,
                data1 = simdata1,
                data2 = simdata2)
#saveRDS(simdata,"simdata10000.rds")

### And one dataset with all 3 items misfitting
# a matrix to specify which dimension each item loads on
wmat3 <- matrix(nrow = 20,
               ncol = 2)
# set all items to load on dimension one
wmat3[1:20,1] <- 1
wmat3[1:20,2] <- 0
# item 9, with good targeting in this item set (closest to sample theta mean of 0), is chosen to belong to a second dimension
wmat3[c(9,13,18),1] <- 0
wmat3[c(9,13,18),2] <- 1

# generate dichotomous data
simdata3 <- eRm::sim.xdim(10000, items1, cutpoint = "randomized", Sigma = sigma, weightmat = wmat3)

#saveRDS(simdata3,"data/simdata3.rds")
```

```{r}
#| eval: false
# simulate data using 40 items
# read pre-generated item locations for reproducibility
items2 <- read_csv("data/rm_items40.csv") %>%
  pull(location)

# a matrix to specify which dimension each item loads on
wmat4 <- matrix(nrow = 40,
               ncol = 2)
# set all items to load on dimension one
wmat4[1:40,1] <- 1
wmat4[1:40,2] <- 0

wmat4[c(9,13,18),1] <- 0
wmat4[c(9,13,18),2] <- 1

# generate dichotomous data
simdata40items <- eRm::sim.xdim(10000, items2, cutpoint = "randomized", Sigma = sigma, weightmat = wmat4)

# check targeting
#RItargeting(as.data.frame(simdata40items), model = "RM")
#RItargeting(simdata3[,1:10], model = "RM")

#saveRDS(simdata40items,"data/simdata40items.rds")
```

```{r}
simdata <- readRDS("data/simdata10000.rds")

simdata3 <- readRDS("data/simdata3.rds") %>% 
  as.data.frame()

demodata <- simdata[[1]] %>% 
  as.data.frame() %>% 
  slice(1:400) %>% 
  select(V1,V11,V3,V12)

simdata40items <- readRDS("data/simdata40items.rds") %>% 
  as.data.frame()
```

# Study 1: Item infit and outfit

Assessing item fit to the Rasch model using infit and outfit was first proposed by Wright and Panchapakesan [-@wright_procedure_1969]. For a historical perspective, see Smith and colleagues [-@smith_using_1998]. Item mean square standardized residuals are either unweighted, which is referred to as "outfit", or information weighted, also known as "infit" [@ostini_polytomous_2006, pp. 86-87]. Outfit is sensitive to outliers, while infit is much less affected by outliers. Both infit and outfit are based on individual response residuals. Conditional item infit and outfit are expected to be near 1, with higher values indicating an item to be underfitting the Rasch model (often due to multidimensionality issues) and lower values indicating overfit. For a more recent overview on item fit and its interpretation, see the book chapter by Christensen and Kreiner [-@christensen_item_2013].

Unconditional item fit uses person and item parameter estimates, which results in biased residuals [@kreiner_exact_2011]. By using the total person score instead of an estimated person score, this bias can be virtually eliminated. The Rasch model describes the probability of a correct or positive response by a person to an item as:

$$
P(X_{vi} = 1) = \frac{\exp(\theta_v - \beta_i)}{1 + \exp(\theta_v - \beta_i)}
$$ {#eq-rasch}

In @eq-rasch, $X$ denotes a response, persons are indexed by $v$, items by $i$, and $\theta$ and $\beta$ are person and item parameters, respectively. To compute the standardized Rasch model response residuals $Z_{vi}$, we use:

$$
Z_{vi} = \frac{X_{vi} - E(X_{vi})}{\sqrt{\text{Var}(X_{vi})}}
$$ {#eq-residuals}

where $E(X_{vi}) = P(X_{vi} = 1)$. The formula for outfit and infit is:

$$
\text{Outfit}_i = \sum\frac{Z_{vi}^2}{n}, \quad \text{Infit}_i = \frac{\sum_{v}
Z_{v_i}^2 \cdot w_{vi}}{\sum_{v}w_{vi}}
$$ {#eq-inoutfit}

where the weights to calculate infit, $w_{vi} = \text{Var}(X_{vi})$. Since $E(X_{vi})$ and $\text{Var}(X_{vi})$ need to be estimated, these are replaced in the conditional fit statistics, where:

$$
\hat{E}(X_{vi}) = \hat{P}(X_{vi} = 1 \mid R_v = r) = \frac{\exp \left(
-\hat{\beta}_i\right) \gamma_r{-1}\left(\hat{\beta}^\left({i}\right) \right) } {\gamma_r (\hat{\beta)}}
$$ {#eq-condfit}


In @eq-condfit, $\hat{\beta}$ indicates the vector of item parameters and $\gamma_r (\hat{\beta})$ is "the elementary symmetrical function of order r of the item parameter estimates" [@muller_item_2020, p. 4]. For further details on conditional item fit we refer to the previously mentioned paper by Müller [-@muller_item_2020] and the book chapter on item fit by Christensen and Kreiner [-@christensen_item_2013].

The function `RIgetfit()` from the `easyRasch` R package is tested here. Its source code can be accessed on GitHub, see @sec-addmat. The function offers the user a choice of the number of bootstrap iterations to use to determine the critical cutoff values for each item's infit and outfit. Our main interest in this study is two-fold. We want to test variations in the number of iterations used in `RIgetfit()` and evaluate how well the critical values based on the parametric bootstrap procedure detect misfitting items. Additionally, a comparison between infit and outfit statistics in terms of detection rate and false positive rate will be conducted.

20 dichotomous items are used, with one item misfitting. Item locations are the same throughout all studies unless otherwise noted. The location of the misfitting item relative to the sample theta mean was selected to be approximately 0, -1, and -2 logits. Three separate datasets were generated with these variations, each with 10 000 simulated respondents. One dataset with all three misfitting items was also generated, using the same sample size.

Then the `RIitemfit()` function is used to summarize the bootstrap results and also calculates the infit and outfit for each item in the observed data and highlights items with infit/outfit values outside of the cutoff values. `RIitemfit()` has a default (user-modifiable) setting to slightly truncate the distribution of values using `stats::quantile()` at 0.001 and 0.999 to remove extreme values. An example is demonstrated in @tbl-itemfit1, using a subset of the items used in the simulations. Note the variation between items in plausible threshold values for infit and outfit based on the bootstrap procedure compared to Smith's rule-of-thumb regarding infit (1±2/$\sqrt{n}$), which would be 0.9-1.1 for a sample size of 400. @fig-itemfit1 provides a visualization of the distribution of bootstrapped infit and outfit values, together with the infit/outfit values from the observed data illustrated using an orange diamond shape. 

This study was rather computationally demanding since each simulation run entailed 100-400 underlying bootstrap iterations. The sample sizes used were 150, 250, 500, and 1000. The number of iterations to determine cutoff values were 100, 200, and 400. Sample size and iteration conditions were fully crossed with each other and the three different targeting variations of the one misfitting item, resulting in 4 * 3 * 3 = 36 conditions. Each combination used 500 replications.

```{r}
#| label: tbl-itemfit1
#| tbl-cap: Conditional item fit with simulation-based cutoff values
#| cache: true

simfit <- RIgetfit(demodata, iterations = 200, cpu = 8)
RIitemfit(demodata, simfit, output = "quarto")
```

```{r}
#| label: fig-itemfit1
#| fig-cap: Distribution of simulation-based item fit and estimated item fit from observed data
RIgetfitPlot(simfit, demodata, output = "both")
```


```{r}
# define function to run simulations for item infit/outfit cutoff values
itemfitboot <- function(dat, iterations, samplesize) {
  
  fit <- list()
  fit <- foreach(i = 1:iterations) %do% {
    data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
      as.data.frame()
    
    # check data for responses in all cells
    n_resp <-
      data %>%
      as.matrix() %>%
      colSums2() %>%
      t() %>%
      as.vector()
    
    if (min(n_resp, na.rm = TRUE) < 11) {
      data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
        as.data.frame()
    } 
    
    # get simulation based cutoff values
    sfit100 <- RIgetfit(data,100,9)
    sfit200 <- RIgetfit(data,200,9)
    sfit400 <- RIgetfit(data,400,9)

    # apply cutoffs and store results
    rfit100 <- RIitemfit(data,sfit100, output = "dataframe") %>% 
      select(infit_msq,outfit_msq,infit_diff,outfit_diff) %>% 
      add_column(item = 1:ncol(data),
                 sims = 100,
                 iteration = i,
                 samplesize = samplesize)
    
    rfit200 <- RIitemfit(data,sfit200, output = "dataframe") %>% 
      select(infit_msq,outfit_msq,infit_diff,outfit_diff) %>% 
      add_column(item = 1:ncol(data),
                 sims = 200,
                 iteration = i,
                 samplesize = samplesize)
    
    rfit400 <- RIitemfit(data,sfit400, output = "dataframe") %>% 
      select(infit_msq,outfit_msq,infit_diff,outfit_diff) %>% 
      add_column(item = 1:ncol(data),
                 sims = 400,
                 iteration = i,
                 samplesize = samplesize)
    
    # combine output 
    fit <- rbind(rfit100,rfit200,rfit400)
  }
  return(fit)
}

```

```{r}
#| eval: false
samplesizes <- c(150,250,500,1000)

ifb0 <- list()
ifb1 <- list()
ifb2 <- list()

#library(tictoc)
#tic()
ifb0 <- map(samplesizes, ~ itemfitboot(simdata[[1]], iterations = 500, samplesize = .x))
#toc() # 14257.358 sec elapsed

ifb1 <- map(samplesizes, ~ itemfitboot(simdata[[2]], iterations = 500, samplesize = .x))
ifb2 <- map(samplesizes, ~ itemfitboot(simdata[[3]], iterations = 500, samplesize = .x))

# saveRDS(ifb0, "data/ifb0_200.rds")
# saveRDS(ifb2, "data/ifb2_200.rds")
# saveRDS(ifb1, "data/ifb1_200.rds")

ifb0_df <- map_dfr(1:4, ~ do.call("rbind", ifb0[[.x]])) %>% 
  add_column(targeting = 0)

ifb1_df <- map_dfr(1:4, ~ do.call("rbind", ifb1[[.x]])) %>% 
  add_column(targeting = 1)

ifb2_df <- map_dfr(1:4, ~ do.call("rbind", ifb2[[.x]])) %>% 
  add_column(targeting = 2)

ifb <- rbind(ifb0_df,ifb1_df,ifb2_df)

write_parquet(ifb,"data/ifb.parquet") 
# This code chunk was processed on three separate occasions (200+200+100 iterations) due to the time-consuming simulation.
# This resulted in three parquet-files with results, files named ifb_part2 and ifb_part3 in addition to the filename two lines above.
```


## Results

```{r}
ifb <- read_parquet("data/ifb.parquet")
ifb_part2 <- read_parquet("data/ifb_part2.parquet") %>% 
  mutate(iteration = iteration + 200)
ifb_part3 <- read_parquet("data/ifb_part3.parquet") %>% 
  mutate(iteration = iteration + 400)

ifb <- rbind(ifb,ifb_part2,ifb_part3)
```

Figures show the percent of simulation runs that have identified an item as misfitting. Items with more than 5% are colored in light red. A number representing the detection rate is shown adjacent to the bar representing the misfitting item. The figure grid columns are labeled with the number of iterations used by `RIgetfit()` to determine cutoff values and grid rows are labeled with the sample size.

Tables describing the results are available under Additional Materials, see @sec-tblinfitresults.

### Infit

```{r}
#| label: fig-ifb0
#| fig-cap: Conditional infit detection rate (misfit item at 0 logits)
#| fig-height: 5
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 0) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(sims = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         samplesize = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~sims) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting 0 logits (item 9 misfit). 500 simulations per combination.",
       #title = "Conditional infit detection rate"
       )
```

@fig-ifb0 shows the detection rate when the misfitting item is located at the sample mean. The detection rate is highest for the condition with 100 iterations with sample sizes 150 and 250, but it also shows higher levels of false positives when the sample size increases to 500 or more.

```{r}
#| label: fig-ifb1
#| fig-cap: Conditional infit detection rate (misfit item at -1 logits)
#| fig-height: 5
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 1) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(sims = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         samplesize = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = 1.2, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~sims) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting -1 logits (item 18 misfit). 500 simulations per combination.",
       #title = "Conditional infit detection rate"
       )
```

When the misfitting item is offset in targeting by -1 logits compared to the sample mean (see @fig-ifb1), the smallest sample size has less power to detect misfit compared to the on-target misfitting item. There are lower rates of false positives across all sample sizes and iterations.

```{r}
#| label: fig-ifb2
#| fig-cap: Conditional infit detection rate (misfit item at -2 logits)
#| fig-height: 5
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 2) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(sims = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         samplesize = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = 1.2, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~sims) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting -2 logits (item 13 misfit). 500 simulations per combination.",
       #title = "Conditional infit detection rate"
       )
```

Finally, when the misfitting item is located at -2 logits compared to the sample mean (see @fig-ifb2), we see a stronger reduction in power for sample sizes 150 and 250. No false positives are identified.

### Outfit

```{r}
#| label: fig-ifb0out
#| fig-cap: Conditional outfit detection rate (misfit item at 0 logits)
#| fig-height: 5
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!outfit_diff == "no misfit",
         targeting == 0) %>% 
  dplyr::count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
    ungroup() %>% 
  mutate(sims = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         samplesize = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = 0, vjust = -0.3) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~sims) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting 0 logits (item 9 misfit). 500 simulations per combination.",
       #title = "Conditional outfit detection rate"
       )
```

```{r}
#| label: fig-ifb1out
#| fig-cap: Conditional outfit detection rate (misfit item at -1 logits)
#| fig-height: 5
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!outfit_diff == "no misfit",
         targeting == 1) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
    ungroup() %>% 
  mutate(sims = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         samplesize = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = 1.2, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~sims) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting -1 logits (item 18 misfit). 500 simulations per combination.",
       #title = "Conditional outfit detection rate"
       )
```

```{r}
#| label: fig-ifb2out
#| fig-cap: Conditional outfit detection rate (misfit item at -2 logits)
#| fig-height: 5
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!outfit_diff == "no misfit",
         targeting == 2) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
    ungroup() %>% 
  mutate(sims = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         samplesize = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~sims) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting -2 logits (item 13 misfit). 500 simulations per combination.",
       #title = "Conditional outfit detection rate"
       )
```

As shown in @fig-ifb0out, @fig-ifb1out, and @fig-ifb2out, outfit is performing worse than infit across the board.

## Comments

Based on these simulation, it is highly recommended to use infit over outfit in assessing item fit. The performance of outfit calls to question whether it is useful at all for detecting item misfit.

Regarding infit and the use of parametric bootstrap with the function `RIgetfit()`, it looks like 100 iterations are to recommend to determine cutoff values when the sample size is 250 or lower, while 200 or 400 iterations reduce the risk for false positives at sample sizes of 500 or larger. False positives are found at sample sizes 500 and 1000 only. The risk for false positives is notably higher when the misfitting item is located at the sample mean compared to when the misfitting item is off-target by -1 logits or more.


# Study 2: Item-restscore

Item-restscore is a metric that compares an expected correlation with the observed correlation. The "restscore" is the ordinal sum score from all items except the item under scrutiny. Since item data are ordinal, the marginal association between item-score and restscore is assessed using Goodman and Kruskal’s $\gamma$ [@goodman_measures_1954;@kreiner_note_2011]. Lower observed values than expected indicate that an item is underfit to the Rasch model, while higher values indicate overfit. The item-restscore function used in this simulation is from the `iarm` package [@mueller_iarm_2022] and outputs Benjamini-Hochberg corrected *p*-values [@benjamini_controlling_1995], which are used to determine whether the differences between the observed and expected values are statistically significant, using *p* < .05 as critical value for each item. The data and procedure in this study follow the same structure as Study 1, with the addition of a smaller sample condition with 100 respondents.

```{r}
ir <- function(dat, iterations, samplesize, cpu = 9) {
  
  require(doParallel)
  registerDoParallel(cores = cpu)
  
  fit <- data.frame()
  fit <- foreach(i = 1:iterations, .combine = rbind) %dopar% {
    data <- dat[sample(1:nrow(dat), samplesize), ]
    
    erm_out <- RM(data)
    
    cfit <- out_infit(erm_out)
    cfit_d <- data.frame(infit = cfit$Infit, outfit = cfit$Outfit) %>% 
      round(3)
    
    i1 <- item_restscore(erm_out)
    i1 <- as.data.frame(i1)
    
    i1d <- data.frame("observed" = as.numeric(i1[[1]][1:ncol(data),1]),
                     "expected" = as.numeric(i1[[1]][1:ncol(data),2]),
                     "se" = as.numeric(i1[[1]][1:ncol(data),3]),
                     "p.value" = as.numeric(i1[[1]][1:ncol(data),4]),
                     "p.adj.BH" = as.numeric(i1[[1]][1:ncol(data),5])
    ) %>% 
      mutate(diff_abs = abs(expected - observed),
             diff = expected - observed,
             ir_padj = ifelse(p.adj.BH < .05, "sign. misfit","no misfit")) %>% 
      select(ir_padj, diff, diff_abs) %>% 
      mutate(item = 1:ncol(data)) %>% 
      add_column(iteration = i,
                 samplesize = samplesize)
    cbind(i1d,cfit_d)
  }
  return(fit)
}

```

```{r}
#| eval: false
samplesizes <- c(100,150,250,500,1000)

ir0 <- list()
ir1 <- list()
ir2 <- list()
#library(tictoc)
#tic()
ir0 <- map(samplesizes, ~ ir(simdata[[1]], iterations = 1000, samplesize = .x))
ir1 <- map(samplesizes, ~ ir(simdata[[2]], iterations = 1000, samplesize = .x))
ir2 <- map(samplesizes, ~ ir(simdata[[3]], iterations = 1000, samplesize = .x))
#toc()
# 483.373 sec elapsed

irall <- list(ir0,ir1,ir2)

#saveRDS(irall,"data/item_restscore1000.rds")
```

## Results
```{r}
#| label: fig-itemrestscore1
#| fig-cap: Item-restscore detection rate across targeting and sample size
#| fig-height: 6

ir_all <- readRDS("data/item_restscore1000.rds")
ir_results <- bind_rows(ir_all[[1]],ir_all[[2]],ir_all[[3]]) %>% 
  add_column(targeting = rep(c(0,1,2), each = 100000))

# using patchwork and separate plots to get geom_text working across targeting variations
results0 <- ir_results %>% 
  filter(targeting == 0) %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit" ) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000,
         targeting = factor(targeting, levels = c(0,1,2), labels = c('On target','-1 logit offset','-2 logits offset')),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "lightblue","lightpink"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting 0 logits") +
  theme(strip.text = element_blank())

results1 <- ir_results %>% 
  filter(targeting == 1) %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit" ) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000,
         targeting = factor(targeting, levels = c(0,1,2), labels = c('On target','1 logit offset','2 logits offset')),
         item = factor(item)) %>%  
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "lightblue","lightpink"))) +
  geom_text(data = . %>% filter(item == 18),
            aes(label = n/10), position = position_dodge(width = 0.9),
            hjust = 1.2, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting -1 logits") +
  theme(strip.text = element_blank())

results2 <- ir_results %>% 
  filter(targeting == 2) %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit" ) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000,
         targeting = factor(targeting, levels = c(0,1,2), labels = c('On target','1 logit offset','2 logits offset'))) %>% 
  mutate(samplesize = factor(samplesize, levels = c(100,150,250,500,1000),
                       labels = c("n = 100",
                                  "n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "lightblue","lightpink"))) +
  geom_text(data = . %>% filter(item == 13),
            aes(label = n/10), position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Targeting -2 logits")

results0 + results1 + results2 +
  plot_layout(axes = "collect",
              axis_titles = "collect") +
  plot_annotation(#title = "Item-restscore detection rate across targeting and sample size",
                  subtitle = "1000 simulated datasets for each combination.",
                  theme = theme_rise(fontfamily = "sans"))
```

This simulation includes an additional condition with 100 respondents, which results in significantly lower detection rates compared to n = 150. Compared to infit at 250 respondents, item-restscore has detection rates of 95.2%, 90.9%, and 62.4% for targeting 0, -1, and -2, while infit has 96.6%, 96.2%, and 70.6%. For sample sizes 500 and 1000, the detection rate is similar, including the increased tendency for false positives at n = 1000. The false positive rate is lower for item-restscore than infit for sample sizes below 1000.

Tables describing the results are available under Additional Materials, see @sec-tblirresults.

# Study 3: Comparing infit and item-restscore

We will now compare the performance of infit and item-restscore when all three items are misfitting at the same time. This simulation will also include a condition with 2000 respondents to examine if the false positive rate increases with more respondents. For infit, we will use 100 iterations with `RIgetfit()` for n < 500 and 200 iterations for n >= 500 since this produced the best results in Study 1. Outfit is also included in the study, mostly to check if its' performance is similar to the condition with only one misfitting item.

```{r}
# define function to run simulations for item infit/outfit cutoff values
itemfitboot2 <- function(dat, iterations, samplesize) {
  
  fit <- list()
  fit <- foreach(i = 1:iterations) %do% {
    data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
      as.data.frame()
    
    # check data for responses in all cells
    n_resp <-
      data %>%
      as.matrix() %>%
      colSums2() %>%
      t() %>%
      as.vector()
    
    if (min(n_resp, na.rm = TRUE) < 11) {
      data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
        as.data.frame()
    } 
    
    if (nrow(data < 400)) {
      sfit <- RIgetfit(data,100,9)
    } else if (nrow(data > 400)) {
      sfit <- RIgetfit(data,200,9)
    }
    
    # apply cutoffs and store results
    rfit <- RIitemfit(data, sfit, output = "dataframe") %>% 
      select(infit_msq,outfit_msq,infit_diff,outfit_diff) %>% 
      add_column(item = 1:ncol(data),
                 iteration = i,
                 samplesize = samplesize)

    # combine output 
    fit <- rfit
  }
  return(fit)
}

# using RIitemfit(cutoff = c(005,995))
itemfitboot2b <- function(dat, iterations, samplesize) {
  
  fit <- list()
  fit <- foreach(i = 1:iterations) %do% {
    data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
      as.data.frame()
    
    # check data for responses in all cells
    n_resp <-
      data %>%
      as.matrix() %>%
      colSums2() %>%
      t() %>%
      as.vector()
    
    if (min(n_resp, na.rm = TRUE) < 11) {
      data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
        as.data.frame()
    } 
    
    if (nrow(data < 400)) {
      sfit <- RIgetfit(data,100,9)
    } else if (nrow(data > 400)) {
      sfit <- RIgetfit(data,200,9)
    }
    
    # apply cutoffs and store results
    rfit <- RIitemfit(data, sfit, output = "dataframe", cutoff = c(.005,.995)) %>% 
      select(infit_msq,outfit_msq,infit_diff,outfit_diff) %>% 
      add_column(item = 1:ncol(data),
                 iteration = i,
                 samplesize = samplesize)

    # combine output 
    fit <- rfit
  }
  return(fit)
}

```


```{r}
#| eval: false
samplesizes <- c(150,250,500,1000,2000)
ifb3 <- list()
ir3 <- list()

ifb3 <- map(samplesizes, ~ itemfitboot2(simdata3, iterations = 500, samplesize = .x))
# I initally used 200 iterations for n = 150 and n = 250, which later seemed illogical, so the next line replaces those simulations
ifb3b <- map(c(150,250), ~ itemfitboot2(simdata3, iterations = 500, samplesize = .x))
#saveRDS(ifb3b,"data/ifb3b.rds")

ifb3c <- map(c(150,250), ~ itemfitboot2b(simdata3, iterations = 500, samplesize = .x))
#saveRDS(ifb3c,"data/ifb3c.rds")

ir3 <- map(samplesizes, ~ ir(simdata3, iterations = 500, samplesize = .x))

#saveRDS(ifb3,"data/ifb3.rds")
#saveRDS(ir3,"data/ir3.rds")
```

```{r}
ifb3 <- readRDS("data/ifb3.rds")
ifb3b <- readRDS("data/ifb3b.rds")
ifb3c <- readRDS("data/ifb3c.rds")

ir3 <- readRDS("data/ir3.rds")

ifb3_df <- bind_rows(ifb3)
ifb3b_df <- bind_rows(ifb3b)

# filter out n = 150-250 from the simulation with 200 iterations and replace with 100 iterations
ifb3_df <- ifb3_df %>% 
  filter(!samplesize %in% c(150,250)) %>% 
  select(!sims)

ifb3_df <- rbind(ifb3_df,ifb3b_df)

ir3_df <- bind_rows(ir3)
```

### Results

```{r}
#| label: fig-ifb3out
#| fig-cap: Conditional outfit detection rate with three misfitting items
#| fig-height: 7
ifb3_df %>% 
  group_by(samplesize) %>% 
  filter(!outfit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000")),
         item = factor(item)) %>% 

  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 0) +
    geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
    geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~.) +
  theme_rise(fontfamily = "sans") +
  labs(str_wrap("3 items misfitting. 200 iterations used to determine critial value. 500 simulations per combination.",68),
       #title = "Conditional outfit detection rate"
       )
```

```{r}
#| label: fig-ifb3
#| fig-cap: Conditional infit detection rate with three misfitting items
#| fig-height: 7
ifb3_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  dplyr::count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000")),
         item = factor(item)) %>% 

  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.6, vjust = 1) +
    geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.6, vjust = 1) +
    geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.6, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~.) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = str_wrap("3 items misfitting. 500 simulations per combination.",68),
       #title = "Conditional infit detection rate"
       )
```

Looking at the performance of infit with three misfitting items (@fig-ifb3), we can see that the detection rate is markedly worse for item 13 (targeting -2 logits) in sample sizes 500 and below, compared to when single items were misfitting. The false positive rate has increased for a sample size of 1000 and we can see it increase strongly at n = 2000. Outfit (@fig-ifb3out) again performs worse than infit.

```{r}
#| label: fig-ifb3b
#| fig-cap: Conditional infit detection rate with three misfitting items and increased truncation
#| fig-height: 4

ifb3c_df <- bind_rows(ifb3c)

ifb3c_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit",
         samplesize < 400) %>% 
  dplyr::count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000")),
         item = factor(item)) %>% 

  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
    geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
    geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.8, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize~.) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = str_wrap("3 items misfitting. 500 simulations per combination. Reduced truncation of bootstrap output values.",68),
       #title = "Conditional infit detection rate"
       )
```

In @fig-ifb3b, the minimal truncation used previously to remove extreme values (quantiles .001 and .999) was increased to .005 and .995. This improves the detection rate, particularly for the n = 250 condition and item 13, but also results in an increased false positive rate.

```{r}
#| label: fig-itemrestscore2
#| fig-cap: Item-restscore detection rate with three misfitting items
#| fig-height: 7

ir3_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "lightblue","lightpink"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
    geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
    geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.5, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item',  guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "3 items misfitting. 500 simulations per combination.",
       #title = "Item-restscore detection rate"
       )

```

```{r}
#| label: fig-comp1
#| fig-cap: Detection rate for item-restscore compared to infit
#| fig-height: 5

ifb3_sum <- ifb3_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(percent = n/500*100) %>% 
  ungroup() %>% 
  filter(item %in% c(9,13,18)) %>% 
  select(!c(n)) %>% 
  add_column(Method = "Infit")

ir3_sum <- ir3_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(percent = n*100/500) %>% 
  ungroup() %>% 
  filter(item %in% c(9,13,18)) %>% 
  select(!n) %>% 
  add_column(Method = "Item-restscore")

rbind(ifb3_sum,ir3_sum) %>% 
  mutate(item = factor(item, levels = c(9,18,13), labels = c("Item 9 (location 0)",
                                                          "Item 18 (location -1)",
                                                          "Item 13 (location -2)"),
                       ordered = T)) %>% 
    mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 

  ggplot(aes(x = item, y = percent, fill = Method)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0(percent,"%"), y = 8), color = "white",
            position = position_dodge(width = 0.9),
            size = 3) +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(y = "Detection rate in percent",
       #title = "Detection rate for item-restscore compared to infit"
       ) +
  scale_fill_brewer(type = "qual", palette = "Dark2") +
  theme(axis.title.x = element_blank())

```

Item-restscore (see @fig-itemrestscore2) shows a comparable detection rate to infit and higher levels of false positives. A comparison is made between the two in @fig-comp1, where item-restscore is performing better than infit at detecting the -2 logits off-target item at n = 250, and better across all items for n = 500 and n = 1000. Infit performs better for samples n = 150 and n = 250 (except for the item with location -2 logits).

Tables describing the results are available under Additional Materials, see @sec-infirresults.

```{r}
#| label: tbl-overunder
#| tbl-cap: Item-restscore summary results across all sample sizes
#| tbl-colwidths: [60,40]
ir3_df %>%
  filter(ir_padj == "sign. misfit") %>% 
  mutate(type = ifelse(diff < 0, "overfit","underfit")) %>% 
  count(item,type) %>% 
  mutate(Percent = n*100/2500) %>% 
  arrange(desc(Percent)) %>% 
  select(!n) %>% 
  set_names(c("Item","Type of misfit","Percent")) %>% 
  knitr::kable()
```

Reviewing the type of misfit identified by item-restscore (see @tbl-overunder), the false positives are all overfitting the Rasch model, except for two instances (out of 2500) indicating underfit for item 12. Items 9, 13, and 18, which were simulated to be misfitting due to loading on a separate dimension, are as expected showing underfit to the Rasch model.

# Study 4: Bootstrapped item-restscore

For this set of simulations, we will use a non-parametric bootstrap procedure with item-restscore. The difference from the parametric bootstrap is that the non-parametric bootstrap samples with replacement directly from the observed response data. First, based on the above problematic sample size of 2000 when three items are misfitting, we will use the bootstrap function to sample with replacement using n = 800 and 250 bootstrap samples. The function `RIbootRestscore()` from the `easyRasch` package will be used. 

Second, we will also apply the bootstrapped item-restscore method to sample sizes 150 and 250, using the complete sample for the same bootstrap procedure to see if this produces more useful information than previously tested strategies for identifying misfitting items.

## Results

```{r}
#| eval: false
ir_boot <- function(dat, iterations, samplesize) {
  
  fit <- list()
  fit <- foreach(i = 1:iterations) %do% {
    data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
      as.data.frame()
    
    fit <- RIbootRestscore(data, iterations = 250, samplesize = 800, cpu = 9, output = "dataframe", cutoff = 0) %>% 
      select(item,item_restscore,percent) %>% 
      filter(!item_restscore == "no misfit") %>% 
      add_column(iteration = i)
    
  }
  return(fit)
}

irb0 <- ir_boot(simdata3, 500, 2000)
saveRDS(irb0,"data/irb0.rds")
```

```{r}
#| label: fig-irb0all
#| fig-cap: Item-restscore bootstrap results (n = 2000)
#| fig-height: 5
irb0 <- readRDS("data/irb0.rds")

irb0_df <- bind_rows(irb0)

# irb0_df %>% 
#   mutate(item = factor(item, levels = paste0("V",20:1), ordered = T)) %>% 
#   ggplot(aes(x = percent, y = item, color = item_restscore)) +
#   stat_pointinterval(point_interval	= "median_hdci") +
#   theme_rise(fontfamily = "sans") +
#   labs(caption = str_wrap("Point indicates median value and lines are highest-density continuous interval for .66 and .95.",80),
#        x = "Distribution of percent values for item misfit", y = "Item",
#        title = "Item-restscore bootstrap results",
#        subtitle = str_wrap("250 bootstrap iterations with 800 respondents from a sample of 2000. 500 simulations were used.",70)
#        ) +
#   scale_color_discrete('Type', labels = c("Overfit","Underfit"))

irb0_df %>% 
  mutate(item = factor(item, levels = paste0("V",20:1), ordered = T)) %>% 
  group_by(item,item_restscore) %>% 
  summarise(median = median(percent),
            mean = mean(percent),
            q_lower = quantile(percent, .05),
            q_upper = quantile(percent, .95),
            sd = sd(percent),
            mad = mad(percent),
            se = sd(percent) / sqrt(500),
            ci_upper = median + (1.96*se),
            ci_lower = median - (1.96*se)) %>% 
  
  ggplot(aes(x = median, y = item, color = item_restscore)) +
  geom_point(size = 3.5, shape = 18) +
  geom_segment(aes(x = median - mad, xend = median + mad)) +
  theme_rise(fontfamily = "sans") +
  labs(caption = str_wrap("Point indicates median value and horizontal lines are median absolute deviation.",80),
       x = "Distribution of percent values for item misfit", y = "Item",
       #title = "Item-restscore bootstrap results",
       subtitle = str_wrap("250 bootstrap iterations with 800 respondents from a sample of 2000. 500 simulations were used.",70)
       ) +
  scale_color_discrete('Type', labels = c("Overfit","Underfit"))
```

```{r}
#| label: tbl-irb0mis
#| tbl-cap: Summary statistics for item-restscore bootstrap simulation
#| tbl-subcap: 
#|   - "Misfitting items"
#|   - "False positives"
#| layout-ncol: 1

irb0_df %>% 
  filter(item %in% c("V9","V13","V18")) %>% 
  rename(Item = item) %>% 
  group_by(Item) %>% 
  summarise(Median = median(percent),
            MAD = mad(percent),
            Mean = mean(percent),
            SD = sd(percent),
            `Percentile .05` = quantile(percent, .05)
            ) %>% 
  ungroup() %>% 
  arrange(desc(`Percentile .05`)) %>% 
  mutate_if(is.numeric, round, 1) %>% 
  knitr::kable()

irb0_df %>% 
  filter(!item %in% c("V9","V13","V18"),
         item_restscore == "overfit") %>% 
  rename(Item = item) %>% 
  group_by(Item) %>% 
  summarise(Median = median(percent),
            MAD = mad(percent),
            Mean = mean(percent),
            SD = sd(percent),
            `Percentile .95` = quantile(percent, .95)
            ) %>% 
  ungroup() %>% 
  arrange(desc(`Percentile .95`)) %>% 
  mutate_if(is.numeric, round, 1) %>% 
  knitr::kable()
```

@fig-irb0all shows that there is variation in the false positive rate and it is nearly always indicating overfit, while the misfitting items are only indicated as underfit. The summary statistics in @tbl-irb0mis show that there can be quite a bit of variation for false positives, but the clear majority of results are below 50%. 3 items have 95th percentile values above 50, with the highest at 58.8.

## Small sample (n = 150)

We will use 200 simulations to check the performance of the bootstrapped item-restscore function for sample size 150. As an additional experimental condition, we will use both 250 and 500 bootstrap iterations for item-restscore in each simulation.

```{r}
#| eval: false
ir_boot2 <- function(dat, iterations, samplesize) {
  
  fit <- list()
  fit <- foreach(i = 1:iterations) %do% {
    data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
      as.data.frame()
    
        # check data for responses in all cells
    n_resp <-
      data %>%
      as.matrix() %>%
      colSums2() %>%
      t() %>%
      as.vector()
    
    if (min(n_resp, na.rm = TRUE) < 11) {
      data <- dat[sample(1:nrow(dat), samplesize), ] %>% 
        as.data.frame()
    } 
    
    fit1 <- RIbootRestscore(data, iterations = 250, samplesize = nrow(data), cpu = 9, output = "dataframe", cutoff = 0) %>% 
      select(item,item_restscore,percent) %>% 
      filter(!item_restscore == "no misfit") %>% 
      add_column(iteration = i,
                 bootit = 250)
    
    fit2 <- RIbootRestscore(data, iterations = 500, samplesize = nrow(data), cpu = 9, output = "dataframe", cutoff = 0) %>% 
      select(item,item_restscore,percent) %>% 
      filter(!item_restscore == "no misfit") %>% 
      add_column(iteration = i,
                 bootit = 500)
    
    fit <- rbind(fit1,fit2)
    
  }
  return(fit)
}

irb150 <- ir_boot2(simdata3, 200, 150)
saveRDS(irb150,"data/irb150.rds")

# irb250 <- ir_boot2(simdata3, 200, 250)
# saveRDS(irb250,"data/irb250.rds")
```

```{r}
#| label: fig-irboot150
#| fig-cap: Item-restscore bootstrap results (n = 150)
#| fig-height: 6
irb150 <- readRDS("data/irb150.rds")

irb150_df <- bind_rows(irb150)

irb150_df %>% 
  mutate(item = factor(item, levels = paste0("V",20:1), ordered = T),
         bootit = factor(bootit, levels = c(250,500),
                         labels = c("250 bootstrap iterations",
                                    "500 bootstrap iterations"))) %>% 
  group_by(item,item_restscore,bootit) %>% 
  summarise(median = median(percent),
            mean = mean(percent),
            q_lower = quantile(percent, .05),
            q_upper = quantile(percent, .95),
            sd = sd(percent),
            mad = mad(percent),
            se = sd(percent) / sqrt(500),
            ci_upper = median + (1.96*se),
            ci_lower = median - (1.96*se)) %>% 
  
  ggplot(aes(x = median, y = item, color = item_restscore)) +
  geom_point(size = 3.5, shape = 18) +
  geom_segment(aes(x = median - mad, xend = median + mad)) +
    geom_text(data = . %>% filter(item == "V9",
                                  item_restscore == "underfit"),
            aes(label = median), color = "black",
            hjust = -0.1, vjust = -1) +
    geom_text(data = . %>% filter(item == "V13",
                                  item_restscore == "underfit"),
            aes(label = median), color = "black",
            hjust = -0.1, vjust = -1) +
    geom_text(data = . %>% filter(item == "V18",
                                  item_restscore == "underfit"),
            aes(label = median), color = "black",
            hjust = -0.1, vjust = -1) +
  theme_rise(fontfamily = "sans") +
  labs(caption = str_wrap("Point indicates median value and horizontal lines are median absolute deviation.",80),
       x = "Distribution of percent values for item misfit", y = "Item",
       #title = "Item-restscore bootstrap results",
       subtitle = str_wrap("250 bootstrap iterations with 800 respondents from a sample of 2000. 500 simulations were used.",70)
       ) +
  scale_color_discrete('Type', labels = c("Overfit","Underfit")) +
  facet_wrap(~bootit)

# irb150_df %>% 
#   mutate(bootit = factor(bootit, levels = c(250,500),
#                          labels = c("250 bootstrap iterations",
#                                     "500 bootstrap iterations"))) %>% 
#   ggplot(aes(x = percent, y = item, color = item_restscore)) +
#   stat_pointinterval(point_interval	= "median_hdci") +
#   theme_rise(fontfamily = "sans") +
#   labs(caption = str_wrap("Point indicates median value and lines are highest-density continuous interval for .66 and .95.",80),
#        x = "Distribution of percent values for item misfit", y = "Item",
#        title = "Item-restscore bootstrap results",
#        subtitle = str_wrap("250 and 500 bootstrap iterations with 150 respondents from a sample of 150. 200 simulations were used.",70)
#        ) +
#   scale_color_discrete('Type', labels = c("Overfit","Underfit")) +
#   facet_wrap(~bootit)
```

Item-restscore bootstrapping improves slightly on the single instance of item-restscore for the n = 150 condition (see @fig-irboot150). When comparing to the previous results in @fig-itemrestscore2, where the detection rate for the same sample size was at 49.2%, 14.6%, and 34.6% (for items 9, 13, and 18 respectively), the corresponding median values from the bootstrapped item-restscore with 250 iterations were 52.4%, 19.2%, and 38.4%. Using 500 bootstrap iterations did not result in relevant improvements over 250 iterations (see @tbl-irb150mis). Compared to the results using infit (@fig-ifb3), with detection rates of 59.2%, 19%, and 51.8%, item-restscore is inferior also when bootstrapped.

```{r}
#| label: tbl-irb150mis
#| tbl-cap: Summary statistics for item-restscore bootstrap simulation (n = 150)

irb150_tbl <- irb150_df %>% 
  filter(item %in% c("V9","V13","V18"),
         item_restscore == "underfit") %>% 
  rename(Item = item) %>% 
  group_by(Item,bootit) %>% 
  summarise(Median = median(percent),
            MAD = mad(percent),
            Mean = mean(percent),
            SD = sd(percent)
            ) %>% 
  ungroup() %>% 
  arrange(bootit,desc(Median)) %>% 
  mutate_if(is.numeric, round, 1) %>% 
  mutate(bootit = as.character(bootit)) %>% 
  relocate(bootit, .before = "Item") %>% 
  rename(`Bootstrap iterations` = bootit)  

irb150_tbl[c(1,3,4,6),1] <- ""

knitr::kable(irb150_tbl)
```


# Study 5: Varying the number of items 

For this simulation, we will vary the number of items and the number of misfitting items. First, 40 dichotomous items will be used, adding 20 new item locations to the previously used set, with the same three items misfitting (items 9, 13, and 18). Second, items 1-10 out of the initial 20 items will be used, which means only item 9 will be misfit. We'll again be using sample sizes of 150, 250, 500, and 1000.

Item-restscore and item infit will be compared. The latter will use 100 bootstrap iterations to determine critical values for sample sizes 150 and 250, and 200 bootstrap iterations for n >= 500.

```{r}
#| eval: false

samplesizes <- c(150,250,500,1000)
ifb10 <- list()
ir10 <- list()

#tic()
ifb40_60 <- map(samplesizes, ~ itemfitboot2(simdata40items, iterations = 60, samplesize = .x))
#toc() # 120s for 2 runs
saveRDS(ifb40_60, "data/ifb40_60i.rds")

ifb40_140 <- map(samplesizes, ~ itemfitboot2(simdata40items, iterations = 140, samplesize = .x))
#toc() # 120s for 2 runs
saveRDS(ifb40_140, "data/ifb40_140i.rds")

ifb10 <- map(samplesizes, ~ itemfitboot2(simdata3[,1:10], iterations = 200, samplesize = .x))

saveRDS(ifb10, "data/ifb10.rds")
```


```{r}
#| eval: false
ir40 <- map(samplesizes, ~ ir(simdata40items, iterations = 500, samplesize = .x))
saveRDS(ir40, "data/ir40.rds")

ir10 <- map(samplesizes, ~ ir(simdata3[,1:10], iterations = 500, samplesize = .x))
saveRDS(ir10, "data/ir10.rds")
```

```{r}
ir40 <- readRDS("data/ir40.rds")
ir10 <- readRDS("data/ir10.rds")

ifb40_60 <- readRDS("data/ifb40_60i.rds")
ifb40_140 <- readRDS("data/ifb40_140i.rds")
ifb10 <- readRDS("data/ifb10.rds")


ir40_df <- bind_rows(ir40)
ir10_df <- bind_rows(ir10)

ifb40_60_df <- bind_rows(ifb40_60)
ifb40_140_df <- bind_rows(ifb40_140)
ifb40_df <- rbind(ifb40_60_df,ifb40_140_df)

ifb10_df <- bind_rows(ifb10)

```

## Results 40 items

```{r}
#| label: fig-ifb40
#| fig-cap: Conditional infit detection rate with 40 items (3 misfit)
#| fig-height: 5
ifb40_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = round(n/200*100,1)) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 

  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
    geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
    geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_continuous('Item', limits = c(1,40), breaks = seq(1,40,1), guide = guide_axis(n.dodge = 2), minor_breaks = NULL) +
  guides(fill = "none") +
  facet_grid(samplesize~.) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "3 items misfitting. 200 simulations per combination.",
       #title = "Conditional infit detection rate"
       )
```


```{r}
#| label: fig-itemrestscore40
#| fig-cap: Item-restscore detection rate with 40 items (3 misfit)
#| fig-height: 5

ir40_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "lightblue","lightpink"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
    geom_text(data = . %>% filter(item == 13),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
    geom_text(data = . %>% filter(item == 18),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = -0.3, vjust = 1) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item', 
                     guide = guide_axis(n.dodge = 2)) +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "3 items misfitting. 500 simulations per combination.",
       #title = "Item-restscore detection rate"
       )

```

Infit performs better when the sample size is 150 or 250 (see @fig-ifb40), while performance is slightly better for item-restscore for n >= 500 in terms of lower rates of false positives (see @fig-itemrestscore40).

## Results 10 items

```{r}
#| label: fig-ifb10
#| fig-cap: Conditional infit detection rate (10 items)
#| fig-height: 5
ifb10_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = round(n/200*100,1)) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000"))) %>% 
  mutate(item = factor(item)) %>% 

  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), 
            position = position_dodge(width = 0.9),
            hjust = 0.5, vjust = 1.5) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete('Item') +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  coord_cartesian(clip = "off") +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Item 9 misfitting (targeting = 0). 200 simulations per combination.",
       #title = "Conditional infit detection rate"
       )
```

```{r}
#| label: fig-itemrestscore10
#| fig-cap: Item-restscore detection rate (10 items)
#| fig-height: 5

ir10_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  
  ggplot(aes(x = item, y = Percent)) +
  geom_col(aes(fill = ifelse(Percent > 5, "a","b"))) +
  geom_text(data = . %>% filter(item == 9),
            aes(label = Percent), position = position_dodge(width = 0.9),
            hjust = 0.5, vjust = 1.5) +
  scale_y_continuous(limits = c(0,105), breaks = seq(0,100,25)) +
  scale_x_discrete() +
  guides(fill = "none") +
  facet_grid(samplesize ~ .) +
  theme_rise(fontfamily = "sans") +
  labs(subtitle = "Item 9 misfitting (targeting = 0). 500 simulations per combination.",
       #title = "Item-restscore detection rate",
       x = "Item")
```


## Summary figure

```{r}
ifb20sum <- ifb3_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = round(n/500*100,1)) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000"))) %>% 
  mutate(item = factor(item)) %>% 
  select(!n) %>% 
  add_column(n_items = 20)

ifb40sum <- ifb40_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = round(n/200*100,1)) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000"))) %>% 
  mutate(item = factor(item)) %>% 
  select(!n) %>% 
  add_column(n_items = 40)

ifb_k <- rbind(ifb20sum,ifb40sum) %>% 
  add_column(Method = "Conditional infit")



ir20sum <- ir3_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  select(!n) %>% 
  add_column(n_items = 20)

ir40sum <- ir40_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  select(!n) %>% 
  add_column(n_items = 40)

ir_k <- rbind(ir20sum,ir40sum) %>% 
  add_column(Method = "Item-restscore")

dfk2040 <- rbind(ifb_k,ir_k)
```


```{r}
# 10 and 20 items separately with one misfitting item
ir10sum <- ir10_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000"))) %>% 
  select(!n) %>% 
  add_column(n_items = 10)

ir20sum2 <- ir_results %>% 
  filter(targeting == 0,
         !samplesize == 100) %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit" ) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  select(!c(targeting,n)) %>% 
    mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  add_column(n_items = 20)


ifb10sum <- ifb10_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = round(n/200*100,1)) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000"))) %>% 
  mutate(item = factor(item)) %>% 
  select(!n) %>% 
  add_column(n_items = 10)

ifb20_1 <- ifb %>% 
  filter(!infit_diff == "no misfit",
         targeting == 0) %>% 
  filter(samplesize %in% c(150,250) & sims == 100) %>% 
  group_by(samplesize) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000")),
         item = factor(item)) %>% 
  select(!n) %>% 
  add_column(n_items = 20)

ifb20_2 <- ifb %>% 
  filter(!infit_diff == "no misfit",
         targeting == 0) %>% 
  filter(samplesize %in% c(500,1000) & sims == 200) %>% 
  group_by(samplesize) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000")),
         item = factor(item)) %>% 
  select(!n) %>% 
  add_column(n_items = 20)

dfk10 <- rbind(ifb10sum,ifb20_1,ifb20_2) %>% 
  add_column(Method = "Conditional infit")

dfk10_2 <- rbind(ir10sum,ir20sum2) %>% 
        add_column(Method = "Item-restscore")


dfk10 <- rbind(dfk10,dfk10_2)
```

```{r}
#| label: fig-comp10
#| fig-cap: Detection rate for item-restscore and infit for 10 items
#| fig-height: 5
dfk10 %>% 
  mutate(item = as.numeric(item)) %>% 
  filter(item %in% c(9),
         !samplesize %in% c("n = 2000")) %>% 
  drop_na(samplesize) %>% 
  mutate(item = factor(item, levels = c(9,18,13), labels = c("Item 9 (location 0)",
                                                             "Item 18 (location -1)",
                                                             "Item 13 (location -2)"),
                       ordered = T),
         n_items = factor(n_items, levels = c(10,20), labels = c("10 items","20 items"))) %>% 

  ggplot(aes(x = item, y = Percent, fill = Method)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0(Percent,"%"), y = 8), color = "white",
            position = position_dodge(width = 0.9),
            size = 3.3) +
  facet_grid(samplesize ~ n_items) +
  theme_rise(fontfamily = "sans") +
  labs(y = "Detection rate in percent",
       #title = "Detection rate for item-restscore and infit across sample sizes",
       subtitle = "One item misfitting with location = 0 logits compared to sample mean.") +
  scale_fill_brewer(type = "qual", palette = "Dark2") +
  theme(axis.title.x = element_blank())
  
```

```{r}
#| label: fig-comp2040
#| fig-cap: Detection rate for item-restscore and infit for 20 and 40 items
#| fig-height: 5
#| fig-width: 8
dfk2040 %>% 
  mutate(item = as.numeric(item)) %>% 
  filter(item %in% c(9,18,13),
         !samplesize %in% c("n = 2000")) %>% 
  drop_na(samplesize) %>% 
  mutate(item = factor(item, levels = c(9,18,13), labels = c("Item 9 (location 0)",
                                                          "Item 18 (location -1)",
                                                          "Item 13 (location -2)"),
                       ordered = T),
         n_items = factor(n_items, levels = c(20,40), labels = c("20 items","40 items"))) %>% 
  mutate(Method = factor(Method, labels = c("Conditional\ninfit","Item-\nrestscore"))) %>% 

  ggplot(aes(x = item, y = Percent, fill = Method)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0(Percent,"%"), y = 8), color = "white",
            position = position_dodge(width = 0.9),
            size = 2.7) +
  facet_grid(samplesize ~ n_items) +
  theme_rise(fontfamily = "sans") +
  labs(y = "Detection rate in percent",
       #title = "Detection rate for item-restscore and infit across sample sizes and number of items",
       subtitle = "Three items misfitting.") +
  scale_fill_brewer(type = "qual", palette = "Dark2") +
  scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  theme(axis.title.x = element_blank())
  
```

@fig-comp10 and @fig-comp2040 summarize the findings from the two different comparisons of item numbers.  Adding more items improves the detection rate substantially for both methods, particularly for smaller samples and the off-target items. 40 items compared to 20 items result in a larger improvement for infit over item-restscore for the n = 150 condition, but also the n = 250 and n = 500 conditions for the -2 logits off-target item.

# Study 6: Conditional likelihood ratio test

While global tests of model fit don't provide any information about reasons for misfit when detected, they can be useful together with more specific tests such as those demonstrated in this paper. A commonly used global goodness-of-fit test is the Likelihood Ratio Test [LRT, @andersen_goodness_1973], which is also implemented in the `eRm` package. While the global test is not central to the purposes of this paper (detecting item misfit), it could provide readers with a familiar test as a reference and perhaps aid the item-focused methods in determining misfit. Previous simulation studies evaluating the LRT have found it to not be sensitive to detect multidimensionality [@debelak_evaluation_2019].

For this simulation, sample sizes were set to 150, 250, 500, 1000, and 2000, using 1000 simulations for each condition. Study 6 follows the previous procedures. First, the three datasets with one misfitting item out of 20, varying the misfitting items location (0, -1, and -2 logits). Then the datasets with three misfitting items, and comparing 20 and 40 items. And finally, one dataset with 10 items and one well-targeted misfitting item.

```{r}
lrtsim <- function(dat, iterations, samplesize, cpu = 9) {
  
  require(doParallel)
  registerDoParallel(cores = cpu)
  
  fit <- data.frame()
  fit <- foreach(i = 1:iterations, .combine = rbind) %dopar% {
    data <- dat[sample(1:nrow(dat), samplesize), ]
    
    erm_out <- RM(data)
    lrt <- LRtest(erm_out)
    data.frame(p = lrt[["pvalue"]],
               n = samplesize)
    # same p-values with clr_tests()
    #clr <- clr_tests(dat.items = data, model = "RM")
    #clr[,3]
    
  }
  return(fit)
}
```

```{r}
#| eval: false
samplesizes <- c(150,250,500,1000,2000)

# single misfit item
lrt0 <- map(samplesizes, ~ lrtsim(simdata[[1]], iterations = 1000, samplesize = .x))
lrt1 <- map(samplesizes, ~ lrtsim(simdata[[1]], iterations = 1000, samplesize = .x))
lrt2 <- map(samplesizes, ~ lrtsim(simdata[[1]], iterations = 1000, samplesize = .x))

# one misfit amongst 10
lrt10 <- map(samplesizes, ~ lrtsim(simdata3[,1:10], iterations = 1000, samplesize = .x))
# 3 misfit items
lrt20 <- map(samplesizes, ~ lrtsim(simdata3, iterations = 1000, samplesize = .x))
lrt40 <- map(samplesizes, ~ lrtsim(simdata40items, iterations = 1000, samplesize = .x))

saveRDS(lrt0,"data/lrt0.rds")
saveRDS(lrt1,"data/lrt1.rds")
saveRDS(lrt2,"data/lrt2.rds")
saveRDS(lrt10,"data/lrt10.rds")
saveRDS(lrt20,"data/lrt20.rds")
saveRDS(lrt40,"data/lrt40.rds")
```

## Results

```{r}
lrt0_df <- readRDS("data/lrt0.rds") %>% 
  bind_rows() %>% 
  select(!i) %>% 
  rename(samplesize = n) %>% 
  add_column(targeting = 0)
lrt1_df <- readRDS("data/lrt1.rds") %>% 
  bind_rows() %>% 
  select(!i) %>% 
  rename(samplesize = n) %>% 
  add_column(targeting = 1)
lrt2_df <- readRDS("data/lrt2.rds") %>% 
  bind_rows() %>% 
  select(!i) %>% 
  rename(samplesize = n) %>% 
  add_column(targeting = 2)

lrt10_df <- readRDS("data/lrt10.rds") %>% 
  bind_rows() %>% 
  select(!i) %>% 
  rename(samplesize = n) %>% 
  add_column(n_items = 10)
lrt20_df <- readRDS("data/lrt20.rds") %>% 
  bind_rows() %>% 
  select(!i) %>% 
  rename(samplesize = n) %>% 
  add_column(n_items = 20)
lrt40_df <- readRDS("data/lrt40.rds") %>% 
  bind_rows() %>% 
  select(!i) %>% 
  rename(samplesize = n) %>% 
  add_column(n_items = 40)

lrt_targ <- rbind(lrt0_df,lrt1_df,lrt2_df)
lrt_n <- rbind(lrt10_df,lrt20_df,lrt40_df)
```


```{r}
#| label: fig-lrt1
#| fig-cap: Likelihood ratio test detection rate
#| fig-subcap: 
#|   - "Across sample sizes and location of misfit item"
#|   - "Across sample sizes and number of items"
#| layout-ncol: 1
lrt_targ %>% 
  mutate(p_corr = ifelse(p < .05, "correct","false")) %>% 
  group_by(samplesize,targeting) %>% 
  dplyr::count(p_corr) %>% 
  mutate(Percent = n*100/sum(n) %>% round(1)) %>% 
  ungroup() %>% 
  filter(p_corr == "correct") %>% 
  mutate(Targeting = factor(targeting, levels = c(0,1,2), labels = c('On target','1 logit offset','2 logits offset')),
         samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000",
                                  "n = 2000"))) %>% 
  
  ggplot(aes(x = samplesize, y = Percent, fill = Targeting)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0(Percent,"%")), 
            color = "black",
            position = position_dodge(width = 0.9),
            hjust = 1.2, vjust = 0.5,
            angle = 90
            ) +
  theme_rise(fontfamily = "sans") +
  scale_fill_brewer(palette = "Greens") +
  labs(y = "Detection rate in percent",
       #title = "Detection rate for LRT across sample sizes and location of misfit item",
       subtitle = "One item out of 20 with misfit. 1000 simulations.",
       x = "Sample size")


lrt_n %>% 
  mutate(p_corr = ifelse(p < .05, "correct","false")) %>% 
  group_by(samplesize,n_items) %>% 
  dplyr::count(p_corr) %>% 
  mutate(Percent = n*100/sum(n) %>% round(1)) %>% 
  ungroup() %>% 
  filter(p_corr == "correct") %>% 
  mutate(n_items = factor(n_items, levels = c(10,20,40), labels = c('10 items (1 misfit)','20 items (3 misfit)','40 items (3 misfit)')),
         samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000",
                                  "n = 2000"))) %>% 
  
  ggplot(aes(x = samplesize, y = Percent, fill = n_items)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0(Percent,"%")), 
            color = "black",
            position = position_dodge(width = 0.9),
            hjust = 1.2, vjust = 0.5,
            angle = 90
            ) +
  theme_rise(fontfamily = "sans") +
  scale_fill_brewer('Number of items') +
  labs(y = "Detection rate in percent",
       #title = "Detection rate for LRT across sample sizes and number of items",
       subtitle = "1000 simulations.",
       x = "Sample size")
```

Results for targeting of the misfit item are summarized in @fig-lrt1, where the detection rate is reported based on the proportion of p-values below .05 in each condition.

LRT performs better than item-restscore for the off-target item conditions, especially when targeting was at -2 logits. Otherwise, performance is similar. Infit is better than LRT for n < 500 when targeting = 0, similar at targeting = 1, and LRT is better than infit at targeting -2.

Looking at LRT for different numbers of items and comparing to results for infit and item-restscore at their best detection rates for the corresponding number of items and sample size, it shows much worse performance for 10 items compared to infit and item-restscore at n = 150, but similar performance at larger sample sizes. For 20 and 40 items, LRT is similar to infit at all sample sizes, and slightly better than item-restscore for n < 500.

It should be noted that LRT also is sensitive to larger sample sizes and increased numbers of items, see <https://pgmj.github.io/clrt.html> for a brief example, and is likely to have a higher error rate than the nominal 5% under many conditions.

# Discussion

This paper was created primarily out of a desire to understand the performance of conditional item fit and item-restscore in detecting item misfit. Studies 1-3 and 5 were originally planned, while study 4 was added due to the results of previous studies showing issues with large sample sizes, and study 6 was added mostly for didactic purposes but also out of curiosity of the global fit LRT performance compared to the other methods. Müller's [-@muller_item_2020] paper on conditional item fit was published five years ago with results that should have sparked discussions in the Rasch community about updating methods (and software) and the justification of rule-of-thumb critical values. We hope that this paper can help spur such discussions by the simple methodology and presentation of results used, which could make this paper more accessible for practitioners.

Main results summarized: 

- small samples (n < 250-500) with small numbers of items make it harder to detect misfitting items
- small samples (n < 250-500) should rely primarily on conditional item infit with simulation-based critical values
- large samples (n > 500) should use bootstrapped item-restscore to reduce the risk of falsely identifying misfit
- misfit in an off-target item is harder to detect than well-targeted items and the LRT test can assist in identifying misfit in off-target items
- increased number of items increases the power to detect misfit

Assessment of item fit and dimensionality should always be done using multiple methods. Study 6 showed the benefits of also looking at the likelihood ratio test when the misfitting item is located -2 logits away from the sample mean location. Previous simulation studies concluded that the LRT did not perform well at detecting multidimensionality, which is how the misfitting items were generated in this study. While no method in this study achieved very good results at the lowest sample size (n = 150), LRT did not do much worse than the other methods.

```{r}
#| label: fig-loadloc
#| fig-cap: Standardized loadings on first residual contrast factor and item location
#| fig-height: 5
RIloadLoc(simdata3[1:400,], model = "RM")
```

Item fit and item-restscore are recommended to be used in parallel, while also examining residual patterns by reviewing standardized factor loadings on the first residual contrast as well as Yen's Q3 residual correlations [@christensen_critical_2017]. @fig-loadloc illustrates an example, using n = 400 and 3/20 items misfitting, where it is clearly seen how the misfitting items form a separate cluster on the x axis. Regarding residual correlations and critical values, the `easyRasch` package also contains a function to use bootstrapping (`RIgetResidCor()`), similar to the `RIgetfit()` function, to determine the appropriate cutoff. This is described briefly in a blog post [@johansson_simulation_2024] and a simulation paper is under preparation.

The findings reported here also make a good argument for removing one item at a time when the analysis indicates several misfitting items, starting with the most underfitting item. Such practice would reduce the risk for type-1 errors. This is especially relevant for n >= 500 and when truly misfitting items are located close to the sample mean, when the risk for false positive misfit indication is elevated. In a situation where the analysis of residual patterns using PCA and Yen's Q3 clearly shows clustering of items, it can be justified to remove several items simultaneously.

Item fit in this paper has been assessed using data from individual respondent residuals in a sample. A useful additional method to analyze and understand item fit and misfit is to inspect item characteristic curves where the sample is divided into class intervals based on their total score and the group residuals are inspected [@buchardt_visualizing_2023].

While the simulations in this paper have all used dichotomous data, all functions evaluated from the `easyRasch` package also work with polytomous data using the Rasch Partial Credit Model [PCM, @masters_rasch_1982].

## Limitations

When doing simulation studies there is always a balance to strike between trying to evaluate many scenarios and not having too high complexity. We have been keeping several things constant, such as item locations and number of items, which makes interpretation easier but may limit the applicability of the results.

The total number of items and the proportion of misfit items clearly have effects on detection rate and could have been investigated further using more variations in sample sizes. The Rasch partial credit model for polytomous data would have been useful to include in a comparison study. When testing the bootstrapped item-restscore method, more variations in the number of bootstrap samples (iterations) might have been of interest, although the difference between 250 and 500 was small.

# Conclusion

For sample sizes under 500, it seems best to rely mostly on item infit with simulation-based critical values, using 100 iterations with `RIgetfit()`. For sample sizes closer to 500, item-restscore is recommended as the primary method, either as a single-run test or bootstrapped. With samples larger than 500, bootstrapped item-restscore controls false positive rates well, while showing high rates of misfit detection. Using 250 iterations for the bootstrapped item-restscore seems adequate. In general, both infit and item-restscore are useful in the analysis process if you have a sample size below 1000.

{{< pagebreak >}}


# References {.unnumbered}

:::{#refs}

:::

{{< pagebreak >}}

# Additional materials {#sec-addmat}

- GitHub link for `easyRasch` source code: <https://github.com/pgmj/easyRasch/>
  - Most functions are defined in this file: <https://github.com/pgmj/easyRasch/blob/main/R/easyRasch.R>
  

## Results tables

```{r}
# library(gt)
# 
# # APA table styling from <https://gist.github.com/pdparker/1b61b6d36d09cb295bf286a931990159>
# apa_style <- function(data) {
#   data %>%
#     opt_table_lines(extent = "none") %>%
#     tab_options(
#       heading.border.bottom.width = 2,
#       heading.border.bottom.color = "black",
#       heading.border.bottom.style = "solid",
#       table.border.top.color = "white",
#       table_body.hlines.color = "white",
#       table_body.border.top.color = "black",
#       table_body.border.top.style = "solid",
#       table_body.border.top.width = 1,
#       heading.title.font.size = 12,
#       table.font.size = 12,
#       heading.subtitle.font.size = 12,
#       table_body.border.bottom.color = "black",
#       table_body.border.bottom.width = 1,
#       table_body.border.bottom.style = "solid",
#       column_labels.border.bottom.color = "black",
#       column_labels.border.bottom.style = "solid",
#       column_labels.border.bottom.width = 1
#     ) %>%
#       opt_table_font(font = "times")
# }
# library(flextable)
# set_flextable_defaults(digits = 1)
options(knitr.kable.NA = "0")
```

### Study 1: Item infit and outfit {#sec-tblinfitresults}

#### 100 iterations

```{r}
#| label: tbl-ifbresults1000
#| tbl-cap: Conditional infit detection rate, 100 bootstrap iterations, targeting = 0 logit

ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 0
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset'))
         ) %>% 
  mutate(Targeting = gsub("<","=",Targeting)) %>% 
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "100 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
```


```{r}
#| label: tbl-ifbresults1001
#| tbl-cap: Conditional infit detection rate, 100 bootstrap iterations, targeting = -1 logit

ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 1
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset'))
         ) %>% 
  mutate(Targeting = gsub("<","=",Targeting)) %>% 
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "100 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
```


```{r}
#| label: tbl-ifbresults1002
#| tbl-cap: Conditional infit detection rate, 100 bootstrap iterations, targeting = -2 logits
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 2
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset'))
         ) %>% 
  mutate(Targeting = gsub("<","=",Targeting)) %>% 
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "100 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
```

#### 200 iterations

```{r}
#| label: tbl-ifbresults2000
#| tbl-cap: Conditional infit detection rate, 200 bootstrap iterations, targeting = 0 logit
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 0
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset'))
         ) %>% 
  mutate(Targeting = gsub("<","=",Targeting),
         `Percent misfit` = car::recode(`Percent misfit`,"NA=0")) %>% 
  
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "200 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
  # #group_by(Targeting) %>% 
  # as_grouped_data("Targeting") %>% 
  # #arrange(Item) %>% 
  # flextable() %>%
  # autofit() %>% 
  # theme_apa() %>% 
  # colformat_double()
```

```{r}
#| label: tbl-ifbresults2001
#| tbl-cap: Conditional infit detection rate, 200 bootstrap iterations, targeting = -1 logit
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 1
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset'))
         ) %>% 
  mutate(Targeting = gsub("<","=",Targeting),
         `Percent misfit` = car::recode(`Percent misfit`,"NA=0")) %>% 
  
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "200 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
  # #group_by(Targeting) %>% 
  # as_grouped_data("Targeting") %>% 
  # #arrange(Item) %>% 
  # flextable() %>%
  # autofit() %>% 
  # theme_apa() %>% 
  # colformat_double()
```

```{r}
#| label: tbl-ifbresults2002
#| tbl-cap: Conditional infit detection rate, 200 bootstrap iterations, targeting = -2 logits
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 2
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset'))
         ) %>% 
  mutate(Targeting = gsub("<","=",Targeting),
         `Percent misfit` = car::recode(`Percent misfit`,"NA=0")) %>% 
  
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "200 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
  # #group_by(Targeting) %>% 
  # as_grouped_data("Targeting") %>% 
  # #arrange(Item) %>% 
  # flextable() %>%
  # autofit() %>% 
  # theme_apa() %>% 
  # colformat_double()
```

#### 400 iterations

```{r}
#| label: tbl-ifbresults4000
#| tbl-cap: Conditional infit detection rate, 400 bootstrap iterations, targeting = 0 logit
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 0
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = car::recode(targeting,"0='Misfit item 0 logit offset';1='Misfit item -1 logit offset';2='Misfit item -2 logits offset'", as.factor = TRUE)
         ) %>% 
  #mutate(Targeting = gsub("<","=",Targeting)) %>% 
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "400 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  #group_by(Targeting) %>% 
  #as_grouped_data("Targeting") %>% 
  arrange(Item) %>% 
  kable()
  # flextable() %>%
  # autofit() %>% 
  # theme_apa() %>% 
  # colformat_double() 
```

```{r}
#| label: tbl-ifbresults4001
#| tbl-cap: Conditional infit detection rate, 400 bootstrap iterations, targeting = -1 logit
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 1
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = car::recode(targeting,"0='Misfit item 0 logit offset';1='Misfit item -1 logit offset';2='Misfit item -2 logits offset'", as.factor = TRUE)
         ) %>% 
  #mutate(Targeting = gsub("<","=",Targeting)) %>% 
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "400 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
```

```{r}
#| label: tbl-ifbresults4002
#| tbl-cap: Conditional infit detection rate, 400 bootstrap iterations, targeting = -2 logits
ifb %>% 
  group_by(targeting, samplesize, sims) %>% 
  filter(!infit_diff == "no misfit",
         targeting == 2
         ) %>% 
  count(item, .drop = F) %>% 
  mutate(`Percent misfit` = n/500*100) %>% 
  ungroup() %>% 
  mutate(`Bootstrap iterations` = factor(sims, levels = c(100,200,400),
                       labels = c("100 iterations",
                                  "200 iterations",
                                  "400 iterations")),
         `Sample size` = factor(samplesize, levels = c(150,250,500,1000),
                       labels = c("n = 150",
                                  "n = 250",
                                  "n = 500",
                                  "n = 1000")),
         Item = factor(item, levels = 1:20, ordered = T),
         Targeting = car::recode(targeting,"0='Misfit item 0 logit offset';1='Misfit item -1 logit offset';2='Misfit item -2 logits offset'", as.factor = TRUE)
         ) %>% 
  #mutate(Targeting = gsub("<","=",Targeting)) %>% 
  select(`Sample size`,`Bootstrap iterations`,Item,`Percent misfit`) %>% 
  filter(`Bootstrap iterations` == "400 iterations") %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
```

### Study 2: Item-restscore {#sec-tblirresults}

```{r}
#| label: tbl-irresults0
#| tbl-cap: Item-restscore detection rate, targeting = 0 logit
ir_results %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit",
         targeting == 0) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000,
         targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset')),
         item = factor(item),
         samplesize = car::recode(samplesize,"100='n < 100';150='n < 150';250='n < 250';500='n < 500';1000='n < 1000'")) %>% 
  mutate(samplesize = gsub("<","=",samplesize)) %>% 
  rename(Targeting = targeting,
         `Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()

```

```{r}
#| label: tbl-irresults1
#| tbl-cap: Item-restscore detection rate, targeting = -1 logit
ir_results %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit",
         targeting == 1) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000,
         targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset')),
         item = factor(item),
         samplesize = car::recode(samplesize,"100='n < 100';150='n < 150';250='n < 250';500='n < 500';1000='n < 1000'")) %>% 
  mutate(samplesize = gsub("<","=",samplesize)) %>% 
  rename(Targeting = targeting,
         `Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()

```

```{r}
#| label: tbl-irresults2
#| tbl-cap: Item-restscore detection rate, targeting = -2 logits
ir_results %>% 
  group_by(targeting, samplesize) %>% 
  filter(ir_padj == "sign. misfit",
         targeting == 2) %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/1000,
         targeting = factor(targeting, levels = c(0,1,2), labels = c('Misfit item 0 logit offset','Misfit item -1 logit offset','Misfit item -2 logits offset')),
         item = factor(item),
         samplesize = car::recode(samplesize,"100='n < 100';150='n < 150';250='n < 250';500='n < 500';1000='n < 1000'")) %>% 
  mutate(samplesize = gsub("<","=",samplesize)) %>% 
  rename(Targeting = targeting,
         `Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()

```

### Study 3: Comparing infit and item-restscore {#sec-infirresults}

```{r}
#| label: tbl-ifb3out
#| tbl-cap: Conditional outfit detection rate with three misfitting items
ifb3_df %>% 
  group_by(samplesize) %>% 
  filter(!outfit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000")),
         item = factor(item)) %>% 
  rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = Item, names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  arrange(Item) %>% 
  # gt() %>% 
  # sub_missing(missing_text = "0") %>% 
  # tab_footnote(footnote = "Note. Misfitting items are 9, 13, and 18.") 
  kable()
```

```{r}
#| label: tbl-ifb3
#| tbl-cap: Conditional infit detection rate with three misfitting items
ifb3_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  dplyr::count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000")),
         item = factor(item)) %>% 
    rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = Item, names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  arrange(Item) %>% 
  # gt() %>% 
  # sub_missing(missing_text = "0") %>% 
  # tab_footnote(footnote = "Note. Misfitting items are 9, 13, and 18.") %>% 
  # apa_style()
  kable()
```

```{r}
#| label: tbl-ifb3b
#| tbl-cap: Conditional infit detection rate with three misfitting items and increased truncation

ifb3c_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit",
         samplesize < 400) %>% 
  dplyr::count(item, .drop = F) %>% 
  mutate(Percent = n/500*100) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000")),
         item = factor(item)) %>% 
  rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = Item, names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  arrange(Item) %>% 
  # gt() %>% 
  # sub_missing(missing_text = "0") %>% 
  # tab_footnote(footnote = "Note. Misfitting items are 9, 13, and 18.") %>% 
  # apa_style()
  kable()
```

```{r}
#| label: tbl-itemrestscore2
#| tbl-cap: Item-restscore detection rate with three misfitting items

ir3_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
  mutate(item = factor(item)) %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
    rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:20, ordered = T)) %>% 
  arrange(Item) %>% 
  # gt() %>% 
  # sub_missing(missing_text = "0") %>% 
  # tab_footnote(footnote = "Note. Misfitting items are 9, 13, and 18.") %>% 
  # apa_style()
  kable()
```

```{r}
#| label: tbl-comp1
#| tbl-cap: Detection rate for item-restscore compared to infit

rbind(ifb3_sum,ir3_sum) %>% 
  mutate(item = factor(item, levels = c(9,18,13), labels = c("Item 9 (location 0)",
                                                             " Item 18 (location -1)",
                                                             " Item 13 (location -2)"),
                       ordered = T)) %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = percent) %>% 
  pivot_wider(id_cols = c(Method,Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  arrange(Item) %>% 
  # gt() %>% 
  # sub_missing(missing_text = "") %>% 
  # apa_style()
  kable()
```

### Study 5: Varying the number of items 

```{r}
#| label: tbl-ifb40
#| tbl-cap: Conditional infit detection rate with 40 items (3 misfit)
ifb40_df %>% 
  group_by(samplesize) %>% 
  filter(!infit_diff == "no misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = round(n/200*100,1)) %>% 
  ungroup() %>% 
  mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:40, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
```

```{r}
#| label: tbl-itemrestscore40
#| tbl-cap: Item-restscore detection rate with 40 items (3 misfit)

ir40_df %>% 
  group_by(samplesize) %>% 
  filter(ir_padj == "sign. misfit") %>% 
  count(item, .drop = F) %>% 
  mutate(Percent = n*100/500) %>% 
  ungroup() %>% 
    mutate(samplesize = factor(samplesize, levels = c(150,250,500,1000,2000),
                             labels = c("n = 150",
                                        "n = 250",
                                        "n = 500",
                                        "n = 1000",
                                        "n = 2000"))) %>% 
  rename(`Sample size` = samplesize,
         Item = item,
         `Percent misfit` = Percent) %>% 
  select(!n) %>% 
  pivot_wider(id_cols = c(Item), names_from = c(`Sample size`),
              values_from = `Percent misfit`) %>% 
  mutate(Item = factor(Item, levels = 1:40, ordered = T)) %>% 
  arrange(Item) %>% 
  kable()
  
```


## Session info

This documents the specific R packages and versions used in this study. Note that the simulations were conducted using `easyRasch` version 0.3.3, while the plots and tables generated directly from `easyRasch` were done using version 0.3.3.2.

```{r}
sessionInfo()
```

