10-Item-response-theory.Rmd

Item Response Theory
======


This tutorial is adapted from 2 articles:  

 * [MultiLCIRT: An R package for multidimensional latent class item response models](https://arxiv.org/pdf/1210.5267.pdf)  
 * [Variable weighting via multidimensional IRT models in Composite Indicators construction](https://www.researchgate.net/profile/Simone_Del_Sarto/publication/303864920_Variable_weighting_via_multidimensional_IRT_models_in_Composite_Indicators_construction/links/575964c308aed884620ad4f0/Variable-weighting-via-multidimensional-IRT-models-in-Composite-Indicators-construction.pdf?origin=publication_detail)  
 
 Thanks to [Simone Del Sarto](delsarto@stat.unipg.it) from University of Perugia, Italy for the guidance on this document.
 
## Vulnerability as a "_latent_" variable

A major task in the measurement of complex and latent household vulnerability dimensions phenomena such as deprivation (or _Basic Needs_), wealth (or _Coping capacity_) and ability (or _Well-being_) consists of __summarising information__ available from a set of dichotomous or ordinal questions (also referred as _items_) from an household survey questionnaire. For example:

 * Deprivation (or _Basic Needs_),  is assessed by collecting data on the extent to which households possess certain commodities, engage in certain activities or are subject to financial pressures. The responses of individuals to deprivation item questionnaires constitute the manifest or observed indicators. (cf. [Item response theory and the measurement of deprivation, Szeles and Fusco, 2013](https://www.researchgate.net/profile/Raileanu_Monica/publication/46452347_Item_response_theory_and_the_measurement_of_deprivation_Evidence_from_PSELL-3/links/00b4952a0e898e8b0d000000/Item-response-theory-and-the-measurement-of-deprivation-Evidence-from-PSELL-3.pdf?origin=publication_detail))

 * Wealth (or _Coping capacity_) is a latent concept to be derived from observable assets, and measures of wealth are calculated from a set of assets a household possesses (cf. [Measuring Household Wealth with Latent Trait Modelling: An Application to Malawian DHS Data, Vandemoortele, 2014](https://link.springer.com/article/10.1007/s11205-013-0447-z)).


## Building composite indicators with IRT 

Item Response Theory (IRT) are specific type of statistical modeling that allows to address this complex challenge of building data-driven vulnerability scores without the use of proxy indicators. The assumptions that IRT approaches allows for are closer to the operational reality and needs:

   *  Latent phenomenon → not directly observable through a proxy, meaning that vulnerability is __not directly measurable__ . 
   
   *  Latent variable   →  assumed to have a __discrete distribution__ rather than continuous vulnerability scales.   
   
   *  Phenomenon manifestation → __response pattern__ to the items of a questionnaire rather to profile of respondents.
   
   *  Summary of the phenomenon  →  __composite indicators__ rather than predicted scores.

Item Response Theory (IRT) (also known as _latent trait theory_, _strong true score theory_, or _modern mental test theory_) is a paradigm for the design, analysis, and scoring of tests, questionnaires, and similar instruments measuring abilities, attitudes, or other latent variables. It is a theory of testing based on the relationship between individuals' performances on a test item and the test takers' levels of performance on an overall measure of the ability that item was designed to measure. IRT models the response of each respondant of a given ability to each item in the test. 

The term `item`, here is generic, covering all kinds of informative items. They might be multiple choice questions that have incorrect and correct responses, but are also commonly statements on questionnaires that allow respondents to indicate level of agreement (a rating or Likert scale), or patient symptoms scored as present/absent, or diagnostic information in complex systems. 

The main advantage of the `discrete` approach consists in the possibility of clustering the statistical units into homogeneous groups (latent classes): units (people, refugees and so on) in the same class share very similar characteristics in terms of the investigated (multidimensional) latent trait. Hence, this model would also allow us to classify refugees in similar groups. Such model can support prioritization needs as it assigns a score to each group, so that it is possible to highlight the most/least vulnerable groups of refugees and/or sketch several refugees’ profiles according to the various dimensions of it. 

Composite indicators can be built using IRT using the following __procedure__:

   *  Step 1 : Assess dimensionality within the variables of the data-set → Multidimensional Latent Class IRT model + clustering algorithm
   *  Step 2 : Adoption the best model for assigning weights to the test items. This steps also allows to account for dimensionality (clusters) within respondents.
   *  Step 3 : Item aggregation using the weights obtained at the previous step → construction of a composite indicator for each dimension


### Step 0: Data Preparation

Data are loaded in a data frame and converted in a matrix, corresponding to the unit-by-unit response configurations.

To make the further estimation of the proposed models and clustering of items faster perform the analyses after aggregating the original records which correspond to the same response pattern so as to obtain a matrix with a record for each distinct response configuration (rather than for each statistical unit).

The function function `aggr_data` is then applied to this matrix. Output from function is: 

   *  `data_dis`: matrix of distinct configurations →  S (used after in model estimation);  
   *  `freq` : vector of corresponding frequencies → yv (used after in model estimation);  
   *  `label`: index of each original response configuration among the distinct ones.  

```{r setup, include=FALSE, warning = FALSE, cache=TRUE}
knitr::opts_chunk$set(echo = TRUE)
using <- function(...) {
    libs <- unlist(list(...))
    req <- unlist(lapply(libs,require,character.only = TRUE))
    need <- libs[req == FALSE]
    if (length(need) > 0) { 
        install.packages(need)
        lapply(need,require,character.only = TRUE)
    }
}
using('tidyverse','gganimate','gghighlight','ggpubr', 'dplyr', 'tidyr',  'ggplot2',  'ggalt', 'forcats', 'R.utils', 'png', 'grid', 'ggpubr', 'scales',  'markdown', 'pander', 'stargazer', 'kableExtra', 'MultiLCIRT', 'MLCIRTwithin', 'mirt', 'tabplot')
options(scipen = 999) # turn-off scientific notation like 1e+48
# #install.packages("stargazer")
# library(stargazer)
# library(kableExtra)
# 
# ### install.packages("MultiLCIRT")
# library(MultiLCIRT)
# 
# #install.packages("MLCIRTwithin")
# library(MLCIRTwithin)
# 
# ### install.packages("mirt")
# library(mirt)
# 
# library(ggplot2)
unhcr_style <- function() {
  font <- "Lato"
  ggplot2::theme(
    
#This sets the font, size, type and colour of text for the chart's title
  plot.title = ggplot2::element_text(family = font, size = 16, face = "bold", color = "#222222"),
#This sets the font, size, type and colour of text for the chart's subtitle,  as well as setting a margin between the title and the subtitle
  plot.subtitle = ggplot2::element_text(family = font, size = 12, margin = ggplot2::margin(9,0,9,0)),
  plot.caption = ggplot2::element_blank(),
#This sets the position and alignment of the legend, removes a title and backround for it and sets the requirements for any text within the legend. The legend may often need some more manual tweaking when it comes to its exact position based on the plot coordinates.
  legend.position = "top",
  legend.text.align = 0,
  legend.background = ggplot2::element_blank(),
  legend.title = ggplot2::element_blank(),
  legend.key = ggplot2::element_blank(),
  legend.text = ggplot2::element_text(family = font, size = 11, color = "#222222"),
#This sets the text font, size and colour for the axis test, as well as setting the margins and removes lines and ticks. In some cases, axis lines and axis ticks are things we would want to have in the chart
  axis.title = ggplot2::element_blank(),
  axis.text = ggplot2::element_text(family = font, size = 11, color = "#222222"),
  axis.text.x = ggplot2::element_text(margin = ggplot2::margin(5, b = 10)),
  axis.ticks = ggplot2::element_blank(),
  axis.line = ggplot2::element_blank(),
#This removes all minor gridlines and adds major y gridlines. In many cases you will want to change this to remove y gridlines and add x gridlines. 
  panel.grid.minor = ggplot2::element_blank(),
  panel.grid.major.y = ggplot2::element_line(color = "#cbcbcb"),
  panel.grid.major.x = ggplot2::element_blank(),
#This sets the panel background as blank, removing the standard grey ggplot background colour from the plot
  panel.background = ggplot2::element_blank(),
#This sets the panel background for facet-wrapped plots to white, removing the standard grey ggplot background colour and sets the title size of the facet-wrap title to font size 22
  strip.background = ggplot2::element_rect(fill = "white"),
  strip.text = ggplot2::element_text(size  = 11,  hjust = 0)
  )
}
#Left align text
left_align <- function(plot_name, pieces){
  grob <- ggplot2::ggplotGrob(plot_name)
  n <- length(pieces)
  grob$layout$l[grob$layout$name %in% pieces] <- 2
  return(grob)
}
## a little help function to better format numbers
format_si <- function(...) {
  function(x) {
    limits <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12,
                1e-9,  1e-6,  1e-3,  1e0,   1e3,
                1e6,   1e9,   1e12,  1e15,  1e18,
                1e21,  1e24)
    prefix <- c("y",   "z",   "a",   "f",   "p",
                "n",   "",   "m",   " ",   "k",
                "M",   "G",   "T",   "P",   "E",
                "Z",   "Y")
    # Vector with array indices according to position in intervals
    i <- findInterval(abs(x), limits)
    # Set prefix to " " for very small values < 1e-24
    i <- ifelse(i == 0, which(limits == 1e0), i)
    paste(format(round(x/limits[i], 1),
                 trim = TRUE, scientific = FALSE, ...),
          prefix[i])
  }
}
```

```{r data , warning = FALSE, message = FALSE, cache=TRUE}
## Data from Bock & Lieberman (1970); contains 5 dichotomously scored items obtained from the Law School Admissions Test, section 7.
# data(LSAT7)
# head(LSAT7)
# S <- as.matrix(LSAT7[,1:5])
# yv <- as.vector(LSAT7[,6])

## dataset contains the responses of a sample of 1510 examinees to 12 binary items on Mathematics. It has been extrapolated from a larger dataset collected in 1996 by the Educational Testing Service within the National Assessment of Educational Progress (NAEP) project.
#data(naep)

##  hads -> Dataset about measurement of anxiety and depression in oncological patients. A data frame with 201 observations on 14 items, measuring depression or anxiety
datadescription <- "Dataset about measurement of anxiety and depression in oncological patients"
data(hads)
# Convert to atrix
data <- as.matrix(hads)

# Prepare data
out_aggr <- aggr_data(data)
S <- out_aggr$data_dis
yv <- out_aggr$freq
labeldata <- out_aggr$label

```

The data-set in this reports contains the responses of a sample of `r nrow(data)` / `r nrow(hads)` records to `r ncol(data)` / `r nrow(S)` binary items.


[Tableplot](https://cran.r-project.org/web/packages/tabplot/vignettes/tabplot-vignette.html) is a powerful visualization method to explore and analyse large multivariate datasets. The representation below are used to help in the understanding of the data,  explore the relationships between the variables, discover strange data patterns, and  check the occurrence and selectivity of missing values.

```{r tableplot,  warning = FALSE, message=FALSE, cache=TRUE }

# tableplot(hads,
#           # sortCol = unconditionnal2, 
#            nBins = 200 , 
#            fontsize = 10,
#            legend.lines = 8, 
#            title = "Overview of all indicators", 
#            fontsize.title = 13)

```


###  Step 1: Clustering algorithm / Dimensionality assessment

The first steps aims at grouping variables measuring the same latent construct in the same cluster

Investigating the best number of latent classes (parameter `k`) is done through Hierarchical clustering of items based on the Rasch model. The function `class_item` is used, it can take a long execution time. The classification is based on a sequence of likelihood ratio tests between pairs of multidimensional models suitably formulated. 


```{r clustering, message=TRUE, warning=TRUE, cache=TRUE, results = "hide"}
out <- class_item(S,
                  yv,
                  k = 3, # latentclasses (partial output is omitted)
                  link = 1, # it is the same with link = 2
                  disc = 0)
```

The dendrogram below highlights three dimensions composing the multidimensional construct. Considering the item contents, the assessed dimension can be interpreted from a narrative point of view.

```{r dendogram, warning = FALSE}

plot(out$dend)

```

The following tables allows to interpret the clustering: 

   *  The first two columns (entitled `items`) indicate items or groups collapsed at each step of the clustering procedure, 
   
   *  the third column(`deviance`) reports the corresponding LR statistic, 
   
   *  the fourth column(`df`)reports the number of degrees of freedom, 
   
   *  the fifth (`p-value`) reports the p-value, and 
   
   *  the last column (`newgroup`) contains the new group of items that is formed at each step.


```{r clusteringout, warning = FALSE }
# , results = "asis"
summary(out)

```


### Step 2: Model selection

The second step is a selection procedure based on estimation of ordinal polytomous multidimensional LCIRT model. The formulation of a specific model in the class of multidimensional Latent Class (LC) IRT models univocally depends on: 

 * number of latent classes `k`, 
 * adopted parameterization in terms of `link` function  , 
 * constraints on the item parameters , and 
 * number of latent dimensions (`s`) and the corresponding allocation of items within each dimension

Parameter estimation for multidimensional IRT models based on discreteness of latent traits is performed through function `est_multi_poly` requires the following main input:

   *  `S`: matrix of all response sequences observed at least once in the sample and listed row-by-row. Usually,S is a matrix of type data_dis obtained by applying function aggr_data to the original data. Missing responses are allowed and they are coded as NaN;

   *  `yv`: vector of the frequencies of every response configuration in S corresponding to the output freq of function aggr_data (default value is given by a vector of ones, implying a unit weight for each record in S);

   *  `k` : number of latent classes - defined in the step 1;

   *  `X` : matrix of observed covariates, having the same dimension as S (default value is  NULL, indicating the absence of covariates in the study);

   *  `start`: method of initialization of the algorithm:  0 (  =  default value) for deterministic starting values, 1 for random starting values, and 2 for arguments given bytheuser.  If start = 2 ,we also need to specify as input the initial values of weights, support points, and discriminating an ddifficulty item parameters (using additional input arguments that are set equal to  NULL otherwise);

   *  `link` : type of link function:  0 (= default value) for the standard Latent Class model (i.e., no link function is specified),1 for global logits, and 2 for local logits. In the case of dichotomous responses, it is the same to specify link =  1 or link = 2;

   *  `disc` : indicator of constraints on the discriminating item parameters: 0 (= default value) if γj= 1, j=1 ,..., r, and 1 otherwise;

   *  `difl` : indicator of constraints on the difficulty item parameters: 0 (= default value) if difficulties are free and 1 if βjx= βj+ τx

   *  `multi` : matrix with a number of rows equal to the number of dimensions and elements in each row equal to the indices of the items measuring the dimension corresponding to that row. Cases where dimensions are measured by a different number of items are allowed, and the number of columns of matrix multi corresponds to the number of items in the largest dimension.



Function `est_multi_poly` supplies the following output:

   *  `Piv`: optional object of type matrix containing the estimated weights of thelatent classes subject-by-subject  (the weights may be different across subjects in the presence of covariates);

  *  `Th`: estimated matrix of ability levels (support points) for each dimension (= row of matrix Th) and latent class (= column of matrix Th);

  *  `Bec` : estimated vector of difficulty item parameters (split in two vectors if difl = 1);

  *  `gac` : estimated vector of discriminating item parameters; if disc = 0 (Rasch-type model), all values of vector gac are constrained to 1;

  *  `fv` : vector indicating the reference item chosen for each latent dimension;

  *  `Phi` : optional object of type array containing the conditional response probabilities (see Eq.(1)) for every item and latent class. The array is made of as many matrices as the latent classes; moreover, the j -th column of each of such matrices refers to item j , where as the x-th row of each matrix refers to the x- th responsecategory(x = 0,..., lj − 1) of item j. In the case of items differing in the number of response categories, zeros are included in the corresponding cells;

  *  `Pp` : optional object of type matrix containing the posterior probabilities of belonging to latent class c (column  c of the Pp matrix), given the response configuration (row of the Pp   matrix);

  *  `lk` : log-likelihood at convergence of the EM algorithm;

  *  `np` : number of free parameters;

  *  `aic` : Akaike Information Criterion index (Akaike,1973);

  *  `bic` : Bayesian Information Criterion index (Schwarz,1978)


#### Number latent classes to consider

Checking here the optimal value of parameter `k`

```{r identifyclass, warning = FALSE, message = FALSE, results = "hide", cache=TRUE}

oneclass <- est_multi_poly(S,yv,
                       k = 1, # number of ability levels (or latent classes)
                       start = 0, # method of initialization of the algorithm (0 = deterministic, 1 = random, 2 = arguments given as input)
                       link = 0) # type of link function (0 = no link function, 1 = global logits, 2 = local logits);

twoclass <- est_multi_poly(S,yv,
                       k = 2,
                       start = 0,
                       link = 0)

threeclass <- est_multi_poly(S, yv,
                       k = 3,
                       start = 0,
                       link = 0)

fourclass <- est_multi_poly(S,yv,
                       k = 4,
                       start = 0,
                       link = 0)

```

```{r identifyclass2, warning = FALSE, , results = "asis", cache=TRUE}
stargazer::stargazer(compare_models(oneclass, twoclass, threeclass, fourclass), 
                     title = "Model Comparison",
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )

```

The model with lowest BIC is retained (i.e. `k` = 3).

#### Define parameterization on link function

To define if Global or local logit shall be used, we need to identify the optimal value of parameter `link`  (1= Glocal logit, 2 = Local logit)

```{r identifylogit, warning = FALSE, message = FALSE, results = "hide", cache=TRUE}
### Check if local or global logit to consider
global.logit <- est_multi_poly(S, yv,
                        k = 3,
                        start = 0,
                        link = 1, # Global logit
                        # indicator of constraints on the discriminating indices (0 = all equal to one, 1 = free)
                        disc = 1, 
                        # indicator of constraints on the difficulty levels (0 = free, 1 = rating scale parameterization)
                        difl = 0,
                        multi = cbind(1:ncol(S)))

local.logit <- est_multi_poly(S, yv,
                        k = 3,
                        start = 0,
                        link = 2, # Local logit
                        disc = 1,
                        difl = 0,
                        multi = cbind(1:ncol(S)))



```

```{r identifylogit2, warning = FALSE, results = "asis", cache=TRUE}

## Findout if global or local logit is better based on value of BIC
stargazer::stargazer(compare_models(global.logit, local.logit), 
                     title = "Model Comparison",
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )
```

The model with lowest BIC is retained (i.e. `link` = 1).


#### Constraints on the item parameters: Likelihood ratio 

To test between nested multidimensional LC IRT models, we can compare different multidimensional model: a restricted model against a larger multidimensional model based on a higher number of dimensions.

A typical example is testing a unidimensional model (and then the hypothesis of unidimensionality) against a bidimensional model  through function `test_dim`.


```{r likelihood,  warning = FALSE, message = FALSE, results = "hide", cache=TRUE}

## Set up the restricted model for the rest of analysis
# # Define matrix to allocate each item to one dimension
multi2 <- rbind(c(2,6,7,8,10,11,12),
                c(1,3,4,5,9,13,14))

unidimensional <- test_dim(S, yv,
                 k = 3,
                 
     # type of link function (1 = global logits, 2 = local logits); with global logits the Graded Response model results;
     # with local logits the Partial Credit results (with dichotomous responses, global logits is the same as 
     # using local logits resulting in the Rasch or the 2PL model depending on the value assigned to disc)
                 link = 1, 
     
                 # Indicator of constraints on the discriminating indices (0 = all equal to one, 1 = free)
                 disc = 1, 
     
                 # indicator of constraints on the difficulty levels (0 = free, 1 = rating scale parametrization) 
                 difl = 0,
     
                 # matrix specifying the multidimensional structure of the restricted model
                 multi0 = multi2,  
     
                 # matrix specifying the multidimensional structure of the larger model
                 multi1 = cbind(1:ncol(S)),  
     
                 # Tolerance level for checking convergence of the algorithm as relative
                 # difference between consecutive log-likelihoods
                 tol = 5*10^-4)

bidimensional <- test_dim(S, yv,
                 k = 3,
                 link = 1,
                 disc = 1,
                 difl = 0,
                 multi1 = multi2,
                 tol = 5*10^-4)
```

```{r likelihood2, warning = FALSE, cache=TRUE}

summary(unidimensional)
summary(bidimensional)

```

We will retain here the model with the smaller BIC value - Retained value for `link` = 1


#### Test of unidimensionality

Once the global logit has been chosen as the best link function, we carry on with the test of unidimensionality.  

An LR test is used to comparemodels which differ in terms of the dimensional structure, all other elements being equal (i.e. free item discriminating and difficulty parameters), that is,
 * (i) a graded response model (GRM) with an r-dimensional structure,
 * (ii) a graded response model (GRM)  with a bidimensional structure (i.e., anxiety and depression), as suggested by the structure of the questionnaire, and
 * (iii) a graded response model (GRM)  with a unidimensional structure (i.e., all the items belong to the same dimension).

For this, models are calculated with different parameters:
 * Difficulty levels (Free or Constrained, ie. constant or non-constant) `difl`.
 * Discriminating indices (Free or Constrained, ie. constant or non-constant) `disc`.


```{r unidimensionnality,  warning = FALSE, message = FALSE, results = "hide", cache=TRUE}

# Unidimensional GRM
Unidimensional.GradedResponse <- est_multi_poly(S, yv,
                         k = 3,
                         start = 0,
                         link = 1,
                         # Indicator of constraints on the discriminating indices ( 1 = free)
                         disc = 1,
                         # indicator of constraints on the difficulty levels (0 = free) 
                         difl = 0)

# Unidimensional  RS-GRM rating scale parameterization
Unidimensional.RatingScale.GradedResponse <- est_multi_poly(S, yv,
                          k = 3,
                          start = 0,
                          link = 1,
                          # Indicator of constraints on the discriminating indices ( 1 = free)
                          disc = 1,
                          # indicator of constraints on the difficulty levels (1 = rating scale parametrization)
                          difl = 1) 
# Unidimensional  1P-GRM
Unidimensional.OneParamater.GradedResponse <- est_multi_poly(S, yv,
                          k = 3,
                          start = 0,
                          link = 1,
                          # Indicator of constraints on the discriminating indices (0 = all equal to one)
                          disc = 0, 
                          # indicator of constraints on the difficulty levels (0 = free)
                          difl = 0)

# Unidimensional  1P-RS-GRM
Unidimensional.OneParamater.RatingScale.GradedResponse <- est_multi_poly(S, yv,
                          k = 3,
                          start = 0,
                          link = 1,
                          # Indicator of constraints on the discriminating indices (0 = all equal to one)
                          disc = 0,
                          # indicator of constraints on the difficulty levels (1 = rating scale parametrization)
                          difl = 1)  

```

```{r unidimensionnality2, warning = FALSE, results = "asis", cache=TRUE}

# comparison of RS-GRM and 1P-GRM with GRM
stargazer::stargazer(compare_models(Unidimensional.GradedResponse,
                                    Unidimensional.RatingScale.GradedResponse,
                                    Unidimensional.OneParamater.GradedResponse,
                                    nested = TRUE), 
                     title = "Model Comparison",
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )


stargazer::stargazer(compare_models(Unidimensional.OneParamater.GradedResponse,
                                    Unidimensional.OneParamater.RatingScale.GradedResponse,
                                    nested = TRUE), 
                     title = "Model Comparison",
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )



## Unidimensional.OneParamater.GradedResponse has the lowest BIC
stargazer::stargazer(rbind(Unidimensional.OneParamater.GradedResponse$Th,
                           prob = Unidimensional.OneParamater.GradedResponse$piv), 
                     title = "Model Comparison",
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )


```




### Step 3: Items aggregation using the weights

The discrimination parameters express the informative capacity of a variable and its relative importance within sub-dimensions of the latent variables (as the model is multidimensional): so that the weight of a variable expresses the relative importance of that variable in the dimension it belongs to.

Weights are extracted from discrimination parameters, estimated through a Multidimensional IRT model, using a `2PL` parametrisation


Two composite indicators can be calculated:

 * __Unweighted score__ (raw relative frequency of correct response) 
 * __Weighted score__, using ˆ γ j (and their transformation) for weighting the items
 


```{r finalmodel,  warning = FALSE, message = FALSE, results = "hide", cache=TRUE}
# Final model selected Unidimensional GRM
# Estimate the model (say, a bidimensional GRM with 3 latent classes)
out_GRM <- est_multi_poly(S,
                          yv,
                          link = 1,
                          k = 3,
                          multi = multi2,
                          disc = 1, # Indicator of constraints on the discriminating indices ( 1 = free)
                          difl = 0, # indicator of constraints on the difficulty levels (0 = free)
                          tol = 10^-3 )

```


#### Multidimensional scores 

gamma_hat: $\hat{\gamma}_j$ (weights) from the model selected in the step above.

```{r,  warning = FALSE, message = FALSE, cache=TRUE}
## `gac` : estimated vector of discriminating item parameters; as this is not a Rasch-type model, values are not constrained to 1;
gamma_hat <- out_GRM$gac

```

Since two dimensions are supposed, these weights can be divided into two groups

```{r,  warning = FALSE, message = FALSE, cache=TRUE}
dim1 <- multi2[1, ] #items of the first dimension
gamma_hat1 <- gamma_hat[dim1]

dim2 <- multi2[2, ] #items of the second dimension
gamma_hat2 <- gamma_hat[dim2]
```

A useful thing to do is to rescale the parameters __(within each dimension)__, so that the most discriminating item (for each dimension) has gamma_hat:$\hat{\gamma}_j = 1$ and the other items have $\hat{\gamma}_j < 1$. This is obtained by dividing each vector by its maximum element.

```{r,  warning = FALSE, message = FALSE, cache=TRUE}

# These are the equivalent ones of $\hat{\gamma}_j$ reported in Table 3 of the paper (slide 38 of the presentation "NFER_lect5"). 
gamma_hat1 <- gamma_hat1/max(gamma_hat1)
gamma_hat2 <- gamma_hat2/max(gamma_hat2)

```



We can now use __scaled discriminating item parameters parameters__ as weights. 

Since these two sets of parameters are dimension-specific so they are not are not comparable.

Two different composite indicators can be constructed at this stage: one for the first dimension (`dim1`) and one for the second dimension (`dim2`). 

Note that dividing by the weight sum (`sum(gamma_hat1)`) is not necessary, it depends on the type of indicator you need (in terms of sum or average).

```{r,  warning = FALSE, message = FALSE, cache=TRUE}

CI1w <- (data[, dim1] %*% gamma_hat1) / sum(gamma_hat1)
  
CI2w <- (data[, dim2] %*% gamma_hat2) / sum(gamma_hat2)
```


These indicators can be compared with the unweighted counterpart.

```{r,  warning = FALSE, message = FALSE, cache=TRUE}
#mean response for each unit
CI1unw <- rowSums(data[, dim1])/length(dim1) #=rowMeans(data[,dim1])

CI2unw <- rowSums(data[, dim2])/length(dim1) #=rowMeans(data[,dim2])
```

```{r,  warning = FALSE, message = FALSE,  results = "asis", cache=TRUE}

stargazer::stargazer(head(data.frame(CI1w, CI1unw, CI2w, CI2unw)), 
                     title = "Final Weight",
                     summary = FALSE,
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )

```


Visualisation of resulting composite indicator

```{r echo=FALSE, message = FALSE, warning = FALSE, cache=TRUE}

scores <- data.frame(data, CI1w, CI1unw, CI2w, CI2unw)
scores$record <- row.names(scores)

grouped_bar_df <- scores %>%
    gather(key = record, 
           value = Value,
           CI1unw, CI2unw) 
  
#Make plot
grouped_bars <- ggplot(grouped_bar_df, 
                       aes(x = Value, 
                           fill = as.factor(record))) +
  
  #coord_flip() +
  #geom_bar(stat = "count", position = "dodge") +
  geom_line(stat = "count") +
  #geom_freqpoly() +
  #geom_histogram(position="dodge2") +   
  geom_hline(yintercept = 0, size = 1, colour = "#333333") +
  unhcr_style() +
  #scale_fill_viridis_d() +
  labs(title = "Scores using different dimensions Unweigthed",
       subtitle = "4 options", 
       caption = datadescription) +
  scale_y_continuous( label = format_si()) + ## Format axis number
  theme(panel.grid.major.x = element_line(color = "#cbcbcb"), 
        panel.grid.major.y = element_blank()) ### changing grid line that should appear


ggpubr::ggarrange(left_align(grouped_bars, c("subtitle", "title")), ncol = 1, nrow = 1)

grouped_bar_df <- scores %>%
    gather(key = record, 
           value = Value,
           CI1w, CI2w) 
  
#Make plot
grouped_bars <- ggplot(grouped_bar_df, 
                       aes(x = Value, 
                           fill = as.factor(record))) +
  
  #coord_flip() +
  #geom_bar(stat = "count", position = "dodge") +
  geom_line(stat = "count") +
  #geom_freqpoly() +
  #geom_histogram(position="dodge2") +   
  geom_hline(yintercept = 0, size = 1, colour = "#333333") +
  unhcr_style() +
  #scale_fill_viridis_d() +
  labs(title = "Scores using different dimensions Weigthed",
       subtitle = "options", 
       caption = datadescription) +
  scale_y_continuous( label = format_si()) + ## Format axis number
  theme(panel.grid.major.x = element_line(color = "#cbcbcb"), 
        panel.grid.major.y = element_blank()) ### changing grid line that should appear


ggpubr::ggarrange(left_align(grouped_bars, c("subtitle", "title")), ncol = 1, nrow = 1)

```

#### Sinle composite indicator

One wish to have a unique composite indicator, computed by using the entire response pattern of each unit. 

In this case, the discrimination parameters need to be made comparable across dimensions. Then, we have to compute $\gamma_j^*$.

First, we need to compute gamma_hat: $\hat{\sigma}_d$

```{r,  warning = FALSE, message = FALSE, cache=TRUE}
#`Piv`: optional object of type matrix containing the estimated weights of thelatent classes subject-by-subject  (the weights may be different across subjects in the presence of covariates);
#3 prior probabilities
piv <- out_GRM$piv  

# `Th`: estimated matrix of ability levels (support points) for each dimension (= row of matrix Th) and latent class (= column of matrix Th)
#2 x 3 latent traits
theta <- out_GRM$Th 

mu_hat <- c(theta %*% piv) # mean latent traits
sigma_hat <- rep(NA, 2) #for each dimension

for (d in 1:2) {
  # Latent trait std.dev.
  sigma_hat[d] <- sqrt((theta[d, ] - mu_hat[d])^2 %*% piv) 
}

#sigma_hat

```

Then, we can compute gamma_hat: $\gamma_j^* = \hat{\sigma}_d \hat{\gamma}_j, j \in \mathcal{J}_d$.

```{r,  warning = FALSE, message = FALSE, cache=TRUE}

gamma_hat1_star <- sigma_hat[1] * gamma_hat1

gamma_hat2_star <- sigma_hat[2] * gamma_hat2

```

As the standard deviation of the latent traits is taken into consideration, the item discriminations are now comparable.

For example, we could rank the items according to their discrimination.

```{r,  warning = FALSE, message = FALSE, cache=TRUE}
item_comp <- data.frame(item = c(dim1,dim2),
                        discr = c(gamma_hat1_star, gamma_hat2_star))

# Item 7 is the least discriminating one, while item 5 is the most discriminating one 

allitem <- item_comp[order(item_comp$discr), ] 
```

Hence, a unique composite indicator can be built using these last weights.

```{r,  warning = FALSE, message = FALSE, cache=TRUE}

weights_unique <- item_comp[order(item_comp$item), "discr"]

data$CI <- data %*% weights_unique / sum(weights_unique)

#CI

```

Finally, one could rescale these last parameters (in order to obtain those denoted by $w_j$ in the paper).
```{r,  warning = FALSE, message = FALSE, results = "asis", cache=TRUE}
allitem$unique <- weights_unique/max(weights_unique)


stargazer::stargazer(allitem, 
                     title = "Final Weight",
                     summary = FALSE,
                     type = "html", 
                     column.sep.width = "4pt",
                     align = TRUE, 
                     no.space = TRUE )



```



If the data at issue are multilevel (e.g., refugees nested in different neighborhoods/villages/camps districts), the model can be extended to cluster high-level units (i.e. neighborhoods/villages/camps districts) as well. Further, covariates can be included to check the effects (if any) of individual (or higher level) covariates in the units’ class membership probabilities. 

Additionally, Discrete Multidimensional Two-Tier IRT Model, an extension of the previous models can be explored to account for the double dependency of single items on more than one dimension of the Latent Variable.

In further research, the proposed procedure can therefore be extended through a classification of higher-order units in homogeneous groups - instead of ranking them one by one -  that would take advantage of the latent class part of the described method, and account for the multilevel structure of the data. So basically not only the analysis would provide the vulnerability score of household but also account for a potential geographic effect within their vulnerability profile.

Another element of research includes to define the impact of sampling on the resulting weights.