#' Boschek's test of reliability
#' @description Implementation of the Boschek's test of reliability as described in
#' "A Probabilistic “True-error” Model for the Reliability Analysis of Nominal Classifications".
#' @references Boschek, P., Vranka, M. A., Bartoš, F., A Probabilistic “True-error” Model
#' for the Reliability Analysis of Nominal Classifications (manuscript under preparation)
#' @import nleqslv
#' @param data A matrix or data frame containing observations. Either in long format - with
#' columns corresponding to raters and rows to subjects rated; or in an aggregated manner - with
#' one column labeled "frequency" containing counts of combinations of ratings which are
#' specified in remaining columns. See Examples for more information.
#' @param C A matrix containing the acceptable probabilities of an incorrect (or correct) classification
#' that specifies the null hypothesis.
#' @param exact_prob Whether the exact probability estimates should be computed. Only possible
#' when the number of categories is 2.
#'
#'
#' @examples ### 2 raters, 2 categories and data in long format
#' # load data1
#' data(d1)
#'
#' # create classification matrix
#' C1 <- matrix(c(.9, .1,
#' .1, .9),byrow = TRUE, nrow = 2, ncol = 2)
#'
#' # fit the model
#' m1 <- boschek.test(data = d1,C = C1)
#'
#' # quickly inspect the model
#' m1
#'
#' # for more details use summary function
#' summary(m1)
#'
#' # to print residuals and aggregated data matrix use residuals function
#' residuals(m1)
#'
#' # in order to print all columns and rows set argument print.all = TRUE
#' residuals(m1, print.all = TRUE)
#'
#' # or access the residuals directly from the fitted object
#' m1$model_A$table
#' m1$model_B$table
#'
#' ### ------------------------------------------------- ###
#' ### 3 raters, 3 categories and data in long format
#' # load data1
#' data(d2)
#'
#' # note that column containing the number of rattings is called frequency
#' head(d2)
#'
#' # create classification matrix
#' C2 <- matrix(c(0.85, 0.1, 0.05,
#' 0.1 , 0.8, 0.1 ,
#' 0.05, 0.1, 0.85 ),byrow = TRUE, nrow = 3, ncol = 3)
#'
#' # fit the model
#' m2 <- boschek.test(data = d2,C = C2)
#'
# check input
#' @export boschek.test
boschek.test <- function(data, C, exact_prob = TRUE){
Q <- C
call <- match.call()
if(any(apply(Q,1,sum) != 1) | any(apply(Q,2,sum) != 1))stop("Q misspecified - doesn't sum to 1")
if(any(diag(Q) <= .5))stop("Assumptions not met - diagonal values of Q must be larger than .5")
if(!is.data.frame(data) & !is.matrix(data))stop("Incompatible data type")
if(any(is.na(data)))stop("NA in data is not allowed")
# reshape data
if(any(colnames(data) == "frequency")){
for(i in c(1:ncol(data))[colnames(data) != "frequency"]){
data <- data[order(data[,i]),]
}
N_flatt <- cbind.data.frame(data[,colnames(data) != "frequency"],"frequency" = data[,"frequency"])
combinations <- apply(as.matrix(data[,colnames(data) != "frequency"]),2,as.numeric)
NR <- ncol(data) - 1
NK <- max(unlist(combinations))
if(nrow(data) != NK^NR)stop("Supply frequency for all combinations of possible raters evaluations")
}else{
N_flatt <- as.data.frame(table(as.data.frame(data)), responseName = "frequency")
combinations <- as.matrix(do.call(expand.grid,lapply(dim(table(as.data.frame(data))),seq)))
NR <- ncol(data)
NK <- max(dim(table(data)))
}
# solve T
X <- nleqslv::nleqslv(rep(1/NK,NK-1),.solve_T, N_flatt = N_flatt, combinations = combinations, Q = Q, NK = NK, NR = NR,
control = list(xtol = 1e-15, ftol = 1e-15))
# compute residuals
N_flatt <- cbind(N_flatt,"expected_frequency" = .solve_P(X$x, combinations = combinations, Q = Q, NK = NK, NR = NR)*sum(N_flatt$frequency))
N_flatt$standardized_residuals <- with(N_flatt,(frequency - expected_frequency)/sqrt(expected_frequency))
N <- sum(N_flatt$frequency)
# data description
data_description <- list(
"observations" = N,
"number_of_categories" = NK,
"number_of_raters" = NR
)
# model B
combinations_help <- apply(combinations,1,function(x){paste(table(c(x,1:NK))-1,sep = "", collapse = "|")})
N_flattB <- NULL
symetry_chisq <- NULL
for(com in unique(combinations_help)){
temp_c <- c(1:nrow(N_flatt))[combinations_help == com]
N_flattB <- rbind.data.frame(N_flattB,
cbind.data.frame(com,
"frequency" = sum(N_flatt[temp_c,"frequency"]),
"expected_frequency" = sum(N_flatt[temp_c,"expected_frequency"])))
# test of symetry
if(!grepl(NR, com)){
symetry_chisq <- c(symetry_chisq,(N_flatt[temp_c,"frequency"]-mean(N_flatt[temp_c,"frequency"]))^2/mean(N_flatt[temp_c,"frequency"]))
}
}
colnames(N_flattB)[1] <- paste(c("evaluation (",
paste(1:NK,sep = "", collapse = "|"),
")"),sep = "", collapse = "")
N_flattB$standardized_residuals <- with(N_flattB,(frequency - expected_frequency)/sqrt(expected_frequency))
# test of symetry
test_of_symmetry <- list(
"X^2" = sum(symetry_chisq),
"df" = NK^NR-(factorial(NK+NR-1)/(factorial(NR)*factorial(NK-1))),
"p.value" = stats::pchisq(sum(symetry_chisq),NK^NR-(factorial(NK+NR-1)/(factorial(NR)*factorial(NK-1))),lower.tail = F)
)
### solve exactly in case of two categories
probabilities_estimates <- NULL
if(NK == 2 & exact_prob){
TQ <- nleqslv::nleqslv(c(
N_flattB$frequency[1]/(N_flattB$frequency[1] + N_flattB$frequency[nrow(N_flattB)]),
(sum(N_flattB$frequency) - N_flattB$frequency[1] - N_flattB$frequency[nrow(N_flattB)])/(NR*sum(N_flattB$frequency))),
.solve_TQ, N_flatt = N_flattB, NR = NR, control = list(xtol = 1e-15, ftol = 1e-15))$x
# check that the solution satisfies boundaries
if((TQ[1] >= 1 | TQ[1] <= 0) | TQ[2] >= .5 | TQ[1] <= 0){
# if not, try alternative starting values
TQ_seq <- sapply(seq(0.01,0.49,0.01),function(q){
nleqslv::nleqslv(c(
N_flattB$frequency[1]/(N_flattB$frequency[1] + N_flattB$frequency[nrow(N_flattB)]),q),
.solve_TQ, N_flatt = N_flattB, NR = NR, control = list(xtol = 1e-15, ftol = 1e-15))$x
})
TQ_seq <- TQ_seq[,!(TQ_seq[1,] >= 1 | TQ_seq[1,] <= 0) & !(TQ_seq[2,] >= .5 | TQ_seq[2,] <= 0)]
if(ncol(TQ_seq) == 0)stop("Solution to routine B not found.")
TQ_seq <- TQ_seq[,!duplicated(round(TQ_seq[2,],3))]
if(is.matrix(TQ_seq))stop("Multiple solutions to routine B found.")
}
t1 <- TQ[1]
t2 <- 1 - TQ[1]
# compute SE
t1_se <- .compute_SE(t1,TQ[2], N)
t2_se <- .compute_SE(t2,TQ[2], N)
probabilities_estimates <- rbind(
"T1" = c(
"Estimate" = t1,
"SE" = t1_se,
"95% l.CI" = ifelse(t1-1.96*t1_se < 0, 0, t1-1.96*t1_se),
"95% u.CI" = ifelse(t1+1.96*t1_se > 1, 1, t1+1.96*t1_se)
),
"T2" = c(
"Estimate" = t2,
"SE" = t2_se,
"95% l.CI" = ifelse(t2-1.96*t2_se < 0, 0, t2-1.96*t2_se),
"95% u.CI" = ifelse(t2+1.96*t2_se > 1, 1, t2+1.96*t2_se)
)
)
}
# chi squared test
chi.sq <- list(
"X^2" = sum(with(N_flatt,(frequency - expected_frequency)^2/expected_frequency)),
"df" = NK^NR - NK,
"p.value" = stats::pchisq(sum(with(N_flatt,(frequency - expected_frequency)^2/expected_frequency)),NK^NR - NK, lower.tail = F)
)
chi.sqB <- list(
"X^2" = sum(with(N_flattB,(frequency - expected_frequency)^2/expected_frequency)),
"df" = factorial(NK+NR-1)/(factorial(NR)*factorial(NK-1)) - NK,
"p.value" = stats::pchisq(sum(with(N_flattB,(frequency - expected_frequency)^2/expected_frequency)),
factorial(NK+NR-1)/(factorial(NR)*factorial(NK-1)) - NK, lower.tail = F)
)
colnames(N_flatt)[1:NR] <- sapply(1:NR,function(x)paste(c("R",x),collapse = ""))
colnames(N_flatt)[(NR+1):ncol(N_flatt)] <- c("Ob Freq","Exp Freq","Std Res")
colnames(N_flattB)[1] <- "N Cat"
colnames(N_flattB)[2:ncol(N_flattB)] <- c("Ob Freq","Exp Freq","Std Res")
ret <- list(
"call" = call,
"data_description" = data_description,
"allowed_classifications" = Q,
"exact_prob" = exact_prob,
"model_A" = list(
"table" = N_flatt,
"test_of_model_fit" = chi.sq
),
"model_B" = list(
"table" = N_flattB,
"probabilities_estimates" = probabilities_estimates,
"test_of_model_fit" = chi.sqB,
"test_of_symmetry" = test_of_symmetry
)
)
class(ret) <- "boschek.test"
return(ret)
}
# methods
print.boschek.test <- function(x){
print_model_fit_A <- paste(c(
'X^2(',x$model_A$test_of_model_fit$df,') = ',sprintf("%.2f", x$model_A$test_of_model_fit$`X^2`),', ',
ifelse(x$model_A$test_of_model_fit$p.value < .001, 'p < .001', paste(c('p = ', sprintf("%.3f", x$model_A$test_of_model_fit$p.value)),collapse = "")),
'\n'),sep = "",collapse = "")
print_test_of_symmetry <- paste(c(
'X^2(',x$model_B$test_of_symmetry$df,') = ',sprintf("%.2f", x$model_B$test_of_symmetry$`X^2`),', ',
ifelse(x$model_B$test_of_symmetry$p.value < .001, 'p < .001', paste(c('p = ', sprintf("%.3f", x$model_B$test_of_symmetry$p.value)),collapse = "")),
'\n'),sep = "",collapse = "")
print_model_fit_B <- paste(c(
'X^2(',x$model_B$test_of_model_fit$df,') = ',sprintf("%.2f", x$model_B$test_of_model_fit$`X^2`),', ',
ifelse(x$model_B$test_of_model_fit$p.value < .001, 'p < .001', paste(c('p = ', sprintf("%.3f", x$model_B$test_of_model_fit$p.value)),collapse = "")),
'\n'),sep = "",collapse = "")
cat("Call:\n")
print(x$call)
cat("\nModel fit:\n")
cat(print_model_fit_A)
cat("\nModel fit assuming symmetry:\n")
cat(print_model_fit_B)
if(x$data_description$number_of_categories == 2 & x$exact_prob){
cat("\nEstimated probabilities:\n")
print(x$model_B$probabilities_estimates[,1])
}
cat("\nTest of symmetry:\n")
cat(print_test_of_symmetry)
}
summary.boschek.test <- function(x){
print_model_fit_A <- paste(c(
'X^2(',x$model_A$test_of_model_fit$df,') = ',sprintf("%.2f", x$model_A$test_of_model_fit$`X^2`),', ',
ifelse(x$model_A$test_of_model_fit$p.value < .001, 'p < .001', paste(c('p = ', sprintf("%.3f", x$model_A$test_of_model_fit$p.value)),collapse = "")),
'\n'),sep = "",collapse = "")
print_test_of_symmetry <- paste(c(
'X^2(',x$model_B$test_of_symmetry$df,') = ',sprintf("%.2f", x$model_B$test_of_symmetry$`X^2`),', ',
ifelse(x$model_B$test_of_symmetry$p.value < .001, 'p < .001', paste(c('p = ', sprintf("%.3f", x$model_B$test_of_symmetry$p.value)),collapse = "")),
'\n'),sep = "",collapse = "")
print_model_fit_B <- paste(c(
'X^2(',x$model_B$test_of_model_fit$df,') = ',sprintf("%.2f", x$model_B$test_of_model_fit$`X^2`),', ',
ifelse(x$model_B$test_of_model_fit$p.value < .001, 'p < .001', paste(c('p = ', sprintf("%.3f", x$model_B$test_of_model_fit$p.value)),collapse = "")),
'\n'),sep = "",collapse = "")
cat("Call:\n")
print(x$call)
cat(paste(c(
'\n','Data description:','\n',
'Observations: ',x$data_description$observations,'\n',
'Categories: ',x$data_description$number_of_categories,'\n',
'Raters: ',x$data_description$number_of_raters,'\n'
), collapse = "")
)
cat("\nModel fit:\n")
cat(print_model_fit_A)
cat("\nModel fit assuming symetry:\n")
cat(print_model_fit_B)
if(x$data_description$number_of_categories == 2 & x$exact_prob){
cat("\nEstimated probabilities:\n")
stats::printCoefmat(x$model_B$probabilities_estimates, digits = 2,
cs.ind = c(1:4), tst.ind = integer(), zap.ind = integer())
}
cat("\nTest of symmetry:\n")
cat(print_test_of_symmetry)
}
residuals.boschek.test <- function(x, print.all = F){
table_A <- x$model_A$table
print_data_summary_add_A <- ""
if(!print.all){
print_data_summary_rows <- 1:ifelse(nrow(table_A)<10,nrow(table_A),10)
if(x$data_description$number_of_raters > 7){
print_data_summary_add1 <- paste(c(' and ',ncol(table_A)-10,' columns'),collapse = "")
table_A <- table_A[,c(1:8,(ncol(table_A)-2):ncol(table_A))]
table_A[,8] <- "..."
colnames(table_A)[8] <- "..."
}else{
print_data_summary_add1 <- NULL
}
if(nrow(table_A) - 10 > 0){
print_data_summary_add_A <- paste(c('... (shortened by ',nrow(table_A) - 10,' rows',print_data_summary_add1,')','\n'),collapse = "")
}
table_A <- table_A[print_data_summary_rows,]
}
table_B <- x$model_B$table
print_data_summary_add_B <- ""
if(!print.all){
print_data_summary_rows <- 1:ifelse(nrow(table_B)<10,nrow(table_B),10)
if(nrow(table_B) - 10 > 0){
print_data_summary_add_B <- paste(c('... (shortened by ',nrow(table_B) - 10,' rows)','\n'),collapse = "")
}
table_B <- table_B[print_data_summary_rows,]
}
cat("Model data matrix:\n")
print(table_A)
cat(print_data_summary_add_A)
cat("\n")
cat("Model data matrix assuming symmetry:\n")
print(table_B)
cat(print_data_summary_add_B)
}
# data documentation
#' d1
#'
#' test dataset for ilustration
#'
#' @format A data frame with 53940 rows and 10 variables:
#' \describe{
#' \item{V1}{ratings by rater 1}
#' \item{V2}{ratings by rater 2}
#' \item{V3}{ratings by rater 3}
#' \item{V4}{ratings by rater 4}
#' }
"d1"
# data documentation
#' d2
#'
#' test dataset for ilustration
#'
#' @format A data frame with 53940 rows and 10 variables:
#' \describe{
#' \item{V1}{rating by rater 1}
#' \item{V2}{rating by rater 2}
#' \item{V3}{rating by rater 3}
#' \item{frequency}{number of rating with the particular combination of ratings}
#' }
"d2"
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.