## ----include = FALSE----------------------------------------------------------
is_cran_check <- !isTRUE(as.logical(Sys.getenv("NOT_CRAN", "false")))
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = !is_cran_check
)

## ----load-data----------------------------------------------------------------
# library(mfrmr)
# 
# list_mfrmr_data()
# 
# data("ej2021_study1", package = "mfrmr")
# head(ej2021_study1)
# 
# study1_alt <- load_mfrmr_data("study1")
# identical(names(ej2021_study1), names(study1_alt))

## ----toy-setup----------------------------------------------------------------
# data("mfrmr_example_core", package = "mfrmr")
# toy <- mfrmr_example_core
# 
# fit_toy <- fit_mfrm(
#   data = toy,
#   person = "Person",
#   facets = c("Rater", "Criterion"),
#   score = "Score",
#   method = "JML",
#   model = "RSM",
#   maxit = 30
# )
# diag_toy <- diagnose_mfrm(fit_toy, residual_pca = "none")
# 
# summary(fit_toy)$overview
# summary(diag_toy)$overview
# names(plot(fit_toy, draw = FALSE))

## ----diagnostics-reporting----------------------------------------------------
# t4_toy <- unexpected_response_table(
#   fit_toy,
#   diagnostics = diag_toy,
#   abs_z_min = 1.5,
#   prob_max = 0.4,
#   top_n = 10
# )
# t12_toy <- fair_average_table(fit_toy, diagnostics = diag_toy)
# t13_toy <- bias_interaction_report(
#   estimate_bias(fit_toy, diag_toy,
#                 facet_a = "Rater", facet_b = "Criterion",
#                 max_iter = 2),
#   top_n = 10
# )
# 
# class(summary(t4_toy))
# class(summary(t12_toy))
# class(summary(t13_toy))
# 
# names(plot(t4_toy, draw = FALSE))
# names(plot(t12_toy, draw = FALSE))
# names(plot(t13_toy, draw = FALSE))
# 
# chk_toy <- reporting_checklist(fit_toy, diagnostics = diag_toy)
# subset(
#   chk_toy$checklist,
#   Section == "Visual Displays",
#   c("Item", "DraftReady", "NextAction")
# )

## ----fit-full-----------------------------------------------------------------
# fit <- fit_mfrm(
#   data = ej2021_study1,
#   person = "Person",
#   facets = c("Rater", "Criterion"),
#   score = "Score",
#   method = "MML",
#   model = "RSM",
#   quad_points = 7
# )
# 
# diag <- diagnose_mfrm(
#   fit,
#   residual_pca = "none"
# )
# 
# summary(fit)
# summary(diag)

## ----fit-full-pca-------------------------------------------------------------
# diag_pca <- diagnose_mfrm(
#   fit,
#   residual_pca = "both",
#   pca_max_factors = 6
# )
# 
# summary(diag_pca)

## ----strict-rsm-pcm-----------------------------------------------------------
# fit_rsm_strict <- fit_mfrm(
#   data = toy,
#   person = "Person",
#   facets = c("Rater", "Criterion"),
#   score = "Score",
#   method = "MML",
#   model = "RSM",
#   quad_points = 7,
#   maxit = 30
# )
# diag_rsm_strict <- diagnose_mfrm(
#   fit_rsm_strict,
#   diagnostic_mode = "both",
#   residual_pca = "none"
# )
# 
# fit_pcm_strict <- fit_mfrm(
#   data = toy,
#   person = "Person",
#   facets = c("Rater", "Criterion"),
#   score = "Score",
#   method = "MML",
#   model = "PCM",
#   step_facet = "Criterion",
#   quad_points = 7,
#   maxit = 30
# )
# diag_pcm_strict <- diagnose_mfrm(
#   fit_pcm_strict,
#   diagnostic_mode = "both",
#   residual_pca = "none"
# )
# 
# summary(diag_rsm_strict)$diagnostic_basis[, c("DiagnosticPath", "Status", "Basis")]
# summary(diag_pcm_strict)$diagnostic_basis[, c("DiagnosticPath", "Status", "Basis")]

## ----strict-screening---------------------------------------------------------
# screen_rsm <- evaluate_mfrm_diagnostic_screening(
#   design = list(person = 18, rater = 3, criterion = 3, assignment = 3),
#   reps = 1,
#   scenarios = c("well_specified", "local_dependence"),
#   model = "RSM",
#   maxit = 30,
#   quad_points = 7,
#   seed = 123
# )
# screen_pcm <- evaluate_mfrm_diagnostic_screening(
#   design = list(person = 18, rater = 3, criterion = 3, assignment = 3),
#   reps = 1,
#   scenarios = c("well_specified", "step_structure_misspecification"),
#   model = "PCM",
#   maxit = 30,
#   quad_points = 7,
#   seed = 123
# )
# 
# screen_rsm$performance_summary[, c("Scenario", "EvaluationUse", "LegacyAnyFlagRate", "StrictAnyFlagRate")]
# screen_pcm$performance_summary[, c("Scenario", "EvaluationUse", "LegacySensitivityProxy", "StrictSensitivityProxy", "DeltaStrictMinusLegacyFlagRate")]

## ----strict-reporting-route---------------------------------------------------
# chk_rsm_strict <- reporting_checklist(fit_rsm_strict, diagnostics = diag_rsm_strict)
# subset(
#   chk_rsm_strict$checklist,
#   Section == "Visual Displays" &
#     Item %in% c("QC / facet dashboard", "Strict marginal visuals", "Precision / information curves"),
#   c("Item", "Available", "DraftReady", "NextAction")
# )

## ----residual-pca-------------------------------------------------------------
# pca <- analyze_residual_pca(diag_pca, mode = "both")
# plot_residual_pca(pca, mode = "overall", plot_type = "scree")

## ----bias-apa-----------------------------------------------------------------
# data("mfrmr_example_bias", package = "mfrmr")
# bias_df <- mfrmr_example_bias
# fit_bias <- fit_mfrm(
#   bias_df,
#   person = "Person",
#   facets = c("Rater", "Criterion"),
#   score = "Score",
#   method = "MML",
#   model = "RSM",
#   quad_points = 7
# )
# diag_bias <- diagnose_mfrm(fit_bias, residual_pca = "none")
# bias <- estimate_bias(fit_bias, diag_bias, facet_a = "Rater", facet_b = "Criterion")
# fixed <- build_fixed_reports(bias)
# apa <- build_apa_outputs(fit_bias, diag_bias, bias_results = bias)
# 
# mfrm_threshold_profiles()
# vis <- build_visual_summaries(fit_bias, diag_bias, threshold_profile = "standard")
# vis$warning_map$residual_pca_overall

## ----reporting-api------------------------------------------------------------
# spec <- specifications_report(fit, title = "Study run")
# data_qc <- data_quality_report(
#   fit,
#   data = ej2021_study1,
#   person = "Person",
#   facets = c("Rater", "Criterion"),
#   score = "Score"
# )
# iter <- estimation_iteration_report(fit, max_iter = 8)
# subset_rep <- subset_connectivity_report(fit, diagnostics = diag)
# facet_stats <- facet_statistics_report(fit, diagnostics = diag)
# cat_structure <- category_structure_report(fit, diagnostics = diag)
# cat_curves <- category_curves_report(fit, theta_points = 101)
# bias_rep <- bias_interaction_report(bias, top_n = 20)
# plot_bias_interaction(bias_rep, plot = "scatter")

## ----design-prediction--------------------------------------------------------
# sim_spec <- build_mfrm_sim_spec(
#   n_person = 30,
#   n_rater = 4,
#   n_criterion = 4,
#   raters_per_person = 2,
#   assignment = "rotating"
# )
# 
# recovery <- suppressWarnings(
#   evaluate_mfrm_recovery(
#     sim_spec = sim_spec,
#     reps = 2,
#     maxit = 30,
#     seed = 2
#   )
# )
# 
# summary(recovery)$recovery_summary[, c("ParameterType", "Facet", "RMSE", "Bias")]
# plot(recovery, type = "summary", metric = "rmse", draw = FALSE)$data$plot_table
# 
# recovery_review <- assess_mfrm_recovery(
#   recovery,
#   min_reps = 2,
#   min_se_available = NULL,
#   max_mcse_rmse_ratio = NULL,
#   max_rmse = c(facet = 1, step = 1, default = 1.5),
#   max_abs_bias = c(default = 0.75)
# )
# 
# summary(recovery_review)$checklist[, c("Section", "Item", "Status")]
# status_plot <- plot(recovery_review, type = "status", draw = FALSE)
# status_plot$data$section_status
# status_plot$data$reading_order
# 
# metric_plot <- plot(recovery_review, type = "metrics", metric = "rmse", draw = FALSE)
# metric_plot$data$plot_table
# metric_plot$data$guidance
# 
# recovery_bundle <- build_summary_table_bundle(
#   recovery_review,
#   appendix_preset = "recommended"
# )
# recovery_bundle$table_index[, c("Table", "Rows", "Role")]
# 
# # For a release-scale check, use the optional validation protocol:
# # source(system.file("validation", "recovery-validation.R", package = "mfrmr"))
# # validation <- mfrmr_run_recovery_validation(tier = "core", output_dir = "recovery-validation")
# # summary(validation)
# # validation_summary <- mfrmr_summarize_recovery_validation(validation)
# # validation_summary$topline_release_decision
# # validation_summary$release_decision_table
# # validation_summary$domain_decision_table
# # The top-line table is the release-recovery conclusion. It uses recovery
# # metrics, convergence, and Monte Carlo precision as the primary evidence.
# # The domain table keeps uncertainty status separate, so OverallStatus =
# # "review" is not read as a recovery-metric failure by itself.
# 
# pred_pop <- predict_mfrm_population(
#   sim_spec = sim_spec,
#   reps = 2,
#   maxit = 30,
#   seed = 1
# )
# 
# summary(pred_pop)$forecast[, c("Facet", "MeanSeparation", "McseSeparation")]
# 
# keep_people <- unique(toy$Person)[1:18]
# toy_mml <- suppressWarnings(
#   fit_mfrm(
#     toy[toy$Person %in% keep_people, , drop = FALSE],
#     person = "Person",
#     facets = c("Rater", "Criterion"),
#     score = "Score",
#     method = "MML",
#     quad_points = 5,
#     maxit = 30
#   )
# )
# 
# new_units <- data.frame(
#   Person = c("NEW01", "NEW01"),
#   Rater = unique(toy$Rater)[1],
#   Criterion = unique(toy$Criterion)[1:2],
#   Score = c(2, 3)
# )
# 
# pred_units <- predict_mfrm_units(toy_mml, new_units, n_draws = 0)
# pv_units <- sample_mfrm_plausible_values(toy_mml, new_units, n_draws = 2, seed = 1)
# 
# summary(pred_units)$estimates[, c("Person", "Estimate", "Lower", "Upper")]
# summary(pv_units)$draw_summary[, c("Person", "Draws", "MeanValue")]

## ----recovery-export, eval = FALSE--------------------------------------------
# export_summary_appendix(
#   list(recovery = recovery, recovery_review = recovery_review),
#   output_dir = tempdir(),
#   prefix = "mfrmr_recovery_appendix",
#   preset = "recommended",
#   include_html = FALSE,
#   overwrite = TRUE
# )

