diff --git a/.gitignore b/.gitignore index 55d5eb60..ead34c4c 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ /Rpackage/iMRMC_1.2.2.zip /Rpackage/iMRMC_1.2.3.tar.gz /Rpackage/iMRMC_1.2.3.zip +/Rpackage/iMRMC/man/*.Rd diff --git a/Rpackage/iMRMC/DESCRIPTION b/Rpackage/iMRMC/DESCRIPTION index af04252c..98683742 100644 --- a/Rpackage/iMRMC/DESCRIPTION +++ b/Rpackage/iMRMC/DESCRIPTION @@ -21,9 +21,9 @@ Description: Package development and documentation is at . License: CC0 Encoding: UTF-8 -SystemRequirements: Java JDK 1.7 or higher +SystemRequirements: Java JDK 1.7 or higher LazyData: true Depends: R (>= 3.5.0) -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.1 Imports: stats Suggests: testthat diff --git a/Rpackage/iMRMC/NAMESPACE b/Rpackage/iMRMC/NAMESPACE index aa297347..1c37cc07 100644 --- a/Rpackage/iMRMC/NAMESPACE +++ b/Rpackage/iMRMC/NAMESPACE @@ -1,4 +1,4 @@ -# Generated by roxygen2: do not edit by hand +# Generated by roxygen2: do not edit by hand export(convertDF) export(convertDFtoDesignMatrix) @@ -13,11 +13,15 @@ export(getMRMCscore) export(getWRBM) export(init.lecuyerRNG) export(laBRBM) +export(laBRBM.anova) export(laWRBM) +export(laWRBM.anova) export(renameCol) export(roc2binary) export(sim.gRoeMetz) export(sim.gRoeMetz.config) +export(sim.new.Hierarchical) +export(sim.new.Hierarchical.config) export(simMRMC) export(simRoeMetz.example) export(successDFtoROCdf) @@ -25,4 +29,9 @@ export(uStat11.conditionalD) export(uStat11.jointD) export(undoIMRMCdf) import(parallel) +importFrom(stats,aov) importFrom(stats,qnorm) +importFrom(stats,rbeta) +importFrom(stats,rgamma) +importFrom(stats,rnorm) +importFrom(stats,var) diff --git a/Rpackage/iMRMC/R/la.anova.R b/Rpackage/iMRMC/R/la.anova.R new file mode 100644 index 00000000..be276dd0 --- /dev/null +++ b/Rpackage/iMRMC/R/la.anova.R @@ -0,0 +1,256 @@ +#' @title MRMC Analysis of Limits of Agreement using ANOVA +#' @description +#' These two functions calculate two types of Limits of Agreement using ANOVA: Within-Reader Between-Modality(WRBM) +#' and Between-Reader Between-Modality(BRBM). The 95\% confidence interval of the mean difference is also provided. +#' The ANOVA method are realized either by applying stats::aov or by matrix multiplication. See more details below +#' about the model structure. +#' +#' @details +#' Suppose the score from reader j for case k under modality i is\eqn{X_{ijk}}, then the difference score from the +#' same reader for the same cases under two different modalities is \eqn{Y_{jk} = X_{1jk} - X_{2jk}}. +#' \itemize{ +#' \item\code{laWRBM} use two-way random effect ANOVA to analyze the difference scores \eqn{Y_{jk}}. The model +#' is \eqn{Y_{jk}=\mu + R_j + C_k + \epsilon_{jk}}, where \eqn{R_j} and \eqn{C_k} are random effects for readers +#' and cases. The variance of mean and individual observation is expressed as the linear combination of the MS +#' given by ANOVA. +#' \item\code{laBRBM} use three-way mixed effect ANOVA to analyze the scores \eqn{X_{ijk}}. The model is given by +#' \eqn{X_{ijk}=\mu + R_j + C_k + m_i + RC_{jk} + mR_{ij} + mC_{ik} + \epsilon_{ijk}}, where \eqn{R_j} and +#' \eqn{C_k} are random effects for readers and cases and \eqn{m_i} is a fixed effect for modality. The variance +#' of mean and individual observation is expressed as the linear combination of the MS given by ANOVA. +#' } +#' +#' @name la.anova +#' +#' @param df +#' Data frame of observations, one per row. Columns identify random effects, fixed effects, +#' and the observation. Namely, +#' \describe{ +#' \item{readerID}{The factor corresponding to the different readers in the study. +#' The readerID is treated as a random effect.} +#' \item{caseID}{The factor corresponding to the different cases in the study. +#' The caseID is treated as a random effect.} +#' \item{modalityID}{The factor corresponding to the different modalities in the study. +#' The modalityID is treated as a fixed effect.} +#' \item{score}{The score given by the reader to the case for the modality indicated.} +#' } +#' +#' @param modalitiesToCompare +#' The factors identifying the modalities to compare. It should be at length 2. Default +#' \code{modalitiesToCompare = c("testA","testB")} +#' +#' @param keyColumns +#' Identify the factors corresponding to the readerID, caseID, modalityID, and score +#' (or alternative random and fixed effects). Default \code{keyColumns = c("readerID", "caseID", +#' "modalityID", "score")} +#' +#' @param if.aov +#' Boolean value to determine whether using aov function to do ANOVA. Default \code{if.aov = TRUE} +#' +#' @return +#' +#' A dataframe with one row. Each column is as following: +#' \describe{ +#' \item{meanDiff}{The mean of difference score.} +#' \item{var.MeanDiff}{The variance of mean difference score} +#' \item{var.1obs}{The variance of a single WRBM/BRBM difference score} +#' \item{ci95meanDiff.bot}{Lower bound of 95\% CI for the mean difference score. \code{meanDiff+ +#' 1.96*sqrt(var.MeanDiff)}} +#' \item{ci95meanDiff.top}{Upper bound of 95\% CI for the mean difference score. \code{meanDiff- +#' 1.96*sqrt(var.MeanDiff)}} +#' \item{la.bot}{Lower bound of WRBM/BRBM Limits of Agreement. \code{meanDiff+2*sqrt(var.1obs)}} +#' \item{la.top}{Upper bound of WRBM/BRBM Limits of Agreement. \code{meanDiff-2*sqrt(var.1obs)}} +#' } +#' +#' The two function shows the same 95\% CI for the mean difference score, but difference Limits of Agreements. +#' +#' @importFrom stats aov var +#' +#' @export +#' +#' @examples +#' library(iMRMC) +#' # Create an MRMC data frame +#' # Refer to Gallas2014_J-Med-Img_v1p031006 +#' simRoeMetz.config <- sim.gRoeMetz.config() +#' +#' # Simulate data +#' df.MRMC <- sim.gRoeMetz(simRoeMetz.config) +#' +#' # Compute Limits of Agreement +#' laWRBM_result <- laWRBM.anova(df.MRMC) +#' print(laWRBM_result) +#' laBRBM_result <- laBRBM.anova(df.MRMC) +#' print(laBRBM_result) +#' +laWRBM.anova <- function(df, modalitiesToCompare = c("testA","testB"), + keyColumns = c("readerID", "caseID", "modalityID", "score"), + if.aov = TRUE +) { + + if (length(modalitiesToCompare) != 2) { + print(paste("length(modalitiesToCompare) =", length(modalitiesToCompare))) + stop("ERROR: modalitiesToCompare should have 2 elements.") + } + + df <- data.frame(readerID = factor(df[[keyColumns[1]]]), + caseID = factor(unclass(df[[keyColumns[2]]])), #unclass for changing Ord.factor to unordered + modalityID = factor(df[[keyColumns[3]]]), + score = df[[keyColumns[4]]]) + + + # Parse out data frames for each modality + df.A <- df[df$modalityID == modalitiesToCompare[1], ] + df.B <- df[df$modalityID == modalitiesToCompare[2], ] + + nReader <- length(unique(df.A$readerID)) + nCase <- length(unique(df.A$caseID)) + + # Compute the difference score between the modalities + diff.df.all <- merge(df.A,df.B, by = c("readerID","caseID")) + diff.df.all$score <- diff.df.all$score.x - diff.df.all$score.y + diff.df <- subset(diff.df.all, select = c("readerID","caseID","score")) + + if(if.aov){ + # Do ANOVA + fit <- aov(score ~ readerID + caseID, data = diff.df) + + # Extract MS + MS <- summary(fit)[[1]]$`Mean Sq` + MSA <- MS[1] + MSB <- MS[2] + sigma2 <- MS[3] + }else{ + diff.mat <- t(convertDFtoScoreMatrix(droplevels(diff.df))) + + MSA <- var(rowMeans(diff.mat)) * nCase + MSB <- var(colMeans(diff.mat)) * nReader + SST <- var(array(diff.mat)) * (nCase*nReader - 1) + sigma2 <- (SST - MSA * (nReader - 1) - MSB * (nCase - 1)) / (nReader - 1) / (nCase - 1) + } + + + # Limit of agreement result + meanDiff <- mean(diff.df$score) + var.MeanDiff <- (MSA + MSB - sigma2)/nReader/nCase + var.1obs <- (nReader * MSA + nCase * MSB + (nReader * nCase - nReader - nCase) * sigma2)/nReader/nCase + # if use 3-way anova to get var.1obs + # var.1obs <- 2*(varRM + varCM + sigma2) + + la.bot <- meanDiff - 2 * sqrt(var.1obs) + la.top <- meanDiff + 2 * sqrt(var.1obs) + + ci95meanDiff.bot <- meanDiff + qnorm(.025) * sqrt(var.MeanDiff) + ci95meanDiff.top <- meanDiff + qnorm(.975) * sqrt(var.MeanDiff) + + result2 <- data.frame( + meanDiff = meanDiff, var.MeanDiff = var.MeanDiff, var.1obs = var.1obs, + ci95meanDiff.bot = ci95meanDiff.bot, ci95meanDiff.top = ci95meanDiff.top, + la.bot = la.bot, la.top = la.top ) + + return(result2) + + +} + + + + +#' @rdname la.anova +#' +#' @export +#' + +laBRBM.anova <- function(df, modalitiesToCompare = c("testA","testB"), + keyColumns = c("readerID", "caseID", "modalityID", "score"), + if.aov = TRUE +) { + + if (length(modalitiesToCompare) != 2) { + print(paste("length(modalitiesToCompare) =", length(modalitiesToCompare))) + stop("ERROR: modalitiesToCompare should have 2 elements.") + } + + df <- data.frame(readerID = factor(df[[keyColumns[1]]]), + caseID = factor(unclass(df[[keyColumns[2]]])), #unclass for changing Ord.factor to unordered + modalityID = factor(df[[keyColumns[3]]]), + score = df[[keyColumns[4]]]) + + + # Parse out data frames for each modality + df.A <- df[df$modalityID == modalitiesToCompare[1], ] + df.B <- df[df$modalityID == modalitiesToCompare[2], ] + + nM <- 2 + nR <- nlevels(droplevels(df.A)$readerID) + nC <- nlevels(droplevels(df.A)$caseID) + + if(if.aov){ + # apply aov function to do 3-way ANOVA + df_2Modality <- rbind(df.A, df.B) + + fit <- aov(score ~ readerID + caseID + modalityID + readerID:caseID + readerID:modalityID + caseID:modalityID, + data = df_2Modality) + # Extract MS + MS <- summary(fit)[[1]]$`Mean Sq` + MSR <- MS[1] + MSC <- MS[2] + MSM <- MS[3] + MSRC <- MS[4] + MSRM <- MS[5] + MSCM <- MS[6] + sigma2 <- MS[7] + + }else{ + # Generate 3-dimentional matrix, reader x case x modality + rcm_mat <-array(0,dim=c(nR, nC, nM)) + + rcm_mat[,,1] <- t(convertDFtoScoreMatrix(droplevels(df.A))) + rcm_mat[,,2] <- t(convertDFtoScoreMatrix(droplevels(df.B))) + + MSR <- var(rowMeans(rcm_mat)) * nC * nM + MSC <- var(rowMeans(colMeans(rcm_mat))) * nR * nM + MSM <- var(colMeans(rcm_mat, dims = 2)) * nR * nC + + SSR <- MSR * (nR - 1) + SSC <- MSC * (nC - 1) + SSM <- MSM * (nM - 1) + + SSRC <- var(array(rowMeans(rcm_mat, dims=2))) * nM * (nR*nC-1) - SSR - SSC + SSRM <- var(array(colMeans(aperm(rcm_mat, c(2,1,3))))) * nC * (nR*nM-1) - SSR - SSM + SSCM <- var(array(colMeans(rcm_mat))) * nR * (nC*nM-1) - SSM - SSC + + MSRC <- SSRC / (nR-1) / (nC-1) + MSRM <- SSRM / (nR-1) / (nM-1) + MSCM <- SSCM / (nC-1) / (nM-1) + + SST <- var(array(rcm_mat)) * (nC*nR*nM - 1) + sigma2 <- (SST - SSR -SSC - SSM - SSRC - SSRM - SSCM) / (nR-1) / (nC-1) / (nM-1) + } + + varR <- (MSR + sigma2 - MSRC - MSRM) / nC / nM + varC <- (MSC + sigma2 - MSRC - MSCM) / nR / nM + #M2 <- (MSM + sigma2 - MSRM - MSCM) / nR / nC *(nM-1) /2 #fixed effect no variance + varRC <- (MSRC - sigma2) / nM + varRM <- (MSRM - sigma2) / nC + varCM <- (MSCM - sigma2) / nR + + # Limit of agreement result + meanDiff <- mean(df.A$score-df.B$score) + var.MeanDiff <- 2*(varRM/nR+varCM/nC+sigma2/nR/nC) + var.1obs <- 2*(varR + varRC + varRM + varCM + sigma2) + + la.bot <- meanDiff - 2 * sqrt(var.1obs) + la.top <- meanDiff + 2 * sqrt(var.1obs) + + ci95meanDiff.bot <- meanDiff + qnorm(.025) * sqrt(var.MeanDiff) + ci95meanDiff.top <- meanDiff + qnorm(.975) * sqrt(var.MeanDiff) + + result2 <- data.frame( + meanDiff = meanDiff, var.MeanDiff = var.MeanDiff, var.1obs = var.1obs, + ci95meanDiff.bot = ci95meanDiff.bot, ci95meanDiff.top = ci95meanDiff.top, + la.bot = la.bot, la.top = la.top ) + + return(result2) +} + + diff --git a/Rpackage/iMRMC/R/sim.new.Hierarchical.R b/Rpackage/iMRMC/R/sim.new.Hierarchical.R new file mode 100644 index 00000000..27136808 --- /dev/null +++ b/Rpackage/iMRMC/R/sim.new.Hierarchical.R @@ -0,0 +1,247 @@ +#' Simulate an MRMC data set comparing two modalities by a hierarchical model +#' +#' @description +#' This procedure simulates an MRMC data set for a MRMC agreement study comparing two +#' modalities.It is a hierarchical model consists of two interaction terms: reader-case +#' interaction and modality-reader-case-replicate interaction. Both the interaction +#' terms are conditional normal distributed, with the case(-related) factor contributing +#' to the conditional mean and the reader(-related) factor contributing to the conditional +#' variance. +#' +#' @details +#' The model is as the following structure: +#' X.ijkl = mu + m.i + RC.jk + mRCE.ijkl +#' \itemize{ +#' \item mu = grand mean +#' \item m.i = modalities (levels: A and B) +#' \item RC.jk|R.j,C.k ~ N(C.k, R.j) reader-case interaction term +#' \item mRCE.ijkl|mR.ij,mC.ik ~ N(mC.ik, mR.ij) modality-reader-case-replicate term +#' \item C.k and mC.ik are Normal/beta distributed +#' \item R.j and mR.ij are Inverse-Gamma distributed +#' } +#' +#' @param config [list] of simulation parameters: +#' \itemize{ +#' \item Experiment labels and size +#' \itemize{ +#' \item modalityID: [vector] label modality A and B. +#' \item nR: [num] number of readers +#' \item nC: [num] number of cases +#' \item C_dist: [chr] distribution of the case. Default \code{C_dist="normal"} +#' } +#' \item Mean and fixed effects: +#' \itemize{ +#' \item mu: [num] grand mean +#' \item tau_A: [num] modality A +#' \item tau_B: [num] modality B +#' } +#' \item Reader-case interaction term +#' \itemize{ +#' \item sigma_C: [num] variance of case factor (if \code{C_dist="normal"}) +#' \item a_C: [num] alpha for distribution of case (if \code{C_dist="beta"}) +#' \item b_C: [num] beta for distribution of case (if \code{C_dist="beta"}) +#' \item alpha_R: [num] shape parameter for reader +#' \item beta_R: [num] scale parameter for reader +#' } +#' \item Modality-reader-case-replicate interaction term for modality A +#' \itemize{ +#' \item sigma_C.A: [num] variance of case factor (if \code{C_dist="normal"}) +#' \item a_C.A: [num] alpha for distribution of case (if \code{C_dist="beta"}) +#' \item b_C.A: [num] beta for distribution of case (if \code{C_dist="beta"}) +#' \item alpha_R.A: [num] shape parameter for reader +#' \item beta_R.A: [num] scale parameter for reader +#' } +#' \item Modality-reader-case-replicate interaction term for modality B +#' \itemize{ +#' \item sigma_C.B: [num] variance of case factor (if \code{C_dist="normal"}) +#' \item a_C.B: [num] alpha for distribution of case (if \code{C_dist="beta"}) +#' \item b_C.B: [num] beta for distribution of case (if \code{C_dist="beta"}) +#' \item alpha_R.B: [num] shape parameter for reader +#' \item beta_R.B: [num] scale parameter for reader +#' } +#' \item Scales for the case related terms and interaction terms +#' \itemize{ +#' \item C_scale: [num] weight for the case factor +#' \item RC_scale: [num] weight for the reader-case interaction term +#' \item tauC_scale: [num] weight for the modality-case term +#' \item tauRCE_scale: [num] weight for the modality-reader-case-replicate interaction term +#' } +#' } +#' @param R [vector] fix the reader factor across different simulation. Default \code{R = NULL} +#' @param AR [vector] fix the modality-reader interaction. Default \code{AR = NULL} +#' @param BR [vector] fix the modality-reader interaction. Default \code{BR = NULL} +#' @param is.within [bol] whether the data are within-modality (AR==BR). Default \code{is.within=FALSE} +#' +#' @return df [data.frame] with nR x nC x 2 rows including +#' \itemize{ +#' \item readerID: [Factor] w/ nR levels "reader1", "reader2", ... +#' \item caseID: [Factor] w/ nC levels "case1", "case2", ... +#' \item modalityID: [Factor] w/ 1 level config$modalityID +#' \item score: [num] reader score +#' } +#' +#' @importFrom stats rbeta rgamma rnorm +#' +#' @export +#' +# @examples +# # Create a sample configuration object +# config <- sim.new.Hierarchical.config() +# # Simulate an MRMC ROC data set +# dFrame <- sim.new.Hierarchical(config) + + +sim.new.Hierarchical = function(config,R = NULL,AR = NULL,BR = NULL,is.within=FALSE) { + + # Initialize ---- + nR = config$nR + nC = config$nC + modalityID = config$modalityID + C_dist = config$C_dist + + mu = config$mu + tau_A = config$tau_A + tau_B = config$tau_B + + alpha_R = config$alpha_R + alpha_R.A = config$alpha_R.A + alpha_R.B = config$alpha_R.B + + beta_R = config$beta_R + beta_R.A = config$beta_R.A + beta_R.B = config$beta_R.B + + C_scale = config$C_scale + RC_scale = config$RC_scale + tauC_scale = config$tauC_scale + tauRCE_scale = config$tauRCE_scale + + + # C and tauC, different case distribution + if(C_dist == 'normal'){ + + C = matrix(rep(rnorm(nC,0,config$sigma_C),nR),nR,nC,byrow =TRUE) + tauC.A = matrix(rep(rnorm(nC,0,config$sigma_C.A),nR),nR,nC,byrow =TRUE) + tauC.B = matrix(rep(rnorm(nC,0,config$sigma_C.B),nR),nR,nC,byrow =TRUE) + + }else if(C_dist == 'beta'){ + + C = matrix(rep(rbeta(nC,config$a_C,config$b_C),nR),nR,nC,byrow =TRUE) + tauC.A = matrix(rep(rbeta(nC,config$a_C.A,config$b_C.A),nR),nR,nC,byrow =TRUE) + tauC.B = matrix(rep(rbeta(nC,config$a_C.B,config$b_C.B),nR),nR,nC,byrow =TRUE) + + } + + # RC + if(is.null(R)){ + R = rgamma(nR, alpha_R, rate = beta_R) + } + var_RC = matrix(rep(R, nC), nR, nC) + + RC = apply(1/sqrt(var_RC), c(1,2), rnorm, n=1, mean=0) # Normal variable with mean 0 and variance inv-Gamma distributed + + + # tauRCE + ## Mod A/ first replicate + if(is.null(AR)){ + AR = rgamma(nR, alpha_R.A, rate = beta_R.A) + } + var_tauRCE.A = matrix(rep(AR, nC), nR, nC) + + tauRCE.A = apply(1/sqrt(var_tauRCE.A), c(1,2), rnorm, n= 1, mean= 0) + + ## Mod B/ second replicate + if(is.within){ + BR = AR + }else{ + if(is.null(BR)){ + BR = rgamma(nR, alpha_R.B, rate = beta_R.B) + } + } + + var_tauRCE.B = matrix(rep(BR, nC), nR, nC) + + tauRCE.B = apply(1/sqrt(var_tauRCE.B), c(1,2), rnorm, n= 1, mean= 0) + + # Aggregate ---- + + modA = mu + tau_A + C * C_scale + RC * RC_scale + tauC.A * tauC_scale + tauRCE.A * tauRCE_scale + modB = mu + tau_B + C * C_scale + RC * RC_scale + tauC.B * tauC_scale + tauRCE.B * tauRCE_scale + + # Five column format ---- + + df = as.data.frame(matrix(0,nrow=nR*nC*2,ncol=4)) + colnames(df) = c("caseID","readerID","modalityID","score") + + df$score = c(as.vector(t(modA)),as.vector(t(modB))) + df$caseID = paste0("Case",rep(1:nC,nR*2)) + df$readerID = as.factor(paste0("reader",rep(rep(1:nR,each=nC),2))) + df$modalityID = as.factor(rep(c(modalityID[1],modalityID[2]),each=nR*nC)) + + + return(df) +} + +#' Create a configuration object for the sim.new.Hierarchical program +#' +#'#' @description +#' This function creates a configuration object for the Hierarchical +#' simulation model to be used as input for the sim.new.Hierarchical program. +#' +#' @details If no arguments, this function returns a default simulation +#' configuration for sim.new.Hierarchical +#' +#' @param nR [num] Number of readers. Default \code{nR = 5} +#' @param nC [num] Number of cases. Default \code{nC = 100} +#' @param modalityID [vector] List of modalityID. Default \code{modalityID = c("testA", "testA*")} +#' @param C_dist [chr] Distribution of the case. Default \code{C_dist="normal"} +#' @param mu [num] grand mean. Default \code{mu = 0} +#' @param tau_A [num] modality A effect. Default \code{tau_A = 0} +#' @param tau_B [num] modality B effect. Default \code{tau_B = 0} +#' @param alpha_R [num] shape parameter for reader. Default \code{alpha_R = 10} +#' @param beta_R [num] scale parameter for reader. Default \code{beta = 1} +#' @param sigma_C [num] variance of case factor (if \code{C_dist="normal"}). Default \code{sigma_C = 1} +#' @param a_C [num] alpha for distribution of case (if \code{C_dist="beta"}). Default \code{a_C = 0.8} +#' @param b_C [num] beta for distribution of case (if \code{C_dist="beta"}). Default \code{b_C = 3} +#' @param sigma_tauC [num] variance of modality-case (if \code{C_dist="normal"}). Default \code{sigma_tauC = 1} +#' @param alpha_tauR [num] shape parameter for modality-reader. Default \code{alpha_tauR = 10} +#' @param beta_tauR [num] scale parameter for modality-reader. Default \code{beta_tauR = 1} +#' @param C_scale [num] weight for the case factor. Default \code{C_scale = 1} +#' @param RC_scale [num] weight for the reader-case interaction term. Default \code{RC_scale = 1} +#' @param tauC_scale [num] weight for the modality-case term. Default \code{tauC_scale = 1} +#' @param tauRCE_scale [num] weight for the modality-reader-case-replicate interaction term. Default \code{tauRCE_scale = 1} +#' +#' @return config [list] Refer to the sim.new.Hierarchical input variable +#' @export +#' + +sim.new.Hierarchical.config = function (nR = 5, nC = 100, modalityID = c("testA", "testA*"), C_dist = 'normal', + mu = 0, tau_A = 0, tau_B = 0, alpha_R = 10, beta_R = 1, + sigma_C = 1, a_C = 0.8, b_C = 3, sigma_tauC = 1, + alpha_tauR = 10, beta_tauR = 1, + C_scale = 1, RC_scale = 1,tauC_scale = 1, tauRCE_scale = 1) +{ + if (C_dist == 'normal'){ + config <- list(modalityID = modalityID, nR = nR, nC = nC, # sample size + C_dist = C_dist, mu = mu, tau_A = tau_A, tau_B = tau_B, # distribution and mean + sigma_C = sigma_C, alpha_R = alpha_R, beta_R = beta_R, # parameter for [RC] + sigma_C.A = sigma_tauC, alpha_R.A = alpha_tauR, beta_R.A = beta_tauR, # parameter for [tauRCE].A + sigma_C.B = sigma_tauC, alpha_R.B = alpha_tauR, beta_R.B = beta_tauR, # parameter for [tauRCE].B + C_scale = C_scale, RC_scale = RC_scale, # scale for [RC] + tauC_scale = tauC_scale, tauRCE_scale = tauRCE_scale) # scale for [tauRCE] + + }else if(C_dist == 'beta'){ + config <- list(modalityID = modalityID, nR = nR, nC = nC, # sample size + C_dist = C_dist, mu = mu, tau_A = tau_A, tau_B = tau_B, # distribution and mean + a_C = a_C, b_C = b_C, alpha_R = alpha_R, beta_R = beta_R, # parameter for [RC] + a_C.A = a_C, b_C.A = b_C, alpha_R.A = alpha_tauR, beta_R.A = beta_tauR,# parameter for [tauRCE].A + a_C.B = a_C, b_C.B = b_C, alpha_R.B = alpha_tauR, beta_R.B = beta_tauR,# parameter for [tauRCE].B + C_scale = C_scale, RC_scale = RC_scale, # scale for [RC] + tauC_scale = tauC_scale, tauRCE_scale = tauRCE_scale) # scale for [tauRCE]) + }else{ + print(paste0("C_dist = ", C_dist)) + stop("ERROR: C_dist should be either normal or beta.") + } + + return(config) +} diff --git a/Rpackage/iMRMC/inst/extra/CompareSimulationModels.Rmd b/Rpackage/iMRMC/inst/extra/CompareSimulationModels.Rmd new file mode 100644 index 00000000..5f6623a8 --- /dev/null +++ b/Rpackage/iMRMC/inst/extra/CompareSimulationModels.Rmd @@ -0,0 +1,70 @@ +--- +title: "Simulate MRMC Data" +output: pdf_document +--- + + +```{r setup, include=FALSE} +library(knitr) +library(iMRMC) +source("../../R/sim.new.Hierarchical.R") +knitr::opts_chunk$set(echo = TRUE) +``` + +There are two functions in this package that can simulate Multi-reader Multi-case (MRMC) data. sim.new.Hierarchical is the function to simulate the MRMC agreement data with no binary truth state. The sim.gRoeMetz is the function to simulate MRMC ROC data. + +## 1. Using the new hierachical model to simulate the MRMC agreement data + +The following is how to use the sim.new.Hierachcial.R to simulate MRMC data with no truth state, the MRMC agreement data + +```{r} + +# configuration +nR = 5 #number of readers +nC = 100 #number of cases + +config <- sim.new.Hierarchical.config(nR=nR, nC=nC, modalityID = c("testA","testB")) + +# simulate MRMC study +set.seed(1, kind = "L'Ecuyer-CMRG") +dFrame.newH <- sim.new.Hierarchical(config) + +# check the first and last few lines of the simulated dataframe +head(dFrame.newH) +tail(dFrame.newH) +``` + +We simulated 5 readers and 100 cases for 2 modalities, so the total number of scores is 5x100x2=1000. + +## 2. Using sim.gRoeMetz in iMRMC package to simulate MRMC ROC data + +The following is how to use the sim.gRoeMetz.R to simulate MRMC data with truth state, the MRMC ROC data + +```{r} + +# configuration +nR = 5 #number of readers +nC.neg = 50 #number of positive cases +nC.pos = 50 #number of negative cases + +config <- sim.gRoeMetz.config(nR = nR, nC.neg = nC.neg, nC.pos = nC.neg) + +# simulate MRMC study +set.seed(1, kind = "L'Ecuyer-CMRG") +dFrame.gRM <- sim.gRoeMetz(config) + +# check the first and last few lines of the simulated dataframe +head(dFrame.gRM) +tail(dFrame.gRM) +``` + +The simulated data starts with truth state of each case, and followed by the reading scores from each of the readers. Since we simulate 50 positive cases and 50 negative cases, there are 100 lines for the truth and 100x5x2 = 1000 lines for the scores from 5 readers for 2 modalities. + +To combine the truth data and the reader scores and change the data to a dataframe with 5 columns: readerID, caseID, modalityID, score, truth, we can use the function undoIMRMCdf in iMRMC pacakge: +```{r} +dFrame.gRM.2 <- undoIMRMCdf(dFrame.gRM) + +# check the first and last few lines of the simulated dataframe +head(dFrame.gRM.2) +tail(dFrame.gRM.2) +``` diff --git a/Rpackage/iMRMC/inst/extra/CompareSimulationModels.pdf b/Rpackage/iMRMC/inst/extra/CompareSimulationModels.pdf new file mode 100644 index 00000000..cc7f73d3 Binary files /dev/null and b/Rpackage/iMRMC/inst/extra/CompareSimulationModels.pdf differ diff --git a/Rpackage/iMRMC/inst/extra/LOA_ANOVA.Rmd b/Rpackage/iMRMC/inst/extra/LOA_ANOVA.Rmd new file mode 100644 index 00000000..6b004969 --- /dev/null +++ b/Rpackage/iMRMC/inst/extra/LOA_ANOVA.Rmd @@ -0,0 +1,231 @@ +--- +title: "Using ANOVA to Estimate Limits of Agreement for MRMC study" +output: pdf_document +--- + +```{r setup, include=FALSE} +library(knitr) +library(iMRMC) +source("../../R/la.anova.R") +source('../../R/sim.new.Hierarchical.R') +knitr::opts_chunk$set(echo = TRUE) +``` + +## 1. Definition + +Suppose $X_{ijk}$ denotes the score for case $k$ $(k=1,...,K)$ from the reader $j$ $(j=1,..,J)$ under modality $i$ $(i=1,2)$ in a Multi-reader Multi-case study comparing two different modalities, where $i=1$ and $i=2$ indicate test modality and reference modality respectively. The difference score between the two modalities is given by $D_{jj'k}^{12}=X_{1jk}-X_{2j'k}$. If $j=j'$, the difference score $D_{jj,k}^{12}$, or simply $Y_{jk}$, denotes the within-reader between-modality (WRBM) difference.That is, the difference score from the same reader under different modalities. Given the mean difference $\overline{D_{WR}^{12}}=E[Y_{jk}]$ and the variance of the difference $V_{WR}^{12}=Var[Y_{jk}]$, the WRBM limits of agreement is defined as: +\[ +\overline{D_{WR}^{12}}\pm2\sqrt{V_{WR}^{12}} +\] +When $j\neq j'$, the difference score $D_{jj'k}^{12}$ denotes the between-reader between-modality (BRBM) difference. Similarly, we have the mean difference $\overline{D_{BR}^{12}}=E[D_{jj'k}^{12}]$ and the variance of the difference $V_{BR}^{12}=Var[D_{jj'k}^{12}]$. Thus, the BRBM limits of agreement is defined as: +\[ +\overline{D_{BR}^{12}}\pm2\sqrt{V_{BR}^{12}} +\] + +To construct the WRBM and BRBM limits of agreement, we need to calculate $\overline{D_{WR}^{12}}$, $V_{WR}^{12}$, $\overline{D_{BR}^{12}}$, and $V_{BR}^{12}$. The two mean differences are easy to estimate. +\[ +\hat{\overline{D_{WR}^{12}}}=\frac{1}{JK}\sum_j\sum_kD_{jj,k}^{12}=\frac{1}{JK}\sum_j\sum_k(X_{1jk}-X_{2jk})=\overline{X_1..}-\overline{X_2..} +\] +where $\overline{X_i..}=\frac{1}{JK}\sum_j\sum_kX_{ijk}$ $i=1,2$ denotes the average score across all the readers and cases for a single modality. Similarly, +\[ +\hat{\overline{D_{BR}^{12}}}=\frac{1}{J(J-1)K}\sum_j\sum_{j\neq j'}\sum_kD_{jj'k}^{12}=\frac{1}{J(J-1)K}\sum_j\sum_{j\neq j'}\sum_k(X_{1jk}-X_{2j'k})=\overline{X_1..}-\overline{X_2..} +\] +Therefore, $\hat{\overline{D_{WR}^{12}}}=\hat{\overline{D_{BR}^{12}}}=\overline{X_1..}-\overline{X_2..}$. The WRBM and BRBM limits of agreement will be different only by $V_{WR}^{12}$ and $V_{BR}^{12}$. In the following tow sections, we will discuss how to use two-way random effect ANOVA to estimate $V_{WR}^{12}$ and use three-way mixed effect ANOVA to estimate $V_{BR}^{12}$ + +## 2. Using two-way random effect ANOVA to estimate $V_{WRBM}$ + +To estimate $V_{WR}^{12}=Var[Y_{jk}]$, we build up a two-way random effect model for the WRBM difference $Y_{jk}$ +\[ +Y_{jk} = \mu+R_j+C_k+\varepsilon_{jk} +\] +where $R_j\sim N(0,\sigma_R^2)$, $C_k\sim N(0,\sigma_C^2)$, and $\varepsilon_{jk}\sim N(0,\sigma_\varepsilon^2)$ are independent random variables. Then, the variance of $Y_{jk}$ can be expressed as +\[ +V_{WR}^{12}=Var[Y_{jk}]=Var(R_j+C_k+\varepsilon_{jk})=\sigma_R^2+\sigma_C^2+\sigma_\varepsilon^2 +\] + +The two-way random effect ANOVA table is given by + +Source DF Sum of Square (SS) Mean Square (MS) E(MS) +------- ----------- ---------------------------- ---------------------- -------- +Reader $J-1$ $SSR=K\sum_j(\overline{Y_{j.}}-\overline{Y_{..}})^2$ $MSR=SSR/(J-1)$ $\sigma_\varepsilon^2+K\sigma_R^2$ +Case $K-1$ $SSC=J\sum_k(\overline{Y_{.k}}-\overline{Y_{..}})^2$ $MSR=SSC/(K-1)$ $\sigma_\varepsilon^2+J\sigma_C^2$ +Error $(J-1)(K-1)$ $SSE=SST-SSR-SSC$ $MSE=SSE/(J-1)(K-1)$ $\sigma_\varepsilon^2$ +Total $JK-1$ $SST=\sum_j\sum_k(Y_{jk}-\overline{Y_{..}})^2$ +------- ----------- ---------------------------- ---------------------- -------- + +In the table above, $\overline{Y_{j.}}=\frac{1}{K}\sum_kY_{jk}$, $\overline{Y_{.k}}=\frac{1}{J}\sum_jY_{jk}$, $\overline{Y_{..}}=\overline{D_{WR}^{12}}$ are the marginal and overall mean of difference score. Hence, the sum of squares (SS) and mean squares (MS) can be calculated from the data. From the last column of the ANOVA table, we can find the relationship between the variance components ($\sigma_R^2$, $\sigma_C^2$, $\sigma_\varepsilon^2$) and the mean squares. So the unbiased estimation for the variance components are +\[ +\begin{aligned} +\hat\sigma_\varepsilon^2&=MSE\\ +\hat\sigma_R^2&=\frac{MSR-MSE}{K}\\ +\hat\sigma_C^2&=\frac{MSC-MSE}{J} +\end{aligned} +\] + +Therefore, the estimation of variance of $Y_{jk}$, +\[ +\hat V_{WR}^{12}=\hat Var[Y_{jk}]=\hat\sigma_R^2+\hat\sigma_C^2+\hat\sigma_\varepsilon^2=\frac{1}{JK}(J*MSR+K*MSC+(JK-J-K)*MSE) +\] + +Following is an example of using function laWRBM.anova to compute the WRBM limits of agreement for the simulated MRMC data +```{r simulation & laWRBM.anova} +# Simulate MRMC data +config <- sim.new.Hierarchical.config(modalityID = c("testA","testB")) +set.seed(1) +data.sim <- sim.new.Hierarchical(config) + +# Using ANOVA to calculate WRBM limits of agreement +laWRBM.anova_result <- laWRBM.anova(data.sim) +print(laWRBM.anova_result) +``` + +The var.1obs in the result is the estimation of variance of WRBM difference, $\hat V_{WR}^{12}$.The limits of agreement is given by [la.bot, la.top]. Since both the ANOVA method and U-statistics method provide unbiased estimation to the variance components, the above result is the same as the one calculated by using U-statistics method. laWRBM is a function in iMRMC package that using U-statistics method to construct WRBM limits of agreement. + +```{r laWRBM} +# Compare the result with laWRBR in iMRMC package +print(laWRBM(data.sim,modalitiesToCompare = c("testA","testB"))) +``` + +## 3. Using three-way mixed effect ANOVA to calculate $V_{WRBM}$ + +To estimate $V_{BR}^{12}=Var[D_{jj'k}^{12}]$, we build up a three-way mixed effect model for the score $X_{ijk}$ +\[ +X_{ijk} = \mu+m_i+R_j+C_k+RC_{jk}+mR_{ij}+mC_{ik}+\varepsilon_{ijk} +\] +where $m_i$ denotes the fixed effect for modality and $\sum_im_i=0$,and the other variables are independently normal distributed: $R_j\sim N(0,\sigma_R^2)$, $C_k\sim N(0,\sigma_C^2)$, $RC_{jk}\sim N(0,\sigma_{RC}^2)$, $mR_{ij}\sim N(0,\sigma_{mR}^2)$, $mC_{ik}\sim N(0,\sigma_{mC}^2)$, and $\varepsilon_{ijk}\sim N(0,\sigma_\varepsilon^2)$ . + +The mixed effect model we applied is the unrestrict mixed model. Compared to the restricted model, the interaction term of the mixed effect and random effect in the unrestricted model will not be correlated at the same random level. This makes the model easier to work with. + +Then, the BRBM difference $D_{jj'k}^{12}$ can be expressed as +\[ +D_{jj'k}^{12}=X_{1jk}-X_{2j'k}=m_1-m_2+R_j-R_j'+RC_{jk}-RC_{j'k}+mR_{1j}-mR_{2j'}+mC_{1k}-mC_{2k}+\varepsilon_{1jk}-\varepsilon_{2j'k} +\] +The variance of $D_{jj'k}^{12}$ is as following +\[ +V_{BR}^{12}=Var[D_{jj'k}^{12}]=2\sigma_R^2+2\sigma_{RC}^2+2\sigma_{mR}^2+2\sigma_{mC}^2+2\sigma_\varepsilon^2 +\] + +The three-way mixed effect ANOVA table is given by + +Source DF Sum of Square (SS) E(MS) +------- ----------- ---------------------------- ----------------------- +Modality $I-1$ $SSM=JK\sum_i(\overline{X_{i..}}-\overline{X_{...}})^2$ $\sigma_\varepsilon^2+K\sigma_{mR}^2+J\sigma_{mC}^2+\frac{JK}{I-1}\sum_im_i^2$ +Reader $J-1$ $SSR=IK\sum_j(\overline{X_{.j.}}-\overline{X_{...}})^2$ $\sigma_\varepsilon^2+I\sigma_{RC}^2+K\sigma_{mR}^2+IK\sigma_R^2$ +Case $K-1$ $SSC=IJ\sum_k(\overline{X_{..k}}-\overline{X_{...}})^2$ $\sigma_\varepsilon^2+I\sigma_{RC}^2+J\sigma_{mC}^2+IJ\sigma_C^2$ +Reader:Case $(J-1)(K-1)$ $SSRC=I\sum_j\sum_k(\overline{X_{.jk}}-\overline{X_{.j.}}-\overline{X_{..k}}+\overline{X_{...}})^2$ $\sigma_\varepsilon^2+I\sigma_{RC}^2$ +Reader:Modality $(J-1)(I-1)$ $SSMR=K\sum_i\sum_j(\overline{X_{ij.}}-\overline{X_{i..}}-\overline{X_{.j.}}+\overline{X_{...}})^2$ $\sigma_\varepsilon^2+K\sigma_{mR}^2$ +Case:Modality $(K-1)(I-1)$ $SSMC=J\sum_i\sum_k(\overline{X_{i.k}}-\overline{X_{i..}}-\overline{X_{..k}}+\overline{X_{...}})^2$ $\sigma_\varepsilon^2+J\sigma_{mC}^2$ +Error $df_{E}$ $SSE=SST-other SS$ $\sigma_\varepsilon^2$ +Total $IJK-1$ $SST=\sum_i\sum_j\sum_k(X_{ijk}-\overline{X_{...}})^2$ +------- ----------- ---------------------------- ----------------------- + +where $\overline{X_{i..}}$, $\overline{X_{.j.}}$, $\overline{X_{..k}}$, $\overline{X_{ij.}}$, $\overline{X_{i.k}}$, $\overline{X_{.jk}}$, and $\overline{X_{...}}$ are marginal or overall mean of the score $X_{ijk}$. The $df_{E}$ denotes the degree of freedom for the error, $df_{E}=IJK-IJ-JK-IK+I+J+K-1$. Similar to the two-way ANOVA table, the last column shows the relationship between the variance components with the mean squares. So the unbiased estimation for the variance components are +\[ +\begin{aligned} +\hat\sigma_\varepsilon^2&=MSE\\ +\hat\sigma_{RC}^2&=\frac{MSRC-MSE}{I}\\ +\hat\sigma_{mC}^2&=\frac{MSMC-MSE}{J}\\ +\hat\sigma_{mR}^2&=\frac{MSMR-MSE}{K}\\ +\hat\sigma_{R}^2&=\frac{MSR-MSRC-MSMR+MSE}{IK}\\ +\hat\sigma_{C}^2&=\frac{MSC-MSRC-MSMC+MSE}{IJ}\\ +\end{aligned} +\] + +Therefore, the estimation of variance of $Y_{jj'k}$ is +\[ +\begin{aligned} +\hat V_{BR}^{12}&=\hat Var[D_{jj'k}^{12}]=2\hat\sigma_R^2+2\hat\sigma_{RC}^2+2\hat\sigma_{mR}^2+2\hat\sigma_{mC}^2+2\hat\sigma_\varepsilon^2\\ +&=\frac{2}{IJK}(J*MSR+J(K-1)*MSRC+J(I-1)*MSMR\\ +&\ \ \ \ +IK*MSMC+(IJK-IJ-IK-JK+J)*MSE) +\end{aligned} +\] + +Following is an example of using function laBRBM.anova to compute the BRBM limits of agreement for the simulated MRMC data +```{r laBRBM.anova} +# Using ANOVA to calculate BRBM limits of agreement +laBRBM.anova_result <- laBRBM.anova(data.sim) +print(laBRBM.anova_result) +``` + +We can compare the above result with the result from laBRBM function in iMRMC package. Again, the ANOVA method and the U-statistics method shows the same result. +```{r laBRBM} +# Compare the result with laWRBR in iMRMC package +print(laBRBM(data.sim,modalitiesToCompare = c("testA","testB"))) +``` + +## 4. Relationship between the two-way random effect ANOVA and the three-way mixed effect ANOVA + +Since the WRBM difference score, $Y_{jk}=X_{1jk}-X_{2jk}$, is a linear combination of the individual score, the variance of $Y_{jk}$ can also be expressed by the mean squares in the three-way mixed effect ANOVA. First, we put the three-way ANOVA model into the WRBM difference score definition. +\[ +Y_{jk}=X_{1jk}-X_{2jk}=m_1-m_2+mR_{1j}-mR_{2j}+mC_{1k}-mC_{2k}+\varepsilon_{1jk}-\varepsilon_{2jk} +\] +Then, the variance of $Y_{jk}$ is as following: +\[ +V_{WR}^{12}=Var[Y_{jk}]=2\sigma_{mR}^2+2\sigma_{mC}^2+2\sigma_\varepsilon^2 +\] +Therefore, by plug in the unbiased estimation of the variance components, we get the estimation of $V_{WRBM}$ +\[ +\hat V_{WR}^{12}=2\hat\sigma_{mR}^2+2\hat\sigma_{mC}^2+2\hat\sigma_\varepsilon^2 +=\frac{2}{JK}(J*MSMR+K*MSMC+(JK-J-K)*MSE) +\] +By comparing this result with the one from the two-way random effect ANOVA, we notice that there is a linear relationship between the sum of squares in the two ANOVA models. In the following we use subscirpt $2w$ to denote the MS or SS for the two-way ANOVA. +\[ +\begin{aligned} +SSR_{2w}&=K\sum_j(\overline{Y_{j.}}-\overline{Y_{..}})^2\\ +&=K\sum_j\left(\overline{X_{1j.}}-\overline{X_{2j.}}-\overline{X_{1..}}+\overline{X_{2..}}\right)^2\\ +&=K\sum_j[(\overline{X_{1j.}}-\overline{X_{1..}}-\overline{X_{.j.}}+\overline{X_{...}})^2+(\overline{X_{2j.}}-\overline{X_{2..}}-\overline{X_{.j.}}+\overline{X_{...}})^2\\ +&\ \ \ \ \ -2(\overline{X_{1j.}}-\overline{X_{1..}}-\overline{X_{.j.}}+\overline{X_{...}})(\overline{X_{2j.}}-\overline{X_{2..}}-\overline{X_{.j.}}+\overline{X_{...}})]\\ +&=K\sum_j\left[\sum_i(\overline{X_{ij.}}-\overline{X_{i..}}-\overline{X_{.j.}}+\overline{X_{...}})^2 +2\cdot\frac{1}{2}\cdot\frac{1}{2}(\overline{Y_{j.}}-\overline{Y_{..}})^2\right]\\ +&=SSMR_{3w}+\frac{1}{2}SSR_{2w} +\end{aligned} +\] +Thus, +\[ +SSR_{2w}=2SSMR_{3w} +\] +Similarly, +\[ +SSC_{2w}=2SSMC_{3w} +\] +\[ +SST_{2w}=2\sum_i\sum_j\sum_k(X_{ijk}-\overline{X_{i..}}-\overline{X_{.jk}}+\overline{X_{...}})^2 +\] +For the total sum of square in the three-way ANOVA, +\[ +\begin{aligned} +SST_{3w}&=\sum_i\sum_j\sum_k(X_{ijk}-\overline{X_{...}})^2\\ +&=\sum_i\sum_j\sum_k\left[(X_{ijk}-\overline{X_{i..}}-\overline{X_{.jk}}+\overline{X_{...}})+(\overline{X_{i..}}-\overline{X_{...}})+(\overline{X_{.jk}}-\overline{X_{...}})\right]^2\\ +&=\sum_i\sum_j\sum_k(X_{ijk}-\overline{X_{i..}}-\overline{X_{.jk}}+\overline{X_{...}})^2+SSM_{3w}+I\sum_j\sum_k(\overline{X_{.jk}}-\overline{X_{...}})^2 +\end{aligned} +\] +The last term on the right hand side of the above fomular can continue to be decomposed as +\[ +\begin{aligned} +&\ \ \ \ \ I\sum_j\sum_k(\overline{X_{.jk}}-\overline{X_{...}})^2\\ +&=I\sum_j\sum_k\left[(\overline{X_{.jk}}-\overline{X_{.j.}}-\overline{X_{..k}}+\overline{X_{...}})+(\overline{X_{.j.}}-\overline{X_{...}})+(\overline{X_{..k}}-\overline{X_{...}})\right]^2\\ +&=SSRC_{3w}+SSR_{3w}+SSC_{3w} +\end{aligned} +\] +Thus, +\[ +\begin{aligned} +SSE_{2w} &= SST_{2w}-SSR_{2w}-SSC_{2w}\\ +&=2(SST_{3w}-SSM_{3w}-SSR_{3w}-SSC_{3w}-SSRC_{3w})-2SSMR_{3w}-2SSMC_{3w}=2SSE_{3w} +\end{aligned} +\] +Since $I=2$, the degree of freedom for $SSMR_{3w}$ is $(J-)(I-1)=J-1$ and for $SSMC_{3w}$ is $(K-)(I-1)=K-1$. $df_E=IJK-IJ-JK-IK+I+J+K-1=JK-J-K+1$ is the degree of freedom for $SSE_{3w}$. We have the same mean square relationship as that for the sum of squares +\[ +\begin{aligned} +MSR_{2w}&=2MSMR_{3w}\\ +MSC_{2w}&=2MSMR_{3w}\\ +MSE_{2w}&=2MSE_{3w}\\ +\end{aligned} +\] +Therefore, +\[ +\begin{aligned} +\hat V_{WRBM}&=\frac{2}{JK}(J*MSMR_{3w}+K*MSMC_{3w}+(JK-J-K)*MSE_{3w})\\ +&=\frac{1}{JK}(J*MSR_{2w}+K*MSC_{2w}+(JK-J-K)*MSE_{2w}) +\end{aligned} +\] + diff --git a/Rpackage/iMRMC/inst/extra/LOA_ANOVA.pdf b/Rpackage/iMRMC/inst/extra/LOA_ANOVA.pdf new file mode 100644 index 00000000..9ab04621 Binary files /dev/null and b/Rpackage/iMRMC/inst/extra/LOA_ANOVA.pdf differ diff --git a/Rpackage/iMRMC/inst/extra/man/la.anova.html b/Rpackage/iMRMC/inst/extra/man/la.anova.html new file mode 100644 index 00000000..852ae445 --- /dev/null +++ b/Rpackage/iMRMC/inst/extra/man/la.anova.html @@ -0,0 +1,149 @@ +R: MRMC Analysis of Limits of Agreement using ANOVA + + + + +
la.anovaR Documentation
+ +

MRMC Analysis of Limits of Agreement using ANOVA

+ +

Description

+ +

These two functions calculate two types of Limits of Agreement using ANOVA: Within-Reader Between-Modality(WRBM) +and Between-Reader Between-Modality(BRBM). The 95% confidence interval of the mean difference is also provided. +The ANOVA method are realized either by applying stats::aov or by matrix multiplication. See more details below +about the model structure. +

+ + +

Usage

+ +
+laWRBM.anova(
+  df,
+  modalitiesToCompare = c("testA", "testB"),
+  keyColumns = c("readerID", "caseID", "modalityID", "score"),
+  if.aov = TRUE
+)
+
+laBRBM.anova(
+  df,
+  modalitiesToCompare = c("testA", "testB"),
+  keyColumns = c("readerID", "caseID", "modalityID", "score"),
+  if.aov = TRUE
+)
+
+ + +

Arguments

+ + + + + + + + + + +
df +

Data frame of observations, one per row. Columns identify random effects, fixed effects, +and the observation. Namely, +

+ +
+
readerID

The factor corresponding to the different readers in the study. +The readerID is treated as a random effect.

+
+
caseID

The factor corresponding to the different cases in the study. +The caseID is treated as a random effect.

+
+
modalityID

The factor corresponding to the different modalities in the study. +The modalityID is treated as a fixed effect.

+
+
score

The score given by the reader to the case for the modality indicated.

+
+
+
modalitiesToCompare +

The factors identifying the modalities to compare. It should be at length 2. Default +modalitiesToCompare = c("testA","testB")

+
keyColumns +

Identify the factors corresponding to the readerID, caseID, modalityID, and score +(or alternative random and fixed effects). Default keyColumns = c("readerID", "caseID", +"modalityID", "score")

+
if.aov +

Boolean value to determine whether using aov function to do ANOVA. Default if.aov = TRUE

+
+ + +

Details

+ +

Suppose the score from reader j for case k under modality i isX_{ijk}, then the difference score from the +same reader for the same cases under two different modalities is Y_{jk} = X_{1jk} - X_{2jk}. +

+ + + + + +

Value

+ +

A dataframe with one row. Each column is as following: +

+ +
+
meanDiff

The mean of difference score.

+
+
var.MeanDiff

The variance of mean difference score

+
+
var.1obs

The variance of a single WRBM/BRBM difference score

+
+
ci95meanDiff.bot

Lower bound of 95% CI for the mean difference score. meanDiff+ + 1.96*sqrt(var.MeanDiff)

+
+
ci95meanDiff.top

Upper bound of 95% CI for the mean difference score. meanDiff- + 1.96*sqrt(var.MeanDiff)

+
+
la.bot

Lower bound of WRBM/BRBM Limits of Agreement. meanDiff+2*sqrt(var.1obs)

+
+
la.top

Upper bound of WRBM/BRBM Limits of Agreement. meanDiff-2*sqrt(var.1obs)

+
+
+ +

The two function shows the same 95% CI for the mean difference score, but difference Limits of Agreements. +

+ + +

Examples

+ +
+library(iMRMC)
+# Create an MRMC data frame
+# Refer to Gallas2014_J-Med-Img_v1p031006
+simRoeMetz.config <- sim.gRoeMetz.config()
+
+# Simulate data
+df.MRMC <- sim.gRoeMetz(simRoeMetz.config)
+
+# Compute Limits of Agreement
+laWRBM_result <- laWRBM.anova(df.MRMC)
+print(laWRBM_result)
+laBRBM_result <- laBRBM.anova(df.MRMC)
+print(laBRBM_result)
+
+
+ + + diff --git a/Rpackage/iMRMC/inst/extra/man/sim.new.Hierarchical.config.html b/Rpackage/iMRMC/inst/extra/man/sim.new.Hierarchical.config.html new file mode 100644 index 00000000..bfe4bfc0 --- /dev/null +++ b/Rpackage/iMRMC/inst/extra/man/sim.new.Hierarchical.config.html @@ -0,0 +1,140 @@ +R: Create a configuration object for the sim.new.Hierarchical... + + + + +
sim.new.Hierarchical.configR Documentation
+ +

Create a configuration object for the sim.new.Hierarchical program

+ +

Description

+ +

#' @description +This function creates a configuration object for the Hierarchical +simulation model to be used as input for the sim.new.Hierarchical program. +

+ + +

Usage

+ +
+sim.new.Hierarchical.config(
+  nR = 5,
+  nC = 100,
+  modalityID = c("testA", "testA*"),
+  C_dist = "normal",
+  mu = 0,
+  tau_A = 0,
+  tau_B = 0,
+  alpha_R = 10,
+  beta_R = 1,
+  sigma_C = 1,
+  a_C = 0.8,
+  b_C = 3,
+  sigma_tauC = 1,
+  alpha_tauR = 10,
+  beta_tauR = 1,
+  C_scale = 1,
+  RC_scale = 1,
+  tauC_scale = 1,
+  tauRCE_scale = 1
+)
+
+ + +

Arguments

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
nR +

[num] Number of readers. Default nR = 5

+
nC +

[num] Number of cases. Default nC = 100

+
modalityID +

[vector] List of modalityID. Default modalityID = c("testA", "testA*")

+
C_dist +

[chr] Distribution of the case. Default C_dist="normal"

+
mu +

[num] grand mean. Default mu = 0

+
tau_A +

[num] modality A effect. Default tau_A = 0

+
tau_B +

[num] modality B effect. Default tau_B = 0

+
alpha_R +

[num] shape parameter for reader. Default alpha_R = 10

+
beta_R +

[num] scale parameter for reader. Default beta = 1

+
sigma_C +

[num] variance of case factor (if C_dist="normal"). Default sigma_C = 1

+
a_C +

[num] alpha for distribution of case (if C_dist="beta"). Default a_C = 0.8

+
b_C +

[num] beta for distribution of case (if C_dist="beta"). Default b_C = 3

+
sigma_tauC +

[num] variance of modality-case (if C_dist="normal"). Default sigma_tauC = 1

+
alpha_tauR +

[num] shape parameter for modality-reader. Default alpha_tauR = 10

+
beta_tauR +

[num] scale parameter for modality-reader. Default beta_tauR = 1

+
C_scale +

[num] weight for the case factor. Default C_scale = 1

+
RC_scale +

[num] weight for the reader-case interaction term. Default RC_scale = 1

+
tauC_scale +

[num] weight for the modality-case term. Default tauC_scale = 1

+
tauRCE_scale +

[num] weight for the modality-reader-case-replicate interaction term. Default tauRCE_scale = 1

+
+ + +

Details

+ +

If no arguments, this function returns a default simulation +configuration for sim.new.Hierarchical +

+ + +

Value

+ +

config [list] Refer to the sim.new.Hierarchical input variable +

+ + + diff --git a/Rpackage/iMRMC/inst/extra/man/sim.new.Hierarchical.html b/Rpackage/iMRMC/inst/extra/man/sim.new.Hierarchical.html new file mode 100644 index 00000000..e000cfeb --- /dev/null +++ b/Rpackage/iMRMC/inst/extra/man/sim.new.Hierarchical.html @@ -0,0 +1,220 @@ +R: Simulate an MRMC data set comparing two modalities by a... + + + + +
sim.new.HierarchicalR Documentation
+ +

Simulate an MRMC data set comparing two modalities by a hierarchical model

+ +

Description

+ +

This procedure simulates an MRMC data set for a MRMC agreement study comparing two +modalities.It is a hierarchical model consists of two interaction terms: reader-case +interaction and modality-reader-case-replicate interaction. Both the interaction +terms are conditional normal distributed, with the case(-related) factor contributing +to the conditional mean and the reader(-related) factor contributing to the conditional +variance. +

+ + +

Usage

+ +
+sim.new.Hierarchical(config, R = NULL, AR = NULL, BR = NULL, is.within = FALSE)
+
+ + +

Arguments

+ + + + + + + + + + + + +
config +

[list] of simulation parameters: +

+ +
    +

  • Experiment labels and size +

    + +
      +

    • modalityID: [vector] label modality A and B. +

      +
      • +

    • nR: [num] number of readers +

      +
      • +

    • nC: [num] number of cases +

      +
      • +

    • C_dist: [chr] distribution of the case. Default C_dist="normal" +

      +
    + +
    • +

  • Mean and fixed effects: +

    + +
      +

    • mu: [num] grand mean +

      +
      • +

    • tau_A: [num] modality A +

      +
      • +

    • tau_B: [num] modality B +

      +
    + +
    • +

  • Reader-case interaction term +

    + +
      +

    • sigma_C: [num] variance of case factor (if C_dist="normal") +

      +
      • +

    • a_C: [num] alpha for distribution of case (if C_dist="beta") +

      +
      • +

    • b_C: [num] beta for distribution of case (if C_dist="beta") +

      +
      • +

    • alpha_R: [num] shape parameter for reader +

      +
      • +

    • beta_R: [num] scale parameter for reader +

      +
    + +
    • +

  • Modality-reader-case-replicate interaction term for modality A +

    + +
      +

    • sigma_C.A: [num] variance of case factor (if C_dist="normal") +

      +
      • +

    • a_C.A: [num] alpha for distribution of case (if C_dist="beta") +

      +
      • +

    • b_C.A: [num] beta for distribution of case (if C_dist="beta") +

      +
      • +

    • alpha_R.A: [num] shape parameter for reader +

      +
      • +

    • beta_R.A: [num] scale parameter for reader +

      +
    + +
    • +

  • Modality-reader-case-replicate interaction term for modality B +

    + +
      +

    • sigma_C.B: [num] variance of case factor (if C_dist="normal") +

      +
      • +

    • a_C.B: [num] alpha for distribution of case (if C_dist="beta") +

      +
      • +

    • b_C.B: [num] beta for distribution of case (if C_dist="beta") +

      +
      • +

    • alpha_R.B: [num] shape parameter for reader +

      +
      • +

    • beta_R.B: [num] scale parameter for reader +

      +
    + +
    • +

  • Scales for the case related terms and interaction terms +

    + +
      +

    • C_scale: [num] weight for the case factor +

      +
      • +

    • RC_scale: [num] weight for the reader-case interaction term +

      +
      • +

    • tauC_scale: [num] weight for the modality-case term +

      +
      • +

    • tauRCE_scale: [num] weight for the modality-reader-case-replicate interaction term +

      +
    + +
+
R +

[vector] fix the reader factor across different simulation. Default R = NULL

+
AR +

[vector] fix the modality-reader interaction. Default AR = NULL

+
BR +

[vector] fix the modality-reader interaction. Default BR = NULL

+
is.within +

[bol] whether the data are within-modality (AR==BR). Default is.within=FALSE

+
+ + +

Details

+ +

The model is as the following structure: +X.ijkl = mu + m.i + RC.jk + mRCE.ijkl +

+ + + + + +

Value

+ +

df [data.frame] with nR x nC x 2 rows including +

+ + + + + +