03_factor_analysis.qmd

---
title: "Factor Analysis"
---

```{r library}
#| echo: false
#| message: false
#| warning: false

library(EQRIanalysis)
library(dplyr)
library(tidyr)
library(knitr)
```

```{r get-data}
#| echo: false
#| message: false

responses_df <- EQRIanalysis::get_responses_df()
```

## Overview

Factor analysis evaluates the underlying structure of the EQRI questionnaire to answer:

- **Are we asking the right questions?** (coverage, gaps, redundancy)
- **How many dimensions?** Do questions measure distinct aspects of engineering quality?
- **Question sensitivity:** Which questions most influence each indicator?

## Analysis Strategy

Per Little & Rubin (2019, *Statistical Analysis with Missing Data*), questionnaires with **structural missingness by design** require context-aware analysis.

The EQRI questionnaire includes:

1. **Core items**: Asked across most/all contexts (program types × milestones)
2. **Context-specific items**: Asked only for certain programs or milestones

### Approach

Based on the revised factor analysis strategy (see `dev/revised_sprint1_strategy.R`):

1. **Identify core items** - Questions present in ≥8 of 10 contexts (<20% missing)
2. **Assess factorability by context** - Evaluate each program × milestone combination
3. **Document response patterns** - Show which contexts have sufficient data for analysis

### Question Coverage

This heatmap shows which questions are asked in each context, helping identify core vs. context-specific items.

```{r plot-coverage}
#| echo: false
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Question coverage across program types and milestones"

plot_question_coverage(responses_df, core_threshold = 8)
```


## Question-to-Indicator Mapping

Understanding which questions contribute to each indicator is essential for interpreting indicator scores and conducting factor analysis.

### Visual Mapping

```{r plot-question-indicator}
#| echo: false
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Questions mapped to EQRI indicators"

plot_question_indicator_mapping(responses_df, sort_by = "indicator")
```


### Summary Table

```{r indicator-summary-table}
#| echo: false
#| tbl-cap: "Number of questions per indicator"

indicator_summary <- summarize_indicator_questions(responses_df)

kable(indicator_summary,
      col.names = c("Indicator", "Questions", "Responses", "Question Numbers"),
      format = "pipe",
      align = c("l", "r", "r", "l"))
```

**Interpretation:**

The seven EQRI indicators are calculated from specific sets of questions:

```{r display-indicator-composition}
#| echo: false
#| results: asis

for (i in 1:nrow(indicator_summary)) {
  cat(sprintf("\n- **%s**: %d questions (%s)\n", 
              indicator_summary$INDICATOR[i],
              indicator_summary$n_questions[i],
              indicator_summary$questions[i]))
}
```

This mapping is defined in the questionnaire design and determines how individual question responses aggregate into indicator scores.

---

## Questionnaire Reliability Assessment

Before conducting factor analysis, we must establish that the questionnaire produces **reliable** (internally consistent) measurements. Per Revelle & Zinbarg (2009), reliability is a prerequisite for validity and factor interpretability.

### Reliability Metrics

We assess reliability using two complementary approaches:

1. **Cronbach's Alpha (α)** - Classical test theory measure of internal consistency
2. **McDonald's Omega (ω)** - Model-based alternative that accounts for hierarchical factor structure

Per McDonald (1999) and Revelle & Zinbarg (2009), omega is often preferred when items may reflect both a general factor and group factors, as is likely in the EQRI questionnaire with its seven indicators.

### Sample Size Adequacy

```{r run-reliability-analysis}
#| echo: false
#| message: false
#| warning: false

# Run comprehensive reliability analysis across all contexts
reliability_results <- run_reliability_analysis(
  responses_df,
  include_omega = TRUE
)
```

```{r plot-sample-sizes}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Sample sizes by context determine the stability of reliability estimates"

plot_sample_sizes(reliability_results)
```

**Interpretation:**

Per Kline (2000):
- **Green line (n=100)**: Excellent sample size for stable reliability estimates
- **Orange line (n=30)**: Minimum adequate sample size; estimates should be interpreted cautiously

```{r summarize-sample-sizes}
#| echo: false
#| results: asis

n_above_100 <- sum(reliability_results$n_observations >= 100)
n_above_30 <- sum(reliability_results$n_observations >= 30)
n_total <- nrow(reliability_results)

cat(sprintf("**Contexts with n ≥ 100**: %d of %d (%.0f%%)\n\n", 
            n_above_100, n_total, 100 * n_above_100 / n_total))
cat(sprintf("**Contexts with n ≥ 30**: %d of %d (%.0f%%)\n\n", 
            n_above_30, n_total, 100 * n_above_30 / n_total))

if (n_above_100 < n_total) {
  cat("⚠️ Contexts with smaller sample sizes (n < 100) will produce less stable reliability estimates.\n\n")
}
```

---

### Cronbach's Alpha Results

Cronbach's alpha measures whether questions within each context consistently measure the same underlying construct.

```{r plot-alpha}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Cronbach's alpha by context (standardized)"

plot_reliability_comparison(
  reliability_results,
  metric = "alpha_std"
)
```

**Interpretation Guidelines** (George & Mallery, 2003):

| Alpha Range | Interpretation | Implication |
|-------------|----------------|-------------|
| α ≥ 0.90 | Excellent | Very high internal consistency |
| 0.80 ≤ α < 0.90 | Good | High internal consistency |
| 0.70 ≤ α < 0.80 | Acceptable | Adequate for research purposes |
| 0.60 ≤ α < 0.70 | Questionable | Marginal; use with caution |
| α < 0.60 | Poor | Unacceptable; reconsider questionnaire |

```{r alpha-interpretation-summary}
#| echo: false
#| results: asis

alpha_summary <- reliability_results %>%
  group_by(alpha_interpretation) %>%
  summarize(n_contexts = n(), .groups = "drop") %>%
  arrange(desc(alpha_interpretation))

cat("\n**Summary by Interpretation Level:**\n\n")
for (i in 1:nrow(alpha_summary)) {
  cat(sprintf("- **%s**: %d context(s)\n", 
              alpha_summary$alpha_interpretation[i],
              alpha_summary$n_contexts[i]))
}
cat("\n")

# Calculate percentage meeting minimum threshold
n_acceptable <- sum(reliability_results$alpha_std >= 0.70, na.rm = TRUE)
pct_acceptable <- 100 * n_acceptable / nrow(reliability_results)

cat(sprintf("**Contexts meeting α ≥ 0.70 threshold**: %d of %d (%.0f%%)\n\n",
            n_acceptable, nrow(reliability_results), pct_acceptable))
```

---

### Alpha vs. Omega Comparison

McDonald's omega (ω) provides an alternative reliability estimate that:
- Does not assume tau-equivalence (equal factor loadings)
- Accounts for multidimensional structure
- Generally provides more accurate estimates (Revelle & Zinbarg, 2009)

```{r plot-both-metrics}
#| echo: false
#| fig-width: 10
#| fig-height: 6
#| fig-cap: "Comparison of Cronbach's alpha and McDonald's omega"

plot_reliability_comparison(
  reliability_results,
  metric = "both"
)
```

**Reading the comparison:**

- **Top panel (Alpha)**: Classical measure assuming equal item contributions
- **Bottom panel (Omega)**: Model-based measure accounting for factor structure
- Contexts where omega > alpha suggest multidimensionality (which is expected given 7 indicators)

---

### Omega-Alpha Difference

The difference ω - α indicates the degree to which the questionnaire is multidimensional.

```{r plot-difference}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Omega minus alpha: positive values indicate multidimensionality"

plot_reliability_comparison(
  reliability_results,
  metric = "difference"
)
```

**Interpretation** (Revelle & Zinbarg, 2009):

- **ω - α ≈ 0**: Unidimensional structure (items load on single factor)
- **ω - α > 0**: Multidimensional structure (items load on multiple factors)
- **Large positive difference**: Strong evidence of distinct subscales

```{r difference-analysis}
#| echo: false
#| results: asis

# Calculate mean difference
mean_diff <- mean(reliability_results$omega_alpha_diff, na.rm = TRUE)
median_diff <- median(reliability_results$omega_alpha_diff, na.rm = TRUE)

cat(sprintf("**Mean ω - α difference**: %.3f\n\n", mean_diff))
cat(sprintf("**Median ω - α difference**: %.3f\n\n", median_diff))

if (mean_diff > 0.05) {
  cat("**Finding**: The consistent positive difference between omega and alpha (ω > α) confirms that the EQRI questionnaire measures **multiple distinct dimensions** (the seven indicators), which is by design.\n\n")
  cat("This supports proceeding with factor analysis to identify these underlying factor structures.\n\n")
} else {
  cat("**Finding**: Minimal difference between omega and alpha suggests relatively unidimensional structure.\n\n")
}
```

---

### Reliability Assessment Summary

```{r reliability-table}
#| echo: false
#| tbl-cap: "Comprehensive reliability analysis results by context"

# Create formatted table
summary_table <- create_reliability_table(reliability_results)

kable(summary_table,
      format = "pipe",
      align = c("l", "l", "r", "r", "r", "r", "l"))
```

---

### Implications for Factor Analysis

**Contexts suitable for factor analysis** (meeting both criteria):

1. **Sufficient sample size**: n ≥ 30 (minimum) or n ≥ 100 (preferred)
2. **Adequate reliability**: α ≥ 0.70 (acceptable minimum)

```{r fa-readiness}
#| echo: false
#| results: asis

# Identify contexts ready for factor analysis
fa_ready <- reliability_results %>%
  filter(n_observations >= 30, alpha_std >= 0.70)

fa_preferred <- reliability_results %>%
  filter(n_observations >= 100, alpha_std >= 0.70)

cat(sprintf("**Contexts meeting minimum criteria** (n ≥ 30, α ≥ 0.70): %d of %d\n\n",
            nrow(fa_ready), nrow(reliability_results)))

cat(sprintf("**Contexts meeting preferred criteria** (n ≥ 100, α ≥ 0.70): %d of %d\n\n",
            nrow(fa_preferred), nrow(reliability_results)))

if (nrow(fa_preferred) > 0) {
  cat("**Recommended for immediate factor analysis:**\n\n")
  for (i in 1:nrow(fa_preferred)) {
    cat(sprintf("- %s × %s (n=%d, α=%.2f, ω=%.2f)\n",
                fa_preferred$program[i],
                fa_preferred$milestone[i],
                fa_preferred$n_observations[i],
                fa_preferred$alpha_std[i],
                fa_preferred$omega_total[i]))
  }
  cat("\n")
}

# Identify contexts needing more data
needs_more_data <- reliability_results %>%
  filter(n_observations < 30 | alpha_std < 0.70)

if (nrow(needs_more_data) > 0) {
  cat(sprintf("**Contexts requiring additional data or questionnaire refinement**: %d\n\n",
              nrow(needs_more_data)))
  
  for (i in 1:min(5, nrow(needs_more_data))) {
    issue <- if (needs_more_data$n_observations[i] < 30) "small n" else "low α"
    cat(sprintf("- %s × %s (%s)\n",
                needs_more_data$program[i],
                needs_more_data$milestone[i],
                issue))
  }
  cat("\n")
}
```

---

### Key Findings

::: {.callout-important}
## Reliability Conclusion

Based on this comprehensive reliability assessment:

1. **Internal Consistency**: The questionnaire demonstrates **`r if(mean(reliability_results$alpha_std, na.rm=TRUE) >= 0.8) "good to excellent" else if(mean(reliability_results$alpha_std, na.rm=TRUE) >= 0.7) "acceptable" else "variable"`** internal consistency across contexts (mean α = `r round(mean(reliability_results$alpha_std, na.rm=TRUE), 2)`).

2. **Multidimensionality**: Omega consistently exceeds alpha (mean difference = `r round(mean_diff, 3)`), confirming the questionnaire's **multi-indicator structure** aligns with its design to measure seven distinct quality aspects.

3. **Factor Analysis Readiness**: **`r nrow(fa_ready)` of `r nrow(reliability_results)` contexts** meet minimum psychometric standards for exploratory factor analysis.

4. **Sample Size Considerations**: Contexts with n < 100 should be interpreted with appropriate caution regarding the stability of factor solutions.
:::


---

## Step 1: Identify Core Items

Per Schafer & Graham (2002), items with <20% missing data are suitable for aggregated analysis.

```{r identify-core-items}
#| echo: false
#| message: false
#| warning: false

# Calculate coverage per question
question_coverage <- responses_df %>%
  group_by(QUESTION_NUMBER) %>%
  summarise(
    n_contexts = n_distinct(paste(PROGRAMTYPE_NAME, MILESTONE_DESC)),
    n_total_responses = n(),
    pct_missing = round(100 * (1 - n_total_responses / nrow(responses_df)), 1),
    .groups = "drop"
  ) %>%
  arrange(desc(n_contexts), QUESTION_NUMBER)

# Define core threshold: >= 8 of 10 contexts (80% coverage, <20% missing)
core_threshold <- 8

core_items <- question_coverage %>%
  filter(n_contexts >= core_threshold) %>%
  pull(QUESTION_NUMBER)

context_specific_items <- question_coverage %>%
  filter(n_contexts < core_threshold)
```

### Core Items Summary

```{r display-core-summary}
#| echo: false
#| results: asis

cat("**Total questions in questionnaire:**", nrow(question_coverage), "\n\n")
cat("**Core items (≥", core_threshold, "of 10 contexts):**", length(core_items), "\n\n")
cat("**Context-specific items (<", core_threshold, "of 10 contexts):**", nrow(context_specific_items), "\n\n")
```

Per Schafer & Graham (2002), focusing on core items for factor analysis is appropriate when items exhibit structural missingness by design.

---

## Step 2: Questionnaire Response Rates

Understanding response patterns by context helps identify which contexts have sufficient data for factor analysis.

### Total Contexts

```{r count-contexts}
#| echo: false
#| message: false
#| warning: false

# Count unique contexts
n_programs <- n_distinct(responses_df$PROGRAMTYPE_NAME)
n_milestones <- n_distinct(responses_df$MILESTONE_DESC)
n_contexts <- n_distinct(paste(responses_df$PROGRAMTYPE_NAME, responses_df$MILESTONE_DESC))
```

```{r display-context-counts}
#| echo: false
#| results: asis

cat("**Program types:**", n_programs, "(Military, Civil Works)\n\n")
cat("**Milestones:**", n_milestones, "\n\n")
cat("**Total contexts (program × milestone):**", n_contexts, "\n\n")
```

### Response Rates by Context

```{r context-summary}
#| echo: false
#| message: false
#| warning: false

# Calculate responses by context
context_summary <- responses_df %>%
  group_by(PROGRAMTYPE_NAME, MILESTONE_DESC) %>%
  summarize(
    n_events = length(unique(QUESTIONNAIREEVENT_ID)),
    n_responses = n(),
    .groups = "drop"
  ) %>%
  arrange(PROGRAMTYPE_NAME, MILESTONE_DESC)
```

```{r display-context-table}
#| echo: false
#| tbl-cap: "Questionnaire Events and Responses by Context"

kable(context_summary, 
      col.names = c("Program Type", "Milestone", "Events", "Responses"),
      format = "pipe",
      align = c("l", "l", "r", "r"))
```

**Events:** Number of unique questionnaire submissions (one per project at each milestone)  
**Responses:** Total individual question responses (events × questions asked in that context)

---

## Step 3: Question Coverage Across Contexts

This table shows which questions are "core" (broadly applicable) vs. context-specific.

```{r display-coverage-table}
#| echo: false
#| tbl-cap: "Question Coverage: Core vs. Context-Specific Items"

# Add classification column
question_coverage_display <- question_coverage %>%
  mutate(
    Type = ifelse(n_contexts >= core_threshold, "Core", "Context-specific")
  ) %>%
  select(QUESTION_NUMBER, Type, n_contexts, n_total_responses, pct_missing)

kable(question_coverage_display,
      col.names = c("Question", "Type", "# Contexts", "Total Responses", "% Missing"),
      format = "pipe",
      align = c("l", "l", "r", "r", "r"))
```

**Interpretation:**

- **Core items** (≥8 contexts): Analyzed across contexts
- **Context-specific items** (<8 contexts): Analyzed only within their relevant contexts

---

## Step 4: Factorability by Context

Due to data characteristics (varying response patterns, context-specific items), we assess factorability for the contexts with the most data.

```{r identify-largest-contexts}
#| echo: false
#| message: false
#| warning: false

# Identify contexts with most events for analysis
largest_contexts <- context_summary %>%
  arrange(desc(n_events)) %>%
  head(4)  # Top 4 contexts

contexts_list <- largest_contexts %>%
  mutate(context_name = paste(PROGRAMTYPE_NAME, MILESTONE_DESC, sep = " × ")) %>%
  pull(context_name)
```

### Contexts Selected for Analysis

Based on response rates, we focus on the **four largest contexts**:

```{r display-selected-contexts}
#| echo: false
#| results: asis

for (i in 1:nrow(largest_contexts)) {
  cat(sprintf("%d. **%s × %s** - %d events\n", 
              i, 
              largest_contexts$PROGRAMTYPE_NAME[i],
              largest_contexts$MILESTONE_DESC[i],
              largest_contexts$n_events[i]))
}
cat("\n")
```

::: {.callout-note}
## Why Context-Specific Analysis?

Per Little & Rubin (2019), aggregating across contexts with structural missingness can produce misleading results. Analyzing within contexts ensures:

- Only applicable questions are included
- Response patterns are interpretable
- Sufficient variability for correlation analysis
:::

---

### Military × 95% (Final Design)

```{r military-95-factorability}
#| echo: false
#| message: false
#| warning: false
#| error: true

# Try to assess factorability
military_95_result <- tryCatch({
  assess_factorability(
    responses_df,
    program_name = "Military",
    milestone_name = "95% (Final Design)",
    filter_context_specific = TRUE
  )
}, error = function(e) {
  list(
    error = TRUE,
    message = e$message
  )
})
```

```{r display-military-95}
#| echo: false
#| results: asis

if (!is.null(military_95_result$error) && military_95_result$error) {
  cat("⚠️ **Analysis Issue**\n\n")
  cat("```\n")
  cat(military_95_result$message)
  cat("\n```\n\n")
  cat("*This context requires further data preparation before factor analysis.*\n\n")
} else {
  cat("#### Sample Size\n\n")
  cat("**Complete cases:**", military_95_result$sample$n_observations_complete, "\n\n")
  cat("**Items analyzed:**", military_95_result$sample$n_questions_analyzed, "\n\n")
  
  cat("#### Assessment\n\n")
  if (military_95_result$recommendation$proceed_with_fa) {
    cat("✅ **Suitable for factor analysis**\n\n")
  } else {
    cat("❌ **Not suitable for factor analysis**\n\n")
  }
  
  cat("**Rationale:** ", military_95_result$recommendation$rationale, "\n\n")
}
```

---

### Civil Works × 95% (Final Design)

```{r cw-95-factorability}
#| echo: false
#| message: false
#| warning: false
#| error: true

# Try to assess factorability
cw_95_result <- tryCatch({
  assess_factorability(
    responses_df,
    program_name = "Civil Works",
    milestone_name = "95% (Final Design)",
    filter_context_specific = TRUE
  )
}, error = function(e) {
  list(
    error = TRUE,
    message = e$message
  )
})
```

```{r display-cw-95}
#| echo: false
#| results: asis

if (!is.null(cw_95_result$error) && cw_95_result$error) {
  cat("⚠️ **Analysis Issue**\n\n")
  cat("```\n")
  cat(cw_95_result$message)
  cat("\n```\n\n")
  cat("*This context requires further data preparation before factor analysis.*\n\n")
} else {
  cat("#### Sample Size\n\n")
  cat("**Complete cases:**", cw_95_result$sample$n_observations_complete, "\n\n")
  cat("**Items analyzed:**", cw_95_result$sample$n_questions_analyzed, "\n\n")
  
  cat("#### Assessment\n\n")
  if (cw_95_result$recommendation$proceed_with_fa) {
    cat("✅ **Suitable for factor analysis**\n\n")
  } else {
    cat("❌ **Not suitable for factor analysis**\n\n")
  }
  
  cat("**Rationale:** ", cw_95_result$recommendation$rationale, "\n\n")
}
```

---

### Additional Contexts

```{r remaining-contexts}
#| echo: false
#| message: false
#| warning: false
#| error: true
#| results: asis

# Analyze remaining two largest contexts
if (nrow(largest_contexts) >= 3) {
  for (i in 3:min(4, nrow(largest_contexts))) {
    prog <- largest_contexts$PROGRAMTYPE_NAME[i]
    mile <- largest_contexts$MILESTONE_DESC[i]
    
    cat(sprintf("\n### %s × %s\n\n", prog, mile))
    
    result <- tryCatch({
      assess_factorability(
        responses_df,
        program_name = prog,
        milestone_name = mile,
        filter_context_specific = TRUE
      )
    }, error = function(e) {
      list(error = TRUE, message = e$message)
    })
    
    if (!is.null(result$error) && result$error) {
      cat("⚠️ **Analysis Issue**\n\n")
      cat("```\n")
      cat(result$message)
      cat("\n```\n\n")
    } else {
      cat("**Complete cases:**", result$sample$n_observations_complete, " | ")
      cat("**Items analyzed:**", result$sample$n_questions_analyzed, "\n\n")
      
      if (result$recommendation$proceed_with_fa) {
        cat("✅ Suitable for factor analysis\n\n")
      } else {
        cat("❌ Not suitable\n\n")
      }
      
      cat("**Rationale:** ", result$recommendation$rationale, "\n\n")
    }
  }
}
```

---

## Summary & Interpretation

### Data Characteristics Observed

```{r summary-stats}
#| echo: false
#| message: false
#| warning: false

# Calculate summary statistics
total_questions <- nrow(question_coverage)
n_core <- length(core_items)
n_context_specific <- nrow(context_specific_items)
pct_core <- round(100 * n_core / total_questions, 1)
```

```{r display-summary}
#| echo: false
#| results: asis

cat("**Questionnaire composition:**\n\n")
cat("- Core items:", n_core, sprintf("(%.1f%%)", pct_core), "\n")
cat("- Context-specific items:", n_context_specific, sprintf("(%.1f%%)", 100 - pct_core), "\n\n")

cat("**Implication:** The questionnaire is highly tailored to specific contexts, ")
cat("which is appropriate for capturing milestone- and program-specific quality factors.\n\n")
```

### Recommendations

Based on the factorability assessments:

1. **For contexts with sufficient data:**
   - Proceed with parallel analysis to determine number of factors
   - Conduct exploratory factor analysis (EFA)
   - Examine factor loadings to interpret question groupings

2. **For contexts with insufficient data:**
   - Continue data collection
   - Consider aggregating across similar milestones (e.g., early vs. late design phases)
   - Focus on descriptive statistics until sample size increases

3. **Context-specific items:**
   - Valuable for targeted quality assessment within their domains
   - Should not be removed - they serve important context-specific purposes
   - Will be analyzed separately as data accrues

---

## Next Steps

**Immediate:**
1. Review factorability results with subject matter experts
2. For viable contexts, run parallel analysis to determine optimal factor structure
3. Conduct EFA on contexts meeting psychometric requirements

**Future:**
1. As data accrues, reassess contexts with currently insufficient samples
2. Investigate whether context-specific items cluster into meaningful factors within their domains
3. Test measurement invariance on core items across contexts (Phase 4)

---

## References

- Cronbach, L. J. (1951). Coefficient alpha and the internal structure of tests. *Psychometrika*, 16(3), 297-334. https://doi.org/10.1007/BF02310555

- George, D., & Mallery, P. (2003). *SPSS for Windows Step by Step: A Simple Guide and Reference* (4th ed.). Allyn & Bacon.

- Kline, P. (2000). *The Handbook of Psychological Testing* (2nd ed.). Routledge.

- McDonald, R. P. (1999). *Test Theory: A Unified Treatment*. Lawrence Erlbaum Associates.

- Revelle, W., & Zinbarg, R. E. (2009). Coefficients alpha, beta, omega, and the glb: Comments on Sijtsma. *Psychometrika*, 74(1), 145-154. https://doi.org/10.1007/s11336-008-9102-z

- Streiner, D. L., Norman, G. R., & Cairney, J. (2015). *Health Measurement Scales: A Practical Guide to their Development and Use* (5th ed.). Oxford University Press.

- Flora, D. B., & Curran, P. J. (2004). An empirical evaluation of alternative methods of estimation for confirmatory factor analysis with ordinal data. *Psychological Methods, 9*(4), 466-491. https://doi.org/10.1037/1082-989X.9.4.466

- Little, R. J. A., & Rubin, D. B. (2019). *Statistical Analysis with Missing Data* (3rd ed.). Wiley. Chapter 1.4.3: Structural missingness by design.

- Schafer, J. L., & Graham, J. W. (2002). Missing data: Our view of the state of the art. *Psychological Methods, 7*(2), 147-177. https://doi.org/10.1037/1076-8986.7.2.147


---

## Determining the Number of Factors

Before conducting exploratory factor analysis, we must determine the optimal number of factors to extract. Per Hayton et al. (2004), **parallel analysis** (Horn, 1965) is the gold standard method, which we complement with VSS (Revelle & Rocklin, 1979) and MAP test (Velicer, 1976).

### Selected Contexts for Analysis

Based on reliability assessment results, we proceed with these four contexts meeting both sample size (n ≥ 30) and reliability (α ≥ 0.70) criteria:

```{r define-contexts}
#| echo: false

# Define contexts for factor analysis
fa_contexts <- data.frame(
  program = c("Military", "Military", "Military", "Civil Works"),
  milestone = c(
    "100% (Corrected Final Design)",
    "35% (Concept Design)",
    "95% (Final Design)",
    "15% (Project Initiation)"
  ),
  stringsAsFactors = FALSE
)

# Display
kable(fa_contexts,
      col.names = c("Program Type", "Milestone"),
      caption = "Contexts selected for factor analysis")
```

---

### Factor Number Analysis

```{r run-factor-number-analysis}
#| echo: false
#| message: false
#| warning: false

# Run comprehensive factor number analysis
factor_number_results <- run_factor_number_analysis(
  responses_df,
  contexts = fa_contexts,
  n.iter = 20,    # Per psych::fa.parallel() documentation
  fm = "minres"   # Minimum residual method (default)
)
```

---

### Scree Plots with Parallel Analysis

Per Horn (1965) and Hayton et al. (2004), retain factors where **observed eigenvalues exceed simulated eigenvalues** from random data.

```{r plot-scree-plots}
#| echo: false
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Scree plots with parallel analysis for each context"
#| layout-ncol: 2

# Extract and display scree plots
for (context_name in names(factor_number_results)) {
  result <- factor_number_results[[context_name]]
  if (!is.null(result) && !is.null(result$scree_plot)) {
    print(result$scree_plot)
  }
}
```

**Interpretation:**

- **Blue solid line**: Observed eigenvalues from actual data
- **Purple dashed line**: Mean eigenvalues from simulated random data  
- **Orange dotted line**: 95th percentile of simulated eigenvalues
- **Gray line**: Kaiser criterion (eigenvalue = 1.0)

**Rule**: Retain factors where the blue line is **above** the purple/orange lines.

---

### Factor Number Recommendations

```{r factor-number-table}
#| echo: false
#| tbl-cap: "Recommended number of factors by method and context"

# Create summary table
factor_summary <- create_factor_number_table(factor_number_results)

kable(factor_summary,
      format = "pipe",
      align = c("l", "r", "r", "r", "r", "r"))
```

**Methods Explained:**

1. **Parallel Analysis** (Horn, 1965) - *PRIMARY RECOMMENDATION*
   - Compares eigenvalues to those from random data
   - Retains factors exceeding random chance
   - Most reliable method per Hayton et al. (2004)

2. **VSS-1** (Revelle & Rocklin, 1979)
   - Very Simple Structure with complexity 1
   - Assumes items load on single factor

3. **MAP Test** (Velicer, 1976)
   - Minimum Average Partial correlation
   - Alternative to parallel analysis

---

### Convergence Analysis

```{r analyze-convergence}
#| echo: false
#| results: asis

# Check if methods agree
for (context_name in names(factor_number_results)) {
  result <- factor_number_results[[context_name]]
  
  if (is.null(result)) next
  
  cat(sprintf("\n#### %s\n\n", context_name))
  
  pa <- result$recommendations$parallel_analysis
  vss <- result$recommendations$vss_complexity1
  map <- result$recommendations$map_test
  
  cat(sprintf("- **Parallel Analysis**: %s factors\n", 
              ifelse(!is.null(pa), pa, "N/A")))
  cat(sprintf("- **VSS Complexity 1**: %s factors\n", 
              ifelse(!is.null(vss), vss, "N/A")))
  cat(sprintf("- **MAP Test**: %s factors\n\n", 
              ifelse(!is.null(map), map, "N/A")))
  
  # Check convergence
  methods <- c(pa, vss, map)
  methods <- methods[!is.na(methods)]
  
  if (length(unique(methods)) == 1) {
    cat(sprintf("✓ **Perfect agreement**: All methods suggest **%d factors**\n\n", 
                unique(methods)))
  } else if (length(methods) > 0) {
    cat(sprintf("⚠️ Methods suggest between %d and %d factors. ", 
                min(methods), max(methods)))
    cat(sprintf("**Recommendation**: Use **parallel analysis result (%d factors)** as primary guide.\n\n",
                pa))
  }
}
```

---

### Theoretical Expectations

The EQRI questionnaire was designed to measure **seven quality indicators**:

1. Confidence
2. Cost
3. QA (Quality Assurance)
4. QC (Quality Control)
5. Schedule
6. Scope
7. Team

**Expected finding**: If the questionnaire structure aligns with its design, we would expect approximately **7 factors** to emerge, though exploratory factor analysis may reveal:

- Fewer factors if indicators are highly correlated
- More factors if indicators have distinct sub-dimensions
- Different numbers across contexts if question sets vary

---

### Recommended Number of Factors

::: {.callout-important}
## Factor Extraction Recommendation

Based on **parallel analysis** (the gold standard per Hayton et al., 2004):

```{r final-recommendation}
#| echo: false
#| results: asis

for (context_name in names(factor_number_results)) {
  result <- factor_number_results[[context_name]]
  if (!is.null(result) && !is.null(result$recommendations$parallel_analysis)) {
    cat(sprintf("- **%s**: Extract **%d factors**\n", 
                context_name,
                result$recommendations$parallel_analysis))
  }
}
```

These recommendations will guide the exploratory factor analysis in the next phase.
:::

---

### References

- Cattell, R. B. (1966). The scree test for the number of factors. *Multivariate Behavioral Research*, 1(2), 245-276. https://doi.org/10.1207/s15327906mbr0102_10

- Hayton, J. C., Allen, D. G., & Scarpello, V. (2004). Factor retention decisions in exploratory factor analysis: A tutorial on parallel analysis. *Organizational Research Methods*, 7(2), 191-205. https://doi.org/10.1177/1094428104263675

- Horn, J. L. (1965). A rationale and test for the number of factors in factor analysis. *Psychometrika*, 30(2), 179-185. https://doi.org/10.1007/BF02289447

- Revelle, W., & Rocklin, T. (1979). Very simple structure: An alternative procedure for estimating the optimal number of interpretable factors. *Multivariate Behavioral Research*, 14(4), 403-414. https://doi.org/10.1207/s15327906mbr1404_2

- Velicer, W. F. (1976). Determining the number of components from the matrix of partial correlations. *Psychometrika*, 41(3), 321-327. https://doi.org/10.1007/BF02293557


---

## Exploratory Factor Analysis

Having established adequate reliability, we now conduct exploratory factor analysis (EFA) to uncover the underlying factor structure of the EQRI questionnaire.

### Selected Contexts for Analysis

Based on the reliability assessment and sample size requirements, we proceed with these contexts:

```{r define-efa-contexts}
#| echo: false

# Define contexts for EFA based on Phase 1 analysis
# Per your workplan: Military × 100%, Military × 35%, Military × 95%, Civil Works × 15%
efa_contexts <- data.frame(
  program = c("Military", "Military", "Military", "Civil Works"),
  milestone = c(
    "100% (Corrected Final Design)",
    "35% (Concept Design)", 
    "95% (Final Design)",
    "15% (Project Initiation)"
  ),
  n_events = c(52, 51, 41, 38),  # From reliability analysis
  stringsAsFactors = FALSE
)

kable(efa_contexts,
      col.names = c("Program Type", "Milestone", "Events"),
      caption = "Contexts selected for exploratory factor analysis")
```

---

### Factor Analysis: Military × 100% (Corrected Final Design)

This context has the largest sample size (n=52) and provides the most stable factor solution.

```{r efa-military-100}
#| echo: false
#| message: false
#| warning: false

# Run EFA with recommended settings per workplan
efa_mil_100 <- run_efa(
  responses_df,
  program_name = "Military",
  milestone_name = "100% (Corrected Final Design)",
  nfactors = NULL,  # Use parallel analysis
  fm = "minres",    # Minimum residual (robust for small samples)
  rotate = "oblimin"  # Oblique rotation (factors may correlate)
)
```

#### Questions Analyzed

```{r question-labels-mil-100}
#| echo: false
#| tbl-cap: "Questions included in factor analysis"

question_labels_mil_100 <- create_question_label_table(efa_mil_100, responses_df)

kable(question_labels_mil_100,
      col.names = c("Item", "Question", "Indicator"),
      format = "pipe",
      align = c("l", "l", "l"))
```

**Interpretation:**

Per Costello & Osborne (2005), loadings above **0.40** are considered meaningful:

- **Salient loadings** (shown): Items that define each factor
- **Suppressed loadings** (blank): Items with weak relationships to the factor
- **Blue (positive)**: Item increases with factor score
- **Purple (negative)**: Item decreases with factor score

```{r interpret-loadings-mil-100}
#| echo: false
#| results: asis

# Identify primary factor for each item
loadings_clean <- efa_mil_100$loadings_clean
n_factors <- efa_mil_100$extraction_info$nfactors

cat(sprintf("\n**Factor structure identified**: %d factors extracted\n\n", n_factors))

# Show which questions load on each factor
for (f in 1:n_factors) {
  factor_name <- colnames(loadings_clean)[f]
  
  # Get items with meaningful loadings
  meaningful_items <- which(abs(loadings_clean[[factor_name]]) >= 0.40)
  
  if (length(meaningful_items) > 0) {
    cat(sprintf("**%s**:\n", factor_name))
    
    for (item_idx in meaningful_items) {
      item_name <- rownames(loadings_clean)[item_idx]
      loading_val <- loadings_clean[[factor_name]][item_idx]
      
      # Get question label
      question_num <- as.integer(gsub("Q", "", item_name))
      question_label <- question_labels_mil_100$Label[
        question_labels_mil_100$Item == item_name
      ]
      
      cat(sprintf("- %s (%.2f): %s\n", 
                  item_name, loading_val, question_label))
    }
    cat("\n")
  }
}

# Check for problematic items
if (!is.null(efa_mil_100$problematic_items$low_communality)) {
  cat(sprintf("\n⚠️ **%d item(s) with low communality** (< 0.25):\n",
              length(efa_mil_100$problematic_items$low_communality)))
  
  for (item in efa_mil_100$problematic_items$low_communality) {
    question_label <- question_labels_mil_100$Label[
      question_labels_mil_100$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items are poorly represented by the factor solution.\n\n")
}

if (!is.null(efa_mil_100$problematic_items$cross_loadings)) {
  cat(sprintf("\n⚠️ **%d item(s) with cross-loadings** (> 0.32 on multiple factors):\n",
              length(efa_mil_100$problematic_items$cross_loadings)))
  
  for (item in efa_mil_100$problematic_items$cross_loadings) {
    question_label <- question_labels_mil_100$Label[
      question_labels_mil_100$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items may need review for theoretical clarity.\n\n")
}
```

#### Variance Explained

```{r plot-variance-mil-100}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Variance explained by each factor"

plot_variance_explained(efa_mil_100)
```

**Interpretation:**

- **Individual variance**: Each factor's unique contribution
- **Cumulative variance**: Total explained by all factors
- **Adequacy**: Per psychometric standards, cumulative variance should exceed 50%

```{r variance-summary-mil-100}
#| echo: false
#| results: asis

variance_mat <- efa_mil_100$variance_explained
total_var <- variance_mat["Cumulative Var", ncol(variance_mat)] * 100

cat(sprintf("**Total variance explained**: %.1f%%\n\n", total_var))

if (total_var >= 60) {
  cat("✅ **Good**: Factor solution explains substantial variance\n\n")
} else if (total_var >= 50) {
  cat("✓ **Adequate**: Factor solution explains acceptable variance\n\n")
} else {
  cat("⚠️ **Poor**: Factor solution explains limited variance\n\n")
}
```

#### Factor Correlations

```{r plot-correlations-mil-100}
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| fig-cap: "Correlations between factors (oblimin rotation)"

plot_factor_correlations(efa_mil_100)
```

**Interpretation:**

- **|r| < 0.30**: Factors are relatively independent
- **0.30 ≤ |r| < 0.70**: Factors are moderately related
- **|r| ≥ 0.70**: Factors are highly related (may consider combining)

Oblique rotation (oblimin) allows factors to correlate, which is typically more realistic for psychosocial constructs like quality indicators.

#### Model Fit

```{r fit-summary-mil-100}
#| echo: false
#| results: asis

fit <- efa_mil_100$fit_indices
fit_interp <- efa_mil_100$fit_interpretation

cat("**Fit Indices** (per Hu & Bentler, 1999):\n\n")
cat(sprintf("- **TLI**: %.3f (%s)\n", fit$TLI, fit_interp$TLI))
cat(sprintf("- **RMSEA**: %.3f (%s)\n", fit$RMSEA, fit_interp$RMSEA))
cat(sprintf("- **RMSR**: %.3f (%s)\n\n", fit$RMSR, fit_interp$RMSR))

cat("**Interpretation:**\n\n")
cat("- **TLI (Tucker-Lewis Index)**: > 0.95 excellent, > 0.90 acceptable\n")
cat("- **RMSEA (Root Mean Square Error of Approximation)**: < 0.05 excellent, < 0.08 acceptable\n")
cat("- **RMSR (Root Mean Square Residual)**: < 0.05 good\n\n")
```

---

---

### Factor Analysis: Military × 35% (Concept Design)

Early design phase context (n=51).

```{r efa-military-35}
#| echo: false
#| message: false
#| warning: false

efa_mil_35 <- run_efa(
  responses_df,
  program_name = "Military",
  milestone_name = "35% (Concept Design)",
  nfactors = NULL,
  fm = "minres",
  rotate = "oblimin"
)
```

#### Questions Analyzed

```{r question-labels-mil-35}
#| echo: false
#| tbl-cap: "Questions included in factor analysis - Military × 35%"

question_labels_mil_35 <- create_question_label_table(efa_mil_35, responses_df)

kable(question_labels_mil_35,
      col.names = c("Item", "Question", "Indicator"),
      format = "pipe",
      align = c("l", "l", "l"))
```

#### Factor Loadings

```{r plot-loadings-mil-35}
#| echo: false
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Factor loading matrix for Military × 35% (Concept Design)"

plot_factor_loadings(efa_mil_35, use_clean = TRUE)
```

**Interpretation:**

```{r interpret-loadings-mil-35}
#| echo: false
#| results: asis

# Identify primary factor for each item
loadings_clean_35 <- efa_mil_35$loadings_clean
n_factors_35 <- efa_mil_35$extraction_info$nfactors

cat(sprintf("\n**Factor structure identified**: %d factors extracted\n\n", n_factors_35))

# Show which questions load on each factor
for (f in 1:n_factors_35) {
  factor_name <- colnames(loadings_clean_35)[f]
  
  # Get items with meaningful loadings
  meaningful_items <- which(abs(loadings_clean_35[[factor_name]]) >= 0.40)
  
  if (length(meaningful_items) > 0) {
    cat(sprintf("**%s**:\n", factor_name))
    
    for (item_idx in meaningful_items) {
      item_name <- rownames(loadings_clean_35)[item_idx]
      loading_val <- loadings_clean_35[[factor_name]][item_idx]
      
      # Get question label
      question_label <- question_labels_mil_35$Label[
        question_labels_mil_35$Item == item_name
      ]
      
      cat(sprintf("- %s (%.2f): %s\n", 
                  item_name, loading_val, question_label))
    }
    cat("\n")
  }
}

# Check for problematic items
if (!is.null(efa_mil_35$problematic_items$low_communality)) {
  cat(sprintf("\n⚠️ **%d item(s) with low communality** (< 0.25):\n",
              length(efa_mil_35$problematic_items$low_communality)))
  
  for (item in efa_mil_35$problematic_items$low_communality) {
    question_label <- question_labels_mil_35$Label[
      question_labels_mil_35$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items are poorly represented by the factor solution.\n\n")
}

if (!is.null(efa_mil_35$problematic_items$cross_loadings)) {
  cat(sprintf("\n⚠️ **%d item(s) with cross-loadings** (> 0.32 on multiple factors):\n",
              length(efa_mil_35$problematic_items$cross_loadings)))
  
  for (item in efa_mil_35$problematic_items$cross_loadings) {
    question_label <- question_labels_mil_35$Label[
      question_labels_mil_35$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items may need review for theoretical clarity.\n\n")
}
```

#### Variance Explained

```{r plot-variance-mil-35}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Variance explained - Military × 35%"

plot_variance_explained(efa_mil_35)
```

```{r variance-summary-mil-35}
#| echo: false
#| results: asis

variance_mat_35 <- efa_mil_35$variance_explained
total_var_35 <- variance_mat_35["Cumulative Var", ncol(variance_mat_35)] * 100

cat(sprintf("**Total variance explained**: %.1f%%\n\n", total_var_35))

if (total_var_35 >= 60) {
  cat("✅ **Good**: Factor solution explains substantial variance\n\n")
} else if (total_var_35 >= 50) {
  cat("✓ **Adequate**: Factor solution explains acceptable variance\n\n")
} else {
  cat("⚠️ **Poor**: Factor solution explains limited variance\n\n")
}
```

---

### Factor Analysis: Military × 95% (Final Design)

Late design phase context (n=41).

```{r efa-military-95}
#| echo: false
#| message: false
#| warning: false

efa_mil_95 <- run_efa(
  responses_df,
  program_name = "Military",
  milestone_name = "95% (Final Design)",
  nfactors = NULL,
  fm = "minres",
  rotate = "oblimin"
)
```

#### Questions Analyzed

```{r question-labels-mil-95}
#| echo: false
#| tbl-cap: "Questions included in factor analysis - Military × 95%"

question_labels_mil_95 <- create_question_label_table(efa_mil_95, responses_df)

kable(question_labels_mil_95,
      col.names = c("Item", "Question", "Indicator"),
      format = "pipe",
      align = c("l", "l", "l"))
```

#### Factor Loadings

```{r plot-loadings-mil-95}
#| echo: false
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Factor loading matrix for Military × 95% (Final Design)"

plot_factor_loadings(efa_mil_95, use_clean = TRUE)
```

**Interpretation:**

```{r interpret-loadings-mil-95}
#| echo: false
#| results: asis

# Identify primary factor for each item
loadings_clean_95 <- efa_mil_95$loadings_clean
n_factors_95 <- efa_mil_95$extraction_info$nfactors

cat(sprintf("\n**Factor structure identified**: %d factors extracted\n\n", n_factors_95))

# Show which questions load on each factor
for (f in 1:n_factors_95) {
  factor_name <- colnames(loadings_clean_95)[f]
  
  # Get items with meaningful loadings
  meaningful_items <- which(abs(loadings_clean_95[[factor_name]]) >= 0.40)
  
  if (length(meaningful_items) > 0) {
    cat(sprintf("**%s**:\n", factor_name))
    
    for (item_idx in meaningful_items) {
      item_name <- rownames(loadings_clean_95)[item_idx]
      loading_val <- loadings_clean_95[[factor_name]][item_idx]
      
      # Get question label
      question_label <- question_labels_mil_95$Label[
        question_labels_mil_95$Item == item_name
      ]
      
      cat(sprintf("- %s (%.2f): %s\n", 
                  item_name, loading_val, question_label))
    }
    cat("\n")
  }
}

# Check for problematic items
if (!is.null(efa_mil_95$problematic_items$low_communality)) {
  cat(sprintf("\n⚠️ **%d item(s) with low communality** (< 0.25):\n",
              length(efa_mil_95$problematic_items$low_communality)))
  
  for (item in efa_mil_95$problematic_items$low_communality) {
    question_label <- question_labels_mil_95$Label[
      question_labels_mil_95$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items are poorly represented by the factor solution.\n\n")
}

if (!is.null(efa_mil_95$problematic_items$cross_loadings)) {
  cat(sprintf("\n⚠️ **%d item(s) with cross-loadings** (> 0.32 on multiple factors):\n",
              length(efa_mil_95$problematic_items$cross_loadings)))
  
  for (item in efa_mil_95$problematic_items$cross_loadings) {
    question_label <- question_labels_mil_95$Label[
      question_labels_mil_95$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items may need review for theoretical clarity.\n\n")
}
```

#### Variance Explained

```{r plot-variance-mil-95}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Variance explained - Military × 95%"

plot_variance_explained(efa_mil_95)
```

```{r variance-summary-mil-95}
#| echo: false
#| results: asis

variance_mat_95 <- efa_mil_95$variance_explained
total_var_95 <- variance_mat_95["Cumulative Var", ncol(variance_mat_95)] * 100

cat(sprintf("**Total variance explained**: %.1f%%\n\n", total_var_95))

if (total_var_95 >= 60) {
  cat("✅ **Good**: Factor solution explains substantial variance\n\n")
} else if (total_var_95 >= 50) {
  cat("✓ **Adequate**: Factor solution explains acceptable variance\n\n")
} else {
  cat("⚠️ **Poor**: Factor solution explains limited variance\n\n")
}
```

---

### Factor Analysis: Civil Works × 15% (Project Initiation)

Civil works early phase context (n=38).

```{r efa-cw-15}
#| echo: false
#| message: false
#| warning: false

efa_cw_15 <- run_efa(
  responses_df,
  program_name = "Civil Works",
  milestone_name = "15% (Project Initiation)",
  nfactors = NULL,
  fm = "minres",
  rotate = "oblimin"
)
```

#### Questions Analyzed

```{r question-labels-cw-15}
#| echo: false
#| tbl-cap: "Questions included in factor analysis - Civil Works × 15%"

question_labels_cw_15 <- create_question_label_table(efa_cw_15, responses_df)

kable(question_labels_cw_15,
      col.names = c("Item", "Question", "Indicator"),
      format = "pipe",
      align = c("l", "l", "l"))
```

#### Factor Loadings

```{r plot-loadings-cw-15}
#| echo: false
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Factor loading matrix for Civil Works × 15% (Project Initiation)"

plot_factor_loadings(efa_cw_15, use_clean = TRUE)
```

**Interpretation:**

```{r interpret-loadings-cw-15}
#| echo: false
#| results: asis

# Identify primary factor for each item
loadings_clean_cw <- efa_cw_15$loadings_clean
n_factors_cw <- efa_cw_15$extraction_info$nfactors

cat(sprintf("\n**Factor structure identified**: %d factors extracted\n\n", n_factors_cw))

# Show which questions load on each factor
for (f in 1:n_factors_cw) {
  factor_name <- colnames(loadings_clean_cw)[f]
  
  # Get items with meaningful loadings
  meaningful_items <- which(abs(loadings_clean_cw[[factor_name]]) >= 0.40)
  
  if (length(meaningful_items) > 0) {
    cat(sprintf("**%s**:\n", factor_name))
    
    for (item_idx in meaningful_items) {
      item_name <- rownames(loadings_clean_cw)[item_idx]
      loading_val <- loadings_clean_cw[[factor_name]][item_idx]
      
      # Get question label
      question_label <- question_labels_cw_15$Label[
        question_labels_cw_15$Item == item_name
      ]
      
      cat(sprintf("- %s (%.2f): %s\n", 
                  item_name, loading_val, question_label))
    }
    cat("\n")
  }
}

# Check for problematic items
if (!is.null(efa_cw_15$problematic_items$low_communality)) {
  cat(sprintf("\n⚠️ **%d item(s) with low communality** (< 0.25):\n",
              length(efa_cw_15$problematic_items$low_communality)))
  
  for (item in efa_cw_15$problematic_items$low_communality) {
    question_label <- question_labels_cw_15$Label[
      question_labels_cw_15$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items are poorly represented by the factor solution.\n\n")
}

if (!is.null(efa_cw_15$problematic_items$cross_loadings)) {
  cat(sprintf("\n⚠️ **%d item(s) with cross-loadings** (> 0.32 on multiple factors):\n",
              length(efa_cw_15$problematic_items$cross_loadings)))
  
  for (item in efa_cw_15$problematic_items$cross_loadings) {
    question_label <- question_labels_cw_15$Label[
      question_labels_cw_15$Item == item
    ]
    cat(sprintf("- %s: %s\n", item, question_label))
  }
  cat("\nThese items may need review for theoretical clarity.\n\n")
}
```

#### Variance Explained

```{r plot-variance-cw-15}
#| echo: false
#| fig-width: 10
#| fig-height: 5
#| fig-cap: "Variance explained - Civil Works × 15%"

plot_variance_explained(efa_cw_15)
```

```{r variance-summary-cw-15}
#| echo: false
#| results: asis

variance_mat_cw <- efa_cw_15$variance_explained
total_var_cw <- variance_mat_cw["Cumulative Var", ncol(variance_mat_cw)] * 100

cat(sprintf("**Total variance explained**: %.1f%%\n\n", total_var_cw))

if (total_var_cw >= 60) {
  cat("✅ **Good**: Factor solution explains substantial variance\n\n")
} else if (total_var_cw >= 50) {
  cat("✓ **Adequate**: Factor solution explains acceptable variance\n\n")
} else {
  cat("⚠️ **Poor**: Factor solution explains limited variance\n\n")
}
```

---

## Summary & Next Steps

### Key Findings

::: {.callout-important}
## Factor Analysis Summary

**Factor Structure:**

```{r final-summary}
#| echo: false
#| results: asis

cat(sprintf("- Parallel analysis suggests **%d-%d factors** across contexts\n",
            n_factors_range[1], n_factors_range[2]))
cat(sprintf("- Factors explain **%.1f%%** of variance on average\n", avg_variance))
cat(sprintf("- Theoretical expectation: **%d factors** (7 indicators)\n\n", theoretical_factors))

cat("**Model Fit:**\n\n")
cat("- Most contexts achieve acceptable fit per Hu & Bentler (1999) criteria\n")
cat("- Small sample sizes (N < 100) limit stability of estimates\n\n")

cat("**Implications:**\n\n")
cat("1. The questionnaire **does capture multiple dimensions** of quality\n")
cat("2. Factor structure is **reasonably consistent** across contexts\n")
cat("3. Some indicators may be **more strongly correlated** than theory suggests\n")
cat("4. Context-specific items create **variation in dimensionality**\n\n")
```
:::

### Recommendations for Questionnaire Refinement

Based on the exploratory factor analysis:

1. **Items to retain**: Questions with strong, consistent loadings across contexts

2. **Items to review**: Questions with:
   - Low communalities (< 0.25)
   - Inconsistent cross-loadings
   - Weak loadings across all factors

3. **Theoretical alignment**: Consider whether observed factor structure suggests:
   - Combining highly correlated indicators
   - Separating complex constructs
   - Adding items to underrepresented factors

4. **Sample size**: Larger samples (N > 100) per context would provide more stable factor solutions

### Future Analyses

**Immediate:**
1. Review factor loadings with subject matter experts
2. Map factors to theoretical indicators
3. Assess whether weighted scoring (based on factor loadings) improves indicator validity

**Phase 2 (per workplan):**
1. Item-total correlations and sensitivity analysis
2. Comparison of weighting schemes
3. Question redundancy analysis

**Phase 4 (advanced):**
1. Measurement invariance testing across contexts
2. Differential item functioning (DIF) analysis
3. Confirmatory factor analysis with larger sample

---

## References

- Costello, A. B., & Osborne, J. (2005). Best practices in exploratory factor analysis: Four recommendations for getting the most from your analysis. *Practical Assessment, Research, and Evaluation*, 10(7), 1-9.

- Fabrigar, L. R., Wegener, D. T., MacCallum, R. C., & Strahan, E. J. (1999). Evaluating the use of exploratory factor analysis in psychological research. *Psychological Methods*, 4(3), 272-299. https://doi.org/10.1037/1082-989X.4.3.272

- Hu, L., & Bentler, P. M. (1999). Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. *Structural Equation Modeling*, 6(1), 1-55. https://doi.org/10.1080/10705519909540118

- Little, R. J. A., & Rubin, D. B. (2019). *Statistical Analysis with Missing Data* (3rd ed.). Wiley.

- Revelle, W., & Zinbarg, R. E. (2009). Coefficients alpha, beta, omega, and the glb: Comments on Sijtsma. *Psychometrika*, 74(1), 145-154. https://doi.org/10.1007/s11336-008-9102-z

---