In rickeycarter/mrmctools: Statistical Tools to Support Multi-Reader, Multi-Case (MRMC) Study Designs

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## Load some packages to tidy up the markdown file. This isn't related to the package.
library(knitr)
library(kableExtra)

Part 1: Data Exploration and Introduction

The package \code{mrmctools} includes a suite of tools to aid in the processing and presentation of Multi-Reader, Multi-Case (MRMC) studies. While the main MRMC analysis is performed in the RJafroc package, there is often significant preprocessing and exploratory analyses that need to be accomplished prior to the finalization of the analysis. This is where this package comes in. This vignette is a viewed as a working draft of an analysis start to finish.

Brief MRMC Study Overview

Briefly, a MRMR is a study design frequently used to evaluate a new imaging strategy's performance at some defined radiologic task. For example, one might be interested in determining how well radiologists read lower dose CT scans while searching for liver cancer (hepatic metastases). A typical study design would be to acquire (or artificially generate) a series of CT scans for a series of patients. Each scan within the series would be at a different dose. The study design is normally fully crossed in that each patient has each dose available. The series of images are randomized and radiologists will read all of the images in a counter-balanced manner (e.g., mixing of doses in a reading day such that only one of hte patient's exams is reviewed within a session).

Reading an image is written generically for accessibility. What occurs is the readers will review each scan and localize (draw a region of interest, ROI) around suspicious features in the dataset. In our studies, readers will assign a confidence score for a detection (confidence some anonomly is present) and the confidence score that they believe the suspicious region is related to the primary task. This latter confidence score is known as the primary task confidence. The primary task concept is an important aspect to these studies. The primary task, which is detection of hepatic metastases in the example that will be illustrated. Decision points around what happens if a reader sees another defect not related to the primary task in the case and notes it need to be addressed. For example, what if liver cirrhosis is identified? Should this ROI be counted in the analysis? Despite one's natural tendancy to say "no, this shouldn't be counted", this is a difficult question to answer. Part of this package provides a standardized framework to account for and document these decisions in the analysis.

The system we utilize to capture the study data utilizes two files. These are described below. To facilitate exploration of this package, we include data from a recently completed MRMC study to demonstrate the packages functions.

Data preparation steps

Our team currently utilizes an in house built reader workstation to record the reader marks. The first step in the MRMC analysis is to prepare the workstation files for analysis. Two excel or CSV files are needed for the analysis. The first is the anonymization file.

1. Anonymization file

This file should contain columns headings for the subject ID number followed by columns of each imaging strategy. This file is often manually edited so that the imaging strategies have clear labels that confirm for standard variable names. The program allows for formatting of the variable names for plots, so essentially, short unique variables names that are informative should be used. Unnecessary columns and rows should be removed removed. An example of a properly formatted anonymization file is provided with the package and examined below.

1. Workstation output file

This is the file that contains all of the reader marks matched, when possible, to the reference marks. Generally, no interaction with this file is required as the generic analysis functions first normalize this file’s column headings programmatically. If should be noted that the format of the workstation file is largely dependent on our workstation file. The file represents a full join of true reference lesions with reader marks after co-localization. While our workstation natively produces this file, this input datafile could be prepared by other means.

The package includes a dataset pairing. The anonymization file is demoanon and the workstation file is demodata. The following code sets these package provided datafiles to a local object for the processing.

library(mrmctools)
data(demodata)
rawdata_t1 <- demodata
data(demoanon)
rawdeidentify <- demoanon

# Levels and Labels to be used for formatting the doses / modalities
# The order should generally be routine/standard down to lowest dose or experimental setting
# These same labels will be used later for additional plotting purposes
# change the order of the dataset. This study starts with routine dose first. I defined these once and 
# used the name and label below
flevels=c("Routine", "Low1", "Low2", "Low3", "Low4")
flabels=c("Routine Dose", "Lower Dose 1", "Lower Dose 2", "Lower Dose 3", "Lower Dose 4")
refdosename <- flevels[1]
refdoselabel <- flabels[1]

This is a listing of the first few rows in each data file.

Reader markings

kable(head(rawdata_t1))

As you can see in the table above, some of the data is customized to an individual survey. There is also a numeric values represented as text columns. The next section will demonstrate the preprocessing steps.

Anonymization File

kable(head(rawdeidentify))

For the anonymization file, you will see a unique subject identifier. The programs that follows will provide a further de-identification of the data.

Preprocessing and normalization

In the preprocessing process, two main functions will be used: readerworkstationfile and applyflowchart.

Function readworkstationfile()

Once the raw anonymization file and workstation file are available, the first package function to be run is `readworkstationfile. The function has parameters that map the workstation variable names to the standardized names. There are also two additional parameters:anonymizereaderandexpandcases. Both of these should generally beTRUE. The reader’s usernames are recorded by the workstation by default.anonymizereaderreplaces these with a randomly assigned integer.expandcases`` provides some error control in case there are some missed combinations of readers and lesions (i.e., a lesion that was never identified by any of the readers on any of the imaging strategies).

raw1 <- readworkstationfile(workstationdata=rawdata_t1, 
                  deidentificationdata=rawdeidentify,
                  refname="Routine",  
                  caseindex="orig_number", 
                  RefID = "REF_SurveyResponseID", 
                  RefDataID = "REF_CaseName",
                  RefDetConf = "REF_DetectionConfidence",
                  RefCode = "REF_DiagnosisCode",
                  RefPTC = "REF_PTC",
                  datasetname = "OBS_CaseName",
                  ObserverPTC = "OBS_PTC",
                  ObserverCode = "OBS_DiagnosisCode",
                  ObserverROIID = "OBS_SurveyResponseID",
                  ReaderFullID = "OBS_UserName",
                  anonymizereader = T, expandcases = T)

This is the first few rows of the normalized data file

kable(head(raw1))

Function applyflowchart()

The applyflowchart function applies the logic in the following figure. This algorithm bins each of the reader ROIs with the primary task confidence (PTC) threshold >= delta into one of 6 bins (Denoted A-F) to facilitate the free response analysis. Bin G (“ignore”, top) does not apply if delta=0. Selecting the PTC threshold >> 0 will generally affect estimated figures of merit (FOM).

The code below is an example function call. The primary task (PT) code in the workstation is 91 for this study. The vector of codes can be easily expanded as the study dictates. However, for the purpose of JAFROC analyses, a mark is either supportive for the primary task ("lesion localization (LL)") or not ("non-lesion localization (NLL)"). This is primarly based on the co-localization and the primary task confidence assigned. It should be noted that as is indicated in the flowchart, the reader-assigned classification (e.g., code = 91) does not affect the calculation of the FOM if all bins are included in the analysis. Only by excluding reader ROIs (e.g., excluding bin A) would the FOM be affected. The NLL and LL indicate whether the lesion counts against or for the FOM, respectively. Informally, a NLL is a false positive and a LL is a true positive.

The function call below also uses a PTCthreshold of 0. This retain all reader ROIs in the dataset. Setting the value at 1 would cause all reader ROIs with a confidence score of 0 for the primary task to be excluded from all analyses. Undermost circumstances, this confidence limit can be 0 or 1 and the analysis will yield essentially the same result. The reason for this is that non-localizations that have zero confidence will always be lower than lesion localizations with non-zero confidence. If a threshold much greater than 1 is used, the results may change significantly. For the purpose of illustrating the package, the confidence threshold has been set to 0.

raw2 <- applyflowchart(raw1, PTcodes = c("91"), PTCthreshold = 0, includeA=T,includeC=T)

A subset of the function's results is in the table below. Note that many new fields have been derived, and there's been additional standardization applied to the file. Some function parameters have been stored in the data to make future processing possible.

kable(head(raw2))

Function plotPTCbins()

A plotting function has been created to examine the confidence score usage of the final flow chart data file (raw2 above). The function will jitter the ROI PTC by bin. The mean and count will be displayed on the figure by default. There is currently only one parameter option for the figure--a descriptive name for the primary task can be supplied. The default simply using the abbreviation PT. The plot objects are standard ggplot2 objects and can be stored as objects and files.

plot_ptcbins <- plotPTCbins(raw2)
plot_ptcbins

Function figuresave()

A small wrapper function has been created to allow for saving of the figures produced by the package. This was created to provide some uniform scaling of the saved images. It's not essential to use.

figuresave("PTCbyBins.PNG", plot_ptcbins)

Prepare JAFROC File

The current implementation of the RJafroc package allows for the use of an Excel data file. There are three required sheets: a sheet for the lesion identification (and lesion weights), lesion localizations (TP) and non-localizations (FP). The makeJAFROCfile function takes care of this for you. In addition, the function exports a fourth sheet that contains the linkage of the reference dataset set to the original study ID and the JAFROC de-identified ID. To de-identify the JAFROC file, this fourth sheet may need to be deleted.

Function makeJAFROCfile()

jaffile <- "demoJAFROC.xlsx"  # it's useful to use a string to store the file name 

makeJAFROCfile(flowchartdata=raw2,outputfile=jaffile)

TIP: In case there is some debate as to what the minimum confidence score for inclusion in the analysis or whether or not Bins A and C should be included, simply run the applyflowchart() and makeJAFROCfile() functions under various settings. You can specify unique file names for the JAFROC excel files and present the analysis in multiple ways.

Function plotPTC()

Once the JAFROC file is completed, it can be used to produce an additional confidence score visulization. This figure pools the reader data to show how the confidence varies between lesion localizations and non-lesion localizations. Normally, one sees clear separation in the score usage. If this is not occurring, the FOMs will be attenuated. This might be an indication the wrong codes were used above when specifying the primary task lesion ID. Or, it could be a hard reading task.

plotPTC(jaffile)

Traditional Sensitivity and Specificity Analyses

In general, the derivation of sensitivity and specifity for a free response study is not well defined. We have developed adaptations that allow for the reporting of sensitivity and specificity.

Function derivesensspec() -- Patient-level Analysis

1. Sensitivity: At least one lesion localization in a case with at least one target to be detected. The primary task confidence for the reader detection is variable and can be used to determine the minimum confidence level for the detection to be considered. The denominator is the number of cases with at least one target lesion.

2. Specificity: No non-lesion localizations above the threshold in a case that does not have any targets. The denominator is the number of cases without any target lesions.

These definitions are readily applied by the derivesensspec function. There are parameters for the JAFROC file to be used along with the confidence thresholds. The function computes the reader-specific performance along with pooled (GEE) estimates. Below is an example of function in use. The resulting dataframe is frequently saved to an Excel format using the WriteXLS function.

specex<-derivesensspec(jaffile,sensitivitythreshold=10, specificitythreshold=10)
kable(specex)

TIP: It has been noted that this function uses thresholds differently than the applyflowchart function (>= vs. >). This may change in future releases. For now, note how the function writes out the logic rule to the datatable. You may adjust the threshold selection as needed.

Function derivelesionsensitivty() -- Lesion-level Sensitivity Analysis

1. Lesion Sensitivity: Lesion detected with confidence above the threshold. The denominator is the number of target lesions in the study.

2. Specificity: Not applicable

The function has the same general call as the previous function. In this case, the sensitivity threshold has been set to -1 so that any co-localization, regardless of classification or confidence is considered. It should be noted, that if the applyflowcart delta has been specified with a value >0, there may be an undesired interaction with the sensitivity and specificity functions.

lesionsens<-derivelesionsens(jaffile, sensitivitythreshold=-1)
kable(lesionsens)

Plotting Functions -- plotthreshold_reader() and plotthreshold_modality()

To support additional examination of the lesion detection data and the impact of the threshold selection, additional plotting functions are available.

The plotthreshold_reader function can plot, by default, the GEE estimates for lesion sensitivity by modality. Or, one can specify an individual reader number to see how an individual reader performances across modalities.

plotthreshold_reader(jaffile)

If you desire to relabel / reorder the modality names, a call such as this can be considered. For this example, reader 1 data is selected. Note, the function calls the derivelesionsensitivty function above. This function formats the reader ID as a text field with a leading zero to ensure ordering of tables. As a result, reader "1" needs to be entered as readerid="01".

plotthreshold_reader(jaffile, readerid="01", addFactor=T, flevels=flevels, flabels =flabels)

A second plot function looks at reader performance for a particular modality. The plotthreshold_modality function will plot all readers along with an optional GEE line. There are some parameter options that will allow changing of line colors if desired. This example requests the study's routine clinical dose. The GEE estimate is included by default.

Both threshold plots allow for a logical parameter getminmax that will search the TP ratings and determine the minimum of the maximum ratings for each reader. This is to mitigate cases where some, or all readers, do not detect at least one lesion with confidence of 100. This can introduce an all zeros condition and cause a function error. If an error is reported, try setting getminmax = F. This will censor the figure at the min-max for the readers. The figure below shows this behavior.

plotthreshold_modality(jaffile, modality=refdosename, modalitylabel=refdoselabel,getminmax = T)

Function essentiallesions()

A final set of functions for lesion-level analyses exists to derive some of the summaries used for analyses. The first concept is look at is the set of essential lesions. An essential lesion is a lesions that was detected by the majority of readers at the routine dose. The function requires the specification of the routine dose along with an export file to store the list in. You can save the resulting data frame for working in memory too.

el<-essentiallesions(jaffile, refdosename,"essentiallesions.xlsx")

kable(head(el,n=20))

Function derivedections()

For a slightly more general look at the detections, the derivedections function can be used. This function utilizes the primary JAFROC file without any additional specification of a threshold (i.e., reader ROIs that are above the PTC threshold for the study). The resulting data frame has the reference lesion ID along with the number of detections by modality.

detections <- derivedetections(jaffile)


kable(head(detections, n=20))

Part 2: JAFROC Analysis

Conducting the JAFROC Analysis

The analysis requires the RJafroc package. The logic and use of that package is left to that package's maintainer. In this section, a few additional plotting and results summaries are highlighted. There are opportunities to do JAFROC and JAFROC1 FOMs. The difference is the "1" version includes non-localization in cases with target lesions in the calculation of the FOMs; otherwise they are excluded. The makeJAFROCfile also creates the default weights such that each lesion is weighted based on the number in the case. For example, if a case has two lesions, the weights for each lesion is 1/2. Different weights can be assigned but sometimes the justification is hard to describe. In general, the weighted analysis is preferred to avoid one case with a large number of lesions having a unreasonable influence on the analysis. The package RJafroc has a vignette that contains many more details about the particular aspects of the analysis that is being conducted here.

The RJafroc package has a way to save a text file out with the classically formatted JAFROC results format. This functionality is not shown activated in this vignette but the code is provided below.

#### Conduct JAFROC analysis
library(RJafroc)
task.mrmc<-RJafroc::DfReadDataFile(fileName=jaffile)
jafroc_results<-RJafroc::StSignificanceTesting(task.mrmc, FOM="wJAFROC", method="DBMH")

## output the text report: commented out to not run here
# RJafroc::UtilOutputReport(method="DBMH", FOM="wJAFROC", dataset=task.mrmc, showWarnings=FALSE)

The resulting jafroc_results object has several elements that need to be parsed for graphical purposes. One of the things that needs to be done is labeling and ordering the modalities. The following code summarizes some approaches. The factor labels get reversed based on the order you want them to show in the final plot.

Function prepareestimates()

The prepareestimates function is a helper function. It is used to derive the raw estimates along with the deltas. These are stored as a list. Also note, that the RJafroc package produces estimates for a variety of analysis considerations by default. The most commonly used options will either be the "Fixed Reader, Random Case" (FRRC, small number of readers) or the "Random Reader, Random Case" (RRRC, larger number of readers). For this study, we will show both.

## From the jafroc results above, go process them for subsequent plotting. 
estimates<- prepareestimates(jafroc_results, refdosename)
## this is a list of two elements: the estimates [[1]] and the deltas to the reference [[2]]
est<- estimates[[1]]
# Format the results for easy printing
## Note: variable strings flevels and flabels are defined above
est$Treatment <- factor(est$Treatment, levels=flevels, labels=flabels)
## For GGPLOT, we are going to reverse the order so the reference is at the top -- this is of course, optional
est$rTreatment <- with(est, factor(Treatment, levels = rev(levels(Treatment))))

Function plotFOM()

Plot the figure using the regular ordering of the factor. This figure requests the fixed reader, random case modelling assumption.

plotFOM(est,"FRRC","Treatment")

Now plot the reverse ordered treatments

plotFOM(est, "FRRC", "rTreatment")

This is a repeat of the above figure only this time electing to use the random reader, random case model.

plotFOM(est, "RRRC", "rTreatment")

As with the previous plotting functions, the plot functions can be saved to an object to facilitate saving using figuresave().

Note: The figures for RRRC and FRRC will look very similar in most cases. The CIs will be slightly different as they are based on different variance estimation approaches. It is a good idea to use descriptive file names to help keep the various versions straight, particularly if multiple jafroc files are generated based on different threshold values.

Function plotDeltaFOM()

Repeat the above now with delta FOM (i.e., change in FOM estimates relative to the reference dose)

del<- estimates[[2]]
del$Treatment <- with(del, factor(Treatment, levels=flevels, labels=flabels))
del$rTreatment <- with(del, factor(Treatment, levels = rev(levels(Treatment))))
### Notice the delta FOM figure has some additional options to control the formatting
delta_fom_rrrc<-plotDeltaFOM(del,"RRRC",refdoselabel,-.05,thetreatment = "rTreatment")

delta_fom_rrrc

Next Steps

At this point, the primary analysis is complete. One could further modify the plot objects, which are standard ggplot2 objects, to refine the layout. Or, conduct sensitivity analyses based on various modeling assumptions. Future package updates may include additional graphical features.

Please post comments and questions to github @rickeycarter

rickeycarter/mrmctools documentation built on May 6, 2019, 8:59 p.m.