Readme.Rmd

---
output: github_document
---

# Exametrika

The Exametrika package is designed for test data engineering and corresponds to the text by Shojima (2022). Test data engineering involves analyzing test response patterns to classify item difficulty and respondent ranking. This analysis includes the following test theory models.

+ Classical Test Theory
+ Item response theory: 2PL, 3PL, 4PL model
+ Latent Class Analysis
+ Latent Rank Analysis
+ Biclustering and Ranklustering
+ Infinite Relational Model for optimal number of classes and fields
+ Bayesian Network Analysis
+ Structure Learning for Bayesian Network Analysis by Genetic Algorithm
+ Local Dependence Latent Rank Analysis
+ Structure Learning for LDLRA by PBIL
+ Local Dependence Biclustering
+ Biclister Network Model

Exametrika is originally implemented and published as a Mathematica and Excel Add-in. Please refer to the following website for more information.

[Accompanying Website for Test Data Engineering](http://shojima.starfree.jp/tde/)

[news](NEWS.md)

## Installation

You can install the development version of Exametrika from [GitHub](https://github.com/) with:

``` r
# install.packages("devtools")
devtools::install_github("kosugitti/Exametrika")
```

## Usage

```{r library}
library(Exametrika)
```

## Example and Sample Data

This package includes the same sample data that is distributed on the original site. The number of test-takers is represented by \( S \), and the number of test items is represented by \( J \). The data is named in a format like JxxSxxx.

#### Exametrika Data Format

Exametrika conducts analysis on data matrices composed solely of 0 or 1. In this matrix, 0 represents an incorrect answer, and 1 indicates a correct answer.For more details, refer to Shojima (2022) regarding the format.

#### Data Matrix Format

- The data matrix must be in a matrix or data frame format.
- It's permissible to include NA as missing values.
- Specific values (for instance, -99) can be designated as missing values.
- Along with the data, you may also provide a missing value index matrix that discerns the presence or absence of missing values.

#### Examinee ID Column

It is possible to incorporate a column of examinee IDs into the data matrix you provide. By default, the first column is presumed to be the examinee ID column.

#### Item Weight Vector

Item weights can be specified by the item weight vector, w. If not given, all elements' weights are set to 1.

#### Item label Vector

Data column names (colnames) are available as item labels. If not specified, a sequential number is assigned.

#### Data Formatting Function

Before any analysis, the dataFormat function decomposes the provided data into an ID vector,Item label vector, data matrix U, missing value index matrix Z, and item weight vector w.

## Test Statistics

```{r Ts}
TestStatistics(J15S500)
```

## ItemStatistics

```{r Is}
ItemStatistics(J15S500)
```

## CTT Example

```{r CTT}
CTT(J15S500)
```

## IRT Example

The IRT function estimates the number of parameters using a logistic model, which can be specified using the `model` option. It supports 2PL, 3PL, and 4PL models.

```{r}
result.IRT <- IRT(J15S500, model = 3)
result.IRT
```

The estimated population of subjects is included in the returned object.
```{r}
head(result.IRT$ability)
```


The plots offer options for Item Characteristic Curves (ICC), Item Information Curves (IIC), and Test Information Curves (TIC), which can be specified through options. Items can be specified using the `items` argument, and if not specified, plots will be drawn for all items. The number of rows and columns for dividing the plotting area can be specified using `nr` and `nc`, respectively.

```{r}
plot(result.IRT, type = "ICC", items = 1:6, nc = 2, nr = 3)
plot(result.IRT, type = "IIC", items = 1:6, nc = 2, nr = 3)
plot(result.IRT, type = "TIC")
```

## LCA Example

Latent Class Analysis requires specifying the dataset and the number of classes.

```{r}
LCA(J15S500, ncls = 5)
```

The returned object contains the Class Membership Matrix, which indicates which latent class each subject belongs to. The Estimate includes the one with the highest membership probability.

```{r}
result.LCA <- LCA(J15S500, ncls = 5)
head(result.LCA$Students)
```

The plots offer options for IRP, CMP, TRP, and LCD. For more details on each, please refer to Shojima (2022).

```{r}
plot(result.LCA, type = "IRP", items = 1:6, nc = 2, nr = 3)
plot(result.LCA, type = "CMP", students = 1:9, nc = 3, nr = 3)
plot(result.LCA, type = "TRP")
plot(result.LCA, type = "LCD")
```

## LRA Example

Latent Class Analysis requires specifying the dataset and the number of classes.

```{r}
LRA(J15S500, nrank = 6)
```

The estimated subject rank membership probabilities and plots are almost the same as those in LCA (Latent Class Analysis). Since a ranking is assumed for the latent classes, rank-up odds and rank-down odds are calculated.

```{r}
result.LRA <- LRA(J15S500, nrank = 6)
head(result.LRA$Students)
```

```{r}
plot(result.LRA, type = "IRP", items = 1:6, nc = 2, nr = 3)
plot(result.LRA, type = "RMP", students = 1:9, nc = 3, nr = 3)
plot(result.LRA, type = "TRP")
plot(result.LRA, type = "LRD")
```

## Biclustering Example

Biclustering and Ranklustering algorithms are almost the same, differing only in whether they include a filtering matrix or not. The difference is specified using the `method` option in the `Biclustering()` function. For more details, please refer to the help documentation.

```{r}
Biclustering(J35S515, nfld = 5, ncls = 6, method = "B")
```


```{r}
result.Ranklusteing <- Biclustering(J35S515, nfld = 5, ncls = 6, method = "R")
plot(result.Ranklusteing, type = "Array")
plot(result.Ranklusteing, type = "FRP", nc = 2, nr = 3)
plot(result.Ranklusteing, type = "RMP", students = 1:9, nc = 3, nr = 3)
plot(result.Ranklusteing, type = "LRD")
```

To find the optimal number of classes and the optimal number of fields, the Infinite Relational Model is available.

```{r}
result.IRM <- IRM(J35S515, gamma_c = 1, gamma_f = 1, verbose = TRUE)
plot(result.IRM, type = "Array")
plot(result.IRM, type = "FRP", nc = 3)
plot(result.IRM, type = "TRP")
```

Additionally, supplementary notes on the derivation of the Infinite Relational Model with Chinese restaurant process is [here](https://kosugitti.github.io/kosugitti10/notes/IRM_memo.pdf).

## Bayesian Network Model

The Bayesian network model is a model that represents the conditional probabilities between items in a network format based on the pass rates of the items. By providing a Directed Acyclic Graph (DAG) between items externally, it calculates the conditional probabilities based on the specified graph. The igraph package is used for the analysis and representation of the network.

There are three ways to specify the graph. You can either pass a matrix-type DAG to the argument adj_matrix, pass a DAG described in a CSV file to the argument adj_file, or pass a graph-type object g used in the igraph package to the argument g.

The methods to create the matrix-type adj_matrix and the graph object g are as follows:

```{r igraph}
library(igraph)
DAG <-
  matrix(
    c(
      "Item01", "Item02",
      "Item02", "Item03",
      "Item02", "Item04",
      "Item03", "Item05",
      "Item04", "Item05"
    ),
    ncol = 2, byrow = T
  )
## graph object
g <- igraph::graph_from_data_frame(DAG)
g
## Adj mmatrix
adj_mat <- as.matrix(igraph::get.adjacency(g))
adj_mat
```

A CSV file with the same information as the graph above in the following format. The first line contains column names (headers) and will not be read as data.

```{r echo=FALSE}
# save csv file
write.csv(DAG, "develop/DAG_file.csv", row.names = FALSE, quote = TRUE)

lines <- readLines("develop/DAG_file.csv")
for (line in lines) {
  cat(line)
  cat("\n")
}
```

While only one specification is sufficient, if multiple specifications are provided, they will be prioritized in the order of file, matrix, and graph object.

An example of executing BNM by providing a graph structure (DAG) is as follows:

```{r BNM}
result.BNM <- BNM(J5S10, adj_matrix = adj_mat)
result.BNM
```

### Structure Larning for Bayesian network with GA

The function searches for a DAG suitable for the data using a genetic algorithm. A best DAG is not necessarily identified. Instead of exploring all combinations of nodes and edges, only the space topologically sorted by the pass rate, namely the upper triangular matrix of the adjacency matrix, is explored. For interpretability, the number of parent nodes should be limited. A null model is not proposed. Utilize the content of the items and the experience of the questioner to aid in interpreting the results. For more details, please refer to Section 8.5 of the text(Shojima,2022).

Please note that the GA may take a considerable amount of time, depending on the number of items and the size of the population.

```{r GAsimple}
StrLearningGA_BNM(J5S10,
  population = 20, Rs = 0.5, Rm = 0.002, maxParents = 2,
  maxGeneration = 100, crossover = 2, elitism = 2
)
```

The method of Population-based incremental learning proposed by Fukuda (2014) can also be used for learning. This method has several variations for estimating the optimal adjacency matrix at the end, which can be specified as options. See help or text Section 8.5.2.

```{r PBIL}
StrLearningPBIL_BNM(J5S10,
  population = 20, Rs = 0.5, Rm = 0.005, maxParents = 2,
  alpha = 0.05, estimate = 4
)
```

## Local Dependent Latent Rank Analysis

LD-LRA is an analysis that combines LRA and BNM, and it is used to analyze the network structure among items in the latent rank. In this function, structural learning is not performed, so you need to provide item graphs for each rank as separate files.

For each class, it is necessary to specify a graph, and there are three ways to do so. You can either pass a matrix-type DAG for each class or a list of graph-type objects used in the igraph package to the arguments adj_list or g_list, respectively, or you can provide a DAG described in a CSV file. The way to specify it in a CSV file is as follows.

```{r read DAG file}
DAG_dat <- matrix(c(
  "From", "To", "Rank",
  "Item01", "Item02", 1,
  "Item04", "Item05", 1,
  "Item01", "Item02", 2,
  "Item02", "Item03", 2,
  "Item04", "Item05", 2,
  "Item08", "Item09", 2,
  "Item08", "Item10", 2,
  "Item09", "Item10", 2,
  "Item08", "Item11", 2,
  "Item01", "Item02", 3,
  "Item02", "Item03", 3,
  "Item04", "Item05", 3,
  "Item08", "Item09", 3,
  "Item08", "Item10", 3,
  "Item09", "Item10", 3,
  "Item08", "Item11", 3,
  "Item02", "Item03", 4,
  "Item04", "Item06", 4,
  "Item04", "Item07", 4,
  "Item05", "Item06", 4,
  "Item05", "Item07", 4,
  "Item08", "Item10", 4,
  "Item08", "Item11", 4,
  "Item09", "Item11", 4,
  "Item02", "Item03", 5,
  "Item04", "Item06", 5,
  "Item04", "Item07", 5,
  "Item05", "Item06", 5,
  "Item05", "Item07", 5,
  "Item09", "Item11", 5,
  "Item10", "Item11", 5,
  "Item10", "Item12", 5
), ncol = 3, byrow = TRUE)

# save csv file
write.csv(DAG_dat, "develop/DAG_file.csv", row.names = FALSE, quote = TRUE)
```

Here, it is shown an example of specifying with matrix-type and graph objects using the aforementioned CSV file. While only one specification is sufficient, if multiple specifications are provided, they will be prioritized in the order of file, matrix, and graph object.

```{r convert_style}
g_csv <- read.csv("develop/DAG_file.csv")
colnames(g_csv) <- c("From", "To", "Rank")
adj_list <- list()
g_list <- list()
for (i in 1:5) {
  adj_R <- g_csv[g_csv$Rank == i, 1:2]
  g_tmp <- igraph::graph_from_data_frame(adj_R)
  adj_tmp <- igraph::get.adjacency(g_tmp)
  g_list[[i]] <- g_tmp
  adj_list[[i]] <- adj_tmp
}
## Example of graph list
g_list
```
```{r example of adj list}
### Example of adj list
adj_list
```

The example of running the LDLRA function using this CSV file would look like this.

```{r LDLRA}
result.LDLRA <- LDLRA(J12S5000,
  ncls = 5,
  adj_file = "develop/DAG_file.csv"
)
result.LDLRA
```

Of course, it also supports various types of plots.
```{r LDLRA plot}
plot(result.LDLRA, type = "IRP", nc = 4, nr = 3)
plot(result.LDLRA, type = "TRP")
plot(result.LDLRA, type = "LRD")
```


### Structure Larning for LDLRA with GA(PBIL)

You can learn item-interaction graphs for each rank using the PBIL algorithm. In addition to various options, the learning process requires a very long computation time. It's also important to note that the result is merely one of the feasible solutions, and it's not necessarily the optimal solution.

```{r LDLRA_PBIL}
result.LDLRA.PBIL <- StrLearningPBIL_LDLRA(J35S515,
  seed = 123,
  ncls = 5,
  method = "R",
  elitism = 1,
  successiveLimit = 15
)
result.LDLRA.PBIL
```

## Local Dependence Biclustring

Latent dependece Biclustering, which incorporates biclustering and a Bayesian network model.

To execute this, you need to determine the number of ranks, the correspondence between items and fields, and the network structure for each rank in advance. For confirmatory Rankclustering, the correspondence between items and fields is provided in the form of a matrix or vector.

```{r settings for LDB, include=FALSE}
fieldFile <- "develop/mtmk15forVer13/FixFieldLDB.csv"
ncls <- 5
edgeFile <- "develop/mtmk15forVer13/EdgesLDB.csv"
FieldData <- read.csv(fieldFile)
conf <- FieldData[, 2]
adj_file <- edgeFile
```

```{r LDB conf}
lines <- readLines(fieldFile)
for (line in lines) {
  cat(line)
  cat("\n")
}
```
```{r}
FieldData <- read.csv(fieldFile)
conf <- FieldData[, 2]
conf
```

The network structure for each rank can be specified using a matrix, a graph object, or a CSV file. Here is an example using a CSV file for specification.

```{r LDB Graph }
lines <- readLines(edgeFile)
for (line in lines) {
  cat(line)
  cat("\n")
}
```

Additionally, as mentioned in the text (Shojima, 2022), it is often the case that seeking the network structure exploratively does not yield appropriate results, so it has not been implemented.
```{r}
result.LDB <- LDB(U = J35S515, ncls = 5, conf = conf, adj_file = edgeFile)
result.LDB
```

Of course, it also supports various types of plots.

```{r LDB plot}
plot(result.LDB, type = "Array")
plot(result.LDB, type = "TRP")
plot(result.LDB, type = "LRD")
plot(result.LDB, type = "RMP", students = 1:9, nc = 3, nr = 3)
plot(result.LDB, type = "FRP", nc = 3, nr = 2)
```

In this model, you can draw a Field PIRP Profile that visualizes the correct answer count for each rank and each field.

```{r LDB FieldPIRP}
plot(result.LDB, type = "FieldPIRP")
```

## Bicluster Network Model

Bicluster Network Model: BINET is a model that combines the Bayesian 
network model and Biclustering. BINET is very similar to LDB and LDR.

The most significant difference is that in LDB, the nodes represent the fields, whereas in BINET, they represent the class. BINET explores the local dependency structure among latent classes at each latent field, where each field is a locus.

To execute this analysis, in addition to the dataset, the same field correspondence file used during exploratory Biclustering is required, as well as an adjacency matrix between classes. 

```{r BINET setting}
fieldFile <- "develop/mtmk15forVer13/FixFieldBINET.csv"
edgeFile <- "develop/mtmk15forVer13/EdgesBINET.csv"
FieldData <- read.csv(fieldFile)
conf <- FieldData[, 2]
```

```{r}
conf <- FieldData[, 2]
conf
```

The adjacency matrix between classes can be provided in a CSV file like as follow:

```{r BINET conf}
lines <- readLines(edgeFile)
for (line in lines) {
  cat(line)
  cat("\n")
}
```

For verification purposes, you also need to specify the number of classes and fields. Once you input these along with the dataset into the function, the analysis will be executed.


```{r BINET}
result.BINET <- BINET(
  U = J35S515,
  ncls = 13, nfld = 12,
  conf = conf, adj_file = edgeFile
)
result.BINET
```

Of course, it also supports various types of plots.

```{r BINET plot}
plot(result.BINET, type = "Array")
plot(result.BINET, type = "TRP")
plot(result.BINET, type = "LRD")
plot(result.BINET, type = "RMP", students = 1:9, nc = 3, nr = 3)
plot(result.BINET, type = "FRP", nc = 3, nr = 2)
```

LDPSR plot shows all Passing Student Rates for all locallyu dependent classes compared with their respective parents.


```{r LDPSRplot}
plot(result.BINET, type = "LDPSR", nc = 3, nr = 2)
```


## Table of Model and Plotting Option Correspondence

| model/type     | IIC | ICC | TIC | IRP | FRP | TRP | LCD/LRD | CMP/RMP | ARRAY | FieldPIRP | LDPSR |
|----------------|-----|-----|-----|-----|-----|-----|---------|---------|-------|-----------|-------|
| IRT            |  ◯  |  ◯  |  ◯  |     |     |     |         |         |       |           |       |
| LCA            |     |     |     |  ◯  |  ◯  |  ◯  |    ◯    |    ◯    |       |           |       |
| LRA            |     |     |     |  ◯  |  ◯  |  ◯  |    ◯    |    ◯    |       |           |       |
| Biclustering   |     |     |     |  ◯  |  ◯  |  ◯  |    ◯    |    ◯    |   ◯   |           |       |
| IRM            |     |     |     |     |  ◯  |  ◯  |         |         |   ◯   |           |       |
| LDLRA          |     |     |     |  ◯  |     |     |    ◯    |    ◯    |       |           |       |
| LDB            |     |     |     |     |  ◯  |  ◯  |    ◯    |    ◯    |   ◯   |     ◯     |       |
| BINET          |     |     |     |     |  ◯  |  ◯  |    ◯    |    ◯    |   ◯   |           |   ◯   |


## Reference

Shojima, Kojiro (2022) Test Data Engineering: Latent Rank Analysis, Biclustering, and Bayesian Network (Behaviormetrics: Quantitative Approaches to Human Behavior, 13),Springer, [https://amzn.to/42eNArd](https://amzn.to/42eNArd)


## Citation

[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.11529926.svg)](https://doi.org/10.5281/zenodo.11529926)