Exploring person, rater, and task attributes

An assessment team wants to describe its raters’ backgrounds before planning training. Do groups summarize different combinations of experience, workload, and training? Would those descriptions change with the number of groups, feature weights, or plausible replacements for unrecorded attributes?

This example follows that question from a feature table to a comparison of partitions. It uses 120 fictional raters, not observations from an assessment program. No true groups are planted. Successful execution does not establish that these features, imputation models, or groups are appropriate for your raters. The groups describe backgrounds; they do not measure rater quality, severity, or ability, and do not identify which training will be effective. Later sections add separate person and task attributes, link all three classifications to a sparse assignment roster, and explain the roles of feature selection and PCA.

These APIs are available from version 0.2.4. Install the optional cluster and mice packages to run the entire example.

Choose attributes and their meaning

Use one row per rater, rather than repeating attributes for every rating. Keep identifiers separate from features. Select variables for the question before examining which choices produce attractive groups.

Column	Meaning and coding	Role here
`ExperienceYears`	Completed years of rating experience; numeric	Feature
`AnnualRatings`	Number of responses rated in the previous year; numeric	Feature
`Specialty`	Main subject area; unordered factor	Feature
`TrainingLevel`	Highest completed training stage; ordered factor	Feature
`Certified`	Holds the relevant certification; logical	Feature
`WorkshopHours`	Hours of workshops attended in the previous year; numeric	Auxiliary imputation predictor
`MentoringYears`	Years in a formal mentoring role; undefined without that role	Reviewed, then excluded from the all-rater feature set

library(mfrmr)
set.seed(20260921)
n <- 120L
raters <- data.frame(
  Rater = sprintf("R%03d", seq_len(n)),
  ExperienceYears = sample(0:30, n, replace = TRUE),
  AnnualRatings = round(exp(rnorm(n, log(250), 0.6))),
  Specialty = factor(sample(c("Language", "Science", "Social science"),
                            n, replace = TRUE)),
  TrainingLevel = ordered(sample(c("Introductory", "Intermediate", "Advanced"),
                                 n, replace = TRUE),
                          levels = c("Introductory", "Intermediate", "Advanced")),
  Certified = sample(c(FALSE, TRUE), n, replace = TRUE),
  WorkshopHours = sample(0:80, n, replace = TRUE)
)
has_mentoring_role <- sample(c(FALSE, TRUE), n, replace = TRUE, prob = c(.7, .3))
raters$MentoringYears <- ifelse(has_mentoring_role,
  pmin(raters$ExperienceYears, sample(0:10, n, replace = TRUE)), NA_real_)

# Introduce unrecorded values only after generating the fictional attributes.
feature_names <- c("ExperienceYears", "AnnualRatings", "Specialty",
                   "TrainingLevel", "Certified")
missing_counts <- c(12L, 18L, 6L, 10L, 10L)
for (j in seq_along(feature_names)) {
  raters[[feature_names[j]]][sample.int(n, missing_counts[j])] <- NA
}
head(raters)

The Gower distance calculation accommodates these mixed types. Numeric features use their observed ranges. Ordered levels use ordered integer codes, so the declared order and the spacing implicit in that coding matter. Logical features treat both matches (FALSE/FALSE and TRUE/TRUE) alike. Certification is therefore not coded as a special presence-only attribute. Equal feature weights are a starting choice; correlated or redundant attributes can still receive too much combined influence. A highly unusual numeric value can also change its feature’s range and hence other raters’ distances.

Distinguish unrecorded and undefined values

Here a missing mentoring duration means that the attribute does not apply. It is not a duration waiting to be estimated. A duration of zero would instead mean an applicable mentoring role with no completed years.

all_features <- c(feature_names, "MentoringYears")
initial <- mfrm_features(raters, "Rater", all_features)
reasons <- initial$missing
reasons$Reason <- ifelse(reasons$Feature == "MentoringYears",
                        "Not applicable: no mentoring role", "Not recorded")
review <- mfrm_features(raters, "Rater", all_features, reasons)
summary(review)
table(review$missing$Reason)

The target is a description of all raters. Requiring a mentoring duration would exclude nonmentors. We therefore leave it out of the selected features and the imputation predictors, retaining the original review. A separate analysis confined to mentors would address a different population.

selected_reasons <- subset(review$missing, Feature %in% feature_names)
features <- mfrm_features(raters, "Rater", feature_names, selected_reasons)
features$feature_summary
table(features$row_summary$Complete)

mfrm_cluster_pam(features, k = 3) would stop on the remaining missing values. Explicit missing = "omit" is available, but selects complete raters and changes the population being described. It does not correct missing-data bias.

Fit an explicit imputation model

Only cells reviewed as unrecorded are eligible here. The fictional omissions were introduced at random; in an actual program, an explanation based on the data collection process is needed. Missingness reasons alone do not identify a statistical missingness mechanism. Unobserved-value-dependent nonresponse requires additional sensitivity assumptions.

The mice documentation describes the methods and predictor controls used below. WorkshopHours can help predict unrecorded attributes without becoming a clustering feature. The identifier is neither imputed nor used to predict other variables.

eligible <- subset(features$missing, Reason == "Not recorded")
imputation_data <- raters[c("Rater", feature_names, "WorkshopHours")]
where <- matrix(FALSE, nrow(imputation_data), ncol(imputation_data),
                dimnames = dimnames(imputation_data))
where[cbind(match(eligible$ID, imputation_data$Rater),
            match(eligible$Feature, names(imputation_data)))] <- TRUE
methods <- mice::make.method(imputation_data, where = where)
methods[feature_names] <- c("pmm", "pmm", "polyreg", "polr", "logreg")
methods["Rater"] <- ""
predictors <- mice::make.predictorMatrix(imputation_data)
predictors[, "Rater"] <- 0
predictors["Rater", ] <- 0
model <- mice::mice(imputation_data, m = 5, maxit = 5,
  method = methods, predictorMatrix = predictors, where = where,
  seed = 724, printFlag = FALSE)
model$loggedEvents
plot(model, c("ExperienceYears", "AnnualRatings"), layout = c(2, 2))

Five imputations and five iterations keep this executable example short; they are not an adequacy recommendation. The plot shows the mean and standard deviation of imputed values along each chain for two numeric features. Persistent drift or differences between chains call for investigation; overlapping lines in this short example do not establish convergence. Inspect logged events, chain behavior, observed versus imputed distributions, and substantive plausibility. A lack of logged events does not establish convergence or a suitable model. mfrm_cluster_imputed() checks consistency of the data, not these assumptions.

Describe groups and compare choices

Reuse the same fitted imputation model for every setting. Changing the imputations as well as the weights would change two things at once.

experience_weights <- setNames(rep(1, length(feature_names)), feature_names)
experience_weights["ExperienceYears"] <- 3
analyses <- list(
  TwoGroups = mfrm_cluster_imputed(features, model, eligible, k = 2),
  ThreeGroups = mfrm_cluster_imputed(features, model, eligible, k = 3),
  ExperienceWeighted = mfrm_cluster_imputed(features, model, eligible, k = 3,
                                          weights = experience_weights)
)
comparison <- mfrm_cluster_compare(analyses)
comparison$analysis_summary
summary(comparison)
weight_effect <- subset(comparison$comparison_summary,
  First == "ThreeGroups" & Second == "ExperienceWeighted")

MeanChangedFraction describes the fraction of rater pairs whose together/apart status changes between settings, averaged over the supplied imputations. The minimum and maximum show its range. The individual rows in comparison$comparisons retain every imputation, including split/join counts and adjusted Rand indices. An adjusted Rand index of one indicates identical partitions, allowing arbitrary group numbering. Neither high agreement nor a large silhouette proves that the groups are meaningful. Silhouettes under different weights use different distances and should not automatically select the weights.

In this fictional example, keeping three groups and raising the experience weight from one to three changes the together/apart status of the fraction shown in the table of rater pairs on average across imputations. Thus, the selected weights materially affect this particular description; presenting a single partition would conceal that dependence. This percentage concerns pairs, not the percentage of raters assigned to a different numbered group. It is not evidence about an actual assessment program.

Also inspect MinGroupSize. A rare category can produce a very small group, and repeated profiles can give high separation values. A high overall silhouette can coexist with these conditions.

Compare feature selections

Does including annual rating workload change the background groups? Hold the raters, group count, method, and imputation model fixed, and remove AnnualRatings from the clustering features. Select the corresponding eligible cells rather than passing cells for an unselected feature.

background_names <- setdiff(feature_names, "AnnualRatings")
background_features <- mfrm_features(raters, "Rater", background_names,
  subset(selected_reasons, Feature %in% background_names))
background_only <- mfrm_cluster_imputed(background_features, model,
  subset(eligible, Feature %in% background_names), k = 3)
feature_comparison <- mfrm_cluster_compare(list(
  WithWorkload = analyses$ThreeGroups, WithoutWorkload = background_only))
feature_comparison$analysis_summary
feature_comparison$weights
summary(feature_comparison)

Here, removing workload changes the together/apart status of the fraction shown in the table of rater pairs on average across the five imputations. This measures the description’s dependence on including workload, not whether workload should be excluded. The Features column records each analysis’s feature count; weights lists only its selected features. Silhouettes now use different distances and cannot automatically choose a feature set.

Workload remains in the imputation model, so this does not examine removing it as a predictor of missing attributes. Different feature selections must reuse the same fitted mids object, and stored completed values are checked against the corresponding imputation. A selection containing only complete attributes can use its empty missing table as impute; it retains a partition for every completion for comparison purposes.

The included entities must match across all results. With missing = "omit", removing an incomplete feature may admit additional raters, in which case the comparison stops. Define and report a common cohort before fitting if that is the intended question; the API does not silently intersect samples.

Interpret the retained profiles

Inspect the actual profiles before attaching interpretations to the groups:

one_completion <- analyses$ThreeGroups$analyses[[1]]
one_completion$profiles$numeric
subset(one_completion$profiles$categorical, Feature == "TrainingLevel")
head(one_completion$membership)
analyses$ThreeGroups$co_membership[1:5, 1:5]

These profiles describe one completion, not an average cluster or a final assignment. Group numbers can change across completions. The co_membership matrix instead records how often each pair is grouped together across all imputations. For example, 0.8 means four of the five completions, conditional on this model and these settings. It is not an 80% probability of a true shared class. Numeric ranges are recalculated within each completion; their changes also contribute to sensitivity.

The workflow answers whether proposed background descriptions depend strongly on the supplied choices. It supplies no automatic preferred group count, consensus assignment, resampling stability, Rubin-pooled inference, or imputation of rating responses. Before using groups to make decisions about raters, establish that the attributes and proposed use address the intended question and assess the consequences with appropriate evidence.

Read separation, feature profiles, and imputation sensitivity

The silhouette plot shows which raters fit their assigned group less clearly under the selected distances. Negative widths warrant examining the rater’s profile and nearby groups; high widths do not establish useful or stable types. The dashed line is the overall mean. Raters excluded from clustering have no width; their IDs remain in the returned plot data.

plot(one_completion)

Inspect one feature at a time so that its units and categories remain clear. Numeric profiles show means and medians in original units, without confidence intervals. Categorical profiles show proportions within each group, preserving the order of training levels. These are descriptions of this completion only.

plot(one_completion, type = "profile", feature = "ExperienceYears")
plot(one_completion, type = "profile", feature = "TrainingLevel")

The development version also offers dedicated ggplot conversion for these saved summaries. For example, remove the heading while keeping the original units, group sizes and interpretation note:

if (requireNamespace("ggplot2", quietly = TRUE)) {
  profile_figure <- as_ggplot(one_completion, type = "profile",
    feature = "ExperienceYears", preset = "monochrome") +
    ggplot2::labs(title = NULL)
  print(profile_figure)
}

Use as_ggplot(one_completion) for the silhouette view, retaining negative widths and the saved overall-mean line. Use the profile call with feature = "TrainingLevel" for a categorical heatmap, including unused levels in their original order. The default and component = "table" keep the full view; categorical profiles also accept component = "matrix". Means and medians are descriptive, not effect estimates or interval bounds. The slight vertical offsets improve visibility without changing the values. Saved summaries and excluded IDs remain available through plot_data().

The next heatmap asks which pairs remain together across imputations. For readability, this view explicitly selects the first 20 IDs in the input table; it does not recluster them, change the denominator, or represent the remaining pairs. Omit ids to display all raters. By default, ID labels are hidden above 50 displayed entities; the data are never sampled automatically.

plot(analyses$ThreeGroups, ids = raters$Rater[1:20])

The colour scale always runs from zero to one. Grey denotes unavailable pairs involving excluded raters, not pairs that never share a group. Rows retain the requested ID order. The map does not infer a consensus partition or a hierarchy, and it does not measure stability under new samples of raters.

For a custom figure, extract the same values without drawing or refitting:

view <- plot(analyses$ThreeGroups, ids = raters$Rater[1:5], draw = FALSE)
plot_data(view)$matrix

The development version also converts this same five-ID view to ggplot:

if (requireNamespace("ggplot2", quietly = TRUE)) {
  co_membership_figure <- as_ggplot(view, component = "matrix") +
    ggplot2::labs(title = NULL)
  print(co_membership_figure)
}

The default and component = "matrix" retain the same complete heatmap. For a new selection, use as_ggplot(analyses$ThreeGroups, ids = ...) or create another saved plot with plot(..., draw = FALSE). A saved view keeps its preset and label choice even after session options change. Grey cells also carry crosses, so unavailable pairs remain distinct from zero in monochrome. plot_data(co_membership_figure) retains the matrix, selected IDs, excluded IDs and number of imputations. Hiding headings with ggplot2::labs() does not remove those records or change the meaning of the fractions.

Custom distribution plots

Profile means do not show individual variation. Retrieve one original-unit value per rater and feature, with group membership and inclusion status:

values <- plot_data(one_completion, "observations", type = "profile",
  feature = "ExperienceYears")
head(values)

Value is the selected completion’s experience in years. Included = FALSE retains omitted raters, even when this particular feature is observed. In this completed example, missing values have been imputed; retain the original features$missing record to distinguish them from observed attributes. A completion’s values are not repeated independent observations, and its arbitrary group labels must not be averaged across imputations.

The following raincloud combines a half-density, a box and individual points. All plotting choices remain editable in ordinary ggplot2. The density is a visual smoothing of a discrete years variable; for small groups or ordinal rubric scores prefer points/counts. Group differences are partly built into the clustering and are not independent tests validating the groups.

if (requireNamespace("ggplot2", quietly = TRUE)) {
  displayed <- subset(values, Included & is.finite(Value))
  displayed$Group <- factor(displayed$Cluster)
  displayed$Position <- as.integer(displayed$Group)
  clouds <- lapply(split(displayed, displayed$Group), function(d) {
    if (nrow(d) < 5L || length(unique(d$Value)) < 2L) return(NULL)
    density <- stats::density(d$Value, from = min(d$Value), to = max(d$Value))
    center <- d$Position[1]
    data.frame(X = c(center, center + .38*density$y/max(density$y), center),
      Value = c(density$x[1], density$x, tail(density$x, 1)),
      Group = d$Group[1])
  })
  clouds <- do.call(rbind, clouds)
  raincloud <- ggplot2::ggplot(displayed,
    ggplot2::aes(x = Position, y = Value, group = Group))
  if (!is.null(clouds)) raincloud <- raincloud +
    ggplot2::geom_polygon(data = clouds,
      ggplot2::aes(x = X, y = Value, group = Group, fill = Group),
      inherit.aes = FALSE, alpha = .35, colour = NA)
  raincloud <- raincloud +
    ggplot2::geom_boxplot(width = .08, outlier.shape = NA) +
    ggplot2::geom_point(ggplot2::aes(x = Position - .15), alpha = .5, size = 1.2,
      position = ggplot2::position_jitter(width = .045, height = 0, seed = 1)) +
    ggplot2::scale_x_continuous(breaks = seq_along(levels(displayed$Group)),
      labels = paste("Group", levels(displayed$Group))) +
    ggplot2::labs(x = NULL, y = "Rating experience (years)") +
    ggplot2::theme_minimal() + ggplot2::theme(legend.position = "none")
  print(raincloud)
}

write.csv(values, "rater-experience-long.csv", row.names = FALSE)
saveRDS(one_completion, "feature-groups.rds")
ggplot2::ggsave("experience-raincloud.pdf", raincloud, width = 7, height = 4.5)

For a selected wide plot table, plot_data(..., format = "long", value_cols = ...) gathers explicitly chosen quantities; other columns remain identifiers or metadata. Matrices retain row/column indices and labels, including missing cells. The network/reporting lesson shows how to export node, edge and coordinate tables without adopting a fixed figure style.

Compare PAM with a hierarchy

A hierarchy asks which groups merge as dissimilarity increases. Fit it separately from PAM: mfrm_cluster_hierarchical() uses average linkage (the mean of all cross-group Gower dissimilarities) or complete linkage (their maximum). The same feature types, weights, and missingness rules apply. Choose the group count explicitly; a tree does not automatically identify it.

For multiply imputed features, reuse the same model and eligible cells. This holds the completed tables fixed while changing the grouping method.

average <- mfrm_cluster_imputed(features, model, eligible, k = 3,
  method = "hierarchical", linkage = "average")
complete <- mfrm_cluster_imputed(features, model, eligible, k = 3,
  method = "hierarchical", linkage = "complete")
method_comparison <- mfrm_cluster_compare(list(
  PAM = analyses$ThreeGroups, Average = average, Complete = complete))
method_comparison$analysis_summary
summary(method_comparison)

The comparison identifies methods and linkages, alongside group sizes and separation. As before, pair-change summaries describe sensitivity to the supplied choices; they do not select a preferred algorithm or establish true rater types. Substantively different profiles or very small groups deserve inspection even when the overall pair agreement is high.

plot(average$analyses[[1]])

This is the full tree for one completion, with boxes marking its three groups. The 120 raters are all included; ID labels are hidden by default above 50 entities. For a smaller cohort, labels appear automatically. Use plot_data(plot(average$analyses[[1]], draw = FALSE))$table to inspect IDs and memberships in leaf order. Excluded raters retain unavailable memberships but have no leaves. No entities are sampled to draw the tree.

In the development version, as_ggplot(average$analyses[[1]]) converts this same saved tree with its group boxes, leaf order and merge heights. The default and component = "tree" retain the full view; use plot_data() for the membership table. Conversion does not refit, recut or combine trees across completions. Add ggplot2::labs(title = NULL, subtitle = NULL) to omit headings, while retaining the caption about what heights mean.

Merge heights are linkage dissimilarities, not significance levels or branch support. Tied distances can produce alternative trees; the group cut uses merge order, so a unique horizontal cutting height may not exist. Inspect each completion’s profiles with type = "profile" or silhouettes with type = "silhouette". plot(average) displays pairwise co-membership across imputations; it does not combine the trees into a consensus hierarchy.

How this relates to MFRM rater-type research

Clustering has been used alongside MFRM, with both hierarchical and nonhierarchical methods. In writing assessment, Eckes (2012), Operational Rater Types in Writing Assessment, clustered MFRM rater-by-criterion bias estimates using hierarchical Ward clustering, also checking complete and average linkage solutions. In music performance assessment, Wesolowski (2019), Predicting Operational Rater-Type Classifications Using Rasch Measurement Theory and Random Forests, reported nonhierarchical k-means clustering of 29 differential rater functioning indices, followed by random-forest prediction of the resulting labels. The forest is supervised prediction, not the clustering step. Avoid using labels discovered from the full sample as if out-of-bag classification alone validated those groups. For residual profiles, use mfrm_residual_network() then mfrm_report(); its help explains its different features, common-person correlation distances, validation requirements and RF follow-up limitations.

Those analyses ask about patterns of rating behaviour. This tutorial instead groups external backgrounds such as experience and specialty. Its Gower distances accommodate the selected mix of variable types, and PAM produces a partition rather than a hierarchy. The separate hierarchical API uses average or complete linkage. Ward’s variance criterion requires appropriate Euclidean geometry; this Gower interface does not define a Euclidean conversion and does not offer Ward clustering (see Murtagh and Legendre, 2014). The cited studies do not validate the present attributes, weights, imputation model, or PAM groups, and do not establish that one method is preferable for all uses. Grouping estimated MFRM effects would additionally require consideration of their measurement uncertainty and comparability. The present external-feature API does not propagate that uncertainty.

Link separately grouped persons, raters, and tasks

An assessment team may also ask whether different kinds of persons and tasks are represented across its rater groups. Keep three feature tables, with one row per person, rater, or task, and a separate rating table containing their IDs. The tables can have different numbers and types of features. Clustering the joined rating rows would repeat attributes according to assignment frequency and answer a different question about rating events.

Table	Example features	Unit classified
Persons	Years of language study, preparation hours, instruction mode	One person
Raters	Experience, workload, specialty, training, certification	One rater
Tasks	Time limit, word limit, response format	One task

The following extends the fictional example with 240 persons and 12 tasks. It uses the existing rater partition from one completed table only to illustrate joining by ID. The other two tables have no missing attributes. The group counts are illustrative choices, not estimated numbers of types.

set.seed(7241)
persons <- data.frame(
  Person = sprintf("P%03d", 1:240),
  LanguageStudyYears = sample(0:15, 240, replace = TRUE),
  PreparationHours = round(exp(rnorm(240, log(20), 0.5))),
  InstructionMode = factor(sample(c("Classroom", "Online", "Hybrid"),
                                  240, replace = TRUE)))
tasks <- data.frame(
  Task = sprintf("T%02d", 1:12),
  TimeLimitMinutes = sample(c(10, 20, 30), 12, replace = TRUE),
  WordLimit = sample(c(100, 200, 400), 12, replace = TRUE),
  ResponseFormat = factor(rep(c("Summary", "Argument", "Explanation"), 4)))
person_groups <- mfrm_cluster_pam(mfrm_features(persons, "Person",
  c("LanguageStudyYears", "PreparationHours", "InstructionMode")), k = 3)
task_groups <- mfrm_cluster_pam(mfrm_features(tasks, "Task",
  c("TimeLimitMinutes", "WordLimit", "ResponseFormat")), k = 3)

Record the planned assignments before collecting scores. Here each person receives three distinct tasks and each performance is assigned two distinct raters. The full Person-by-Rater-by-Task cross-product is not the assignment roster. Random scores below are placeholders for demonstrating coverage; they provide no evidence about achievement, rater effects, or model fit.

assignments <- do.call(rbind, lapply(persons$Person, function(id) {
  do.call(rbind, lapply(sample(tasks$Task, 3), function(task) {
    data.frame(Person = id, Rater = sample(raters$Rater, 2), Task = task)
  }))
}))
ratings <- assignments
ratings$Score <- sample(1:4, nrow(ratings), replace = TRUE)
ratings$Score[sample(nrow(ratings), 72)] <- NA_integer_
rating_review <- describe_mfrm_data(ratings, person = "Person",
  facets = c("Rater", "Task"), score = "Score",
  rating_min = 1, rating_max = 4, keep_original = TRUE,
  include_agreement = FALSE, expected_design = assignments)
rating_review$structural_missingness$summary
rating_review$design_connectivity

There are 1,440 assigned cells out of 345,600 possible combinations, and 72 assigned scores are absent. Assignment coverage is therefore 95%, even though fewer than 0.5% of all combinations were assigned. Neither a zero score nor an imputed score is inserted into unassigned cells. With actual data, supply the original roster: reconstructing it from observed scores would conceal missing assigned ratings.

Start from the assignment roster so that even an entirely absent rating row remains in the denominator. This example has one score per Person-Rater-Task key; check uniqueness before joining. For repeated occasions, include the occasion in the key. Attach the three independently obtained group labels by ID with match(); row position in a feature table is not an ID. The check stops this illustration if an ID is unknown or an entity has no classification. In an analysis allowing omitted entities, report their unclassified ratings separately rather than silently dropping them.

rating_key <- c("Person", "Rater", "Task")
stopifnot(!anyDuplicated(assignments[rating_key]),
          !anyDuplicated(ratings[rating_key]))
labeled <- merge(assignments, ratings, by = rating_key,
                 all.x = TRUE, sort = FALSE)
labeled$PersonGroup <- person_groups$membership$Cluster[
  match(labeled$Person, person_groups$membership$ID)]
labeled$RaterGroup <- one_completion$membership$Cluster[
  match(labeled$Rater, one_completion$membership$ID)]
labeled$TaskGroup <- task_groups$membership$Cluster[
  match(labeled$Task, task_groups$membership$ID)]
group_columns <- c("PersonGroup", "RaterGroup", "TaskGroup")
stopifnot(!anyNA(labeled[group_columns]))
for (column in group_columns) labeled[[column]] <- factor(labeled[[column]], 1:3)
labeled$Assigned <- 1L
labeled$Observed <- as.integer(!is.na(labeled$Score))
planned_counts <- xtabs(Assigned ~ PersonGroup + RaterGroup + TaskGroup,
                        data = labeled)
observed_counts <- xtabs(Observed ~ PersonGroup + RaterGroup + TaskGroup,
                         data = labeled)
group_coverage <- as.data.frame(planned_counts, responseName = "Assigned")
group_coverage$Observed <- as.vector(observed_counts)
group_coverage$Coverage <- ifelse(group_coverage$Assigned > 0,
  group_coverage$Observed / group_coverage$Assigned, NA_real_)
group_coverage

This table answers which group combinations were assigned and how many of their scores were recorded. A combination with no assignments has unavailable coverage, not zero coverage. Counts reflect the assignment design; they do not demonstrate associations between traits or differences in rating quality. Ratings share persons, raters, and tasks, so an ordinary independent-row test of group differences would not account for this dependence. Raw score means also confound person composition, rater severity, and task difficulty. Such comparisons need an appropriate measurement model and sufficient overlap, rather than treating the group labels as established latent classes.

If attributes are imputed for several facets, each imputation model must respect its unit of analysis and relevant dependencies. The current API analyzes one entity table at a time; it does not fit a joint cross-facet imputation model or pool group-by-group rating effects. Pairing the first completion of independently fitted person and rater models does not establish a valid joint imputation. Nor can numeric group labels be matched across completions without addressing label changes and membership uncertainty.

Distinguish sparse assignment from missing attributes

Four situations require different decisions:

Situation	Example	What to retain or review
Unassigned rating	A rater was never scheduled to score a performance	Original assignment roster; no invented observation
Missing assigned rating	A scheduled score was not submitted	Assigned denominator, absence reason, observed overlap
Unrecorded external attribute	Training hours were not recorded	Eligible cells, predictors, imputation assumptions
Undefined external attribute	Mentoring years for someone with no mentoring role	Meaning of the attribute; not an ordinary imputation target

Sparseness describes how much of a design is observed and how those observations connect. MCAR means that missingness does not depend on the observed or missing values. MAR allows dependence on observed values, but not on missing values after conditioning on the observed information. MNAR allows additional dependence on the missing values. Planned allocation does not automatically imply MCAR: assignment may depend on recorded proficiency or on information absent from the analysis. A missingness rate or diagnostic test alone cannot establish MAR. See Rubin (1976), Inference and Missing Data for the conditions for ignoring a missingness process in inference.

Inspect both planned and observed connectivity, per-level support, and coverage within relevant subgroups. Connected Person-facet graphs do not by themselves establish identification of a complete model, especially with interactions, or adequate precision. In a simulation of sparse MFRM designs, Wind and Ge (2021) found that detecting differential rater functioning depended on the linking design, not merely on having incomplete ratings. Their findings do not give a universal acceptable missingness percentage. Imputing external attributes does not add observed links or resolve confounding between disconnected rating blocks. Grouping noisy estimated effects from sparse data may also reflect differences in precision, which this feature API does not model.

Many features, k-means, and PCA

More columns are not necessarily more useful information. Repeated measures of one attribute can dominate a distance through their combined weight; unrelated columns can obscure the differences of interest. Steinley and Brusco (2008) examined variable selection for clustering and found that non-informative variables and their distributions affected performance. Choose features and their substantive domains before choosing weights. For example, equal total weight for a domain containing ten features and one containing two requires different per-feature weights. This is a substantive choice, not an automatic correction for correlations. The comparison API supports both positive feature weights and different selections on the same entities, as illustrated above. Choosing features to maximize an attractive grouping still requires separate validation; the comparison does not automate feature selection.

PCA and k-means answer different questions. PCA represents variation in numeric features with orthogonal continuous components. K-means partitions entities to minimize within-group squared Euclidean distances to centroids. Ding and He (2004), Principal Component Analysis and Effective K-means Clustering describe their relationship through a continuous relaxation of discrete cluster membership. This does not make ordinary PCA and discrete clustering equivalent, or establish that the first two PCs recover useful groups.

For a fixed numeric input, retaining all nonzero PC scores without further rescaling preserves its Euclidean distances. Truncating components changes those distances and may remove low-variance separation as well as noise. Standardizing features before PCA, or whitening the scores afterwards, also changes the chosen geometry. An explained-variance threshold alone is not evidence that a clustering structure was preserved. A two-component plot is a projection, not a test for the existence of clusters.

Ordinary PCA and k-means should not receive arbitrary numeric codes for nominal categories such as specialty. Keep Gower/PAM or average/complete linkage for mixed features. For selected numeric attributes, mfrm_pca() and mfrm_cluster_kmeans() use R’s prcomp() and kmeans(). They center features and, by default, divide by sample standard deviations. Positive feature weights apply to squared Euclidean distances. The fitted transformation, original units, entity IDs and missingness accounting remain available. No automatic feature, component or group-count selection is performed.

Here the question is whether grouping experience, annual workload and workshop hours changes when those three attributes are represented by two components. The same already fitted imputation model supplies both analyses. Workshop hours remain a predictor in that model whether or not they are selected for clustering. These numeric selections answer a different question from the mixed-background groups above; specialty and training level are not numeric inputs to this PCA.

numeric_names <- c("ExperienceYears", "AnnualRatings", "WorkshopHours")
numeric_review <- mfrm_features(raters, "Rater", numeric_names)
numeric_direct <- mfrm_cluster_imputed(numeric_review, model,
  numeric_review$missing, k = 3, method = "kmeans", seed = 42)
numeric_reduced <- mfrm_cluster_imputed(numeric_review, model,
  numeric_review$missing, k = 3, method = "kmeans", components = 2, seed = 42)
numeric_comparison <- mfrm_cluster_compare(list(
  StandardizedFeatures = numeric_direct, TwoComponents = numeric_reduced))
numeric_comparison$analysis_summary
summary(numeric_comparison)

Read Space, Scaling and Components before comparing the partitions. ChangedFraction measures changed together/apart pairs, allowing arbitrary renumbering of clusters. Different component counts change distances, so silhouette widths do not provide a common-scale test that one analysis is better. The default 25 random starts and fixed seed make this execution reproducible, but do not prove a global optimum. Assess another seed separately if initialization sensitivity matters.

first_numeric <- numeric_reduced$analyses[[1]]
summary(first_numeric$pca)
first_numeric$pca$loadings
plot(first_numeric$pca)  # The full explained-variance table, including omitted PCs.
plot(first_numeric$pca, type = "scores", groups = first_numeric, labels = FALSE)
plot(first_numeric$pca, type = "loadings", components = 1)
plot(first_numeric, type = "profile", feature = "ExperienceYears")

The first-completion score plot annotates a two-dimensional projection with that completion’s group labels. Colour and point shape distinguish groups; preset = "monochrome" retains the shapes. Shapes repeat after six groups, so inspect plot_data() for IDs and encodings in crowded views. Filled/open scree points distinguish retained/omitted components. Loadings are coefficients in the transformed feature space, not correlations in original units. Their signs are arbitrary. Inspect the original-unit profiles to explain what distinguishes the groups. No view establishes ability, severity, rater quality or an effect of training.

Every completion retains its own centering, standard deviations and PCA basis. The comparison pairs the same completions, while allowing these transformations to change as imputed values change. It compares partitions; it does not average PCA scores, loadings or arbitrary group numbers, and it does not Rubin-pool inference. These functions impute no rating scores.

For an already complete numeric table, use mfrm_cluster_kmeans(review, k=3) directly. To inspect and retain a specific reduction, fit pca <- mfrm_pca(review, components=2) and then mfrm_cluster_kmeans(pca, k=3). Passing a PCA object uses its scores exactly, without further standardization, range scaling or whitening. Save that PCA and its grouping together. plot() and mfrm_cluster_compare() reuse saved results; changing features, component counts, weights or seeds requires the corresponding new analysis. mfrmr_output_guide("features") lists the dedicated output routes. Use summary() for partition comparisons, and plot() or plot_data() for fitted PCA and grouping views.

In the development version, as_ggplot() has dedicated conversions for PCA scree, scores and loadings. The scores view keeps the selected components, equal axis units and the saved group colours and shapes. Only included entities are drawn; excluded IDs and the feature transformation remain in plot_data(). This example shows the same first completion used above:

if (requireNamespace("ggplot2", quietly = TRUE)) {
  numeric_figure <- as_ggplot(first_numeric$pca, type = "scores",
    groups = first_numeric, labels = FALSE, preset = "monochrome") +
    ggplot2::labs(title = NULL, subtitle = NULL)
  print(numeric_figure)
  plot_data(numeric_figure)$components
}

Use type = "scree" for explained-variance percentages and type = "loadings" for coefficients in the transformed feature space. The default and explicit component = "table" keep the same complete view. For already saved plot payloads, as_ggplot() uses their saved preset and label choice; set those choices when making the payload. A change of heading does not change the analysis. Saved dendrograms also have dedicated conversion as described above. Imputation co-membership also has dedicated conversion as described above. Silhouettes and numeric/categorical profiles also have dedicated conversion. All these views retain their saved summaries; use plot_data() for further custom graphics.

Reopen an analysis without repeating imputation or clustering

Keep the comparison result to retain both analyses, every completion, its transformation and partition, the fitted imputation model, and the co-membership summaries. A saved plot view also retains the selected display settings. The same saving route applies to comparison, feature_comparison and method_comparison above. The example uses temporary files; choose permanent paths to continue in a later session.

feature_file <- tempfile(fileext = ".rds")
feature_view <- plot(first_numeric$pca, type = "scores", groups = first_numeric,
  labels = FALSE, preset = "monochrome", draw = FALSE)
saveRDS(list(comparison = numeric_comparison, view = feature_view), feature_file)
saved_features <- readRDS(feature_file)
summary(saved_features$comparison)
saved_features$comparison$analysis_summary  # Each setting and completion.
plot_data(saved_features$view)    # Inspect this first-completion view.
feature_summary_file <- tempfile(fileext = ".csv")
write.csv(summary(saved_features$comparison), feature_summary_file, row.names = FALSE)
if (requireNamespace("ggplot2", quietly = TRUE)) {
  saved_figure <- as_ggplot(saved_features$view)
}

Changing a session plotting option does not restyle a saved view. Create a new view from the saved analysis to change its preset; changing the analysis settings requires a new analysis. The CSV contains only the comparison summary; keep the RDS for its analyses, transformations, omitted IDs and assumptions. These feature results use their own summary(), plot() and plot_data() methods. They are separate from mfrm_results() and from the pooled response-score intervals described in vignette("mfrmr-response-imputation"). Saved feature files retain identifiers and original attributes, so review what they contain before sharing them.