---
title: |
  <div><b>Test and Item Analysis Report</b><img src="hexbin_html.png" width="160"></img></div>

  **************
  
output:
    html_document:
       number_sections: false
       toc: true
author: |
  `r ifelse(params$dataset == "", "", params$dataset)`
  
  `r ifelse(params$author == "", "", params$author)`
  
  This report was created on `r format(Sys.Date(), "%d.%m.%Y")`

  &copy; 2020 ShinyItemAnalysis: Test and item analysis, version `r packageVersion("ShinyItemAnalysis")`
    
  **************
  
params:
  author: NA
  dataset: NA
  itemNames: NA
  a: NA
  k: NA
  gr: NA
  results: NA
  cutScore: NA
  histogram_totalscores: NA
  standardscores_table: NA
  corr_plot: NA
  corr_plot_numclust: NA
  corr_plot_clustmethod: NA
  corr_type: NA
  scree_plot: NA
  isCriterionPresent: NA
  validity_check: NA
  validity_plot: NA
  validity_table: NA
  DDplot: NA
  DDplotRange1: NA
  DDplotRange2: NA
  DDplotNumGroups: NA
  DDplotDiscType: NA
  cronbachs_alpha_table: NA
  itemexam: NA
  distractor_plot: NA
  type_distractor_plot: NA
  distractor_plot_legend_length: NA
  logreg: NA
  zlogreg: NA
  zlogreg_irt: NA
  nlsplot: NA
  multiplot: NA
  wrightMap: NA
  irt_type: NA
  irt: NA
  irtiic: NA
  irttif: NA
  irtcoef: NA
  irtfactor: NA
  irtability: NA
  isGroupPresent: NA
  histCheck: NA
  DIF_total_table: NA
  DIF_total_hist: NA
  DIF_total_ttest: NA
  deltaplotCheck: NA
  DIF_deltaplot: NA
  DIF_deltaplot_text: NA
  MHCheck: NA
  DIF_MH_print: NA
  logregCheck: NA
  DIF_logistic_print: NA
  DIF_logistic_plot: NA
  multiCheck: NA
  DDF_multinomial_print: NA
  DDF_multinomial_plot: NA
  sessionInfo: NA
---

<style>
body {text-align: justify}
</style>
**************

# Introduction

This report was created by by R version `r getRversion()` and its package ShinyItemAnalysis version `r packageVersion("ShinyItemAnalysis")`. ShinyItemAnalysis provides test and item analysis and it is available on [CRAN](https://CRAN.R-project.org/package=ShinyItemAnalysis) and also [online](https://shiny.cs.cas.cz/ShinyItemAnalysis/). 

To cite ShinyItemAnalysis application in publications, please, use:


  Martinkova P., \& Drabinova A.
  (2018)
  ShinyItemAnalysis for teaching psychometrics and to enforce routine analysis of educational tests. 
  The R Journal, 10(2), 503-515.
  [doi: 10.32614/RJ-2018-074](https://doi.org/10.32614/RJ-2018-074)
  
  
This program is free software and you can redistribute it and or modify it under the terms of the GNU GPL 3 as published by the Free Software Foundation. This program is distributed in the hope that it will be useful, but without any warranty; without even the implied warranty of merchantability of fitness for a particular purpose. Visit [www.ShinyItemAnalysis.org](www.ShinyItemAnalysis.org) to learn more!

This project was supported by Czech Science Foundation under grant number [GJ15-15856Y](http://www.cs.cas.cz/martinkova/psychometrics.html).

**************

# Summary

## Total scores

### Summary table of total scores
The table below summarizes basic characteristics of total scores including minimum and maximum, mean, median, standard deviation, skewness and kurtosis. The kurtosis here is estimated by a sample kurtosis $\frac{m_4}{m_2^2}$, where $m_4$ is the fourth central moment and $m_2$ is the second central moment. The skewness is estimated by a sample skewness $\frac{m_3}{m_2^{3/2}}$, where $m_3$ is the third central moment. The kurtosis for normally distributed scores is near the value of 3 and the skewness is near the value of 0. 

```{r total-scores-table, echo = FALSE}
library(knitr)
table = params$results
colnames(table) = NULL

kable(t(table), digits = c(0, 0, 2, 2, 2, 2, 2), align = 'c')
```

### Histogram of total scores
For cut-score `r params$cutScore`, the blue part of a histogram shows respondents with a total score above the cut-score, the grey column shows respondents with a total score equal to the cut-score and the red part of a histogram shows respondents below the cut-score.

```{r total-scores-histogram, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300}
params$histogram_totalscores + 
  ggtitle('Histogram of total scores') + 
  theme(text = element_text(size = 11), 
        plot.title = element_text(face = "bold")) + 
  ggtitle("")
```

**************

# Scoring

### Summary table of standard scores
The total score, also known as the raw score, is the total number of correct answers. It can be used to compare an individual score to a norm group, e.g. if the mean is 12, then an individual score can be compared to see if it is below or above this average. The percentile indicates the value below which a percentage of observations falls, e.g. an individual score at the 80th percentile means that the individual score is the same or higher than the scores of 80% of all respondents. The success rate is the percentage of correct answers, e.g. if the maximum points of a test is equal to 20 and an individual score is 12 then the success rate is 12/20 = 0.6, i.e. 60%. The Z-score, or the standardized score, is a linear transformation of the total score with a mean of 0 and with a variance of 1. If X is the total score, M is its mean and SD is its standard deviation then Z-score = (X - M) / SD. The T-score is the transformed Z-score with a mean of 50 and a standard deviation of 10. If Z is Z-score then T-score = (Z * 10) + 50.

```{r standard-scores-table, echo = FALSE}
table = params$standardscores_table
kable(table, align = 'c')
```

`r if (params$irt_type != "none") {"More complex estimates of ability are provided in the IRT section."}`

`r if (any(params$corr_plot != "") | params$validity_check) {"**************"}`

`r if (any(params$corr_plot != "") | params$validity_check) {"# Validity"}`

`r if (any(params$corr_plot != "")) {"## Correlation structure"}`

`r if (any(params$corr_plot != "")) {"### Correlation heat map "}`
`r if (any(params$corr_plot != "")) {if (params$corr_type == "pearson") {"A correlation heat map displays Pearson correlations of items. Pearson correlation coefficient describes linear correlation between two random variables X and Y."}}`

`r if (any(params$corr_plot != "")) {if (params$corr_type == "polychoric") {"A correlation heat map displays polychoric correlations of items. Polychoric/tetrachoric correlation between two ordinal/binary variables is calculated from their contingency table, under the assumption that the ordinal variables dissect continuous latent variables that are bivariate normal."}}`

`r if (any(params$corr_plot != "")) {if (params$corr_type == "spearman") {"A correlation heat map displays Spearman correlations of items. Spearman's rank correlation coefficient describes strength and direction of monotonic relationship between random variables X and Y, i.e. dependence between the rankings of two variables."}}`

`r clustmethod_names <- c("None" = "none",
                          "Ward's"  = "ward.D",
                          "Ward's n. 2" = "ward.D2",
                          "single" = "single",
                          "complete" = "complete",
                          "average" = "average",
                          "McQuitty" = "mcquitty",
                          "median" = "median",
                          "centroid" = "centroid")`
                          
`r if (any(params$corr_plot != "") & params$corr_plot_clustmethod != "none") {paste("A correlation heat map is reordered using hierarchical clustering with", names(clustmethod_names)[clustmethod_names == params$corr_plot_clustmethod], "linkage method. The number of highlighted clusters is", params$corr_plot_numclust, ". ")}`

```{r corr-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300}
if (any(params$corr_plot != "")) {
  corP <- params$corr_plot
  numclust <- params$corr_plot_numclust
  clustmethod <- params$corr_plot_clustmethod
  
  tlcex <- ifelse(ncol(corP) > 30, 0.7, 1)
  
  if (clustmethod == "none"){
    corrplot(corP, tl.cex = tlcex)
  } else {
    corrplot(corP, tl.cex = tlcex, order = "hclust", 
             hclust.method = clustmethod, addrect = numclust)
  }
}
```

`r if (any(params$corr_plot != "")) {"### Scree plot"}` 

`r if (any(params$corr_plot != "")) {"A scree plot displays the eigenvalues associated with a component or a factor in descending order."}`

```{r scree-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300}
if (any(params$corr_plot != "")) {
  params$scree_plot + 
    ggtitle("") + 
    theme(text = element_text(size = 11), 
          plot.title = element_text(face = "bold"))
}
```

`r if (params$validity_check) {"## Predictive validity"}`

`r if(!(params$isCriterionPresent) & params$validity_check){"No criterion variable vector is present. The predictive validity analysis could not have been generated."}`

```{r validity-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300, message = FALSE}
if(params$isCriterionPresent & params$validity_check){
  params$validity_plot + 
    theme(text = element_text(size = 11),
          legend.box.just = "top",
          legend.justification = c("left", "top"), 
          legend.position = c(0, 1), 
          legend.box = "horizontal", 
          legend.box.margin = margin(3, 3, 3, 3))
}
```


```{r validity-table, echo = FALSE, results = 'hide'}
if(params$isCriterionPresent & params$validity_check){
  tab <- params$validity_table
  p.val <- tab["p-value"]
  rho <- tab[1]

  txt1 <- ifelse(p.val < 0.05, "less than", "larger than")
  txt2 <- ifelse(p.val < 0.05, "reject", "do not reject")
  txt3 <- ifelse(p.val < 0.05, 
                 paste("total score and the criterion variable are", ifelse(rho > 0, "positively", "negatively"),
                       "correlated."),
                 paste("we cannot conclude that a significant correlation between total score
                        and criterion variable exists."))
  txt4 <- paste("A test for association between the total score and the criterion variable is based on Spearman`s ",
                expression(rho), 
                ". The null hypothesis is that the correlation is 0. ", sep = "") 
  txt5 <- paste("Results: ", expression(rho), " = ", rho, " (p-value ", p.val ,").", sep = "") 
  txt6 <- paste("Interpretation: The p-value is", txt1, "0.05, thus we", txt2, "the null hypothesis;", txt3)
  
  txt <- paste(txt4, txt5, txt6, collapse = "/n")
}

```

`r if(params$isCriterionPresent & params$validity_check){txt}`

**************

# Traditional item analysis

## Item analysis

### Difficulty/Discrimination plot

```{r DD-plot-interpretation, echo = FALSE}
txt1 <- "Difficulty (red) of items is estimated as a percent of the respondents who answered correctly to that item. "

range1 <- params$DDplotRange1
range2 <- params$DDplotRange2
numgroups <- params$DDplotNumGroups

txt2 <- switch(params$DDplotDiscType,
               "RIR" = "Discrimination (blue) is described by Pearson correlation between an item and the rest of the items.",
               "RIT" = "Discrimination (blue) is described by Pearson correlation between an item and the total score.",
               "none" = "",
               "ULI" = if (any(range1 != 1, range2 != 3, numgroups != 3)) {
                 paste0("Discrimination (blue) is described by the difference of the percent correct in the ",
                        "<b>", range1, "</b>",
                        ifelse(range1 >= 4, "-th", switch(range1, "1" = "-st", "2" = "-nd", "3" = "-rd")),
                        " and <b>", range2, "</b>",
                        ifelse(range2 >= 4, "-th", switch(range2, "1" = "-st", "2" = "-nd", "3" = "-rd")),
                        " group of respondents out of a total number of ",
                        "<b>", numgroups, "</b>",
                        " groups. ")
                 } else {
                   "Discrimination (blue) is described by the difference of the percent correct in the upper and lower third of respondents (Upper-Lower Index, ULI). By a rule of thumb, it should not be lower than 0.2 (borderline in the plot), except for very easy or very difficult items."
                   })

HTML(paste0(txt1, txt2))
```

```{r DD-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300}
params$DDplot + 
  ggtitle("") + 
  theme(text = element_text(size = 11), 
        plot.title = element_text(face = "bold")) 
```

### Cronbach's alpha

Cronbach's alpha is an estimate of internal consistency of a psychometric test. It is a function of the number of items in a test, the average covariance between item-pairs, and the variance of the total score (Cronbach, 1951).

```{r cronbachs-alpha, echo = FALSE}
table = params$cronbachs_alpha_table

kable(table, digits = c(0, 3, 0), align = 'c')
```

### Traditional item analysis table

```{r item-table-text, echo = FALSE}
range1 <- params$DDplotRange1
range2 <- params$DDplotRange2
numgroups <- params$DDplotNumGroups

txt1 <- "Diff. - difficulty of the item is estimated as its average score divided by its range, Avg. score - average score of the item, SD - standard deviation, ULI - Upper-Lower Index calculated as the difference in the percent correct between the upper and lower third of respondents, RIT - Pearson correlation between an item and the total score, RIR - Pearson correlation between an item and the rest of the items, Alpha Drop - Cronbach's alpha of a test without a given item"

if (params$DDplotNumGroups == 3 & params$DDplotRange1 == 1 & params$DDplotRange2 == 3){
  txt2 <- ". "
} else {
  txt2 <- paste0(", gULI - generalized ULI calculated as the difference between the difficulty recorded in the ", 
                 "<b>", range1, "</b>",
                 ifelse(range1 >= 4, "-th", switch(range1, "1" = "-st", "2" = "-nd", "3" = "-rd")),
                 " and <b>", range2, "</b>",
                 ifelse(range2 >= 4, "-th", switch(range2, "1" = "-st", "2" = "-nd", "3" = "-rd")),
                 " group of respondents out of a total number of ",
                 "<b>", numgroups, "</b>",
                 " groups. ")
}

HTML(paste0(txt1, txt2))
```

```{r item-table, echo = FALSE}

if (params$DDplotNumGroups == 3 & params$DDplotRange1 == 1 & params$DDplotRange2 == 3){
  table = params$itemexam
  tab = format(round(table[, -1], 2), nsmall = 2)
  table = data.frame(table[, 1], tab[, -8])
  colnames(table) = c("Item", "Diff.", "Avg. score", "SD", "ULI", "RIT", "RIR", "Alpha Drop")
  rownames(table) = NULL

  kable(table, digits = c(0, 2, 2, 2, 2, 2, 2), align = 'c')
} else {
  table = params$itemexam
  tab = format(round(table[, -1], 2), nsmall = 2)
  table = data.frame(table[, 1], tab)
  colnames(table) = c("Item", "Diff.", "Avg. score", "SD", "ULI", "RIT", "RIR", "Alpha Drop",   "gULI")
  rownames(table) = NULL

  kable(table, digits = c(0, 2, 2, 2, 2, 2, 2, 2), align = 'c')
}
```

## Distractor analysis

### Distractor plot

Respondents are divided into a selected number of groups by their total score. Subsequently, the percentage of respondents in each group who selected a given answer (correct answer or distractor) is displayed. The correct answer should be selected more often by the respondents with a higher total score than by those with a lower total score, i.e. the solid line should be increasing. The distractor should work in the opposite direction, i.e. the dotted lines should be decreasing.

### Multinomial plot

In the multinomial plot, points represent proportion of selected option with respect to total score. Their size is determined by count of respondents who achieved given level of total score and who selected given option.

```{r distractor-multi-plot-short, echo = FALSE, fig.height = 4.5, fig.width = 11}
if (params$distractor_plot_legend_length <= 6){
  if (params$type_distractor_plot == "Combinations"){
    for (i in 1:length(params$distractor_plot)){
      g1 <- params$distractor_plot[[i]]
      g2 <- params$multiplot[[i]]
      if (is.null(g2)) {
        g2 <- textGrob("")
        g1 <- g1 + theme(legend.position = c(0.14, 0.8),
                         legend.margin = margin(0, 0, 0, 0, "cm"))
      } else {
        g2 <- g2 + theme(legend.position = c(0.04, 0.97),
                         legend.margin = margin(0, 0, 0, 0, "cm"),
                         legend.spacing.x = unit(0.25, "cm"))
      }
      grid.arrange(g1, g2, ncol = 2)
    }
  } else {
    for (i in 1:length(params$distractor_plot)){
      g1 <- params$distractor_plot[[i]] + theme(legend.position = c(0.14, 0.8),
                                                legend.margin = margin(0, 0, 0, 0, "cm"))
      g2 <- params$multiplot[[i]] + theme(legend.position = c(0.04, 0.97),
                                          legend.margin = margin(0, 0, 0, 0, "cm"),
                                          legend.spacing.x = unit(0.25, "cm"))
      if (is.null(g2)) g2 <- textGrob("")
      grid.arrange(g1, g2, ncol = 2)
    }
  }
}
```

```{r distractor-multi-plot-long, echo = FALSE, fig.height = 4.5, fig.width = 11}
if (params$distractor_plot_legend_length > 6){
  if (params$type_distractor_plot == "Combinations"){
    for (i in 1:length(params$distractor_plot)){
      g1 <- params$distractor_plot[[i]]
      g2 <- params$multiplot[[i]] 
      if (is.null(g2)) {
        g2 <- textGrob("")
        tmp <- ggplotGrob(g1 + theme(legend.position = "right",
                                     legend.text = element_text(size = 10),
                                     legend.key.size = unit(0.8, 'lines')))$grobs
        legend <- tmp[[which(sapply(tmp, function(x) x$name) == "guide-box")]]
      } else {
        tmp <- ggplotGrob(g2 + theme(legend.position = "right", 
                                     legend.box = "vertical",
                                     legend.spacing = unit(-0.2, "cm"),
                                     legend.text = element_text(size = 10),
                                     legend.key.size = unit(0.8, 'lines')))$grobs
        legend <- tmp[[which(sapply(tmp, function(x) x$name) == "guide-box")]]
        g2 <- g2 + theme(legend.position = "none")
      }
      grid.arrange(g1, g2, legend, ncol = 3, widths = c(0.45, 0.45, 0.1))
    }
  } else {
    for (i in 1:length(params$distractor_plot)){
      g1 <- params$distractor_plot[[i]] + theme(legend.position = c(0.18, 0.75), 
                                                legend.margin = margin(0, 0, 0, 0, "cm"))
      g2 <- params$multiplot[[i]] 
        
      if (is.null(g2)) {
        g2 <- textGrob("")
        legend <- textGrob("")
      } else {
        tmp <- ggplotGrob(g2 + theme(legend.position = "right", 
                                     legend.box = "vertical",
                                     legend.spacing = unit(-0.2, "cm"),
                                     legend.text = element_text(size = 10),
                                     legend.key.size = unit(0.8, 'lines')))$grobs
        legend <- tmp[[which(sapply(tmp, function(x) x$name) == "guide-box")]]
        g2 <- g2 + theme(legend.position = "none")
      }
      grid.arrange(g1, g2, legend, ncol = 3, widths = c(0.45, 0.45, 0.1))
    }  
  }
}
```

<!-- IRT SECTION -->

```{r IRT-section-settings, echo = FALSE}
show_IRT <- (params$irt_type != "none")
```

```{r IRT-section, eval = show_IRT, echo = FALSE, results = "asis"}
cat("
************** 
# IRT models

Item Response Theory (IRT) models are mixed-effect regression models in which the respondent's ability $\\theta$ is assumed to be a latent and is estimated together with item parameters.")
```

```{r IRT-wright-map-text, echo = FALSE, results = "asis", eval = show_IRT}
cat("
## Wright (item-person) map using 1PL IRT model
    
A Wright map, also called an item-person map, is a graphical tool to display person estimates and item parameters. The person side (left) represents a histogram of estimated knowledge of the  respondents. The item side (right) displays estimates of the difficulty of particular items.
")
```

```{r IRT-wright-map-plot, echo = FALSE, eval = show_IRT, fig.height = 4, fig.width = 11, fig.align = 'center', dpi = 300}
params$wrightMap + theme(text = element_text(size = 11))
```

```{r IRT-equation, echo = FALSE, results = "hide", eval = show_IRT}
IRT_selected_model <- switch(params$irt_type,
                             "rasch" = "Rasch IRT model",
                             "1pl"   = "1PL IRT model",
                             "2pl"   = "2PL IRT model",
                             "3pl"   = "3PL IRT model",
                             "4pl"   = "4PL IRT model")

IRT_equation <- switch(params$irt_type,
                             "rasch" = "$${P}(Y_{ij} = 1|\\theta_{i}, a, b_j) = \\frac{e^{(\\theta_i - b_j)}}{1 + e^{(\\theta_i - b_j)}}$$",
                             "1pl"   = "$${P}(Y_{ij} = 1|\\theta_{i}, a, b_j) = \\frac{e^{a(\\theta_i - b_j)}}{1 + e^{a(\\theta_i - b_j)}}$$",
                             "2pl"   = "$${P}(Y_{ij} = 1|\\theta_{i}, a_j, b_j) = \\frac{e^{a_j(\\theta_i - b_j)}}{1 + e^{a_j(\\theta_i - b_j)}}$$",
                             "3pl"   = "$${P}(Y_{ij} = 1|\\theta_{i}, a_j, b_j, c_j) = c_j + (1 + c_j) \\frac{e^{a_j(\\theta_i - b_j)}}{1 + e^{a_j(\\theta_i - b_j)}}$$",
                             "4pl"   = "$${P}(Y_{ij} = 1|\\theta_{i}, a_j, b_j, c_j, d_j) = c_j + (d_j + c_j) \\frac{e^{a_j(\\theta_i - b_j)}}{1 + e^{a_j(\\theta_i - b_j)}}$$")
```

```{r IRT-equation-text, echo = FALSE, results = "asis", eval = show_IRT}
txt <- ifelse(params$irt_type == "rasch", 
              "Model parameters are estimated using a marginal maximum likelihood method. Ability $\\theta$ is assumed to follow normal distribution with a freely estimated variance.",
              "Model parameters are estimated using a  marginal maximum likelihood method. Ability $\\theta$ is assumed to follow standard normal distribution." )


cat("
## Equation \n\n",

paste0("All subsequent analyses are based on the selected ", IRT_selected_model, ":"),
IRT_equation,
txt
)
```

```{r IRT-ability-estimates-text, echo = FALSE, results = "asis", eval = show_IRT}
cat("
## Ability estimates \n\n",

paste0("Ability is estimated using an EAP algorithm and a ", IRT_selected_model, ". The relationship between ability estimates (factor scores, F-scores) and standardized total test scores (Z-scores) can be seen on the scatter plot below. A table with ability estimates for all respondents can be downloaded from the application. ")
)
```

```{r IRT-ability-estimates-table, echo = FALSE, results = "asis", eval = show_IRT}
table = round(params$irtability, 2)
table = cbind(c("Total Scores", "Z-Scores", "F-Scores"), table)
colnames(table)[1] = "Score type"
  
kable(table, align = c('l', rep('c', 7)))
```

```{r IRT-ability-scatterplot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = "center", dpi = 300, eval = show_IRT}
params$irtfactor + 
  ggtitle('') +
  theme(text = element_text(size = 11))
```

```{r IRT-curves-text, echo = FALSE, results = "asis", eval = show_IRT}
cat("
## Item characteristic and information curves
")
```

```{r IRT-curves-plots-short, echo = FALSE, results = "hide", fig.height = 14.5, fig.width = 11, eval = show_IRT}
if (ncol(params$a) < 50){
  g1 <- params$irt + 
    ggtitle("Item characteristic curves") + 
    xlab(expression(theta)) + 
    ylab(expression("P("*theta*")")) + 
    guides(colour = guide_legend(ncol = max(ncol(params$a) %/% 30, 1))) + 
    theme(legend.position = "right", 
          legend.text = element_text(size = 10))
  g2 <- params$irtiic + 
      ggtitle("Item information curves") + 
      xlab(expression(theta)) + 
      ylab(expression("I("*theta*")")) + 
      guides(colour = guide_legend(ncol = max(ncol(params$a) %/% 30, 1))) + 
      theme(legend.position = "right", 
            legend.text = element_text(size = 10))
  g3 <- params$irttif + 
      ggtitle("Test information function") + 
      xlab(expression(theta)) + 
      ylab(expression("I("*theta*")")) + 
      theme(legend.position = c(0.9, 0.8), 
            legend.text = element_text(size = 10))
  
  tmp <- ggplotGrob(g1)$grobs
  legend <- tmp[[which(sapply(tmp, function(x) x$name) == "guide-box")]]

  grid.arrange(arrangeGrob(arrangeGrob(g1 + theme(legend.position = "none"), 
                                       g2 + theme(legend.position = "none"), ncol = 1), 
                           legend, ncol = 2, widths = c(0.85, 0.15)), 
               arrangeGrob(g3, textGrob(""), ncol = 2, widths = c(0.9, 0.1)), 
               ncol = 1, heights = c(0.66, 0.33))
}
```

```{r IRT-curves-plots-long, echo = FALSE, results = "hide", fig.height = 16, fig.width = 11, eval = show_IRT}
if (ncol(params$a) >= 50){
    g1 <- params$irt + 
      ggtitle("Item characteristic curves") + 
      xlab(expression(theta)) + 
      ylab(expression("P("*theta*")")) + 
      guides(colour = guide_legend(nrow = max(ncol(params$a) %/% 12, 1))) + 
      theme(legend.position = "bottom", 
            plot.margin = unit(c(0, 0, 0, 0), "cm"),
            legend.margin = margin(t = -5, unit = 'cm'),
            legend.text = element_text(size = 10),
            legend.key.size = unit(0.8, 'lines'))
    g2 <- params$irtiic + 
      ggtitle("Item information curves") + 
      xlab(expression(theta)) + 
      ylab(expression("I("*theta*")")) + 
      guides(colour = guide_legend(nrow = max(ncol(params$a) %/% 12, 1))) + 
      theme(legend.position = "none", 
            plot.margin = unit(c(0, 0, 0, 0), "cm"))
    g3 <- params$irttif + 
      ggtitle("Test information function") + 
      xlab(expression(theta)) + 
      ylab(expression("I("*theta*")")) + 
      theme(legend.position = c(0.9, 0.8), 
            legend.text = element_text(size = 10), 
            plot.margin = unit(c(-4, 0, 4, 0), "cm"))
    
    tmp <- ggplotGrob(g1)$grobs
    legend <- tmp[[which(sapply(tmp, function(x) x$name) == "guide-box")]]

    grid.arrange(g1 + theme(legend.position = "none"), 
                 g2,
                 legend,
                 g3,
                 ncol = 1)
    
}
```

```{r IRT-parameter-estimates-text, echo = FALSE, results = "asis", eval = show_IRT}
cat("
## Parameter estimates and item fit

Estimates of parameters are completed by their standard errors (SE) and by signed Chi-squared statistics S-X2 (see Orlando and Thissen, 2000). P-values lower than 0.05 indicate suspicious items not fitting the selected IRT model. S-X2 is computed only when no missing data are present.
")
```

```{r IRT-parameter-estimates-table, echo = FALSE, eval = show_IRT}
table = params$irtcoef
colnames(table) = c("a", "SE(a)", "b", "SE(b)", "c", "SE(c)", "d", "SE(d)", "S-X2", "df", "p-value")
  
kable(table, align = c('l', rep('c', 11)))
```

<!-- DIF SECTION -->

```{r DIF-section-settings, echo = FALSE}
show_DIF <- any(c(params$histCheck, params$deltaplotCheck, params$MHCheck, params$logregCheck, params$multiCheck))
check_group <- !params$isGroupPresent
```

```{r DIF-section, eval = show_DIF, echo = FALSE, results = 'asis'}
txt <- ifelse(check_group, "No group vector is present. DIF and DDF analyses could have not been generated!", "")
cat(paste("
************** 
# DIF/Fairness analysis", 
txt))
```

<!-- DIF SECTION - histograms -->

```{r DIF-section-histograms-settings, echo = FALSE, results = 'asis'}
eval_histograms <- (params$histCheck & params$isGroupPresent)
```

```{r DIF-section-histograms-part1, eval = eval_histograms, echo = FALSE, results = 'asis'}
cat("
## Total scores by group
DIF analysis may come to a different conclusion than test of group differences in total scores. Two groups may have the same distribution of total scores, yet, some item may function differently for the two groups. Also, one of the groups may have a signifficantly lower total score, yet, it may happen that there is no DIF item [(see Martinkova et al., 2017)](https://doi.org/10.1187/cbe.16-10-0307).

### Summary table of total scores in a reference and focal group
")
```

```{r DIF-section-histograms-table, echo = FALSE, results = "asis", eval = eval_histograms}
kable(params$DIF_total_table, 
          digits = c(0, 0, 0, 3, 3, 3, 3), 
          align = c('l', rep('c', 6)))
```

```{r DIF-section-histograms-part2, eval = eval_histograms, echo = FALSE, results = 'asis'}
cat(paste("### Histograms of total scores by group \n"))
```

```{r DIF-section-histograms-plot, eval = eval_histograms, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300}
params$DIF_total_hist +
  theme(text = element_text(size = 11), 
        plot.title = element_text(face = "bold"),
        legend.position = c(0.1, 0.87)) 
```

```{r DIF-section-histograms-part3, eval = eval_histograms, echo = FALSE, results = 'asis'}
cat("
### Comparison of total scores 

Test for difference in total scores between reference and focal group is based od Welch two sample t-test. Explanation: Diff. (CI) - difference in means of total scores with 95% confidence interval, t-value - test statistic, df - degrees of freedom, p-value - if it is lower than 0.05, it means significant difference in total scores.
")
```

```{r DIF-section-histograms-ttest, echo = FALSE, eval = eval_histograms}
kable(params$DIF_total_ttest,
      align = rep('c', 4), row.names = F)
```

<!-- DIF SECTION - delta plot -->

```{r DIF-section-deltaplot-settings, echo = FALSE, results = 'asis'}
eval_deltaplot <- (params$deltaplotCheck & params$isGroupPresent)
if (eval_deltaplot) res_deltaplot <- params$DIF_deltaplot_text
```

```{r DIF-section-deltaplot-part1, eval = eval_deltaplot, echo = FALSE, results = 'asis'}
txt1 <- ifelse(res_deltaplot$purify, 
                     "Item purification was used.", 
                     "Item purification was not applied.")
txt2 <- paste0("The detection threshold is ", round(res_deltaplot$thr[length(res_deltaplot$thr)], 2), ". ")
txt3 <- ifelse(any(res_deltaplot$DIFitems == "no DIF item detected"), 
              "**No DIF item was detected. **", 
              paste0("**Items detected as DIF: ", 
                     paste(params$itemNames[res_deltaplot$DIFitems], collapse = ", "), ". ** "))

cat(paste("
## Delta plot method
A delta plot compares the proportions of correct answers per item in the two groups. It displays non-linear transformation of these proportions using quantiles of standard normal distributions (so called delta scores) of each item for the two groups in a scatterplot called diagonal plot or delta plot. An item is under suspicion of DIF if the delta point departs considerably from the major axis. \n

### Summary table \n", 
txt1, txt2, "\n\n", 
txt3))
```

```{r DIF-section-deltaplot-table, echo = FALSE, eval = eval_deltaplot}
res <- params$DIF_deltaplot_text
tab <- cbind(res$Props, res$Deltas, res$Dist[, ncol(res$Dist)])
colnames(tab) <- c("Prop. Ref", "Prop. Foc", "Delta Ref", "Delta Foc", "Dist.")
rownames(tab) <- params$itemNames
    
kable(tab, digits = c(0, 0, 3, 3, 3, 3, 3), 
      align = c('l', rep('c', 6)))
```

```{r DIF-section-deltaplot-part2, eval = eval_deltaplot, echo = FALSE, results = 'asis'}
cat("
### Delta plot
Points represents delta scores. Black solid is a major axis. Red dashed lines stand for detection threshold.
")
```

```{r DIF-section-deltaplot-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300, eval = eval_deltaplot}
print(params$DIF_deltaplot +
        geom_text(size = 0.1) +
        theme(text = element_text(size = 11), 
              plot.title = element_text(face = "bold")) + 
        ggtitle(""))
```

<!-- DIF SECTION - Mantel-Haenszel -->

```{r DIF-section-MH-settings, echo = FALSE, results = 'asis'}
eval_MH <- (params$MHCheck & params$isGroupPresent)
if (eval_MH) res_MH <- params$DIF_MH_print
```

```{r DIF-section-MH-part1, eval = eval_MH, echo = FALSE, results = 'asis'}
txt1 <- paste(ifelse(res_MH$purification, 
                     "Item purification was used.", 
                     "Item purification was not applied."),
              ifelse(is.null(res_MH$p.adjust.method), 
                     "No p-value adjustment for multiple comparisons was used.",
                     paste(switch(res_MH$p.adjust.method, 
                                  "holm" = "Holm's", 
                                  "hochberg" = "Hochberg's", 
                                  "hommel" = "Hommel's", 
                                  "bonferroni" = "Bonferroni", 
                                  "BH" = "Benjamini-Hochberg", 
                                  "BY" = "Benjamini-Yekutieli", 
                                  "fdr" = "Benjamini-Hochberg", 
                                  "none" = "No"), 
                           "p-value adjustment for multiple comparisons was used.")))
txt2 <- paste0("The detection threshold is ", round(res_MH$thr[length(res_MH$thr)], 2), ". ")
txt3 <- ifelse(any(res_MH$DIFitems == "No DIF item detected"), 
                      "**No DIF item was detected. **", 
                      paste0("**Items detected as DIF: ", 
                             paste(params$itemNames[res_MH$DIFitems], collapse = ", "), ".**  "))

cat(paste("
## Mantel-Haenszel test
 \n

### Summary table \n", 
txt1, txt2, "\n\n",
txt3))
```

```{r DIF-section-MH-table, echo = FALSE, eval = eval_MH}
if (res_MH$p.adjust.method == "none"){
  star <- symnum(round(res_MH$p.value, 4), 
               c(0, 0.001, 0.01, 0.05, 0.1, 1), 
               symbols = c("***", "**", "*", ".", ""))
  pval <- paste0(sprintf("%.3f", res_MH$p.value), star)
} else {
  star <- symnum(round(res_MH$adjusted.p, 4), 
               c(0, 0.001, 0.01, 0.05, 0.1, 1), 
               symbols = c("***", "**", "*", ".", ""))
  pval <- paste0(sprintf("%.3f", res_MH$adjusted.p), star)
  pval <- cbind(sprintf("%.3f", res_MH$p.value),
                paste0(sprintf("%.3f", res_MH$adjusted.p), star))
}

deltaMH <- round(-2.35 * log(res_MH$alphaMH), 4)
esize <- symnum(abs(deltaMH), 
                c(0, 1, 1.5, Inf), 
                symbols = c("A", "B", "C"))  
deltaMH <- paste(sprintf("%.3f", deltaMH), esize)

tab <- cbind(sprintf("%.3f", res_MH$MH), pval, 
             sprintf("%.3f", res_MH$alphaMH), deltaMH)
rownames(tab) <- params$itemNames

if (dim(tab)[2] == 5){
  colnames(tab) <- c("Stat.", "P-value", "Adj. p-value", "alpha_MH", "delta_MH")
  kable(tab, digits = c(0, 3, 3, 0, 0, 0), 
        align = 'c')
} else {
  colnames(tab) <- c("Stat.", "P-value", "alpha_MH", "delta_MH")
  kable(tab, digits = c(0, 3, 0, 0, 0, 0), 
        align = 'c')
}
```

```{r DIF-section-MH-part2, eval = eval_MH, echo = FALSE, results = 'asis'}
cat("Signif. codes: 0 '&#42;&#42;&#42;' 0.001 '&#42;&#42;' 0.01 '&#42;' 0.05 '.' 0.1 '&nbsp;' 1  

Effect size (ETS Delta scale): 

 'A': negligible effect 
 
 'B': moderate effect 
 
 'C': large effect  

Effect size codes: 0 'A' 1.0 'B' 1.5 'C' (for absolute values of delta_MH) 
")
```

<!-- DIF SECTION - logistic regression -->

```{r DIF-section-logistic-settings, echo = FALSE, results = 'asis'}
eval_logistic <- (params$logregCheck & params$isGroupPresent)
eval_logistic_plot <- eval_logistic & !is.null(params$DIF_logistic_plot[[1]])

if (eval_logistic) res_logistic <- params$DIF_logistic_print
```

```{r DIF-section-logistic-part1, eval = eval_logistic, echo = FALSE, results = 'asis'}

txt1 <- switch(res_logistic$type, 
               "both" = "Both types of DIF were tested.", 
               "udif" = "Uniform DIF was tested only.", 
               "nudif" = "Non-uniform DIF was tested only.")
txt2 <- paste(ifelse(res_logistic$purification, 
                     "Item purification was used.", 
                     "Item purification was not applied."),
              ifelse(is.null(res_logistic$p.adjust.method), 
                     "No p-value adjustment for multiple comparisons was used.",
                     paste(switch(res_logistic$p.adjust.method, 
                                  "holm" = "Holm's", 
                                  "hochberg" = "Hochberg's", 
                                  "hommel" = "Hommel's", 
                                  "bonferroni" = "Bonferroni", 
                                  "BH" = "Benjamini-Hochberg", 
                                  "BY" = "Benjamini-Yekutieli", 
                                  "fdr" = "Benjamini-Hochberg", 
                                  "none" = "No"), 
                           "p-value adjustment for multiple comparisons was used.")))
txt3 <- paste0("The detection threshold is ", round(res_logistic$thr[length(res_logistic$thr)], 2), ". ")
txt4 <- ifelse(any(res_logistic$DIFitems == "No DIF item detected"), 
                      "\\textbf{No DIF item was detected. }  ", 
                      paste0("\\textbf{Items detected as DIF: ", 
                             paste(params$itemNames[res_logistic$DIFitems], collapse = ", "), ".}  "))

cat(paste("
## Logistic regression method
Logistic regression allows for detection of uniform and non-uniform DIF by adding a group specific intercept (uniform DIF) and group specific interaction (non-uniform DIF) into the model and by testing for their significance. \n

### Summary table \n", 
txt1, txt2, txt3, "\n\n",
txt4, "\\vspace{1em}"))
```

```{r DIF-section-logistic-table, echo = FALSE, eval = eval_logistic}
if (res_logistic$p.adjust.method == "none"){
  star <- symnum(round(res_logistic$p.value, 4), 
               c(0, 0.001, 0.01, 0.05, 0.1, 1), 
               symbols = c("***", "**", "*", ".", ""))
  pval <- paste0(sprintf("%.3f", res_logistic$p.value), star)
} else {
  pval <- res_logistic$adjusted.p
  star <- symnum(round(res_logistic$adjusted.p, 4), 
               c(0, 0.001, 0.01, 0.05, 0.1, 1), 
               symbols = c("***", "**", "*", ".", ""))
  pval <- cbind(sprintf("%.3f", res_logistic$p.value),
                paste0(sprintf("%.3f", res_logistic$adjusted.p), star))
}

ZT <- symnum(round(res_logistic$deltaR2, 4), 
             c(0, 0.13, 0.26, 1), 
             symbols = c("A", "B", "C"))
JG <- symnum(res_logistic$deltaR2, 
             c(0, 0.035, 0.07, 1), 
             symbols = c("A","B", "C"))
  
tab <- cbind(sprintf("%.3f", res_logistic$Logistik), pval, 
             sprintf("%.3f", res_logistic$deltaR2), ZT, JG)
rownames(tab) <- params$itemNames

if (dim(tab)[2] == 6){
  colnames(tab) <- c("Stat.", "P-value", "Adj. p-value", "R^2", "ZT", "JG")

  kable(tab, 
        digits = c(0, 3, 3, 0, 0, 0, 0), 
        align = 'c')
} else {
  colnames(tab) <- c("Stat.", "P-value", "R^2", "ZT", "JG")

  kable(tab, 
        digits = c(0, 3, 0, 0, 0, 0),
        align = 'c')
}
```

```{r DIF-section-logistic-part2, eval = eval_logistic, echo = FALSE, results = 'asis'}
cat("Signif. codes: 0 '&#42;&#42;&#42;' 0.001 '&#42;&#42;' 0.01 '&#42;' 0.05 '.' 0.1 '&nbsp;' 1 

Effect size is based on Nagelkerke's R^2. 

'A' means negligible, 'B' moderate and 'C' large effect size 

The thresholds are: 

Zumbo & Thomas (ZT): 0 'A' 0.13 'B' 0.26 'C' 1 

Jodoin & Gierl (JG): 0 'A' 0.035 'B' 0.07 'C' 1.
")
```

```{r DIF-section-logistic-part3, eval = eval_logistic_plot, echo = FALSE, results = 'asis'}
cat("
### Characteristic curves of DIF items

Plots are based on DIF logistic procedure without any correction method.
Points represent a proportion of a correct answer with respect to standardized total score. Their size is determined by count of the respondents who achieved a given level of standardized total score.
")
```

```{r DIF-section-logistic-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300, eval = eval_logistic_plot}
for (i in 1:length(params$DIF_logistic_print$DIFitems)){
  print(params$DIF_logistic_plot[[i]] + 
          theme(text = element_text(size = 16), 
                plot.title = element_text(size = 16, face = "bold")) + 
          theme(legend.box.just = "top",
                legend.justification = c("left", "top"),
                legend.position = c(0, 1), 
                legend.box = "horizontal", 
                legend.box.margin = margin(3, 3, 3, 3)))
}
```

<!-- DIF SECTION - DDF -->

```{r DIF-section-DDF-settings, echo = FALSE, results = 'asis'}
eval_DDF <- (params$multiCheck & params$isGroupPresent)
eval_DDF_plot <- eval_DDF & !is.null(params$DDF_multinomial_plot[[1]])

if (eval_DDF) res_DDF <- params$DDF_multinomial_print
```

```{r DIF-section-DDF-part1, eval = eval_DDF, echo = FALSE, results = 'asis'}
txt1 <- switch(res_DDF$type, 
               "both" = "Both types were of DDF tested.", 
               "udif" = "Uniform DDF was tested only.", 
               "nudif" = "Non-uniform DDF was tested only.")
txt2 <- ifelse(res_DDF$purification, 
                     "Item purification was used.", 
                     "Item purification was not applied.")
txt3 <- ifelse(is.null(res_DDF$p.adjust.method), 
               "No p-value adjustment for multiple comparisons was used.", 
               paste(switch(res_DDF$p.adjust.method, 
                            "holm" = "Holm's", 
                            "hochberg" = "Hochberg's", 
                            "hommel" = "Hommel's", 
                            "bonferroni" = "Bonferroni", 
                            "BH" = "Benjamini-Hochberg", 
                            "BY" = "Benjamini-Yekutieli", 
                            "fdr" = "Benjamini-Hochberg", 
                            "none" = "No"), 
                     "p-value adjustment for multiple comparisons was used."))
txt4 <- ifelse(any(res_DDF$DDFitems == "No DDF item detected"), 
               "**No DDF item was detected. **  ", 
               paste0("**Items detected as DDF: ", paste(params$itemNames[res_DDF$DDFitems], collapse = ", "), ". **  "))

cat(paste("
## DDF detection using multinomial regression method

Differential Distractor Functioning (DDF) occurs when respondents from different groups but with the same knowledge have a different probability of selecting at least one distractor choice. DDF is examined here by a  multinomial log-linear regression model with Z-score and group membership as covariates. \n

### Summary table \n",

txt1, txt2, txt3, "\n\n",
txt4))
```

```{r DIF-section-DDF-table, echo = FALSE, eval = eval_DDF}
if (res_DDF$p.adjust.method == "none"){
      pval <- res_DDF$pval
      star <- ifelse(pval < 0.001, "***", 
                 ifelse(pval < 0.01, "**",
                        ifelse(pval < 0.05, "*",
                               ifelse(pval < 0.1, ".", ""))))
      pval <- paste(sprintf("%.3f", res_DDF$pval), star, sep = "")
    } else {
      star <- ifelse(res_DDF$adj.pval < 0.001, "***", 
                 ifelse(res_DDF$adj.pval < 0.01, "**",
                        ifelse(res_DDF$adj.pval < 0.05, "*",
                               ifelse(res_DDF$adj.pval < 0.1, ".", ""))))
      apval <- paste(sprintf("%.3f", res_DDF$adj.pval), star, sep = "")
      pval <- cbind(sprintf("%.3f", res_DDF$pval), apval)
    }
    
tab <- cbind(sprintf("%.3f", res_DDF$Sval), pval)

if (ncol(tab) == 2){
      colnames(tab) <- c("Stat.", "P-value")
      rownames(tab) <- params$itemNames
      kable(tab, align = 'c')
} else {
      colnames(tab) <- c("Stat.", "P-value", "Adj. p-value")
      rownames(tab) <- params$itemNames
      kable(tab, align = 'c')
}
```

```{r DIF-section-DDF-part2, eval = eval_DDF, echo = FALSE, results = 'asis'}
cat("
Signif. codes: 0 '&#42;&#42;&#42;' 0.001 '&#42;&#42;' 0.01 '&#42;' 0.05 '.' 0.1 '&nbsp;' 1  
")
```

```{r DIF-section-DDF-part3, eval = eval_DDF_plot, echo = FALSE, results = 'asis'}
cat(
"
### Characteristic curves of DDF items

Points represent a proportion of a selected answer with respect to standardized total score. Their size is determined by count of the respondents who achieved a given level of standardized total score and who selected a given option with respect to the group membership. 
")
```

```{r DIF-section-DDF-plot, echo = FALSE, fig.width = 11, fig.height = 4, fig.align = 'center', dpi = 300, eval = eval_DDF_plot}
for (i in 1:length(params$DDF_multinomial_print$DDFitems)){
  print(params$DDF_multinomial_plot[[i]] + 
          theme(text = element_text(size = 16),
                plot.title = element_text(size = 16, face = "bold")))
}
```

```{r DIF-section-closing, eval = show_DIF, echo = FALSE, results = 'asis'}
cat("**************")
```

<!-- SESSION INFO -->

```{r session-info-settings, echo = FALSE, results = 'asis'}
show_sessioninfo <- (params$sessionInfo == 'yes')
```

```{r session-info-part1, echo = FALSE, results = 'asis', eval = show_sessioninfo}
cat("
# Session info

Session info provides information about settings of the R console and used packages and their versions which were used for the analysis.
")
```

```{r session-info-print, echo = FALSE, results = "markup", eval = show_sessioninfo}
sessionInfo()  
```

```{r session-info-closing, eval = show_sessioninfo, echo = FALSE, results = 'asis'}
cat("**************")
```
