## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ## ----sensitivity-data--------------------------------------------------------- utils::data("data_complete_30", package = "dyadicMarkov") utils::data("data_complete_60", package = "dyadicMarkov") utils::data("data_complete_90", package = "dyadicMarkov") utils::data("data_complete_180", package = "dyadicMarkov") utils::data("data_complete_720", package = "dyadicMarkov") ## ----extract-chain------------------------------------------------------------ simulation_variables <- setNames( rev(sort(unique(as.character(data_complete_30$variable)))), c("V1", "V2") ) simulation_chain <- function(x, n, dyad, variable, member) { measurement_columns <- paste0("TM", seq_len(n)) variable_code <- unname(simulation_variables[[variable]]) keep <- x$dyad == as.character(dyad) & x$variable == variable_code & x$members == as.character(member) as.integer(unlist( x[keep, measurement_columns, drop = FALSE], use.names = FALSE )) } simulation_pattern <- function(x, dyad, variable, member) { variable_code <- unname(simulation_variables[[variable]]) x[ x$dyad == as.character(dyad) & x$variable == variable_code & x$members == as.character(member), "caseSimulated", drop = TRUE ] } ## ----univariate-chains-------------------------------------------------------- FM_30 <- simulation_chain(data_complete_30, 30L, 241L, "V2", 1L) SM_30 <- simulation_chain(data_complete_30, 30L, 241L, "V2", 2L) FM_90 <- simulation_chain(data_complete_90, 90L, 241L, "V2", 1L) SM_90 <- simulation_chain(data_complete_90, 90L, 241L, "V2", 2L) ## ----univariate-counts-------------------------------------------------------- emp_uni_30 <- dyadicMarkov::countEmp(FM_30, SM_30, states = 2L) emp_uni_90 <- dyadicMarkov::countEmp(FM_90, SM_90, states = 2L) rbind( `L = 30` = rowSums(emp_uni_30), `L = 90` = rowSums(emp_uni_90) ) ## ----univariate-patterns------------------------------------------------------ uni_30 <- dyadicMarkov::univariatePattern( chainFM = FM_30, chainSM = SM_30, states = 2L, alpha = 0.05 ) uni_90 <- dyadicMarkov::univariatePattern( chainFM = FM_90, chainSM = SM_90, states = 2L, alpha = 0.05 ) c( `L = 30` = uni_30$pattern, `L = 90` = uni_90$pattern ) c( `L = 30` = simulation_pattern( data_complete_30, 241L, "V2", 1L ), `L = 90` = simulation_pattern( data_complete_90, 241L, "V2", 1L ) ) ## ----bivariate-chains--------------------------------------------------------- FM_V1_30 <- simulation_chain(data_complete_30, 30L, 241L, "V1", 1L) SM_V1_30 <- simulation_chain(data_complete_30, 30L, 241L, "V1", 2L) FM_V2_30 <- FM_30 SM_V2_30 <- SM_30 FM_V1_90 <- simulation_chain(data_complete_90, 90L, 241L, "V1", 1L) SM_V1_90 <- simulation_chain(data_complete_90, 90L, 241L, "V1", 2L) FM_V2_90 <- FM_90 SM_V2_90 <- SM_90 emp_bi_30 <- dyadicMarkov::countEmpBivariate( chainFM_V1 = FM_V1_30, chainSM_V1 = SM_V1_30, chainFM_V2 = FM_V2_30, chainSM_V2 = SM_V2_30, states = 2L ) emp_bi_90 <- dyadicMarkov::countEmpBivariate( chainFM_V1 = FM_V1_90, chainSM_V1 = SM_V1_90, chainFM_V2 = FM_V2_90, chainSM_V2 = SM_V2_90, states = 2L ) ## ----bivariate-cases---------------------------------------------------------- case_bi_30 <- dyadicMarkov::bivariateCase(emp_bi_30, alpha = 0.05) case_bi_90 <- dyadicMarkov::bivariateCase(emp_bi_90, alpha = 0.05) c( `L = 30` = case_bi_30$case, `L = 90` = case_bi_90$case ) c( `L = 30` = simulation_pattern( data_complete_30, 241L, "V1", 1L ), `L = 90` = simulation_pattern( data_complete_90, 241L, "V1", 1L ) ) ## ----bivariate-local-pattern-------------------------------------------------- complete_bi_90 <- dyadicMarkov::completePattern(emp_bi_90) complete_bi_90$pattern ## ----full-sensitivity-reproduction, eval=FALSE-------------------------------- # # ------------------------------------------------------------------------- # # Helper 1: recover the A-code from a univariate result. # # # # univariatePattern() returns labels such as: # # "IM (A0)" # # "APM (A1)" # # "AM (A2)" # # "PM (A3)" # # ------------------------------------------------------------------------- # univariate_code <- function(x) { # code <- sub( # "^.*\\(([A-Z][0-9]+)\\)$", # "\\1", # x$pattern # ) # # if (identical(code, x$pattern)) { # stop("The univariate pattern code could not be extracted.") # } # # code # } # # # # ------------------------------------------------------------------------- # # Helper 2: recover the current pattern code selected by an AIC procedure. # # # # partialPattern() and completePattern() return: # # - $pattern: the selected pattern written out in full; # # - $aic: the candidate table, including the short matrix code. # # # # Matching these two fields gives codes such as B2, D4, or E4. # # ------------------------------------------------------------------------- # selected_matrix_code <- function(x) { # i <- match(x$pattern, x$aic$pattern) # # if (is.na(i)) { # stop("The selected pattern could not be matched to its AIC table.") # } # # x$aic$matrix[i] # } # # # # ------------------------------------------------------------------------- # # Helper 3: classify one selected member-variable combination. # # # # For each classification, the selected member is placed in the # # first-member position for the main variable, while the other member is # # placed in the second-member position. The second variable is supplied # # using the same member order. # # # # The global bivariate classification determines whether the result is: # # - trivial; # # - univariate; # # - partial bivariate; # # - complete bivariate. # # # # The corresponding local identification method is then applied. # # ------------------------------------------------------------------------- # classify_local_pattern <- function(data, n, dyad, variable, member, # alpha = 0.05) { # # other_member <- if (member == 1L) 2L else 1L # second_variable <- if (variable == "V1") "V2" else "V1" # # # Main variable # main_first <- simulation_chain( # data, n, dyad, variable, member # ) # # main_second <- simulation_chain( # data, n, dyad, variable, other_member # ) # # # Second variable, using the same member order # second_first <- simulation_chain( # data, n, dyad, second_variable, member # ) # # second_second <- simulation_chain( # data, n, dyad, second_variable, other_member # ) # # # Empirical 16 x 2 transition-count matrix # empirical <- dyadicMarkov::countEmpBivariate( # chainFM_V1 = main_first, # chainSM_V1 = main_second, # chainFM_V2 = second_first, # chainSM_V2 = second_second, # states = 2L # ) # # # Global bivariate case # global <- dyadicMarkov::bivariateCase( # empirical, # alpha = alpha # ) # # # Proceed to the appropriate local identification # if (global$case == "trivial") { # return("trivial") # } # # if (global$case == "univariate") { # local <- dyadicMarkov::univariatePattern( # main_first, # main_second, # states = 2L, # alpha = alpha # ) # # return(univariate_code(local)) # } # # if (global$case == "partial") { # local <- dyadicMarkov::partialPattern(empirical) # return(selected_matrix_code(local)) # } # # if (global$case == "complete") { # local <- dyadicMarkov::completePattern(empirical) # return(selected_matrix_code(local)) # } # # stop("Unknown bivariate case: ", global$case) # } # # # # ------------------------------------------------------------------------- # # Helper 4: run the classification for all 1,000 dyads at one sequence # # length. # # # # Each dyad contributes four member-variable combinations: # # # # V1, first member # # V1, second member # # V2, first member # # V2, second member # # # # The simulation label is read directly from caseSimulated rather than # # entered manually. # # ------------------------------------------------------------------------- # classify_one_length <- function(data, n, alpha = 0.05) { # # dyads <- sort(unique(as.integer(data$dyad))) # # orientations <- data.frame( # variable = c("V1", "V1", "V2", "V2"), # member = c(1L, 2L, 1L, 2L), # stringsAsFactors = FALSE # ) # # out <- vector( # "list", # length(dyads) * nrow(orientations) # ) # # k <- 1L # # for (i in seq_along(dyads)) { # # d <- dyads[i] # # # Optional progress information when the full analysis is run locally # if (i %% 100L == 0L) { # message( # "Length ", n, ": completed ", # i, " of ", length(dyads), " dyads" # ) # } # # for (j in seq_len(nrow(orientations))) { # # variable <- orientations$variable[j] # member <- orientations$member[j] # # # Pattern used to simulate this member-variable combination # variable_code <- unname(simulation_variables[[variable]]) # # keep <- data$dyad == as.character(d) & # data$variable == variable_code & # data$members == as.character(member) # # simulated <- unique(data$caseSimulated[keep]) # # if (length(simulated) != 1L) { # stop( # "Expected exactly one simulation label for dyad ", # d, ", variable ", variable, # ", member ", member, "." # ) # } # # # Pattern identified by the current dyadicMarkov workflow # identified <- classify_local_pattern( # data = data, # n = n, # dyad = d, # variable = variable, # member = member, # alpha = alpha # ) # # out[[k]] <- data.frame( # length = n, # dyad = d, # variable = variable, # member = member, # simulated = simulated, # identified = identified, # stringsAsFactors = FALSE # ) # # k <- k + 1L # } # } # # do.call(rbind, out) # } # # # # ------------------------------------------------------------------------- # # Run all five sequence lengths. # # # # 5 lengths x 1,000 dyads x 4 member-variable combinations # # = 20,000 classifications. # # ------------------------------------------------------------------------- # simulation_data <- list( # `30` = data_complete_30, # `60` = data_complete_60, # `90` = data_complete_90, # `180` = data_complete_180, # `720` = data_complete_720 # ) # # classification_results <- do.call( # rbind, # lapply(names(simulation_data), function(n) { # classify_one_length( # data = simulation_data[[n]], # n = as.integer(n), # alpha = 0.05 # ) # }) # ) # # # ------------------------------------------------------------------------- # # Calculate sensitivity and specificity for each pattern and sequence length. # # # # Sensitivity = TP / (TP + FN) # # Specificity = TN / (TN + FP) # # # # The simulated pattern is treated as the reference classification. # # ------------------------------------------------------------------------- # patterns <- c("A1", "B2", "D4", "E4") # lengths <- c(30L, 60L, 90L, 180L, 720L) # # sensitivity_results <- data.frame( # Length = lengths # ) # # for (pattern in patterns) { # # sensitivity <- numeric(length(lengths)) # specificity <- numeric(length(lengths)) # # for (i in seq_along(lengths)) { # # current <- classification_results[ # classification_results$length == lengths[i], # ] # # truth <- current$simulated == pattern # identified <- current$identified == pattern # # TP <- sum(truth & identified) # FN <- sum(truth & !identified) # TN <- sum(!truth & !identified) # FP <- sum(!truth & identified) # # sensitivity[i] <- TP / (TP + FN) # specificity[i] <- TN / (TN + FP) # } # # sensitivity_results[[paste(pattern, "Se")]] <- # round(sensitivity, 3) # # sensitivity_results[[paste(pattern, "Sp")]] <- # round(specificity, 3) # } # # sensitivity_results ## ----sensitivity-results------------------------------------------------------ # Results obtained from the full sensitivity-analysis reproduction above. # Values are stored here so that the vignette does not rerun the 20,000 # classifications every time it is built. sensitivity_results <- data.frame( Length = c(30L, 60L, 90L, 180L, 720L), `A1 Se` = c(0.125, 0.836, 0.936, 0.941, 0.946), `A1 Sp` = c(0.992, 0.942, 0.951, 0.997, 1.000), `B2 Se` = c(0.114, 0.680, 0.801, 0.835, 0.850), `B2 Sp` = c(0.994, 0.940, 0.935, 0.991, 1.000), `D4 Se` = c(0.000, 0.119, 0.426, 0.876, 0.951), `D4 Sp` = c(0.999, 0.972, 0.955, 0.952, 0.960), `E4 Se` = c(0.015, 0.358, 0.621, 0.765, 0.808), `E4 Sp` = c(1.000, 0.997, 0.997, 0.999, 0.998), check.names = FALSE ) knitr::kable( sensitivity_results, digits = 3, align = c("r", rep("r", 8)), caption = paste( "Sensitivity (Se) and specificity (Sp) for patterns A1, B2, D4,", "and E4 across the five sequence lengths." ) )