06-rsm.Rmd

# The radiological search model (RSM) {#rsm}


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(exactci)
```


## How much finished 99% {#rsm-how-much-finished}


<!-- Switch to $\mu$, $\lambda$ and $\nu$ and do away with prime parameter? Yes -->



## Introduction {#rsm-intro}

All models of ROC data *that do not incorporate search* involve two fundamental parameters (i.e., not including binning-related threshold parameters). For example, the unequal variance binormal model requires the $a,b$ parameters. Alternative ROC models (e.g., CBM and PROPROC) also require two fundamental parameters.

Two fundamental parameters of ROC models are needed (1) to accommodate the average visibility of lesions in the dataset (e.g., the $a$ or separation parameter) and (ii) the fact that the observed diseased case distribution is usually wider than that of the non-diseased cases (e.g., the $b <1$ parameter). If one assumes same widths for both distributions, so in effect $b=1$ is no longer a free parameter, and one allows a varying number of latent marks on all cases, then it becomes possible that the distribution of the highest rating on diseased cases will have greater width than that on non-diseased cases simply due to the fact that latent NLs on diseased cases will have lower z-samples than latent LLs on diseased cases (i.e., a mix of NL and LLs) while on non-diseased cases there will be only NL z-samples. So the basic idea is to have a visibility parameter, a parameter describing the distribution of the number of latent NLs per case and a parameter describing the distribution of the number of latent LLs per case, i.e., a three-parameter model should suffice. And in fact the RSM contains three fundamental parameters: $\mu$, $\lambda$ and $\nu$. In addition the lowest threshold $\zeta_1$ needs to be included as a parameter as it determines the extent and shape of the RSM predicted operating characteristics. This will become clearer in the next chapter but for now can be illustrated by considering the extreme case $\zeta_1 = \infty$ when the predicted FROC is the single point (0,0). 


## The radiological search model {#rsm-details}

The radiological search model (RSM) for the free-response paradigm is a statistical parameterization of the Nodine-Kundel model. It consists of:

* A *search stage* in which suspicious regions, i.e., the latent marks, are identified via peripheral vision. The total number of latent marks on a case is random non-negative integer and in fact some cases may have zero latent marks, a fact that will turn out to have important consequences for the shapes of all RSM predicted operating characteristics.

* A *decision stage* during which each latent mark is closely examined via foveal scanning, relevant features are extracted and analyzed and the observer calculates a decision variable or z-sample for each latent mark.

* If the z-sample exceeds a pre-selected minimum reporting threshold, denoted $\zeta_1$ the location is marked, i.e., the latent mark becomes an actual mark. 

* Latent marks can be either latent NLs (corresponding to non-diseased regions) or latent LLs (corresponding to lesions). The number of latent NLs or LLs on a case are denoted $l_1, l_2$ respectively. Latent NLs can occur on non-diseased or diseased cases but latent LLs can only occur on diseased cases. Assume that every diseased case has $L$ actual lesions (this will later be extended to arbitrary number of lesions per diseased case). ^[Since the RSM is a parametric model one does not need the four subscript notation needed to account for case and location dependence necessary to describe observed data, as in Chapter \@ref(empirical). This allows for simpler notation, as the reader may have noticed, unencumbered by 4 subscripts as in $z_{k_ttl_ss}$ in Table \@ref(empirical-notation).]


## RSM assumptions {#rsm-assumptions}

**Assumption 1:** The number of latent NLs, $l_1 \geq 0$, is sampled from the Poissson distribution $\text{Pois}()$ with mean $\lambda$:

\begin{equation} 
l_1 \sim \text{Pois}\left ( \lambda \right ) 
(\#eq:rsm-poisson-sampling)
\end{equation}


The probability mass function (pmf) of the Poissson distribution is defined by:

\begin{equation} 
\text{pmf}_{P}\left ( l_1, \lambda \right ) = exp\left ( -\lambda \right ) \frac{{(\lambda)^{l_1}}}{l_1!}
(\#eq:rsm-poisson-pmf)
\end{equation}

**Assumption 2:** The number of latent LLs, $l_2$, where $0 \leq l_2 \leq L$ (since the number of latent LLs cannot exceed the number of lesions) is sampled from the binomial distribution $\text{B}$ with success probability $\nu$ and trial size $L$: 


\begin{equation} 
l_2 \sim \text{B}\left ( L, \nu \right ) 
(\#eq:rsm-binomial-sampling)
\end{equation}

The probability mass function (pmf) of the binomial distribution is defined by:

\begin{equation} 
\text{pmf}_{B}\left ( l_2, L, \nu \right ) = \binom{L}{l_2} \left (\nu  \right )^{l_2} \left (1-\nu  \right )^{L-l_2}
(\#eq:rsm-binomial-pmf)
\end{equation}

> Collectively $\lambda$ and $\nu$ are termed the *search* parameters.

**Assumption 3:** Each latent mark is associated with a z-sample. That for a latent NL is denoted $z_{l_11}$ while that for a latent LL is denoted $z_{l_22}$. Latent NLs can occur on non-diseased and diseased cases while latent LLs can only occur on diseased cases.

**Assumption 4:** For latent NLs the z-samples are obtained by sampling $N \left ( 0, 1 \right )$:

\begin{equation} 
z_{l_11} \sim N \left ( 0, 1 \right )
(\#eq:rsm-sampling-l1)
\end{equation}

**Assumption 5:** For latent LLs the z-samples are obtained by sampling $N \left ( \mu, 1 \right )$:

\begin{equation} 
z_{l_22} \sim N \left ( \mu, 1 \right )
(\#eq:rsm-sampling-l2)
\end{equation}


The probability density function $\phi\left ( z | \mu \right )$ of the normal distribution $N \left ( \mu, 1 \right )$ is defined by:


\begin{equation} 
\phi\left ( z | \mu \right )=\frac{1}{\sqrt{2\pi}}\exp\left ( -\frac{(z-\mu)^2}{2} \right )
(\#eq:rsm-pdf-phi-mu)
\end{equation}

> The parameter $\mu$ is termed the *classification* parameter.


**Bning rule:** In an FROC study with R ratings, the observer adopts $R$ ordered cutoffs $\zeta_r$, where $\left ( r = 1, 2, ..., R \right )$. Defining  $\zeta_0 = -\infty$ and $\zeta_{R+1} = \infty$, then if $\zeta_r \leq z_{l_ss} < \zeta_{r+1}$ the corresponding latent site is marked and rated in bin $r$, and if $z_{l_ss} \leq \zeta_1$ the site is not marked. ($R$ is the number of FROC bins.)

**Mark location:** The location of the mark is assumed to be at the exact center of the latent site that exceeded a cutoff and an infinitely precise proximity criterion is adopted. Consequently, there is no confusing a mark made because of a latent LL z-sample exceeding the cutoff with one made because of a latent NL z-sample exceeding the cutoff. Therefore, any mark made because of a latent NL z-sample that satisfies $\zeta_r \leq z_{l_11} < \zeta_{r+1}$ will be scored as a non-lesion localization (NL) and rated $r$. Likewise, any mark made because of a latent LL z-sample that satisfies $\zeta_r \leq z_{l_22} < \zeta_{r+1}$ will be scored as a lesion-localization (LL) and rated $r$. 

**Rating assigned to unmarked sites:** Unmarked LLs are assigned the zero rating (or any rating lower than the lowest allowed FROC-1 rating). Note that even lesions that were not found by the search stage, and therefore do not qualify as latent LLs, are assigned the zero rating. This is because they represent observable events (and less suspicious than the lowest allowed FROC-1 rating). In contrast, unmarked latent NLs are unobservable events. Unlike lesions there is no a-priori reader-independent list of non-lesion locations; what constitutes a NL is reader dependent, see Fig. \@ref(fig:visual-search-eye-tracking).

By choosing $R$ large enough the preceding discrete rating model is applicable to quasi-continuous z-samples.


## Physical meanings of the RSM parameters {#rsm-parameter-physical-meanings}

The parameters have the following physical meanings:

### The $\mu$ parameter {#rsm-mu-parameter}

The $\mu$ parameter is the lesion *perceptual signal to noise ratio pSNR*, as described in (print book) Chapter 12.5.2, between latent NLs and latent LLs. For white noise background this is similar to the physical SNR [@chakraborty1997computer] after correction for the non-linear response of the visual system to visual stimuli [@siddiqui2005discrete]. For clinical backgrounds pSNR is determined by the competition for the observer's foveal attention from other regions that could be mistaken for lesions.  

The $\mu$ parameter is similar to detectability index $d'$, which is the separation parameter of two unit normal distributions required to achieve the observed probability of correct choice (PC) in a two alternative forced choice task between cued NLs and cued LLs. Individually and for each reader one determines the locations of the latent marks using eye-tracking apparatus and then runs a 2AFC study as follows: pairs of images are shown, each with a cued location, one a latent NL and the other a latent LL, where all locations were recorded in prior eye-tracking sessions for the specific radiologist. The radiologist's task is to pick the image with the latent LL. The probability correct $\text{PC}$ in this task is related to the $\mu$ parameter by:

\begin{equation} 
\mu = \sqrt{2} \Phi^{-1} \left ( \text{PC} \right )
(\#eq:rsm-mu-2afc)
\end{equation}

The radiologist on whom the eye-tracking measurements are performed and the one who performs the two alternative forced choice tasks must be the same, as two radiologists may not agree on latent NL marks. A complication in conducting such a study is that because of memory effects a lesion can only be shown once to each reader: clinical images are distinctive - once a radiologist has found a lesion in a clinical image, that event may become imprinted in long-term memory; one cannot repeatedly compare this lesion to other NLs in the 2AFC task as the radiologist will always pick the remembered lesion. This is a difficult study to conduct as I found out.


### The $\lambda$ parameter {#rsm-summary-lambda-parameter}

The $\lambda$ parameter determines the tendency of the observer to generate latent NLs. The mean number of latent NLs per case is an estimate of $\lambda$. ^[It can be measured via eye-tracking apparatus. This time it is only necessary to cluster the marks and classify each mark as a latent NL or latent LL according to the adopted acceptance radius. An eye-tracking based estimate would be the total number of latent NLs in the dataset divided by the total number of cases.]


I have found it best to illustrate sampling to non-statistics majors with numerical examples. Consider two observers, one with $\lambda = 1$ and the other with $\lambda = 2$. While one cannot predict the exact number of latent NLs on any specific case, the value of $\lambda$ determines the average number of latent NLs. 
 

The following code illustrates Poissson sampling, estimation of the mean and confidence interval for 100 samples from two Poissson distributions. The number of samples has been set to $K_1=100$ (the first argument to `rpois()` is the number of non-diseased cases; the second argument is the value of $\lambda$).


```{r rsm-code-snippet-poisson-samples-ci-100, echo=TRUE}
K1 <- 100
lambda <- c(1,2)
seed <- 1;set.seed(seed);samples1 <- rpois(K1,lambda = lambda[1])
seed <- 1;set.seed(seed);samples2 <- rpois(K1,lambda = lambda[2])

ret11 <- poisson.exact(sum(samples1),K1)
ret21 <- poisson.exact(sum(samples2),K1)
```

```{r, echo=FALSE}
cat ("K1 = ", K1,", lambda 1st reader = ", lambda[1],", lambda 2nd reader = ", lambda[2])
cat("obs. mean, reader 1 = ", mean(samples1), "\n")
cat("obs. mean, reader 2 = ", mean(samples2), "\n")
cat ("Rdr. 1: 95% CI = ", "[", ret11$conf.int[1:2], "]\n")
cat ("Rdr. 2: 95% CI = ", "[", ret21$conf.int[1:2], "]\n")
```


For reader 1 the estimate of the Poissson parameter (the mean parameter of the Poissson distribution is frequently referred to as the Poissson parameter) is 1.01 with 95% confidence interval (0.823, 1.227); for reader 2 the corresponding estimates are 2.02 and (1.751, 2.319). As the number of cases increases, the confidence interval shrinks. For example, with 10000 cases, i.e., 100 times the value in the previous example:

```{r rsm-code-snippet-poisson-samples-ci-10000, echo=FALSE}
K1 <- 10000
lambda <- c(1,2)
seed <- 1;set.seed(seed);samples1 <- rpois(K1,lambda = lambda[1])
seed <- 1;set.seed(seed);samples2 <- rpois(K1,lambda = lambda[2])

ret12 <- poisson.exact(sum(samples1),K1)
ret22 <- poisson.exact(sum(samples2),K1)
```


```{r, echo=FALSE}
cat ("K1 = ", K1,", lambda 1st reader = ", lambda[1],", lambda 2nd reader = ", lambda[2],"\n")
cat("obs. mean, reader 1 = ", mean(samples1), "\n")
cat("obs. mean, reader 2 = ", mean(samples2), "\n")
cat ("Rdr. 1: 95% CI = ", "[", ret12$conf.int[1:2], "]\n")
cat ("Rdr. 2: 95% CI = ", "[", ret22$conf.int[1:2], "]\n")
```

This time for reader 1, the estimate of the Poissson parameter is 1.01 with 95% confidence interval (0.986, 1.025); for reader 2 the corresponding estimate is 2.01 with 95% confidence interval (1.978, 2.034). The width of the confidence interval is inversely proportional to the square root of the number of cases (the example below is for reader 1):


```{r rsm-code-snippet-poisson-samples-ci-100a}
ret11$conf.int[2] - ret11$conf.int[1]
ret12$conf.int[2] - ret12$conf.int[1]
```

Since the number of cases was increased by a factor of 100, the width decreased by a factor of 10, the square-root of the ratio of the numbers of cases.



### The $\nu$ parameter {#rsm-summary-nu-parameter}

The $\nu$ parameter determines the ability of the observer to find lesions. Assuming the same number of lesions per diseased case, the fraction of latent LLs per diseased case is an estimate of $\nu$. ^[It too can be measured via eye-tracking apparatus performed on a radiologist. An eye-tracking based estimate would be the total number of latent LLs in the dataset divided by the total number of lesions.] 


Consider two observers, one with  $\nu = 0.5$ and the other with $\nu = 0.9$. Again, while one cannot predict the number of latent LLs on any specific diseased case, or which lesions will be correctly localized, one can predict the average number of latent LLs per diseased case. 

The following code uses $K_2 = 100$ samples, the number of diseased cases, each with one lesion. The arguments to `rbinom()` - for random binomial samples - are the number of diseased cases, the number of lesions per case and the value of $\nu$.

```{r rsm-code-snippet-binomial-example1, echo=TRUE}
K2 <- 100
nu <- c(0.5, 0.9)
seed <- 1;set.seed(seed);samples1 <- rbinom(K2,1,nu[1])
seed <- 1;set.seed(seed);samples2 <- rbinom(K2,1,nu[2])

ret1 <- binom.exact(sum(samples1),K2)
ret2 <- binom.exact(sum(samples2),K2)
```


```{r, echo=FALSE}
cat ("K2 = ", K2,", nu 1st reader = ", 0.5,", nu 2nd reader = ", 0.9,"\n")
cat("mean, reader 1 = ", mean(samples1), "\n")
cat("mean, reader 2 = ", mean(samples2), "\n")
cat ("Rdr. 1: 95% CI = ", "[", ret1$conf.int[1:2],"]\n")
cat ("Rdr. 2: 95% CI = ", "[", ret2$conf.int[1:2],"]\n")
```

The result shows that for reader 1 the estimate of the binomial success rate parameter is 0.48 with 95% confidence interval (0.379, 0.582). For reader 2 the corresponding estimates are 0.94 and (0.874, 0.978). 


As a more complicated but clinically realistic example, consider a dataset with 100 cases where 97 cases have one lesion per case, two have two lesions per case and one has three lesions per case (these are typical lesion distributions observed in screening mammography). The code follows:


```{r rsm-code-snippet-binomial-example2, echo=TRUE}
K2 <- c(97,2,1);Lk <- c(1,2,3);nu <- c(0.5, 0.9)
samples1 <- array(dim = c(sum(K2),length(K2)))
seed <- 1;set.seed(seed)
# I am using el instead of l as the latter looks like 1
for (el in 1:length(K2)) {
  samples1[1:K2[el],el] <- rbinom(K2[el],Lk[el],nu[1])
}

samples2 <- array(dim = c(sum(K2),length(K2)))
seed <- 1;set.seed(seed)
for (el in 1:length(K2)) {
  samples2[1:K2[el],el] <- rbinom(K2[el],Lk[el],nu[2])
}

ret1 <- binom.exact(sum(samples1[!is.na(samples1)]),sum(K2*Lk))
ret2 <- binom.exact(sum(samples2[!is.na(samples2)]),sum(K2*Lk))
```

```{r, echo=FALSE}
cat("K2[1] =", K2[1],", K2[2] =", K2[2],", K2[3] =", K2[3], ", nu1 =", nu[1], ", nu2 =", nu[2])
cat("obsvd. mean, reader 1 = ", sum(samples1[!is.na(samples1)])/sum(K2*Lk))
cat("obsvd. mean, reader 2 = ", sum(samples2[!is.na(samples2)])/sum(K2*Lk))
cat ("Rdr. 1: 95% CI = ", ret1$conf.int[1:2])
cat ("Rdr. 2: 95% CI = ", ret2$conf.int[1:2])
```


For reader 1, the estimate of the binomial success probability is 0.490 with 95% confidence interval (0.391, 0.590); for reader 2 the corresponding estimates are 0.933 and (0.866, 0.973). 


## Intrinsic RSM parameters {#rsm-intrinsic-parameters}

While the parameters $\lambda$ and $\nu$ are physically meaningful a little thought reveals that they must depend on $\mu$. From the solar-analogy described in Section \@ref(froc-paradigm-solar-analogy) we know that if $\mu = 0$ the lesions have zero contrast and therefore cannot be found by the search mechanism implying $\nu = 0$. Moreover attempting to find these zero contrast lesions must generate a large number of non-lesion localizations implying $\lambda = \infty$. 


The following is a simple model of the $\mu$ dependence of $\lambda$ and $\nu$. The model re-parameterizes the *physical* parameters $\lambda$ and $\nu$ in terms of *intrinsic* parameters $\lambda_i$ and $\nu_i$ that are $\mu$ independent^[The need for the first re-parameterization, involving $\nu$, was foreseen in the original search model papers [@chakraborty2006search; @chakraborty2006roc] but the need for the second re-parameterization (involving $\lambda$) became evident more recently.]:


\begin{equation}
\left. 
\begin{aligned}
\nu =& 1 - \text{exp}\left ( - \mu \nu_i \right ) \\
\lambda =& \frac{\lambda_i}{\mu}
\end{aligned}
\right \}
(\#eq:rsm-transform)
\end{equation}


The inverse transformations are:


\begin{equation}
\left. 
\begin{aligned}
\nu_i =& - \frac{\ln \left ( 1-\nu \right )}{\mu}\\
\lambda_i =& \mu \lambda 
\end{aligned}
\right \}
(\#eq:rsm-inv-transform)
\end{equation}


The intrinsic parameters obey $\lambda_i \ge 0$ and $\nu_i \ge 0$.  

Since it determines $\nu$, the $\nu_i$ parameter can be considered as the intrinsic (i.e., $\mu$-independent) ability to find lesions; specifically, *it is the rate of increase of $\nu$ with $\mu$ at small $\mu$*:


\begin{equation} 
\nu_i = \left (\frac{\partial \nu}{\partial \mu}  \right )_{\mu = 0}
(\#eq:rsm-nup-limit)
\end{equation}


According to Eqn. \@ref(eq:rsm-transform), as $\mu \rightarrow \infty$, $\nu \rightarrow 1$ and conversely, as $\mu \rightarrow 0$, $\nu \rightarrow 0$. The solar analogy in Section \@ref(froc-paradigm-solar-analogy) is instructive. The dependence of $\nu$ on $\mu$ is consistent with the fact that higher contrast lesions are easier to find. A non-expert is expected to find a high contrast lesion whereas a low contrast lesion will be more difficult to find even by an expert observer. 

According to Eqn. \@ref(eq:rsm-transform) the value of $\mu$ also determines $\lambda$: as $\mu \rightarrow \infty$, $\lambda \rightarrow 0$, and conversely, as $\mu \rightarrow 0$, $\lambda \rightarrow \infty$. Here too the solar analogy in Section \@ref(froc-paradigm-solar-analogy) is instructive. Since the sun has very high contrast, there is no reason for the observer to search for other suspicious regions which have no possibility of resembling it. On the other hand, attempting to locate a faint star can generate many false sightings because the expected small contrast from the faint real star could be comparable to that from a number of regions in the nearby background.


## Summary {#rsm-discussion-summary}

This chapter has described a statistical parameterization of the Nodine-Kundel model of visual search. The model accounts for key aspects of the process: 

* Search: finding lesions and finding non-lesions. These are characterized by the two search parameters $\lambda$ and $\nu$. 
* Classification: The ability to correctly rate a lesion higher than a NL is characterized by the third (classification) parameter of the model $\mu$. 

While the 2 search parameters have relatively simple physical meanings they depend on $\mu$. Consequently, it is necessary to introduce intrinsic parameters $\lambda_i$ and $\nu_i$ which are independent of $\mu$. 


The next chapter explores the ROC curve predictions of the radiological search model.