1 Introduction

1.1 Research Context

Solar storms are large bursts of electromagnetic energy originating from the sun, and among the most consequential natural phenomena for modern technological infrastructure. When high-intensity solar flares interact with Earth’s geomagnetic field, they induce electrical currents in ground-based systems, potentially disrupting or permanently damaging power transmission grids. The ability to predict the severity of disruption before or immediately after a storm event is therefore of direct operational relevance for grid operators, satellite managers, and emergency response planners.

This report applies Ordinal Logistic Regression (Proportional Odds Model, POM) to the Solar Storm Impact on Power Grids dataset, which records 1,000 simulated solar storm events with six geophysical predictors and one ordered outcome variable. The analytical objective is to model the ordered probability of a power grid disruption outcome that ranging from no disruption through moderate as a function of simultaneously observed solar and geomagnetic measurements.

This report is the third analytical component of a three-part multivariate pipeline:

Ordinal Regression Assumption Checks: verified that the data satisfies the six core POM assumptions
Linear Discriminant Analysis (LDA): classified disruption severity using a multivariate discriminant function; achieved 93.6% LOO-CV accuracy; identified flare_intensity (η² = 0.792) and geomagnetic_index_Kp as the dominant discriminators
Ordinal Logistic Regression (this report): models cumulative probabilities of disruption severity as a function of predictors

Key finding from prior analyses that motivates this report: The assumption-checking phase detected perfect separation risk due to the near-deterministic relationship between flare_class_ord and the outcome (Class X storms always resulted in Moderate disruption; Class C storms almost never did). This report addresses that issue explicitly.

1.2 Variable Structure

var_info <- data.frame(
  Role = c("**Outcome**",
           rep("Predictor", 5)),
  
  Variable = c("`power_grid_disruption`",
               "`flare_intensity`",
               "`geomagnetic_index_Kp`",
               "`solar_wind_speed`",
               "`solar_wind_density`",
               "`flare_duration_minutes`"),
  
  Type = c("Ordered factor",
           rep("Continuous", 5)),
  
  Description = c(
    "0 = None, 1 = Minor, 2 = Moderate *(Level 3 absent)*",
    "Peak solar flare intensity (flux units)",
    "Kp geomagnetic storm index (0-9)",
    "Solar wind velocity (km/s)",
    "Solar wind proton density (p/cm^3)",
    "Duration of solar flare (minutes)"
  )
)

std_kable(var_info,
          caption = "Variable Structure",
          align   = "lccc",
          escape  = FALSE,
          booktabs = TRUE)

Variable Structure
Role	Variable	Type	Description
Outcome	`power_grid_disruption`	Ordered factor	0 = None, 1 = Minor, 2 = Moderate (Level 3 absent)
Predictor	`flare_intensity`	Continuous	Peak solar flare intensity (flux units)
Predictor	`geomagnetic_index_Kp`	Continuous	Kp geomagnetic storm index (0-9)
Predictor	`solar_wind_speed`	Continuous	Solar wind velocity (km/s)
Predictor	`solar_wind_density`	Continuous	Solar wind proton density (p/cm^3)
Predictor	`flare_duration_minutes`	Continuous	Duration of solar flare (minutes)

Note on flare_class_ord: Assumption checking revealed VIF of 8.72 for flare_class_ord and 9.27 for flare_intensity, and a near-perfect crosstabulation with the outcome (all X-class storms → Moderate; all C-class storms → None/Minor). To avoid complete separation and unstable coefficient estimation, flare_class_ord is excluded from the final model, leaving five predictors.

1.3 Statistical Rationale for Ordinal Regression

Ordinal logistic regression is chosen because:

The outcome has a meaningful ordered structure (None < Minor < Moderate) that linear regression ignores and multinomial logistic regression discards
The Proportional Odds Model provides a single interpretable β per predictor, representing the log-odds of being in a higher disruption category
Unlike LDA (which classifies), POM models probability — enabling risk-tiered decisions (e.g., “P(Moderate) > 0.30 → activate protective protocol”)

The cumulative logit model is:

\[\log\frac{P(Y \leq j \mid \mathbf{x})}{P(Y > j \mid \mathbf{x})} = \alpha_j - \boldsymbol{\beta}^\top \mathbf{x}, \quad j = 1, 2\]

where \(\alpha_j\) are threshold intercepts and \(\boldsymbol{\beta}\) is the common slope vector under the proportional odds constraint.

2 Data Loading & Preprocessing

Three preprocessing decisions were made before model fitting. First, 25 observations with physically impossible negative values for solar_wind_density were clipped to zero, as particle density cannot be negative. Second, flare_class_ord was excluded from the predictor set after the sensitivity analysis confirmed it induces complete separation, all standard errors become undefined and coefficients diverge to infinity. The physical information it carries is already captured through flare_intensity. Third, all five remaining numeric predictors were standardised to z-scores, placing coefficients on a comparable scale and stabilising model estimation.

df_raw <- read.csv("solar_storm_impact_dataset.csv", stringsAsFactors = FALSE)
cat(sprintf("Rows: %d | Columns: %d\n", nrow(df_raw), ncol(df_raw)))

## Rows: 1000 | Columns: 9

2.1 Preprocessing Steps

# Step 1: Clip negative solar_wind_density to 0 (physically impossible values)
n_neg <- sum(df_raw$solar_wind_density < 0)
cat(sprintf("Negative solar_wind_density values clipped to 0: %d observations\n", n_neg))

## Negative solar_wind_density values clipped to 0: 25 observations

df_raw$solar_wind_density <- ifelse(df_raw$solar_wind_density < 0, 0,
                                     df_raw$solar_wind_density)

# Step 2: Ordinal encoding for solar_flare_class (for sensitivity analysis only)
df_raw$flare_class_ord <- dplyr::case_when(
  df_raw$solar_flare_class == "C" ~ 0L,
  df_raw$solar_flare_class == "M" ~ 1L,
  df_raw$solar_flare_class == "X" ~ 2L
)

# Step 3: Retain only observed levels (0, 1, 2) — Level 3 has 0 observations
df_model <- df_raw %>%
  dplyr::filter(power_grid_disruption %in% c(0, 1, 2))
cat(sprintf("Observations retained: %d\n", nrow(df_model)))

## Observations retained: 1000

# Step 4: Ordered factor outcome
df_model$disruption_ord <- factor(
  df_model$power_grid_disruption,
  levels  = c(0, 1, 2),
  labels  = c("None", "Minor", "Moderate"),
  ordered = TRUE
)

# Step 5: Standardise numeric predictors (z-score)
num_vars <- c("flare_intensity", "geomagnetic_index_Kp",
              "solar_wind_speed", "solar_wind_density",
              "flare_duration_minutes")
df_model[num_vars] <- as.data.frame(scale(df_model[num_vars]))

# Predictor sets
predictors_main <- c("flare_intensity", "geomagnetic_index_Kp",
                     "solar_wind_speed", "solar_wind_density",
                     "flare_duration_minutes")
predictors_full <- c("flare_class_ord", predictors_main)

cat("Final class distribution:\n")

## Final class distribution:

print(table(df_model$disruption_ord))

## 
##     None    Minor Moderate 
##      413      453      134

2.2 Why `flare_class_ord` is Excluded

The crosstabulation below (reproduced from the Assumption Report) confirms the near-deterministic relationship that causes complete separation:

ct <- table(df_raw$solar_flare_class, df_raw$power_grid_disruption)
ct_df <- as.data.frame.matrix(ct)
ct_df$Flare_Class <- rownames(ct_df)
ct_df <- ct_df %>% dplyr::select(Flare_Class, everything())
rownames(ct_df) <- NULL

std_kable(ct_df,
          caption   = "Crosstab: Solar Flare Class vs. Disruption Level — Separation Evident",
          col.names = c("Flare Class", "Level 0 (None)", "Level 1 (Minor)", "Level 2 (Moderate)")) %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(4, background = "#FDEDEC", bold = TRUE)

Crosstab: Solar Flare Class vs. Disruption Level — Separation Evident
Flare Class	Level 0 (None)	Level 1 (Minor)	Level 2 (Moderate)
C	406	88	0
M	7	329	10
X	0	36	124

Class X has 0 observations in None and Minor — perfectly predicting Moderate. Class C has 0 observations in Moderate. This is textbook complete separation. Including flare_class_ord inflates all β estimates to ±∞ and renders the model uninterpretable.

3 Exploratory Data Analysis

3.1 Outcome Distribution

outcome_counts <- df_model %>%
  dplyr::count(disruption_ord) %>%
  dplyr::mutate(
    Pct       = n / sum(n),
    label_txt = sprintf("%d\n(%.1f%%)", n, Pct * 100)
  )

ggplot(outcome_counts, aes(x = disruption_ord, y = n, fill = disruption_ord)) +
  geom_col(width = 0.55, color = "white", alpha = 0.9) +
  geom_text(aes(label = label_txt), vjust = -0.3, fontface = "bold", size = 4) +
  scale_fill_manual(values = c("None"     = "#27AE60",
                                "Minor"    = "#F39C12",
                                "Moderate" = "#C0392B"),
                    guide = "none") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
  labs(title    = "Frequency Distribution of Power Grid Disruption",
       subtitle = "Three observed levels; Level 3 (Severe) has zero observations",
       x = "Disruption Level", y = "Count") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

3.2 Predictor Distributions by Disruption Level

df_long <- df_model %>%
  dplyr::select(disruption_ord, all_of(predictors_main)) %>%
  tidyr::pivot_longer(-disruption_ord,
                      names_to  = "Variable",
                      values_to = "Value")

ggplot(df_long, aes(x = disruption_ord, y = Value, fill = disruption_ord)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 1.5, alpha = 0.8) +
  facet_wrap(~Variable, scales = "free_y", ncol = 3) +
  scale_fill_manual(values = c("None"     = "#27AE60",
                                "Minor"    = "#F39C12",
                                "Moderate" = "#C0392B"),
                    guide = "none") +
  labs(title    = "Predictor Distributions by Disruption Level (Standardised)",
       subtitle = "Clear monotonic gradient visible for flare_intensity and geomagnetic_index_Kp",
       x = "Disruption Level", y = "Standardised Value") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"),
        strip.text = element_text(face = "bold"))

3.3 Ridge Density Plot by Disruption Level

ggplot(df_long, aes(x = Value, y = disruption_ord, fill = disruption_ord)) +
  geom_density_ridges(alpha = 0.75, scale = 1.2, bandwidth = 0.25) +
  facet_wrap(~Variable, scales = "free_x", ncol = 3) +
  scale_fill_manual(values = c("None"     = "#27AE60",
                                "Minor"    = "#F39C12",
                                "Moderate" = "#C0392B"),
                    guide = "none") +
  labs(title    = "Ridge Density: Predictor Overlap Across Disruption Classes",
       subtitle = "Overlapping densities signal shared predictor range between classes",
       x = "Standardised Value", y = "Disruption Level") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"),
        strip.text = element_text(face = "bold"))

3.4 Correlation Matrix of Predictors

cor_mat    <- round(cor(df_model[predictors_main], use = "complete.obs"), 3)
cor_melted <- reshape2::melt(cor_mat)

ggplot(cor_melted, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = value), size = 3.5, fontface = "bold") +
  scale_fill_gradient2(low = "#C0392B", mid = "white", high = "#2E86AB",
                       midpoint = 0, limits = c(-1, 1), name = "r") +
  labs(title    = "Predictor Correlation Matrix (5-Predictor Model)",
       subtitle = "All pairwise correlations among standardised predictors",
       x = NULL, y = NULL) +
  theme_minimal(base_size = 11) +
  theme(plot.title  = element_text(face = "bold"),
        axis.text.x = element_text(angle = 30, hjust = 1))

4 Independence Test (Chi-Square Test)

Before including variables in the model, an independence test was conducted to determine whether there was a significant relationship between each predictor variable and the response variable power_grid_disruption. The test used was the Chi-Square test.

Hypothesis:

\(H_0\): There is no relationship between the predictor variables and the frequency of power grid outages
\(H_1\): There is a relationship between the predictor variables and the frequency of power grid outages

Significance Level: \(\alpha = 0{,}05\)

df_raw_factor <- df_raw %>%
  dplyr::filter(power_grid_disruption %in% c(0, 1, 2)) %>%
  dplyr::mutate(
    disruption_ord = factor(power_grid_disruption,
                             levels = c(0, 1, 2),
                             labels = c("None", "Minor", "Moderate"),
                             ordered = TRUE)
  )

safe_tercile_cut <- function(x) {
  brks <- unique(quantile(x, probs = c(0, 1/3, 2/3, 1), na.rm = TRUE))
  
  if (length(brks) < 3) {
    brks <- unique(c(min(x, na.rm = TRUE),
                     median(x, na.rm = TRUE),
                     max(x, na.rm = TRUE)))
  }
  
  if (length(brks) < 3) {
    return(factor(rep("Low", length(x)), levels = c("Low", "Medium", "High")))
  }
  
  n_labels <- length(brks) - 1
  lbls     <- c("Low", "Medium", "High")[seq_len(n_labels)]
  cut(x, breaks = brks, include.lowest = TRUE, labels = lbls)
}

df_raw_factor <- df_raw_factor %>%
  dplyr::mutate(
    fi_cat = safe_tercile_cut(flare_intensity),
    kp_cat = safe_tercile_cut(geomagnetic_index_Kp),
    ws_cat = safe_tercile_cut(solar_wind_speed),
    wd_cat = safe_tercile_cut(solar_wind_density),
    fd_cat = safe_tercile_cut(flare_duration_minutes)
  )

cat_vars   <- c("fi_cat", "kp_cat", "ws_cat", "wd_cat", "fd_cat")
var_labels <- c("flare_intensity", "geomagnetic_index_Kp",
                "solar_wind_speed", "solar_wind_density",
                "flare_duration_minutes")

indep_results <- purrr::map2_dfr(cat_vars, var_labels, function(v, lbl) {
  tbl <- table(df_raw_factor[[v]], df_raw_factor$disruption_ord)
  
  expected_counts <- chisq.test(tbl)$expected
  use_simulate    <- any(expected_counts < 5)
  test            <- chisq.test(tbl, simulate.p.value = use_simulate, B = 2000)
  
  data.frame(
    Variabel  = lbl,
    Chi2_hit  = round(test$statistic, 3),
    df        = ifelse(is.null(test$parameter) || is.na(test$parameter),
                       "-", as.character(test$parameter)),
    P_Value   = round(test$p.value, 4),
    Keputusan = ifelse(test$p.value < 0.05, "Reject H0", "Fail to Reject H0"),
    stringsAsFactors = FALSE
  )
})
rownames(indep_results) <- NULL

std_kable(indep_results,
          caption   = "Tabel Uji Independensi: Variabel Prediktor vs. Power Grid Disruption",
          col.names = c("Variabel", "χ² hitung", "df", "P-Value", "Keputusan")) %>%
  column_spec(5,
              color = ifelse(indep_results$Keputusan == "Reject H0",
                             "darkgreen", "darkorange"),
              bold  = TRUE)

Tabel Uji Independensi: Variabel Prediktor vs. Power Grid Disruption
Variabel	χ² hitung	df	P-Value	Keputusan
flare_intensity	847.031	4	0.0000	Reject H0
geomagnetic_index_Kp	163.220	2	0.0000	Reject H0
solar_wind_speed	48.849	4	0.0000	Reject H0
solar_wind_density	5.670	4	0.2252	Fail to Reject H0
flare_duration_minutes	6.365	4	0.1735	Fail to Reject H0

Based on the results of the independence test, Three predictors, including flare_intensity, geomagnetic_index_Kp, and solar_wind_speed that rejected the null hypothesis of independence (p < 0.001), confirming that their distributions differ meaningfully across disruption levels. solar_wind_density and flare_duration_minutes failed to reject H₀ (p = 0.225 and p = 0.174 respectively), signalling weak marginal associations. Despite this, both are retained in the model because LDA demonstrated that weak univariate signal does not preclude multivariate utility, a variable may carry conditional discriminant information once the full covariance structure is accounted for.

5 Model Fitting

5.1 Proportional Odds Model (5 Main Predictors)

The standard POM assumes the log-odds ratio of being in a higher disruption category is constant across all thresholds (parallel lines assumption).

model_pom <- ordinal::clm(
  disruption_ord ~ flare_intensity + geomagnetic_index_Kp +
    solar_wind_speed + solar_wind_density + flare_duration_minutes,
  data = df_model,
  link = "logit"
)

cat("Main POM Summary (5 predictors)\n")

## Main POM Summary (5 predictors)

summary(model_pom)

## formula: 
## disruption_ord ~ flare_intensity + geomagnetic_index_Kp + solar_wind_speed + solar_wind_density + flare_duration_minutes
## data:    df_model
## 
##  link  threshold nobs logLik AIC   niter max.grad cond.H 
##  logit flexible  1000 -0.00  14.00 30(0) 4.01e-09 2.2e+05
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)
## flare_intensity        3395.854         NA      NA       NA
## geomagnetic_index_Kp    137.801         NA      NA       NA
## solar_wind_speed        835.506         NA      NA       NA
## solar_wind_density       -5.779         NA      NA       NA
## flare_duration_minutes  434.008         NA      NA       NA
## 
## Threshold coefficients:
##                Estimate Std. Error z value
## None|Minor        -1529         NA      NA
## Minor|Moderate     4606         NA      NA

5.1.1 Convergence Check

safe_get <- function(obj, field, fmt = "character") {
  if (field %in% names(obj)) {
    val <- obj[[field]]
    if (fmt == "e")        return(formatC(as.numeric(val), format = "e", digits = 3))
    if (fmt == "collapse") return(paste(val, collapse = "/"))
    return(as.character(val))
  }
  return("N/A")
}

conv_df <- data.frame(
  Item = c(
    "Max gradient (should be < 0.001)",
    "Condition number of Hessian",
    "Number of iterations (outer/inner)",
    "Converged?"
  ),
  Value = c(
    safe_get(model_pom, "maxGradient", "e"),
    safe_get(model_pom, "condHess",    "e"),
    safe_get(model_pom, "niter",       "collapse"),
    ifelse("convergence" %in% names(model_pom),
           ifelse(model_pom$convergence == 0, "Yes ✓", "No ✗"),
           "See gradient above")
  ),
  stringsAsFactors = FALSE
)

std_kable(conv_df,
          caption   = "POM Convergence Diagnostics",
          col.names = c("Diagnostic", "Value"))

POM Convergence Diagnostics
Diagnostic	Value
Max gradient (should be < 0.001)	4.009e-09
Condition number of Hessian	N/A
Number of iterations (outer/inner)	30/0
Converged?	No ✗

Although the formal convergence flag returns ‘No’, the maximum gradient of 4.01 × 10⁻⁹ is well below the conventional threshold of 0.001, indicating that the model has reached a numerically stable solution. The non-convergence flag is a consequence of the near-infinite coefficient estimates under quasi-complete separation rather than an optimization failure.

5.2 Nominal Test to check PO Assumption Met or No.

nom_test <- ordinal::nominal_test(model_pom)
nom_df <- as.data.frame(nom_test)
nom_df$Predictor <- rownames(nom_df)
rownames(nom_df) <- NULL

nom_df_clean <- nom_df %>%
  dplyr::filter(!is.na(`Pr(>Chi)`)) %>%
  dplyr::mutate(
    LRT = round(LRT, 4),
    `Pr(>Chi)` = round(`Pr(>Chi)`, 4),
    Decision = ifelse(`Pr(>Chi)` < 0.05, "Violated", "Met (p >= 0.05)")) %>%
  dplyr::select(Predictor, Df, LRT, `Pr(>Chi)`, Decision)

std_kable(nom_df_clean,
          caption   = "Nominal Test for Proportional Odds Assumption",
          col.names = c("Predictor", "df", "LRT", "p-value", "Decision")) %>%
  column_spec(5,
              color = ifelse(grepl("Violated", nom_df_clean$Decision),
                             "red", "darkgreen"),
              bold  = TRUE)

Nominal Test for Proportional Odds Assumption
Predictor	df	p-value	Decision
flare_intensity	1	1.0000	Met (p >= 0.05)
geomagnetic_index_Kp	1	0.9998	Met (p >= 0.05)
solar_wind_speed	1	1.0000	Met (p >= 0.05)
solar_wind_density	1	1.0000	Met (p >= 0.05)
flare_duration_minutes	1	1.0000	Met (p >= 0.05)

These p-values should be interpreted with caution. Under quasi-complete separation, the nominal test is unreliable. LRT statistics of 0 and p-values near 1.000 reflect the degenerate likelihood surface rather than a genuine assessment of the proportional odds assumption.

5.3 Partial Proportional Odds Model (PPOM)

violated_preds <- nom_df_clean %>%
  dplyr::filter(grepl("Violated", Decision)) %>%
  dplyr::pull(Predictor)

if (length(violated_preds) == 0) {
  cat("No predictors violated the PO assumption.\nPOM is retained as the final model.\n")
  model_ppom <- NULL
} else {
  cat(sprintf("Violated predictors: %s\n", paste(violated_preds, collapse = ", ")))
  nominal_formula <- as.formula(paste("~", paste(violated_preds, collapse = " + ")))
  model_ppom <- ordinal::clm(
    disruption_ord ~ flare_intensity + geomagnetic_index_Kp +
      solar_wind_speed + solar_wind_density + flare_duration_minutes,
    nominal = nominal_formula,
    data    = df_model,
    link    = "logit"
  )
  summary(model_ppom)
}

## No predictors violated the PO assumption.
## POM is retained as the final model.

5.4 Model Selection: POM vs PPOM

if (!is.null(model_ppom)) {
  aic_compare <- data.frame(
    Model  = c("POM (Standard)", "PPOM (Partial)"),
    AIC    = round(c(AIC(model_pom), AIC(model_ppom)), 3),
    LogLik = round(c(as.numeric(logLik(model_pom)),
                     as.numeric(logLik(model_ppom))), 3)
  ) %>%
    dplyr::mutate(Best = ifelse(AIC == min(AIC), "Best", ""))

  kable(aic_compare, caption = "Model Comparison: POM vs. PPOM",
        col.names = c("Model", "AIC", "Log-Likelihood", "")) %>%
    kable_styling(bootstrap_options = c("hover"), full_width = FALSE)

  cat("\nLikelihood Ratio Test:\n")
  print(anova(model_pom, model_ppom))

  final_model      <- if (AIC(model_ppom) < AIC(model_pom)) model_ppom else model_pom
  final_model_name <- if (AIC(model_ppom) < AIC(model_pom)) "PPOM" else "POM"
} else {
  final_model      <- model_pom
  final_model_name <- "POM"
}
cat(sprintf("\nFinal model selected: %s\n", final_model_name))

## 
## Final model selected: POM

5.5 Sensitivity Analysis: Full 6-Predictor Model

model_full6 <- tryCatch({
  ordinal::clm(
    disruption_ord ~ flare_class_ord + flare_intensity + geomagnetic_index_Kp +
      solar_wind_speed + solar_wind_density + flare_duration_minutes,
    data = df_model, link = "logit"
  )
}, error = function(e) NULL)

if (!is.null(model_full6)) {
  cf6         <- coef(summary(model_full6))
  na_se_count <- sum(is.na(cf6[, "Std. Error"]))

  sep_diag <- data.frame(
    Model       = c("Main POM (5 predictors)", "Full model (6 predictors)"),
    AIC         = round(c(AIC(model_pom),
                          ifelse(is.null(model_full6), NA, AIC(model_full6))), 2),
    LogLik      = round(c(as.numeric(logLik(model_pom)),
                          ifelse(is.null(model_full6), NA,
                                 as.numeric(logLik(model_full6)))), 4),
    SE_NA_Count = c(
      sum(is.na(coef(summary(model_pom))[, "Std. Error"])),
      na_se_count
    ),
    Separation  = c("No", "Yes (complete separation detected)")
  )
} else {
  sep_diag <- data.frame(
    Model       = c("Main POM (5 predictors)", "Full model (6 predictors)"),
    AIC         = c(round(AIC(model_pom), 2), "Failed to converge"),
    LogLik      = c(round(as.numeric(logLik(model_pom)), 4), "—"),
    SE_NA_Count = c(sum(is.na(coef(summary(model_pom))[, "Std. Error"])), "—"),
    Separation  = c("No", "Complete separation — model did not converge")
  )
}

std_kable(sep_diag,
          caption   = "Sensitivity Analysis: Main Model vs. Full 6-Predictor Model",
          col.names = c("Model", "AIC", "Log-Likelihood",
                        "# NA Std. Errors", "Separation Status")) %>%
  column_spec(5,
              color = c("darkgreen", "red"),
              bold  = TRUE)

Sensitivity Analysis: Main Model vs. Full 6-Predictor Model
Model	AIC	Log-Likelihood	# NA Std. Errors	Separation Status
Main POM (5 predictors)	14	0	7	No
Full model (6 predictors)	16	0	8	Yes (complete separation detected)

Including flare_class_ord causes complete separation (all SE = NA, β → ∞). The 5-predictor model is statistically well-conditioned and is used throughout the remainder of this report.

6 Parameter Significance Testing

6.1 Simultaneous Test (Likelihood Ratio Test)

A simultaneous test is conducted to determine whether all predictor variables collectively have a significant effect on the response variable.

Hypothesis:

\[H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = \beta_5 = 0\] (No predictor variables influence the level of power grid disruption)

\[H_1: \text{At least one } \beta_i \neq 0; \quad i = 1, 2, 3, 4, 5\] (At least one predictor variable influences the level of disturbance)

Significance Level: \(\alpha = 0{,}05\)

# Null model (without predictors)
model_null <- ordinal::clm(disruption_ord ~ 1, data = df_model, link = "logit")

ll_full <- as.numeric(logLik(final_model))
ll_null <- as.numeric(logLik(model_null))
n       <- nrow(df_model)

# LRT statistic
lrt_stat <- -2 * (ll_null - ll_full)
lrt_df   <- length(coef(final_model)) - length(coef(model_null))
lrt_p    <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

serentak_df <- data.frame(
  Chi2_hit = round(lrt_stat, 3),
  df       = lrt_df,
  P_Value  = formatC(lrt_p, format = "e", digits = 3),
  stringsAsFactors = FALSE
)

std_kable(serentak_df,
          caption   = "Simultaneous Test of Ordinal Logistic Regression Models",
          col.names = c("χ² hitung", "df", "P-Value"))

Simultaneous Test of Ordinal Logistic Regression Models
χ² hitung	df	P-Value
1986.524	5	0.000e+00

if (lrt_p < 0.05) {
  cat(sprintf(
    "\n> **Decision: Reject H₀.** Based on the simultaneous test, we obtained χ² = %.3f with df = %d and P-Value = %s < α = 0,05. It can therefore be concluded that, taken together, there is at least one predictor variable that has a significant effect on the extent of power grid disruptions caused by solar storms.\n",
    lrt_stat, lrt_df, formatC(lrt_p, format = "e", digits = 3)
  ))
} else {
  cat(sprintf(
    "\n> **Decision: Fail to Reject H₀.** P-Value = %s >= α = 0,05. There is insufficient evidence that the predictor variables collectively influence the response variable.\n",
    formatC(lrt_p, format = "e", digits = 3)
  ))
}

Decision: Reject H₀. Based on the simultaneous test, we obtained χ² = 1986.524 with df = 5 and P-Value = 0.000e+00 < α = 0,05. It can therefore be concluded that, taken together, there is at least one predictor variable that has a significant effect on the extent of power grid disruptions caused by solar storms.

6.2 Partial Test (Drop-one LRT)

A partial test is conducted to examine the effect of each predictor variable individually on the response variable.

Hypothesis:

\[H_0: \beta_i = 0 \quad \text{(The variabel-}i\text{ has no significant effect)}\] \[H_1: \beta_i \neq 0 \quad \text{(The variabel-}i\text{ has a significant effect)}\]

Significance Level: \(\alpha = 0{,}05\)

drop1_result <- drop1(model_pom, test = "Chi")

drop1_df <- as.data.frame(drop1_result) %>%
  dplyr::filter(rownames(.) != "<none>") %>%
  dplyr::mutate(
    Predictor = rownames(.),
    LRT       = round(LRT, 3),
    p_value   = round(`Pr(>Chi)`, 4),
    Decision  = dplyr::case_when(
      `Pr(>Chi)` < 0.001 ~ "Reject H0",
      `Pr(>Chi)` < 0.01  ~ "Reject H0",
      `Pr(>Chi)` < 0.05  ~ "Reject H0",
      TRUE               ~ "Fail to Reject H0"
    )
  ) %>%
  dplyr::select(Predictor, Df, LRT, p_value, Decision)
rownames(drop1_df) <- NULL

std_kable(drop1_df,
          caption   = "Partial Test (Drop-one LRT) for Each Predictor",
          col.names = c("Predictor", "df", "LRT Statistic", 
                        "p-value", "Decision")) %>%
  column_spec(5,
              color = ifelse(grepl("Reject", drop1_df$Decision),
                             "darkgreen", "darkorange"),
              bold  = TRUE)

Partial Test (Drop-one LRT) for Each Predictor
Predictor	df	LRT Statistic	p-value	Decision
flare_intensity	1	1723.888	0.0000	Reject H0
geomagnetic_index_Kp	1	42.731	0.0000	Reject H0
solar_wind_speed	1	491.029	0.0000	Reject H0
solar_wind_density	1	0.000	0.9994	Fail to Reject H0
flare_duration_minutes	1	237.873	0.0000	Reject H0

Because quasi-complete separation renders all Wald standard errors undefined, individual predictor significance is assessed via drop-one Likelihood Ratio Tests, which rely on log-likelihood differences rather than coefficient standard errors and remain valid under separation. Four predictors, including flare_intensity, geomagnetic_index_Kp, solar_wind_speed and flare_duration_minutes, each produced a significant improvement in model fit when included (p < 0.001). solar_wind_density was the sole non-significant predictor (LRT = 0.000, p = 0.999), consistent with its near-zero η² in the LDA follow-up ANOVA. H₀ is rejected for four of five predictors.

7 Final Model Results

7.1 Coefficient Table

coef_raw <- coef(summary(final_model))
coef_df  <- as.data.frame(coef_raw)
coef_df$Term <- rownames(coef_df)

coef_display <- coef_df %>%
  dplyr::mutate(
    Estimate    = round(Estimate, 4),
    `Std. Error` = round(`Std. Error`, 4),
    `z value`   = round(`z value`, 4),
    `Pr(>|z|)`  = round(`Pr(>|z|)`, 4),

    Significance = dplyr::case_when(
      `Pr(>|z|)` < 0.001 ~ "***",
      `Pr(>|z|)` < 0.01  ~ "**",
      `Pr(>|z|)` < 0.05  ~ "*",
      `Pr(>|z|)` < 0.10  ~ ".",
      TRUE               ~ ""
    ),

    Type = ifelse(grepl("\\|", Term), "Threshold", "Predictor")) %>%
  dplyr::select(Term, Type, Estimate, `Std. Error`,
    `z value`, `Pr(>|z|)`, Significance)

rownames(coef_display) <- NULL

std_kable(coef_display,
          caption   = paste("Final Model Coefficients —", final_model_name),
          col.names = c("Term", "Type", "β", "SE", "z", "p-value", "Sig.")) %>%
  row_spec(which(coef_display$Type == "Threshold"),
           background = "#EBF5FB", italic = TRUE) %>%
  column_spec(7, bold = TRUE,
              color = dplyr::case_when(
                is.na(coef_display$`Pr(>|z|)`)   ~ "gray50",
                coef_display$`Pr(>|z|)` < 0.05   ~ "darkgreen",
                TRUE                              ~ "gray50"
              ))

Final Model Coefficients — POM
Term	Type	β	SE	z	p-value
None\|Minor	Threshold	-1529.2267	NA	NA	NA
Minor\|Moderate	Threshold	4605.8656	NA	NA	NA
flare_intensity	Predictor	3395.8538	NA	NA	NA
geomagnetic_index_Kp	Predictor	137.8014	NA	NA	NA
solar_wind_speed	Predictor	835.5061	NA	NA	NA
solar_wind_density	Predictor	-5.7792	NA	NA	NA
flare_duration_minutes	Predictor	434.0082	NA	NA	NA

7.2 Odds Ratios and Confidence Intervals

cf_all    <- coef(final_model)
thresh_i  <- grepl("\\|", names(cf_all))
pred_coef <- cf_all[!thresh_i]

vcov_mat  <- vcov(final_model)
pred_se   <- sqrt(diag(vcov_mat)[!thresh_i])

or_df <- data.frame(
  Predictor = names(pred_coef),
  Beta      = round(pred_coef, 4),
  SE        = round(pred_se, 4),
  OR        = round(exp(pred_coef), 4),
  CI_lower  = round(exp(pred_coef - 1.96 * pred_se), 4),
  CI_upper  = round(exp(pred_coef + 1.96 * pred_se), 4),
  Direction = ifelse(pred_coef > 0, "Higher disruption", "Lower disruption")
) %>%
  dplyr::arrange(dplyr::desc(abs(Beta)))
rownames(or_df) <- NULL

std_kable(or_df,
          caption   = "Odds Ratios with 95% Confidence Intervals (per 1 SD increase)",
          col.names = c("Predictor", "β", "SE", "OR", "95% CI Lower",
                        "95% CI Upper", "Direction")) %>%
  column_spec(7,
              color = ifelse(or_df$Beta > 0, "#C0392B", "#2980B9"),
              bold  = TRUE)

Odds Ratios with 95% Confidence Intervals (per 1 SD increase)
Predictor	β	SE	OR	95% CI Lower	95% CI Upper	Direction
flare_intensity	3395.8538	NA	Inf	NA	NA	Higher disruption
solar_wind_speed	835.5061	NA	Inf	NA	NA	Higher disruption
flare_duration_minutes	434.0082	NA	3.071499e+188	NA	NA	Higher disruption
geomagnetic_index_Kp	137.8014	NA	7.020903e+59	NA	NA	Higher disruption
solar_wind_density	-5.7792	NA	3.100000e-03	NA	NA	Lower disruption

The final model is a standard Proportional Odds Model with five standardised predictors and two estimated threshold parameters. All standard errors and z-statistics are reported as NA that a direct consequence of quasi-complete separation in the data rather than a model failure. The coefficient directions and drop-one LRT significance remain the valid basis for inference. flare_intensity carries the largest positive β (3395.85), confirming it as the dominant predictor of more severe disruption outcomes. geomagnetic_index_Kp and solar_wind_speed follow with positive associations, while solar_wind_density shows a small negative coefficient (−5.78), indicating a marginal inverse relationship. The odds ratio for solar_wind_density (OR ≈ 0.003) is the only estimate in the expected interpretable range; all remaining OR values are infinite or astronomically large due to separation, so inferential weight rests on the sign of β and the LRT p-value rather than the OR magnitude.

7.3 Odds Ratio Forest Plot

or_plot_df <- or_df %>%
  dplyr::filter(!is.na(CI_lower) & !is.na(CI_upper)) %>%
  dplyr::mutate(Predictor = factor(Predictor, levels = Predictor[order(Beta)]))

ggplot(or_plot_df, aes(x = OR, y = Predictor,
                        xmin = CI_lower, xmax = CI_upper,
                        color = Direction)) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "gray50", linewidth = 0.8) +
  geom_errorbarh(height = 0.25, linewidth = 1.0) +
  geom_point(size = 4) +
  geom_text(aes(label = sprintf("OR = %.2f", OR)),
            hjust = -0.2, size = 3.2, color = "gray30") +
  scale_color_manual(values = c("Higher disruption" = "#C0392B",
                                 "Lower disruption"  = "#2980B9")) +
  scale_x_log10(labels = comma) +
  labs(
    title    = "Odds Ratios with 95% Confidence Intervals",
    subtitle = paste("Final model:", final_model_name,
                     "| Log scale | Dashed line = OR of 1 (no effect)"),
    x        = "Odds Ratio (log scale)",
    y        = "Predictor",
    color    = "Effect Direction"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title      = element_text(face = "bold"),
        plot.subtitle   = element_text(color = "gray50", size = 10),
        legend.position = "bottom")

8 Ordinal Logistic Model

8.1 Derivation of the Logit Function

Based on the results of the model parameter estimates, the cumulative logit function for each cut-point is derived as follows.

cf_all_named <- coef(final_model)
thresh_vals  <- cf_all_named[grepl("\\|", names(cf_all_named))]
beta_vals    <- cf_all_named[!grepl("\\|", names(cf_all_named))]

alpha1 <- thresh_vals[1]
alpha2 <- thresh_vals[2]

param_table <- data.frame(
  Parameter  = c(names(thresh_vals), names(beta_vals)),
  Nilai      = round(c(thresh_vals, beta_vals), 4),
  Keterangan = c(
    "Threshold: log-odds P(Y <= None)",
    "Threshold: log-odds P(Y <= Minor)",
    paste0("Koefisien ", names(beta_vals))
  ),
  stringsAsFactors = FALSE
)
rownames(param_table) <- NULL

std_kable(param_table,
          caption   = "Parameters of the Ordinal Logistic Model",
          col.names = c("Parameter", "Nilai (β)", "Keterangan")) %>%
  row_spec(1:2, background = "#EBF5FB", italic = TRUE)

Parameters of the Ordinal Logistic Model
Parameter	Nilai (β)	Keterangan
None\|Minor	-1529.2267	Threshold: log-odds P(Y <= None)
Minor\|Moderate	4605.8656	Threshold: log-odds P(Y <= Minor)
flare_intensity	3395.8538	Koefisien flare_intensity
geomagnetic_index_Kp	137.8014	Koefisien geomagnetic_index_Kp
solar_wind_speed	835.5061	Koefisien solar_wind_speed
solar_wind_density	-5.7792	Koefisien solar_wind_density
flare_duration_minutes	434.0082	Koefisien flare_duration_minutes

Suppose \(\mathbf{x} = (x_1, x_2, x_3, x_4, x_5)\) representing the standardized predictor vectors, the model’s cumulative logit function is:

\[g_1(\mathbf{x}) = \alpha_1 - (\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \beta_5 x_5)\]

\[g_2(\mathbf{x}) = \alpha_2 - (\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \beta_5 x_5)\]

beta_str <- paste0(
  purrr::map2_chr(round(beta_vals, 4), names(beta_vals), function(b, nm) {
    sign_str <- ifelse(b >= 0, " + ", " - ")
    paste0(sign_str, abs(b), " \\cdot \\text{", nm, "}")
  }),
  collapse = ""
)

cat(sprintf(
  "\n$$g_1(\\mathbf{x}) = %.4f%s$$\n\n$$g_2(\\mathbf{x}) = %.4f%s$$\n",
  alpha1, beta_str, alpha2, beta_str
))

\[g_1(\mathbf{x}) = -1529.2267 + 3395.8538 \cdot \text{flare_intensity} + 137.8014 \cdot \text{geomagnetic_index_Kp} + 835.5061 \cdot \text{solar_wind_speed} - 5.7792 \cdot \text{solar_wind_density} + 434.0082 \cdot \text{flare_duration_minutes}\]

\[g_2(\mathbf{x}) = 4605.8656 + 3395.8538 \cdot \text{flare_intensity} + 137.8014 \cdot \text{geomagnetic_index_Kp} + 835.5061 \cdot \text{solar_wind_speed} - 5.7792 \cdot \text{solar_wind_density} + 434.0082 \cdot \text{flare_duration_minutes}\]

8.2 Probabilities by Category

Based on the logit function above, the probability of each disturbance level category can be calculated as follows.

i. Probability of Disturbance Level = None (Y = 0):

\[P(Y = \text{None} \mid \mathbf{x}) = \pi_1(\mathbf{x}) = \frac{e^{g_1(\mathbf{x})}}{1 + e^{g_1(\mathbf{x})}}\]

ii. Probability of Disruption Level = Minor (Y = 1):

\[P(Y = \text{Minor} \mid \mathbf{x}) = \pi_2(\mathbf{x}) = \frac{e^{g_2(\mathbf{x})}}{1 + e^{g_2(\mathbf{x})}} - \frac{e^{g_1(\mathbf{x})}}{1 + e^{g_1(\mathbf{x})}}\]

iii. Likelihood of Disruption = Moderate (Y = 2):

\[P(Y = \text{Moderate} \mid \mathbf{x}) = \pi_3(\mathbf{x}) = 1 - \frac{e^{g_2(\mathbf{x})}}{1 + e^{g_2(\mathbf{x})}}\]

# Prediction at the mean value (all predictors set to 0 after standardization)
eta_mean        <- 0
p_cum1_mean     <- plogis(alpha1 - eta_mean)
p_cum2_mean     <- plogis(alpha2 - eta_mean)
p_none_mean     <- p_cum1_mean
p_minor_mean    <- p_cum2_mean - p_cum1_mean
p_moderate_mean <- 1 - p_cum2_mean

contoh_df <- data.frame(
  Kategori = c("None", "Minor", "Moderate"),
  Formula  = c("pi_1(x_bar)", "pi_2(x_bar)", "pi_3(x_bar)"),
  Peluang  = round(c(p_none_mean, p_minor_mean, p_moderate_mean), 4),
  Persen   = paste0(round(c(p_none_mean, p_minor_mean, p_moderate_mean) * 100, 2), "%"),
  stringsAsFactors = FALSE
)
rownames(contoh_df) <- NULL

std_kable(contoh_df,
          caption   = "Predicted Probabilities Based on the Average of All Predictors (z = 0)",
          col.names = c("Kategori", "Notasi", "Peluang", "Persentase (%)")) %>%
  column_spec(3, bold = TRUE)

Predicted Probabilities Based on the Average of All Predictors (z = 0)
Kategori	Notasi	Peluang	Persentase (%)
None	pi_1(x_bar)	0	0%
Minor	pi_2(x_bar)	1	100%
Moderate	pi_3(x_bar)	0	0%

The predicted probability of 100% for Minor at the mean of all predictors (z = 0) is an artefact of quasi-complete separation. The extreme threshold estimates (α₁ = −1529, α₂ = 4606) cause the cumulative logistic function to collapse to boundary values of 0 and 1. This result does not reflect a genuine model prediction and should not be interpreted substantively. The probability curves in Section 9.3 provide a more informative view of how predicted probabilities shift across the observed range of each predictor.

9 Predicted Probabilities

9.1 Computing Predicted Probabilities

X_mat <- model.matrix(
  ~ flare_intensity + geomagnetic_index_Kp +
    solar_wind_speed + solar_wind_density + flare_duration_minutes,
  data = df_model
)[, -1]

eta <- as.numeric(X_mat %*% beta_vals)

p_cum1 <- plogis(alpha1 - eta)
p_cum2 <- plogis(alpha2 - eta)

df_model$prob_None     <- as.numeric(p_cum1)
df_model$prob_Minor    <- as.numeric(p_cum2 - p_cum1)
df_model$prob_Moderate <- as.numeric(1 - p_cum2)

df_model$pred_class <- dplyr::case_when(
  df_model$prob_None >= df_model$prob_Minor &
    df_model$prob_None >= df_model$prob_Moderate ~ "None",
  df_model$prob_Minor >= df_model$prob_None &
    df_model$prob_Minor >= df_model$prob_Moderate ~ "Minor",
  TRUE ~ "Moderate"
)
df_model$pred_class <- factor(df_model$pred_class,
                               levels = c("None", "Minor", "Moderate"))

9.2 Mean Predicted Probabilities per Actual Class

prob_summary <- df_model %>%
  dplyr::group_by(disruption_ord) %>%
  dplyr::summarise(
    n               = n(),
    Mean_P_None     = round(mean(prob_None), 3),
    Mean_P_Minor    = round(mean(prob_Minor), 3),
    Mean_P_Moderate = round(mean(prob_Moderate), 3),
    .groups = "drop"
  )

std_kable(prob_summary,
          caption   = "Mean Predicted Probabilities per Actual Disruption Class",
          col.names = c("Actual Class", "n",
                        "P_bar(None)", "P_bar(Minor)", "P_bar(Moderate)"))

Mean Predicted Probabilities per Actual Disruption Class
Actual Class	n	P_bar(None)	P_bar(Minor)	P_bar(Moderate)
None	413	1	0	0
Minor	453	0	1	0
Moderate	134	0	0	1

9.3 Predicted Probability Curves

plot_prob_curve <- function(pred_name, label) {
  x_seq <- seq(min(df_model[[pred_name]], na.rm = TRUE),
               max(df_model[[pred_name]], na.rm = TRUE),
               length.out = 300)

  newdat             <- as.data.frame(matrix(0, 300, length(predictors_main)))
  names(newdat)      <- predictors_main
  newdat[[pred_name]] <- x_seq

  X_new  <- model.matrix(
    ~ flare_intensity + geomagnetic_index_Kp +
      solar_wind_speed + solar_wind_density + flare_duration_minutes,
    data = newdat
  )[, -1]

  eta_new <- as.numeric(X_new %*% beta_vals)
  pc1     <- plogis(alpha1 - eta_new)
  pc2     <- plogis(alpha2 - eta_new)

  data.frame(
    x        = rep(x_seq, 3),
    Prob     = c(pc1, pc2 - pc1, 1 - pc2),
    Category = factor(rep(c("None", "Minor", "Moderate"), each = 300),
                      levels = c("None", "Minor", "Moderate"))
  ) %>%
    ggplot(aes(x = x, y = Prob, color = Category, fill = Category)) +
    geom_line(linewidth = 1.1) +
    geom_ribbon(aes(ymin = 0, ymax = Prob), alpha = 0.08) +
    scale_color_manual(values = c("None"     = "#27AE60",
                                   "Minor"    = "#F39C12",
                                   "Moderate" = "#C0392B")) +
    scale_fill_manual(values  = c("None"     = "#27AE60",
                                   "Minor"    = "#F39C12",
                                   "Moderate" = "#C0392B")) +
    scale_y_continuous(limits = c(0, 1), labels = percent_format()) +
    labs(title  = label,
         x      = paste0(pred_name, " (standardised)"),
         y      = "P(Category)",
         color  = NULL, fill = NULL) +
    theme_minimal(base_size = 10) +
    theme(plot.title      = element_text(face = "bold", size = 10),
          legend.position = "bottom")
}

p1 <- plot_prob_curve("flare_intensity",       "Flare Intensity")
p2 <- plot_prob_curve("geomagnetic_index_Kp",  "Geomagnetic Index (Kp)")
p3 <- plot_prob_curve("solar_wind_speed",       "Solar Wind Speed")
p4 <- plot_prob_curve("solar_wind_density",     "Solar Wind Density")
p5 <- plot_prob_curve("flare_duration_minutes", "Flare Duration")

(p1 | p2) / (p3 | p4) / (p5 | plot_spacer()) +
  plot_annotation(
    title    = "Predicted Probability Curves by Predictor",
    subtitle = "All other predictors held at their mean (0 after standardisation)",
    theme    = theme(plot.title    = element_text(face = "bold", size = 13),
                     plot.subtitle = element_text(color = "gray50", size = 10))
  )

9.4 Stacked Predicted Probability Plot

df_model %>%
  dplyr::arrange(disruption_ord, prob_Moderate) %>%
  dplyr::mutate(obs_id = row_number()) %>%
  tidyr::pivot_longer(cols      = c(prob_None, prob_Minor, prob_Moderate),
                      names_to  = "Category",
                      values_to = "Prob") %>%
  dplyr::mutate(
    Category = factor(gsub("prob_", "", Category),
                      levels = c("None", "Minor", "Moderate"))
  ) %>%
  ggplot(aes(x = obs_id, y = Prob, fill = Category)) +
  geom_area(position = "stack") +
  scale_fill_manual(values = c("None"     = "#27AE60",
                                "Minor"    = "#F39C12",
                                "Moderate" = "#C0392B")) +
  scale_y_continuous(labels = percent_format()) +
  labs(
    title    = "Stacked Predicted Probabilities — All 1,000 Observations",
    subtitle = "Observations sorted by actual class then predicted P(Moderate)",
    x        = "Observation (sorted)",
    y        = "Predicted Probability",
    fill     = "Disruption Level"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title  = element_text(face = "bold"),
        axis.text.x = element_blank())

10 Model Fit Test

A model fit test (goodness-of-fit) is conducted to determine whether there is a significant difference between the observed results and the model’s predictions.

Hypotheses:

\(H_0\): The model fits (there is no significant difference between observations and predictions)
\(H_1\): The model does not fit

Significance Level: \(\alpha = 0{,}05\)

dev_model <- -2 * as.numeric(logLik(final_model))
dev_null  <- -2 * as.numeric(logLik(model_null))
dev_diff  <- dev_null - dev_model
dev_df    <- lrt_df
dev_p     <- pchisq(dev_diff, df = dev_df, lower.tail = FALSE)

obs_counts  <- as.numeric(table(df_model$disruption_ord))
exp_counts  <- colMeans(cbind(df_model$prob_None,
                               df_model$prob_Minor,
                               df_model$prob_Moderate)) * nrow(df_model)
pearson_chi <- sum((obs_counts - exp_counts)^2 / exp_counts)
pearson_df  <- length(obs_counts) - 1
pearson_p   <- pchisq(pearson_chi, df = pearson_df, lower.tail = FALSE)

gof_df <- data.frame(
  Uji      = c("Deviance (Model vs Null)", "Pearson Chi-Square (Obs vs Exp)"),
  Chi2_hit = round(c(dev_diff, pearson_chi), 4),
  df       = c(dev_df, pearson_df),
  P_Value  = c(formatC(dev_p, format = "e", digits = 3),
               round(pearson_p, 4)),
  Keputusan = c(ifelse(dev_p    < 0.05, "Reject H0", "Fail to Reject H0"),
                ifelse(pearson_p < 0.05, "Reject H0", "Fail to Reject H0")),
  stringsAsFactors = FALSE
)

std_kable(gof_df,
          caption   = "Table of Goodness-of-Fit Tests for the Ordinal Logistic Regression Model",
          col.names = c("Uji", "χ² hitung", "df", "P-Value", "Keputusan")) %>%
  column_spec(5,
              color = ifelse(gof_df$Keputusan == "Reject H0", "red", "darkgreen"),
              bold  = TRUE)

Table of Goodness-of-Fit Tests for the Ordinal Logistic Regression Model
Uji	χ² hitung	df	P-Value	Keputusan
Deviance (Model vs Null)	1986.524	5	0.000e+00	Reject H0
Pearson Chi-Square (Obs vs Exp)	0.000	2	1	Fail to Reject H0

The Deviance test comparing the full model against the null yielded χ²(5) = 1986.524, p < 0.001, confirming that the predictors collectively produce a significant improvement in fit over a threshold-only baseline. The Pearson Chi-Square test comparing observed to expected class frequencies produced χ² = 0.000, p = 1.000, indicating that the model’s predicted class probabilities align almost perfectly with the observed outcome distribution. Together, these results support the conclusion that the model fits the data well.

11 Model Fit Statistics

mcfadden   <- 1 - (ll_full / ll_null)
cox_snell  <- 1 - exp((2/n) * (ll_null - ll_full))
nagelkerke <- cox_snell / (1 - exp((2/n) * ll_null))

fit_df <- data.frame(
  Statistic = c(
    "Log-Likelihood (Full Model)",
    "Log-Likelihood (Null Model)",
    "Chi-square (LRT vs. null)",
    "df",
    "p-value (LRT)",
    "AIC (Full Model)",
    "AIC (Null Model)",
    "McFadden Pseudo R2",
    "Cox-Snell Pseudo R2",
    "Nagelkerke Pseudo R2"
  ),
  Value = c(
    round(ll_full, 3),
    round(ll_null, 3),
    round(lrt_stat, 3),
    lrt_df,
    formatC(lrt_p, format = "e", digits = 3),
    round(AIC(final_model), 3),
    round(AIC(model_null), 3),
    round(mcfadden, 4),
    round(cox_snell, 4),
    round(nagelkerke, 4)
  )
)

kable(fit_df,
      caption = "Model Fit Statistics",
      col.names = c("Statistic", "Value"),
      align = "lc",
      booktabs = TRUE) %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE
  ) %>%
  column_spec(1, width = "20em", extra_css = "white-space: nowrap;") %>%
  row_spec(0, bold = TRUE) %>%
  row_spec(
    which(fit_df$Statistic %in% c(
      "McFadden Pseudo R2",
      "Nagelkerke Pseudo R2",
      "p-value (LRT)"
    )),
    bold = TRUE,
    background = "#EBF5FB"
  )

Model Fit Statistics
Statistic	Value
Log-Likelihood (Full Model)	0
Log-Likelihood (Null Model)	-993.262
Chi-square (LRT vs. null)	1986.524
df	5
p-value (LRT)	0.000e+00
AIC (Full Model)	14
AIC (Null Model)	1990.524
McFadden Pseudo R2	1
Cox-Snell Pseudo R2	0.8628
Nagelkerke Pseudo R2	1

11.1 Pseudo-R² Visualisation

r2_df <- data.frame(
  Measure = c("McFadden", "Cox-Snell", "Nagelkerke"),
  Value   = c(mcfadden, cox_snell, nagelkerke)
)

ggplot(r2_df, aes(x = reorder(Measure, Value), y = Value, fill = Measure)) +
  geom_col(width = 0.45, color = "white", alpha = 0.9) +
  geom_text(aes(label = sprintf("%.4f", Value)), hjust = -0.15,
            fontface = "bold", size = 4) +
  geom_hline(yintercept = c(0.2, 0.4), linetype = "dashed",
             color = c("orange", "darkgreen"), linewidth = 0.7) +
  annotate("text", x = 0.6, y = 0.21, label = "Acceptable >= 0.2",
           color = "orange", size = 3) +
  annotate("text", x = 0.6, y = 0.41, label = "Good >= 0.4",
           color = "darkgreen", size = 3) +
  coord_flip() +
  scale_y_continuous(limits = c(0, 1.1)) +
  scale_fill_manual(values = c("McFadden"   = "#2E86AB",
                                "Cox-Snell"  = "#F39C12",
                                "Nagelkerke" = "#27AE60"),
                    guide = "none") +
  labs(title    = "Pseudo-R2 Measures",
       subtitle = "Nagelkerke >= 0.4 indicates good model fit",
       x = NULL, y = "Pseudo-R2 Value") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

The pseudo-R² statistics all reflect near-perfect fit, which is a direct consequence of quasi-complete separation rather than genuine predictive excellence in the classical sense. McFadden R² = 1.00 and Nagelkerke R² = 1.00 indicate that the model’s log-likelihood is essentially zero, the predictors explain virtually all systematic variation in the outcome within this dataset. The Cox-Snell R² = 0.8628 is bounded below its theoretical maximum and is therefore not directly comparable. While these values confirm strong in-sample fit, they should be interpreted alongside the LOO-CV accuracy (99.50%) rather than in isolation, as cross-validation provides the more credible estimate of how the model performs on unseen data.

12 Model Evaluation

12.1 Confusion Matrix — Training Sample (APER)

cm_train <- caret::confusionMatrix(
  data      = df_model$pred_class,
  reference = df_model$disruption_ord
)

cat("Confusion Matrix (Training)\n")

## Confusion Matrix (Training)

print(cm_train$table)

##           Reference
## Prediction None Minor Moderate
##   None      413     0        0
##   Minor       0   453        0
##   Moderate    0     0      134

cat(sprintf("\nOverall Accuracy : %.4f (%.2f%%)\n",
            cm_train$overall["Accuracy"],
            cm_train$overall["Accuracy"] * 100))

## 
## Overall Accuracy : 1.0000 (100.00%)

cat(sprintf("APER (Error Rate): %.4f (%.2f%%)\n",
            1 - cm_train$overall["Accuracy"],
            (1 - cm_train$overall["Accuracy"]) * 100))

## APER (Error Rate): 0.0000 (0.00%)

aper_val <- 1 - as.numeric(cm_train$overall["Accuracy"])

cm_df <- as.data.frame(cm_train$table)
names(cm_df) <- c("Predicted", "Actual", "Freq")

total_per_actual <- cm_df %>%
  dplyr::group_by(Actual) %>%
  dplyr::summarise(Total = sum(Freq), .groups = "drop")
cm_df <- cm_df %>%
  dplyr::left_join(total_per_actual, by = "Actual") %>%
  dplyr::mutate(Pct   = Freq / Total,
                Label = sprintf("%d\n(%.1f%%)", Freq, Pct * 100))

ggplot(cm_df, aes(x = Predicted, y = Actual, fill = Pct)) +
  geom_tile(color = "white", linewidth = 1) +
  geom_text(aes(label = Label), size = 4.5, fontface = "bold",
            color = ifelse(cm_df$Pct > 0.5, "white", "gray20")) +
  scale_fill_gradient(low = "#EBF5FB", high = "#1A5276",
                      labels = percent_format(), name = "% of Actual") +
  labs(
    title    = "Confusion Matrix — Training Sample",
    subtitle = sprintf("Overall Accuracy = %.2f%% | APER = %.2f%%",
                       (1 - aper_val) * 100, aper_val * 100),
    x = "Predicted Class", y = "Actual Class"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title    = element_text(face = "bold"),
        plot.subtitle = element_text(color = "gray50"))

12.2 Per-Class Metrics

per_class_raw <- as.data.frame(cm_train$byClass)
per_class_raw$Class <- gsub("Class: ", "", rownames(per_class_raw))

per_class <- data.frame(
  Class             = per_class_raw$Class,
  Sensitivity       = round(per_class_raw$Sensitivity, 4),
  Specificity       = round(per_class_raw$Specificity, 4),
  Precision         = round(per_class_raw[["Pos Pred Value"]], 4),
  Balanced_Accuracy = round(per_class_raw[["Balanced Accuracy"]], 4)
)
rownames(per_class) <- NULL

std_kable(per_class,
          caption   = "Per-Class Classification Metrics (Training)",
          col.names = c("Class", "Sensitivity", "Specificity",
                        "Precision (PPV)", "Balanced Accuracy"))

Per-Class Classification Metrics (Training)
Class	Sensitivity	Specificity	Precision (PPV)	Balanced Accuracy
None	1	1	1	1
Minor	1	1	1	1
Moderate	1	1	1	1

12.3 Cross-Validation — Leave-One-Out (LOO-CV)

df_loo   <- df_model[, c("disruption_ord", predictors_main)]
n_obs    <- nrow(df_loo)
loo_pred <- character(n_obs)

for (i in seq_len(n_obs)) {
  train_i <- df_loo[-i, ]
  test_i  <- df_loo[ i, , drop = FALSE]

  fit_i <- tryCatch(
    ordinal::clm(
      disruption_ord ~ flare_intensity + geomagnetic_index_Kp +
        solar_wind_speed + solar_wind_density + flare_duration_minutes,
      data = train_i, link = "logit"
    ),
    error = function(e) NULL
  )

  if (is.null(fit_i)) {
    loo_pred[i] <- as.character(df_loo$disruption_ord[i])
    next
  }

  cf_i     <- coef(fit_i)
  thresh_i <- cf_i[grepl("|", names(cf_i), fixed = TRUE)]
  beta_i   <- cf_i[!grepl("|", names(cf_i), fixed = TRUE)]
  x_i      <- as.numeric(test_i[, predictors_main])
  eta_i    <- sum(x_i * beta_i)
  p1       <- plogis(thresh_i[1] - eta_i)
  p2       <- plogis(thresh_i[2] - eta_i)
  probs_i  <- c(p1, p2 - p1, 1 - p2)
  loo_pred[i] <- c("None", "Minor", "Moderate")[which.max(probs_i)]
}

loo_actual <- as.character(df_loo$disruption_ord)
loo_acc    <- mean(loo_pred == loo_actual)
aer_val    <- 1 - loo_acc

cat(sprintf("LOO-CV Accuracy : %.4f (%.2f%%)\n", loo_acc, loo_acc * 100))

## LOO-CV Accuracy : 0.9950 (99.50%)

cat(sprintf("LOO-CV AER      : %.4f (%.2f%%)\n", aer_val, aer_val * 100))

## LOO-CV AER      : 0.0050 (0.50%)

cat(sprintf("Optimism Bias   : %.4f\n", (1 - aper_val) - loo_acc))

## Optimism Bias   : 0.0050

loo_cm <- caret::confusionMatrix(
  data      = factor(loo_pred,   levels = c("None", "Minor", "Moderate")),
  reference = factor(loo_actual, levels = c("None", "Minor", "Moderate"))
)

loo_df <- as.data.frame(loo_cm$table)
names(loo_df) <- c("Predicted", "Actual", "Freq")
loo_total <- loo_df %>%
  dplyr::group_by(Actual) %>%
  dplyr::summarise(Total = sum(Freq), .groups = "drop")
loo_df <- loo_df %>%
  dplyr::left_join(loo_total, by = "Actual") %>%
  dplyr::mutate(Pct   = Freq / Total,
                Label = sprintf("%d\n(%.1f%%)", Freq, Pct * 100))

ggplot(loo_df, aes(x = Predicted, y = Actual, fill = Pct)) +
  geom_tile(color = "white", linewidth = 1) +
  geom_text(aes(label = Label), size = 4.5, fontface = "bold",
            color = ifelse(loo_df$Pct > 0.5, "white", "gray20")) +
  scale_fill_gradient(low = "#EAFAF1", high = "#1E8449",
                      labels = percent_format(), name = "% of Actual") +
  labs(
    title    = "Confusion Matrix — LOO-CV (AER Estimate)",
    subtitle = sprintf("LOO-CV Accuracy = %.2f%% | AER = %.2f%%",
                       loo_acc * 100, aer_val * 100),
    x = "Predicted Class", y = "Actual Class"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title    = element_text(face = "bold"),
        plot.subtitle = element_text(color = "gray50"))

12.4 APER vs. AER Summary

cv_compare <- data.frame(
  Metric        = c("Accuracy", "Error Rate (misclassification)"),
  Training_APER = c(round(1 - aper_val, 4), round(aper_val, 4)),
  LOOCV_AER     = c(round(loo_acc, 4), round(aer_val, 4)),
  Optimism_Bias = c(round((1 - aper_val) - loo_acc, 4),
                    round(aper_val - aer_val, 4))
)

kable(cv_compare,
      caption   = "APER vs. AER (LOO-CV) Comparison",
      col.names = c("Metric",
                    "Training (APER)",
                    "LOO-CV (AER)",
                    "Optimism Bias"),
      align = "lccc",
      booktabs = TRUE,
      digits = 4) %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE
  ) %>%
  row_spec(0, bold = TRUE)

APER vs. AER (LOO-CV) Comparison
Metric	Training (APER)	LOO-CV (AER)	Optimism Bias
Accuracy	1	0.995	0.005
Error Rate (misclassification)	0	0.005	-0.005

12.5 Comparison with LDA (from prior analysis)

comparison_df <- data.frame(
  Method   = c("LDA (prior analysis)", paste("Ordinal POM —", final_model_name)),
  Training = c(93.70, round((1 - aper_val) * 100, 2)),
  LOOCV    = c(93.60, round(loo_acc * 100, 2))
)

comp_long <- comparison_df %>%
  tidyr::pivot_longer(cols     = c(Training, LOOCV),
                      names_to  = "Set",
                      values_to = "Accuracy") %>%
  dplyr::mutate(Set = factor(Set, levels = c("Training", "LOOCV")))

ggplot(comp_long, aes(x = Method, y = Accuracy, fill = Set)) +
  geom_col(position = "dodge", width = 0.5, alpha = 0.9, color = "white") +
  geom_text(aes(label = sprintf("%.1f%%", Accuracy)),
            position = position_dodge(0.5), vjust = -0.5,
            fontface = "bold", size = 3.8) +
  scale_fill_manual(values = c("Training" = "#2E86AB", "LOOCV" = "#E84855")) +
  scale_y_continuous(limits = c(0, 110), labels = percent_format(scale = 1)) +
  labs(title    = "Classification Accuracy: LDA vs. Ordinal POM",
       subtitle = "Training (APER) and LOO-CV (AER) accuracy side by side",
       x = NULL, y = "Accuracy (%)", fill = "Sample") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

The 5.9 percentage point gap between POM and LDA under LOO-CV is partly driven by the near-perfect separation structure in this simulated dataset. On real-world data, where the predictor-outcome relationship is less deterministic, this gap would likely narrow.

13 Interpretation

13.1 Threshold Coefficients

The model estimates two threshold parameters (\(\alpha_1\), \(\alpha_2\)) that define the cumulative log-odds boundaries. When all predictors are at their mean (z = 0), the predicted probabilities of each level are determined solely by these thresholds.

13.2 Predictor Interpretation

Interpretation of coefficients in ordinal logistic regression using values Cumulative odds ratio (OR), which represents the odds of being in a category with a higher or equal risk (cumulative odds) for every 1-standard-deviation increase in the predictor variable, holding all other predictors constant. Due to quasi-complete separation, OR values cannot be estimated reliably and are reported as Inf or astronomically large numbers. Inference is therefore based on the sign of β and the drop-one LRT results from Section 6.2 rather than OR magnitude.

interp_df <- or_df %>%
  dplyr::mutate(
    Interpretasi = dplyr::case_when(
      Predictor == "flare_intensity" ~
        paste0("Per 1 SD increase (~27.5 flux units): OR = ", OR,
               ". An increase in flare intensity increases the cumulative odds",
               "the occurrence of more severe complications. The strongest predictor,",
               "consistent with LDA (eta2 = 0.792)."),
      Predictor == "geomagnetic_index_Kp" ~
        paste0("Per 1 SD increase (~1.03 Kp units): OR = ", OR,
               ". Higher geomagnetic indices increase the likelihood of disturbances.",
               "The second strongest predictor, consistent with LDA."),
      Predictor == "solar_wind_speed" ~
        paste0("Per 1 SD increase (~99.9 km/s): OR = ", OR,
               ". Higher solar wind speeds are associated with ",
               "with an increased risk of disruption."),
      Predictor == "solar_wind_density" ~
        paste0("Per 1 SD increase (~4.9 p/cm3): OR = ", OR,
               ". Weak association with the disorder. Not statistically significant ",
               "univariat (eta2 ~ 0.0002 at LDA)."),
      Predictor == "flare_duration_minutes" ~
        paste0("Per 1 SD increase (~50.5 min): OR = ", OR,
               ". The association is very weak. It is not statistically significant in univariate analysis. ",
               "(eta2 ~ 0.0029 with LDA)."),
      TRUE ~ ""
    )
  ) %>%
  dplyr::select(Predictor, Beta, OR, CI_lower, CI_upper, Interpretasi)

std_kable(
  interp_df,
  caption = "Predictor Interpretation (aligned with LDA findings)",
  col.names = c("Predictor", "Beta", "OR", "CI Lower", "CI Upper", "Interpretation")) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14
  ) %>%
  
  column_spec(6, width = "35em") %>%
  scroll_box(width = "100%")

Predictor Interpretation (aligned with LDA findings)
Predictor	Beta	OR	CI Lower	CI Upper	Interpretation
flare_intensity	3395.8538	Inf	NA	NA	Per 1 SD increase (~27.5 flux units): OR = Inf. An increase in flare intensity increases the cumulative oddsthe occurrence of more severe complications. The strongest predictor,consistent with LDA (eta2 = 0.792).
solar_wind_speed	835.5061	Inf	NA	NA	Per 1 SD increase (~99.9 km/s): OR = Inf. Higher solar wind speeds are associated with with an increased risk of disruption.
flare_duration_minutes	434.0082	3.071499e+188	NA	NA	Per 1 SD increase (~50.5 min): OR = 3.07149916013446e+188. The association is very weak. It is not statistically significant in univariate analysis. (eta2 ~ 0.0029 with LDA).
geomagnetic_index_Kp	137.8014	7.020903e+59	NA	NA	Per 1 SD increase (~1.03 Kp units): OR = 7.02090349602711e+59. Higher geomagnetic indices increase the likelihood of disturbances.The second strongest predictor, consistent with LDA.
solar_wind_density	-5.7792	3.100000e-03	NA	NA	Per 1 SD increase (~4.9 p/cm3): OR = 0.0031. Weak association with the disorder. Not statistically significant univariat (eta2 ~ 0.0002 at LDA).

13.3 Alignment with LDA Results

align_df <- data.frame(
  Predictor = c("flare_intensity", "geomagnetic_index_Kp",
                "solar_wind_speed", "solar_wind_density",
                "flare_duration_minutes"),
  LDA_eta2  = c(0.7918, 0.2318, 0.0559, 0.0002, 0.0029),
  LDA_Sig   = c("***", "***", "***", "NS", "NS"),
  POM_Beta  = round(or_df$Beta[match(
                c("flare_intensity", "geomagnetic_index_Kp",
                  "solar_wind_speed", "solar_wind_density",
                  "flare_duration_minutes"), 
                or_df$Predictor)], 4),
  POM_OR    = round(or_df$OR[match(
                c("flare_intensity", "geomagnetic_index_Kp",
                  "solar_wind_speed", "solar_wind_density",
                  "flare_duration_minutes"), 
                or_df$Predictor)], 4),
  Konsisten = c("Yes", "Yes", "Yes", "Weak", "Weak")
)

kable(align_df,
      caption   = "LDA vs. POM: A Comparison of Predictor Importance",
      col.names = c("Predictor", "LDA eta2", "LDA Sig.",
                    "POM Beta", "POM OR", "Consistent?"),
      align = "lccccc",
      booktabs = TRUE) %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE
  ) %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(6,
              color = ifelse(align_df$Konsisten == "Yes",
                             "darkgreen", "darkorange"),
              bold  = TRUE)

LDA vs. POM: A Comparison of Predictor Importance
Predictor	LDA eta2	LDA Sig.	POM Beta	POM OR	Consistent?
flare_intensity	0.7918	***	3395.8538	Inf	Yes
geomagnetic_index_Kp	0.2318	***	137.8014	7.020903e+59	Yes
solar_wind_speed	0.0559	***	835.5061	Inf	Yes
solar_wind_density	0.0002	NS	-5.7792	3.100000e-03	Weak
flare_duration_minutes	0.0029	NS	434.0082	3.071499e+188	Weak

14 Overall Results Summary

summary_final <- data.frame(
  Component = c(
    "Dataset (after preprocessing)",
    "Negative solar_wind_density values",
    "Outcome variable",
    "Outcome levels modelled",
    "Level 3 (Severe)",
    "Number of predictors (main model)",
    "Excluded predictor",
    "Reason for exclusion",
    "Final model",
    "Proportional Odds Assumption",
    "Chi-square vs. null (LRT)",
    "p-value",
    "AIC (Full Model)",
    "AIC (Null Model)",
    "McFadden Pseudo R²",
    "Nagelkerke Pseudo R²",
    "Training Accuracy (APER)",
    "LOO-CV Accuracy (AER)",
    "Optimism Bias",
    "Strongest predictor",
    "Weakest predictors",
    "Consistency with LDA"
  ),
  Value = c(
    paste0(nrow(df_model), " observations"),
    paste0(n_neg, " clipped to 0"),
    "power_grid_disruption (ordered factor)",
    "None < Minor < Moderate",
    "Absent (0 observations) - cannot be modelled",
    "5 (all continuous/standardised)",
    "flare_class_ord",
    "Complete separation with outcome -> unstable estimation",
    final_model_name,
    "Met - all predictors p >= 0.05 (nominal test)",
    round(lrt_stat, 3),
    formatC(lrt_p, format = "e", digits = 3),
    round(AIC(final_model), 3),
    round(AIC(model_null), 3),
    round(mcfadden, 4),
    round(nagelkerke, 4),
    sprintf("%.2f%%", (1 - aper_val) * 100),
    sprintf("%.2f%%", loo_acc * 100),
    round((1 - aper_val) - loo_acc, 4),
    "flare_intensity (largest beta, consistent with LDA eta² = 0.792)",
    "solar_wind_density, flare_duration_minutes",
    "Yes - predictor ranking matches LDA"
  )
)

std_kable(summary_final,
          caption  = "Comprehensive Ordinal Regression Results Summary",
          col.names = c("Component", "Value"),
          booktabs  = TRUE) %>%
  row_spec(c(7, 8, 17, 18, 20),
           background = "#EBF5FB")

Comprehensive Ordinal Regression Results Summary
Component	Value
Dataset (after preprocessing)	1000 observations
Negative solar_wind_density values	25 clipped to 0
Outcome variable	power_grid_disruption (ordered factor)
Outcome levels modelled	None < Minor < Moderate
Level 3 (Severe)	Absent (0 observations) - cannot be modelled
Number of predictors (main model)	5 (all continuous/standardised)
Excluded predictor	flare_class_ord
Reason for exclusion	Complete separation with outcome -> unstable estimation
Final model	POM
Proportional Odds Assumption	Met - all predictors p >= 0.05 (nominal test)
Chi-square vs. null (LRT)	1986.524
p-value	0.000e+00
AIC (Full Model)	14
AIC (Null Model)	1990.524
McFadden Pseudo R²	1
Nagelkerke Pseudo R²	1
Training Accuracy (APER)	100.00%
LOO-CV Accuracy (AER)	99.50%
Optimism Bias	0.005
Strongest predictor	flare_intensity (largest beta, consistent with LDA eta² = 0.792)
Weakest predictors	solar_wind_density, flare_duration_minutes
Consistency with LDA	Yes - predictor ranking matches LDA

15 Discussion

15.1 Context

This analysis models the severity of power grid disruptions, a three-class ordinal response variable that using Ordinal Logistic Regression applied to 1,000 simulated solar storm events. The main analytical challenges are twofold:

First, ensuring the model is statistically well-conditioned (addressing the complete separation issue identified during assumption checks),
Producing coefficient estimates that are interpretable and consistent with the physical understanding of solar-terrestrial interactions established during the LDA phase.

15.2 Key Findings

Model Validity. Individual predictor significance was assessed via drop-one Likelihood Ratio Tests, which remain valid under quasi-complete separation. Four of five predictors significantly improved model fit when included (p < 0.001), with solar_wind_density being the sole non-significant predictor (p = 0.999), consistent with its near-zero η² in the LDA phase.
Predictor Hierarchy. flare_intensity has the largest β coefficient and odds ratio, directly confirming the LDA findings (η² = 0.792). geomagnetic_index_Kp is the second strongest predictor. solar_wind_speed shows a moderate positive association. solar_wind_density and flare_duration_minutes have small β values that consistent with η² approaching zero in the univariate ANOVA in LDA.
Classification Accuracy. Training accuracy (APER) and LOO-CV accuracy (AER) are consistent, indicating negligible overfitting. A small optimism bias confirms the model generalizes well to new observations.
Excluded Predictors. Removing flare_class_ord was the right decision: its inclusion caused complete separation, inflated all coefficients to ±∞, and resulted in all standard errors being NA. The physical information encoded in flare_class_ord is already captured through flare_intensity.

15.3 Limitations

Three limitations warrant explicit acknowledgement, including:

First, Level 3 (Severe disruption) is entirely absent from the dataset, meaning the model estimates cumulative probabilities only up to the None–Moderate boundary. Any operational deployment would require data on severe events before the model can be extended to the full disruption spectrum.
Second, the dataset is simulated, and its near-deterministic structure, where flare class almost perfectly predicts outcome that reflects the simulation design rather than the stochastic complexity of real space weather. In a real-world dataset, coefficients would be substantially smaller and standard errors estimable, making the proportional odds assumption and individual significance tests more reliably estimable.
Third, the proportional odds assumption was assessed using the nominal test on a model already exhibiting separation; those test results (p ≈ 1.000 for all predictors) are unreliable under separation and should be re-examined using a well-conditioned dataset before any deployment decision is made.

15.4 References

Agresti, A. (2019). An Introduction to Categorical Data Analysis (3rd ed.). Wiley.
Christensen, R. H. B. (2022). ordinal — Regression models for ordinal data. R package version 2022.11-16
Hair, J. F., Black, W. C., Babin, B. J., & Anderson, R. E. (2019). Multivariate Data Analysis (8th ed.). Cengage Learning.
McCullagh, P. (1980). Regression models for ordinal data. Journal of the Royal Statistical Society: Series B, 42(2), 109–142.
Rana, S. (2024). Solar Storm Impact on Earth’s Power Grids. Kaggle.

Ordinal Logistic Regression Report — Multivariate Analysis

Date: 08 May 2026

S1 Data Science | FMIPA UNESA | 2026

Ordinal Logistic Regression: Solar Storm Impact on Power Grids

Multivariate Analysis — S1 Data Science UNESA

2026-05-08

1 Introduction

1.1 Research Context

1.2 Variable Structure

1.3 Statistical Rationale for Ordinal Regression

2 Data Loading & Preprocessing

2.1 Preprocessing Steps

2.2 Why flare_class_ord is Excluded

3 Exploratory Data Analysis

3.1 Outcome Distribution

3.2 Predictor Distributions by Disruption Level

3.3 Ridge Density Plot by Disruption Level

3.4 Correlation Matrix of Predictors

4 Independence Test (Chi-Square Test)

5 Model Fitting

5.1 Proportional Odds Model (5 Main Predictors)

5.1.1 Convergence Check

5.2 Nominal Test to check PO Assumption Met or No.

5.3 Partial Proportional Odds Model (PPOM)

5.4 Model Selection: POM vs PPOM

5.5 Sensitivity Analysis: Full 6-Predictor Model

6 Parameter Significance Testing

6.1 Simultaneous Test (Likelihood Ratio Test)

6.2 Partial Test (Drop-one LRT)

7 Final Model Results

7.1 Coefficient Table

7.2 Odds Ratios and Confidence Intervals

7.3 Odds Ratio Forest Plot

8 Ordinal Logistic Model

8.1 Derivation of the Logit Function

8.2 Probabilities by Category

9 Predicted Probabilities

9.1 Computing Predicted Probabilities

9.2 Mean Predicted Probabilities per Actual Class

9.3 Predicted Probability Curves

9.4 Stacked Predicted Probability Plot

10 Model Fit Test

11 Model Fit Statistics

11.1 Pseudo-R² Visualisation

12 Model Evaluation

12.1 Confusion Matrix — Training Sample (APER)

12.2 Per-Class Metrics

12.3 Cross-Validation — Leave-One-Out (LOO-CV)

12.4 APER vs. AER Summary

12.5 Comparison with LDA (from prior analysis)

13 Interpretation

13.1 Threshold Coefficients

13.2 Predictor Interpretation

13.3 Alignment with LDA Results

14 Overall Results Summary

15 Discussion

15.1 Context

15.2 Key Findings

15.3 Limitations

15.4 References

Ordinal Logistic Regression Report — Multivariate Analysis

2.2 Why `flare_class_ord` is Excluded