EligibilityScreening_LS.Rmd

---
title: "Eligibility Screening: 5-pt-Likert-Scale"
author: "Jutta Pieper"
date: "28.04.2022"
output: 
  html_document:
    toc: true
    number_sections: true
    toc_float: 
      collapsed: true
    toc_depth: 2
    keep_md: yes
bibliography: references/bibliography.bibtex
csl: references/apa.csl
link-citations: true
---

```{r setup, echo = FALSE, include=TRUE, warning = FALSE, message = FALSE}
library(tidyverse)
library(DT)
options(DT.options = list(dom = 'Blfrtip', pageLength = 5, searching = FALSE, buttons = c( 'csv', 'excel', 'pdf'))) 
library(ggpubr)
library(knitr)
library(kableExtra)
options(knitr.kable.NA = "")
kable <- function(data, ...) {
   knitr::kable(data, digits=3, ...) %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
}
knit_print.data.frame <- function(x, ...) {
  res <- paste(c("", "", kable(x)), collapse = "\n")
  asis_output(res)
}
registerS3method("knit_print", "data.frame", knit_print.data.frame)
```

Read in questionnaire:

```{r read_merge_data}
filler_items = read.csv("./judgments/all/LikertSkala_filler.csv", fileEncoding = "UTF-8")
test_items = read.csv("./judgments/all/LikertSkala_test.csv", fileEncoding = "UTF-8") %>% 
  mutate(ITEM_FUNCTION = "test")
questionnaire = bind_rows(filler_items, test_items) %>%
  mutate_at(.vars = c("CONDITION_NO","KEY_CONDITION", "workerId"), .funs = as_factor)
```


```{r init_worker_profile, echo = FALSE}
## Store all gathered information in data.frame for worker overview and final decision 
worker_profile = data.frame(
  workerId = character(), 
  criterion = character(), 
  score = numeric(), 
  accept = logical())
```

# Progress (incomplete submissions) 

We require that at least 90\% of the questionnaire has been completed. 
(Note: we assume that there are no missing trials due to data loss / technical errors).

```{r available_trials}
available_trials = questionnaire %>%
  group_by(workerId) %>%
  summarise(trials = n(), 
            trials_prop = n()/length(unique(questionnaire$trial_index))) %>%
  mutate(accept = ifelse(trials_prop < 0.9, FALSE, TRUE)) %>%
  arrange(trials_prop)
```

Let's check whether there are participants that left the questionnaire prematurely:

```{r}
available_trials %>% filter(trials_prop < 1) %>% kable() 
```

```{r wp_available, echo = FALSE}
worker_profile = rbind(worker_profile, 
                       available_trials %>% 
                         select(workerId, trials_prop, accept) %>% 
                         rename(score = trials_prop) %>% 
                         mutate(criterion = "progress"))
```


Remove incomplete data from questionnaire:

```{r}
 questionnaire = questionnaire %>%
   filter(!(workerId %in% (
     available_trials %>% filter(accept == FALSE) %>% pull(workerId)
   )))
```


# Latency-based identification 


## Spammers

We will reject workers with *mean*(RT) < 3000 ms (as simple spammers) and those  with *median*(RT) < 3000 ms (as clever spammers):


```{r}
simple_spammer = questionnaire %>%
                         group_by(workerId) %>%
                         summarise(score = mean(rt)) %>%
                         mutate(criterion = "meanRT", 
                                accept = ifelse(score < 3000, FALSE, TRUE))
```

Are there any simple spammers?
```{r}
simple_spammer %>% filter(accept == FALSE) %>% kable() 
```

```{r}
clever_spammer = questionnaire %>%
                         group_by(workerId) %>%
                         summarise(score = median(rt)) %>%
                         mutate(criterion = "medianRT", 
                                accept = ifelse(score < 3000, FALSE, TRUE))
```

Are there any clever spammers?
```{r}
clever_spammer %>% filter(accept == FALSE) %>% kable() 
```

```{r echo = FALSE}
worker_profile = rbind(worker_profile, simple_spammer)
worker_profile = rbind(worker_profile, clever_spammer)
```


## RT distributions

Participants and different (groups of) items presumably exhibit different RT distributions:

```{r rt_density_fc}
rt_density_worker = questionnaire %>%
  ggplot(aes(x=log(rt), group = workerId)) + 
  geom_density(alpha=.3, fill = "#DFF1FF")

rt_density_item_funs = questionnaire %>% 
  ggplot(aes(x = log(rt), fill=ITEM_FUNCTION)) +
  geom_density(alpha=.3) + 
  scale_fill_brewer(palette = "Greys")
```
  
```{r rt_density_fc_plot, echo = FALSE, fig.width = 10, fig.height=5, fig.cap="Figure 1 in @Pieper_et_al_2022",fig.topcaption=TRUE}
ggpubr::annotate_figure(
  ggpubr::ggarrange(
    rt_density_item_funs + 
      theme(legend.position = "None") + 
      ggtitle("per item function") +
      theme_bw(base_size = 16), 
    rt_density_worker  + 
      theme(legend.position = "bottom", legend.title = element_blank()) + 
      ggtitle("per participant") +
      theme_bw(base_size = 16),
    nrow = 1, common.legend = TRUE, legend = "bottom")) +
  theme_minimal(base_size = 20)
```

The more diverse the distribution the stronger is the superiority of ReMFOD (see next chapter) above generic cutoff points for outlier detection. 

## Underperforming 

**Recursive multi-factorial outlier detection (ReMFOD , @Pieper_et_al_2022)**
```{r}
source("./R Sources/ReMFOD.R")
```

ReMFOD (see [source code](./R Sources/ReMFOD.R)) aims at identifying individual trials as genuine intermissions  and rushes. In doing so, ReMFOD accounts for different RT distributions of different participants and item functions, as well as swamping and masking effects. Underpinned by these suspicious individual trials, underperforming participants can be determined by means of proportion of trials not responded to wholeheartedly. We propose to discard participants who have responded genuinely to less than 90~\% of trials because they supposedly did not meet the task with the necessary seriousness.

To account for different RT distributions, ReMFOD compares the RT of each trial to a lower and an upper cutoff point, which each consider two  cutoff criteria, respectively: The first criterion  is computed with respect to the group of trials with the same *item function* (i.e. attention trials only, control trials only, etc.) regardless of the participant responding, the second one is computed with respect to  all trials of the corresponding participant (regardless of the item function). Only if an RT surmounts or  falls below *both*  criteria, it will be designated as a *genuine intermission*
or as a *genuine rush*.[^1]

[^1]: @Miller_1991 proposes the values of 3 (very conservative), 2.5 (moderately conservative) or even 2 (poorly
conservative).  @Haeussler_Juzek_2016 suggest using an asymmetric criterion (using standard deviations) of -1.5 for the lower and +4 for the upper cutoff point.

\begin{align}    \label{eq:cutoff_outlier_rt}
    \textit{cutoff_}&\textit{intermission} = \max \left\{  \right.\\
&\left. \text{median}(\textit{RTs:participant}) + 2.5 \times \text{mad}(\textit{RTs:participant}), \right.\nonumber\\
 &\left. \text{median}(\textit{RTs:item_function}) + 2.5 \times \text{mad}(\textit{RTs:item_function})\right.\nonumber\}
\end{align}

\begin{align}  \label{eq:cutoff_guesses_rt}
    \textit{cutoff_}&\textit{rush} = \min \left\{  \right.\\
&\left. \text{median}(\textit{RTs:participant}) - 1.5 \times \text{mad}(\textit{RTs:participant}), \right.\nonumber\\
 &\left. \text{median}(\textit{RTs:item_function}) - 1.5 \times \text{mad}(\textit{RTs:item_function})\right.\nonumber\}
\end{align}

To account for swamping and masking effects (see @Ben-Gal_2005), 
the process described above will be repeated on a reduced data set (i.e. excluding already detected outliers) until no more outliers  can be found. Therefore, in each iteration step, the cutoff points must be computed afresh. 



### Overview plot for item functions

Different outlier types, computed with respect to different groups, are marked by different shapes: Box-shaped trials are the only RTs we consider as genuine intermissions or rushes. Note that the shapes may overlap as these outliers have been computed by various procedures differing in the groups they included to identify outliers.

```{r ls_remfod_plot_item_function,  fig.width = 10, fig.height=5}
## compute different outlier types based on the whole questionnaire and plot these
remfod_plot = questionnaire  %>% outlier_plots_remfod()
item_plot = remfod_plot # we are going to reuse remfod_plot for workers
## remove data we are currently not interested in from the plot 
item_plot$data = item_plot$data %>% 
  filter(!ITEM_FUNCTION %in% c("calibration", "filler"))
## structure plot as you like
item_plot + facet_wrap(~ ITEM_FUNCTION, nrow = 1)
```


### Performance of participants

We expect that 90 % of the trials are answered without genuine intermission or rushes, i.e. that 90 % of the RTs are *valid*

```{r genuine_outliers_summ, message = FALSE, warning = FALSE}
rt_outlier_count = remfod(questionnaire)  %>% 
  group_by(workerId,  direction) %>% tally() %>%
  spread(key = "direction", value = "n") %>% 
  rename( none = "<NA>") %>%  
  mutate_if(is.numeric, replace_na, 0) %>%
  mutate(trials_total = sum(long, short, none),
         prop_long = long/trials_total, 
         prop_short = short/trials_total,
         prop_out = (short+long)/trials_total,
         prop_valid = none/trials_total,
         accept = ifelse(prop_valid < 0.9, FALSE, TRUE)) %>%
  arrange(prop_valid) %>%
  mutate_if(is.numeric, round, 4)
```

```{r genuine_outliers_summ_table,  echo = FALSE, message = FALSE, warning = FALSE, echo = FALSE}
rt_outlier_count %>%  datatable(rownames = FALSE, extensions = 'Buttons')
```


```{r worker_profile_remfod, echo = FALSE}
worker_profile = rbind(worker_profile, 
  rt_outlier_count %>% 
    select(workerId, prop_valid, accept) %>% 
    rename(score = prop_valid) %>%
    mutate(criterion = "validRTs")) 
```

#### Plots of individual participants

Some underperforming participants

```{r FC_underperforming, echo = FALSE}
worker_plot = remfod_plot ## computed different outlier types based on the whole questionnaire and plot these
## remove data we are currently not interested in from the plot 
worker_plot$data = worker_plot$data %>% 
  filter(workerId %in% (rt_outlier_count %>% filter(accept == FALSE) %>% pull(workerId)))
## structure plot as you like
no_pages = min(9,floor((rt_outlier_count %>% filter(accept == FALSE) %>% nrow()) / 3)) # prevent bug facet_wrap_paginate ~ show only full pages
for(i in c(1:no_pages)){
  plot(worker_plot + ggforce::facet_wrap_paginate(~ workerId, nrow = 1, ncol = 3, page = i))
}
```


```{r FC_exemplative, eval = all(c(161,205,174) %in% unique(remfod_plot$data$workerId)), echo = FALSE, results='asis'}
cat("***\nFurther exemplative workers\n\n")

worker_plot = remfod_plot ## computed different outlier types based on the whole questionnaire and plot these
## remove data we are currently not interested in from the plot
worker_plot$data = worker_plot$data %>%
  filter(workerId %in% c(161,205,174))
## structure plot as you like
worker_plot + facet_wrap(~ workerId, nrow = 1, ncol = 3)
```


# Item-based identification 

```{r attention_trials}
control_trials = questionnaire %>% 
  filter(ITEM_FUNCTION == "control")  %>%
  droplevels()

attention_trials = questionnaire %>% 
  filter(ITEM_FUNCTION == "attention")  %>%
  droplevels()
```

Each control group should provide chances of less than 5% to pass controls by guessing. 

## Guessing probabilities  
```{r}
source("./R Sources/GuessingProbs.R")
```

To compute the probability (by standard binomial expansions) of (*exactly!*) k correct answers out of N trials, where 

- ***k*** amount of trials answered correctly
- ***N*** mount of total trials
- ***p*** the probability of a correct response
- ***q*** the probability of an incorrect response

we can use the formula  (see @Frederick_Speed_2007), implemented by the  function `k_out_of_N`:

\begin{align} 
 \label{eq:probs_binomial}
 \frac{N!}{k!(N - k)!}&p^{k}q^{N-k}\text{, where}   \\
p &= \text{ probability of a correct response}  \nonumber  \\
q &= \text{ probability of an incorrect response} \nonumber 
\end{align}

To compute the probability (by standard binomial expansions) of k *or fewer* correct answers out of N trials, we can use the function `k_out_of_N_cumulative`, which returns the sum of all the probabilities from 0 to k  (see [source code](./R Sources/GuessingProbs.R)).


## Criterion selection

Attention trials only exist in ungrammatical conditions, the analysis hence amounts to counting acceptable response.


Let's look at the chances to answer **at least** k out of N items correct if three of the options are considered correct, and two incorrect in a 5-pt-LS are considered correct:

```{r eval = FALSE}
k_out_of_N_matrix_cumulative(p = 3/5, Ns = seq(4,16,2), ks = c(2:12))
```


```{r echo = FALSE}
k_out_of_N_matrix_cumulative(p = 3/5, Ns = seq(4,16,2), ks = c(2:12)) %>% 
  rename("N\\k" = N) %>%
  mutate_if(is.numeric, round, digits = 3) %>% kable(caption = "Table 3a in @Pieper_et_al_2022") %>%
  column_spec(1,bold=T) 
```

If we assume a probability of less than 5~\% to pass a test by chance and combine it with the qualification that not all trials need to be answered correctly, it follows that nine out of ten trials need to be responded to correctly. 

If we take the neutral point not to be acceptable, this number is reduced to five correct responses out of six trials:

```{r eval = FALSE}
k_out_of_N_matrix_cumulative(p = 2/5, Ns = seq(4,16,2), ks = c(2:12))
```


```{r echo = FALSE}
k_out_of_N_matrix_cumulative(p = 2/5, Ns = seq(4,16,2), ks = c(2:12)) %>% 
  rename("N\\k" = N) %>%
  mutate_if(is.numeric, round, digits = 3) %>%
  kable(caption = "Table 3b in @Pieper_et_al_2022") %>%
  column_spec(1,bold=T) 
```


As we stick to the positional account  for control trials (see @Pieper_et_al_2022), we need to evaluate 
grammatical and ungrammatical trials separately. 

On a positional account, the assessment of participants' performance is carried out separately for grammatical (`CONDTION_NO == 1`) and ungrammatical (`CONDTION_NO == 2`) stimuli,  focusing on the pertinent  side of the scale in each case. We recommend allowing the neutral point as legitimate response to grammatical stimuli but not to ungrammatical stimuli as it does neither reliably indicate the rejection of the grammatical version nor the rejection of the ungrammatical version. The chance to pass control trials is then the joint probability of passing the two groups individually. Regarding our rule of thumb, these shall not exceed 5\%.

Let's have a look at some joint probabilities for *up to eight* trials per group, where group sizes are equal. 

```{r}
probs_joint = do.call(rbind, lapply(1:8,function(i){
  probs_joint_positional(i, 3/5, # N, p grammatical
                         i, 2/5) #N, p ungrammatical
}))
```

To receive more balanced options, we further set the constraints that passing grammatical and ungrammatical trials needs to be below 0.6 and that $k$ needs to be less than $N$ in both conditions.


```{r}
probs_joint %>%
  filter(joint_probs <= 0.05 & gram_probs < 0.6 & ungram_probs < 0.6) %>%
  filter(gram_k < gram_N & ungram_k < ungram_N) %>%
  mutate_if(is.numeric, round, digits = 4) %>%
  unite(N, gram_N, ungram_N, sep = "-") %>% relocate(N) %>%
  arrange(N, ungram_k) %>%
  kable(caption = "Table 4 in @Pieper_et_al_2022")
```


### Evaluation groups

There are different ways to proceed: we could evaluate attention and control items as one group, evaluate them separately, and even evaluate related and unrelated control items separately. 
The best choice may depend on the exact nature of your trials.
We are going to compute evaluations based on all groups mentioned.

To facilitate computation, we combine the different groups (named by `ITEM_FUNCTION`) into a single frame -- whereby adapting `ITEM_FUNCTION` in some cases: 

```{r}
eval_trials = rbind(control_trials, attention_trials) %>% 
  # attenion or control items evaluated as one group
  mutate(ITEM_FUNCTION = "attention|control") %>% 
  bind_rows(attention_trials) %>%
  bind_rows(control_trials) %>% 
  # separate evaluation of related and unrelated control trials
  bind_rows(control_trials %>% 
                     mutate(ITEM_FUNCTION = paste(ITEM_FUNCTION, ITEM_SUBGROUP, sep="_"))
                   ) 
```

Find number of trials in each group (per questionnaire):
```{r message = FALSE}
Ns =   eval_trials %>% 
  select(ITEM_FUNCTION, CONDITION_NO, itemId) %>%
  distinct() %>% 
  group_by(ITEM_FUNCTION, CONDITION_NO) %>% 
  tally() %>%
  spread(key = CONDITION_NO, value = n, sep = "_") %>%
  mutate_if(is.numeric, replace_na, 0)
```

```{r echo = FALSE}
Ns %>% kable()
```

compute joint probabilities (of passing grammatical and ungrammatical trials):

```{r}
probs_joint = mapply(probs_joint_positional,
       Ns$CONDITION_NO_1, 3/5, Ns$CONDITION_NO_2, 2/5) %>%
  t() %>% 
  as.data.frame() %>%
  mutate(ITEM_FUNCTION = Ns$ITEM_FUNCTION) %>% 
  relocate(ITEM_FUNCTION) %>%
  unnest(cols = names(.))
```

Note: we may now have inbalanced group sizes. 
Find required k to each group:

```{r}
eval_criteria = probs_joint %>% 
  filter(joint_probs <= 0.05) %>%
                            ## do not remove attention           
  filter((gram_probs < 0.5||gram_probs == 1) & ungram_probs < 0.5) %>%
                            ## do not remove attention
  filter((gram_k < gram_N || gram_N == 0) & ungram_k < ungram_N) %>% 
  mutate_if(is.numeric, round, digits = 3)  %>%
  group_by(ITEM_FUNCTION) %>%
  slice_max(joint_probs)
```

suggested thresholds (i.e. minimal requirements)

```{r echo = FALSE, fig.align='center'}
eval_criteria %>% arrange(ITEM_FUNCTION) %>% kable()
```
put in long format:

```{r}
eval_criteria_long = rbind(
  eval_criteria %>% 
    select(ITEM_FUNCTION, starts_with("gram")) %>%
    rename_all(~stringr::str_replace(.,"^gram_","")) %>%
    mutate(CONDITION_NO = 1)
  ,
  eval_criteria %>% 
    select(ITEM_FUNCTION, starts_with("ungram")) %>%
    rename_all(~stringr::str_replace(.,"^ungram_","")) %>%
    mutate(CONDITION_NO = 2)
) %>% 
  rowwise() %>%
  mutate(prop_k = k/N) %>%
  relocate(CONDITION_NO, .after = ITEM_FUNCTION) %>%
  arrange(ITEM_FUNCTION)
```


```{r echo = FALSE, fig.align='center'}
eval_criteria_long %>% arrange(ITEM_FUNCTION) %>% kable()
```

count correct responses and compare to number of required correct responses


```{r eval_item_based, warning = FALSE, message= FALSE}
item_based_eval = 
  eval_trials %>% 
  mutate(ANSWER_correct = ifelse(CONDITION_NO == 1, # grammatical
                                  ifelse(ANSWER < 3, "incorrect", "correct"), 
                                 #ungrammatical 
                                  ifelse(ANSWER < 3, "correct", "incorrect"))) %>%
  group_by(workerId, ITEM_FUNCTION, CONDITION_NO, ANSWER_correct) %>% 
  tally() %>% ungroup(workerId) %>%
  spread(key = ANSWER_correct, value = n) %>% 
  mutate_if(is.numeric, replace_na,0) %>% 
  ## we use proportions in order to deal with potentially missing data
  mutate(prop_correct = correct / (correct+incorrect)) %>%
  merge(eval_criteria_long) %>% 
  mutate(accept = ifelse(prop_correct < prop_k, FALSE, TRUE)) %>%
  arrange(prop_correct)
```


```{r echo = FALSE}
item_based_eval  %>% 
  mutate_if(is.numeric, round, digits = 3) %>%
  datatable(rownames = FALSE, extensions = 'Buttons')
```


```{r worker_profile_item_based, echo = FALSE}
worker_profile = rbind(worker_profile, 
  item_based_eval %>% 
    select(workerId, prop_correct, accept, ITEM_FUNCTION, CONDITION_NO) %>% 
    unite(criterion, ITEM_FUNCTION, CONDITION_NO, sep = "_") %>%
    rename(score = prop_correct)
)
```


#  Rejected participants 

Let's have a look at the criteria we used to evaluate participants: 
```{r}
unique(worker_profile$criterion)
```
If we applied all of those criteria, and reject all participants who failed on any of them, how many participants would we be left with, i.e. accept?

```{r}
worker_acceptance = worker_profile %>%
  select(-score) %>%
  group_by(workerId) %>%
  summarise(accept = !any(!accept))

table(worker_acceptance$accept)
```

## Rejection reasons 

For convenience, we have (covertly) stored all information gathered in a table named `worker_profile`. 
Let's have a look at how many participants are rejected the individual criteria: 

```{r}
## individual
rejection_reasons = worker_profile %>% 
  group_by(criterion, accept) %>% 
  tally() %>%
  spread(key = accept, value = n) %>%
  rename(acccept = 'TRUE', reject = 'FALSE')
```
 
```{r echo=FALSE, fig.align = 'center'}
kable(rejection_reasons)
```
As expected, many participants fail frequently on attention trials. As opposed to Forced Choice tasks, these can merely measure attention, but cannot draw participants attention to the area of the manipulation (and such improving performance in the future). 


Let's have a look at how these reasons combine, but, for the sake of simplicity, under restrictions to groups as proposed in @Pieper_et_al_2022, i.e. with separate evaulation of attention and control trials, but joined evaluation of related and unrelated controls. 
```{r}
## Combined Rejection Reasons 
rejection_reasons_combined = worker_profile %>%
  filter(accept == FALSE) %>%
  filter(!grepl("\\|",criterion) & !(grepl("related",criterion))) %>%
  group_by(workerId) %>%
  arrange(criterion) %>% 
  summarise(criteria = paste(criterion, collapse = " + ")) %>%
  group_by(criteria) %>% tally() %>%
  arrange(criteria)
```

```{r echo=FALSE, fig.align = 'center'}
kable(rejection_reasons_combined)
```




## Final decision
If we do not want to apply certain criteria, we can delete them now from the worker profile: 

As we find that enough participants passed attention trials (although it was required that all trials are responded to correctly), we do remove the group `attention|control` as it is hence not needed (and for illustration purposes):

```{r}
worker_profile = worker_profile %>%
  ## restrict to attention and control as groups
  filter(!grepl("\\|",criterion) & !(grepl("related",criterion)))
```


## Participants overview table

we might also check on individual participants: 

```{r}
worker_profile %>%
  mutate(score = ifelse(score > 1, 
                        score/1000, # seconds 
                        score*100) # percent
         ) %>%
  mutate(score = format(round(score, 2), nsmall = 2)) %>%
  select(-accept) %>%
  spread(key = criterion, value = score) %>%
  merge(worker_acceptance) %>%  relocate(accept) %>%
  mutate(accept = ifelse(accept,"yes","no")) %>%
  datatable(rownames = FALSE, extensions = 'Buttons',  options = list(
            columnDefs = list(list(className = 'dt-right',
                                   targets = 1:(1+length(unique(worker_profile$criterion)))))
            ))
```


## Remove ineligible participants

Rejected Workers:
```{r}
reject = worker_profile %>%
  filter(accept == FALSE) %>%
  pull(workerId) %>% unique() %>% sort()

length(reject)
reject
```

Remove their data from fillers:
```{r}
filler_items = filler_items %>%
  filter(! workerId %in% reject)
write.csv(filler_items, "./judgments/eligible/LikertSkala_filler.csv", fileEncoding = "UTF-8", row.names = FALSE)
```


Remove their data from test items:
```{r}
test_items = test_items %>%
  filter(! workerId %in% reject)
write.csv(test_items, "./judgments/eligible/LikertSkala_test.csv", fileEncoding = "UTF-8", row.names = FALSE)
```

# References