REGRESSION ANALYSIS OF SEISMIC NETWORK PARAMETERS AND THEIR INFLUENCE ON EARTHQUAKE DATA QUALITY

REGRESSION ANALYSIS OF SEISMIC NETWORK PARAMETERS AND THEIR INFLUENCE ON EARTHQUAKE DATA QUALITY Dosen Pengampu: Prof. Dr. Anang Kurnia, S.Si., M.Si. Dr. Agus Mohamad Soleh, S.Si., M.T. Disusun oleh: Sanyyatul Rohma Ahfa (G7501242012) PROGRAM STUDI BIOFISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR 2025 REGRESSION ANALYSIS OF SEISMIC NETWORK PARAMETERS AND THEIR INFLUENCE ON EARTHQUAKE DATA QUALITY Sanyyatul Rohma Ahfa1 ¹Department of Biophysics, Faculty of Mathematics and Natural Sciences, IPB University, Bogor, Indonesia ABSTRACT The reliability of earthquake data products fundamentally depends on the configuration and performance of seismic monitoring networks. This study introduces a new composite metric that integrates magnitude consistency, depth reliability, and station coverage to holistically assess the quality of seismic data. We apply multiple linear regression to a global dataset of 51 earthquakes (𝑀 ≥ 6.0) from the revised USGS catalog to quantify the extent to which key network parameters azimuthal gap, distance to the nearest station, root mean square (RMS) travel time residuals, and horizontal location uncertainty explain the variability in this quality value. The resulting model is statistically significant (𝐹(5,45) = 3.074,𝑝 = 0.018) and explains approximately 25.5% of the variance (𝑅² = 0.255). Horizontal location error is identified as the sole significant independent predictor (𝛽 = −2.443,𝑝 = 0.041) , demonstrating that epicentral accuracy is a primary control on overall data reliability. Diagnostic analyzes confirm the statistical robustness of the model. These findings provide quantitative evidence that network performance, particularly location accuracy, is a crucial determinant of data quality. The proposed composite evaluation framework offers a reproducible, data-driven approach to network evaluation and optimization, with direct implications for improving the reliability of seismic catalogs used in hazard analysis and early warning systems. Keywords: composite quality metric, horizontal location uncertainty, global seismicity, network geometry, statistical diagnostics. INTRODUCTION The quality of earthquake data is a fundamental prerequisite for advancing seismological research and supporting operational applications such as seismic hazard assessment, early warning systems, and the design of earthquake-resistant infrastructure (Allen and Melgar 2019; Lagomarsino et al. 2021). Modern global monitoring networks, including those of the United States Geological Survey (USGS), routinely generate large amounts of earthquake observations. However, the reliability of these data products is not only determined by the characteristics of the source, but is intrinsically controlled by the configuration of the seismic network, the distribution of stations, and the measurement precision (Bondár and McLaughlin, 2009; Schorlemmer et al. 2018). Consequently, assessing the quality of earthquake data requires explicit consideration of network-related factors that determine the accuracy of parameter estimation. Previous studies have extensively investigated individual source parameters such as magnitude and depth (Kanamori and Rivera 2008), but the broader concept of seismic data quality which encompasses accuracy, internal consistency, and observational completeness has received relatively less quantitative attention (Wyss et al. 2012). In practice, uncertainties in reported earthquake parameters propagate into downstream scientific analyzes and operational products, highlighting the importance of systematic approaches to evaluate and compare data quality across events and regions. The geometry of the seismic network plays a central role in determining the accuracy of earthquake solutions. Parameters such as azimuthal gap, station density, and the proximity of recording stations to the hypocenter directly influence waveform sampling and the stability of the inversion (van Driel et al. 2012). Dense and well-distributed networks generally reduce uncertainties in epicenters and depths, while sparse or geometrically biased configurations can introduce systematic errors and greater variability in estimated parameters (Husen and Hardebeck, 2010; Mignan et al. 2011). Additionally, quality indicators such as horizontal location uncertainty and root mean square (RMS) residuals provide integrated measures of network performance and inversion reliability (Waldhauser and Schaff 2008). Despite their recognized importance, the combined influence of these network parameters on overall earthquake data quality has not yet been extensively quantified within a uniform statistical framework on a global scale (Roselli et al. 2016). A significant limitation in existing research is the lack of standardized, composite metrics capable of capturing multiple dimensions of seismic data quality simultaneously (Akciz and Rockwell 2022). Most studies focus on individual performance indicators or assess specific components of the earthquake solution separately (Tavera H and Seismology Group 2012; Lyra et al. 2018), making it difficult to compare data quality across different events or evaluate trade-offs between competing network properties. Furthermore, relatively few studies use formal regression-based approaches to quantify the relative contributions of various seismic network parameters to overall data quality, particularly in diverse tectonic environments and with different event characteristics (Dudzisz et al. 2022). This gap limits efforts to develop evidence-based strategies for network optimization and quality control in global earthquake monitoring systems (Liu et al. 2020). To address these challenges, this study introduces a composite earthquake data quality score that integrates magnitude consistency, depth reliability, and station coverage into a single evaluative measure. Next, we apply multiple linear regression analysis to systematically quantify the influence of key seismic network parameters azimuthal gap, distance to the nearest station, root mean square error, and horizontal location uncertainty on observed variations in earthquake data quality. The central research question guiding this analysis is: to what extent do the parameters of seismic networks collectively explain the variability in the quality of earthquake data on a global scale? Using a composite dataset of 51 globally distributed earthquakes (M ≥ 6.0) from the revised USGS catalog (USGS, 2024), we develop and evaluate a regression model that explicitly links network configuration to data quality results. Extensive diagnostic tests are used to assess the regression assumptions, including multicollinearity, normality of residuals, heteroscedasticity, and autocorrelation, ensuring statistical rigor and interpretability. By identifying horizontal location error as a dominant factor determining the quality of earthquake data, this study provides quantitative support for prioritizing location accuracy in the design of seismic networks. More generally, the proposed framework offers a reproducible, data-driven approach for evaluating network performance and improving the reliability of global earthquake catalogs. METHODS Data Collection and Processing Earthquake data were obtained from the publicly accessible earthquake catalog of the United States Geological Survey (USGS) (USGS, 2024). To ensure consistency and reliability, the dataset was limited to events with a moment magnitude 𝑀𝑚≥ 6.0 that occurred between November 2023 and March 2024. This selection yielded 51 globally distributed earthquakes from a wide range of tectonic environments, all of which had been assessed and validated by the USGS. Only events with complete metadata for all network and quality-related parameters were retained, making data imputation unnecessary and minimizing potential biases associated with missing values. To quantitatively represent the quality of earthquake data, a composite seismic quality value was constructed as the primary dependent variable. This metric integrates three complementary dimensions of earthquake characterization: (1) magnitude consistency, (2) depth reliability, and (3) station coverage adequacy. Composite measurements of this type have been shown to provide more stable and interpretable assessments of seismic data quality than single-parameter evaluations, especially when comparing events across heterogeneous network configurations (Wiemer and Wyss 2000; Roselli et al. 2016). The composite score was defined as: 𝑆𝑄𝑆 = 𝑍(𝑀𝑚) + 𝑍( 1 𝐷𝑒𝑝𝑡ℎ + 1) + 𝑍(𝑛𝑠𝑡 100) where 𝑍(⋅) denotes z-score standardization, defined as 𝑍(𝑥)=𝑥 −𝜇𝑥 𝜎𝑥. Standardization ensures equal weighting of individual components and mitigates differences in scale and units between the contributing variables (Rahilly et al. 2021). The term for inverse depth emphasizes the reduced reliability typically associated with deeper or poorly constrained events, while the number of stations serves as a proxy for observational redundancy and network coverage. Seismic Network Parameters and Variable Transformations Independent variables were selected to represent key aspects of the seismic network's geometry and performance, which are known to influence the quality of earthquake solutions (Bondár and McLaughlin 2009; van Driel et al. 2012). The following predictors were included: 1. Azimuthal gap: Defined as the largest angular separation between recording stations around an event. This parameter was normalized by dividing by 100 to improve the numerical stability and interpretability of the coefficients. 2. Distance to the nearest station (𝑑𝑚𝑖𝑛): Transformation as 𝑙𝑜𝑔(𝑑𝑚𝑖𝑛 + 0.1) to reduce positive skewness while preserving zero and near-zero values, consistent with previous seismological studies (Husen and Hardebeck 2010). 3. Root Mean Square Error (RMS): The RMS travel time residuals, which reflect waveform mismatch and inversion performance. This variable was retained in its original units, as its distribution approached normality. 4. Horizontal Location Error: A measure of epicentral uncertainty, normalized by dividing by 10 to ensure comparable scaling with other predictors. 5. Interaction Term (Gap × RMS): An interaction between the azimuthal gap and the RMS error was included to capture potential synergistic effects between the network geometry and waveform misfit on data quality. The interaction term was scaled as 𝑔𝑎𝑝 × 𝑅𝑀𝑆 1000 to maintain numerical stability (Tavera H and Seismology Group 2012). These variables collectively represent both geometric and performance-related characteristics of seismic networks, enabling an integrated assessment of their influence on the quality of earthquake data. Regression Modeling Framework Multiple linear regression was used to quantify the relationships between seismic network parameters and the composite earthquake data quality score. The regression model was specified as follows: 𝑆𝑒𝑖𝑠𝑚𝑖𝑐 𝑄𝑢𝑎𝑙𝑖𝑡𝑦 𝑆𝑐𝑜𝑟𝑒 = 𝛽₀ + 𝛽₁(𝐺𝑎𝑝𝑛𝑜𝑟𝑚) + 𝛽₂ log min 0.1) + 𝛽₃(𝑅𝑀𝑆) + 𝛽₄(𝐻𝑜𝑟𝑖𝑧𝑜𝑛𝑡𝑎𝑙 𝐸𝑟𝑟𝑜𝑟 𝑁𝑜𝑟𝑚) + 𝛽₅(𝐺𝑎𝑝 × 𝑅𝑀𝑆) + 𝜀 where 𝛽0is the intercept, 𝛽1− 𝛽5 are regression coefficients, and 𝜀 represents the error term. This formulation allows for a direct comparison of the relative contributions of network parameters to observed variations in earthquake data quality. All analyzes were performed using R version 4.3.1 (R Core Team 2023). Data processing and visualization were performed using the tidyverse suite (Wickham et al. 2019), model outcomes were summarized with broom (Robinson 2023), multicollinearity diagnostics were assessed with car (Fox and Weisberg 2018), and statistical tests were conducted with lmtest (Zeileis and Hothorn 2002). Model Validation and Diagnostic Procedures To ensure statistical robustness and interpretability, the regression model underwent a comprehensive series of diagnostic evaluations according to best practices in geophysical regression analysis (Dudzisz et al. 2022): 1. Multicollinearity: Assessed using Variance Inflation Factors (VIFs), where values above 5 indicate potentially problematic correlations between predictors (James et al. 2021). 2. Normality of residuals: Evaluated using the Shapiro-Wilk test and visual inspection of quantile-quantile (Q-Q) plots (Royston 1995). 3. Homoscedasticity: Tested using the Breusch-Pagan test to assess the constant variance of the residuals across the predicted values (Breusch dan Pagan 1979). 4. Autocorrelation: Investigated using the Durbin-Watson test to verify the independence of the residuals (Durbin and Watson 1950). 5. Influential observations: identified using Cook's distance, with a threshold value defined as 4 𝑛−𝑘−1, where n is the sample size and k is the number of predictors (Cook 1977). Statistical significance was evaluated at the α = 0.05 level for all hypothesis tests. Together, these procedures ensure that the reported regression results reflect genuine relationships between seismic network parameters and earthquake data quality, rather than artifacts of model specification errors or violations of assumptions. RESULTS AND DISCUSSION Descriptive Statistics and Regression Results The analyzed dataset consists of 51 globally distributed earthquakes with moment magnitudes 𝑀𝑤≥ 6.0. The quality of the earthquake data was quantified using a composite seismic quality value that integrates magnitude consistency, depth reliability, and station coverage, as described in the Methods section. This composite measure serves as a holistic proxy for the reliability of earthquake solutions derived from heterogeneous seismic network configurations. The multiple linear regression model, which includes normalized seismic network parameters and an interaction term, is statistically significant at the model level (𝐹(5, 45) = 3.074,𝑝 = 0.018) . The model explains approximately 25.5% of the variance in the composite seismic quality score (𝑅2= 0.2546; 𝑎𝑑𝑗𝑢𝑠𝑡𝑒𝑑 𝑅2= 0.1717). Although a significant portion of the variability remains unexplained reflecting the inherently complex and multi-factorial nature of global earthquake observations this level of explanatory power is consistent with previous regression-based analyzes of seismic data quality and network performance. Research into the regression coefficients (Table 1) shows that horizontal location error is the only statistically significant independent predictor of seismic data quality (𝛽 = −2.443,𝑝 = 0.041). The negative coefficient indicates that a larger epicentral uncertainty is systematically associated with lower overall data quality. This finding underscores the central role of accurate event location in seismic characterization, as uncertainty in location propagates into subsequent estimates of depth, magnitude, and source parameters, thereby decreasing the reliability of earthquake catalogs. Table 1. Multiple Linear Regression Model for Seismic Data Quality (N=51) Predictor Coefficient (β) Std. Error tvalue pvalue 95% Confidence Interval (Intercept) 5.723 2.738 2.090 0.042 * (0.21, 11.24) Azimuthal Gap (norm.) -6.585 7.430 -0.886 0.380 (-21.55, 8.38) Log(Dmin + 0.1) 0.468 0.290 1.614 0.114 (-0.12, 1.05) RMS Error -5.074 3.312 -1.532 0.133 (-11.74, 1.59) Horizontal Error (norm.) -2.443 1.161 -2.105 0.041 ** (-4.78, -0.11) Gap × RMS Interaction 64.742 94.034 0.688 0.495 (-124.54, 254.02) Note: *p < 0.05, **p < 0.01, **p < 0.001 The multiple linear regression model yielded a statistically significant fit to the data (𝐹(5, 45) = 3.074,𝑝 = 0.018). The model accounted for 25.46% of the variance in the composite seismic data quality scores (𝑅² = 0.2546,𝐴𝑑𝑗𝑢𝑠𝑡𝑒𝑑 𝑅² = 0.1717) , with a residual standard error of 1.499 . This indicates that the selected network parameters collectively provide a meaningful, though partial, explanation of data quality variability in global earthquake monitoring, where numerous unmodeled physical and operational factors also play a role. The relative size and uncertainty of the effect of each predictor are illustrated in Figure 1, which shows regression coefficients with their 95% confidence intervals. Only the interval associated with the horizontal location error does not cross zero, which confirms its statistical significance as a unique contribution to the model. In contrast, the wide confidence intervals for the azimuthal gap, the root mean square error, and the interaction term indicate significant estimation uncertainty. These results suggest that, within the current model specification, the influence of these parameters cannot be determined independently with a high degree of certainty a point that is further investigated thru diagnostic analysis. Figure 1. Regression coefficients with 95% confidence intervals. Partial Effects of Horizontal Location Error The relationship between the most influential predictor horizontal location error and the composite seismic quality value is further investigated using a partial residual plot (Figure 2). After correcting for all other predictors, a clear negative linear trend is visible. As horizontal uncertainty increases, the partial residuals of the quality score systematically decrease, indicating a monotonic deterioration in data quality with decreasing location precision. This pattern provides strong visual support for the regression results and is consistent with theoretical and empirical expectations regarding the sensitivity of earthquake solutions to epicentral accuracy (Husen and Hardebeck 2010). Figure 2. Partial residual plot for normalized horizontal location error Model Diagnostics, Robustness, and Limitations A comprehensive range of diagnostic tests and plots were used to evaluate the validity of the regression assumptions (Fox and Weisberg 2018). The distribution of the model residuals (Figure 3) is approximately symmetrical and bell-shaped, which corresponds to the nonsignificant Shapiro-Wilk test for normality (𝑊 = 0.96,𝑝 = 0.084). This indicates that the assumption of normally distributed errors is reasonably met. Figure 3. Histogram of regression model residuals Further diagnostic plots (Figure 4) provide additional evidence for the robustness of the model. The residuals versus the predicted values show no systematic pattern, supporting the assumptions of linearity and homoscedasticity, which is confirmed by a non-significant Breusch-Pagan test (𝑝 = 0.410). The normal Q-Q plot shows only small deviations in the tails, while the scale-location plot indicates an approximately constant variance across the predicted values. The residuals versus leverage plot does not show any influential observations, as all points lie within the Cook's distance contours. The independence of the residuals is confirmed by the Durbin-Watson test (𝐷𝑊 = 2.12,𝑝 = 0.64). 22 2024-0108 20:48:42 4.9225 126.1575 62.6 6.7 93 km SE of Sarangan i, Philippin es 138 22 23 2024-0101 07:18:41 37.1895 136.8272 10.0 6.2 8 km SW of Anamizu , Japan 185 32 24 2024-0101 07:10:09 37.4874 137.2710 10.0 7.5 2024 Noto Peninsul a, Japan Earthqua ke 282 36 25 2023-1230 17:16:23 -2.9934 139.3720 33.0 6.3 146 km WSW of Abepura, Indonesi a 190 15 26 2023-1228 09:15:16 44.5960 149.0388 31.0 6.5 115 km SE of Kuril’sk, Russia 156 34 27 2023-1222 17:36:32 -52.0879 27.9516 10.0 6.1 South of Africa 284 19 28 2023-1221 14:55:56 51.2097 -175.3313 20.0 6.1 116 km SE of Adak, Alaska 126 52 29 2023-1220 12:11:21 -15.8585 -72.5207 93.0 6.2 11 km E of Iray, Peru 184 70 30 2023-1211 06:33:32 -18.8139 -175.6541 248.0 6.2 174 km NW of Fangale’ ounga, Tonga 233 47 31 2023-1207 12:56:30 -20.6152 169.3089 48.0 7.1 118 km S of Isangel, Vanuatu 146 20 32 2023-1203 19:49:35 8.9705 126.6121 20.0 6.9 34 km ENE of Arasasan, Philippin es 151 25 33 2023-1203 14:35:56 8.7320 126.8230 13.0 6.0 58 km E of Marihata g, Philippin es 88 76 34 2023-1203 10:35:52 8.4880 126.7454 19.0 6.6 47 km ENE of Hinatuan , Philippin es 135 32 35 2023-1202 20:52:15 8.4475 126.7782 9.0 6.0 49 km NE of Barcelon a, Philippin es 176 42 36 2023-1202 18:09:25 8.4402 126.9387 46.4 6.3 63 km ENE of Barcelon a, Philippin es 225 15 37 2023-1202 17:40:15 8.4230 126.7220 36.1 6.1 41 km NE of Barcelon a, Philippin es 275 13 38 2023-1202 16:03:41 8.4377 126.7570 35.0 6.4 47 km NE of Barcelon a, Philippin es 125 65 39 2023-1202 14:37:04 8.5266 126.4161 40.0 7.6 19 km E of Gamut, Philippin es 128 43 40 2023-1127 21:46:42 -3.5605 144.0313 10.0 6.5 44 km E of Wewak, Papua New Guinea 248 15 41 2023-1124 09:05:03 20.1319 145.5195 22.2 6.9 Maug Islands region, 328 33 Northern Mariana Islands 42 2023-1122 04:47:31 -14.9603 167.9718 13.0 6.7 97 km E of PortOlry, Vanuatu 120 23 43 2023-1122 02:48:51 1.7827 127.1887 102.0 6.0 91 km W of Tobelo, Indonesi a 132 23 44 2023-1117 08:14:10 5.5627 125.0013 52.0 6.7 32 km SW of Kablalan , Philippin es 244 21 45 2023-1114 07:00:56 -4.0330 87.1044 9.0 6.1 South Indian Ocean 132 54 46 2023-1113 07:43:34 -3.8855 151.0670 10.0 6.1 126 km WNW of Rabaul, Papua New Guinea 120 50 47 2023-1110 20:45:11 -6.0976 130.0606 10.0 6.1 Banda Sea 186 37 48 2023-1108 13:02:06 -6.1310 129.8738 10.0 6.7 Banda Sea 179 28 49 2023-1108 04:53:49 -6.4160 129.5466 6.0 7.1 Banda Sea 142 53 50 2023-1108 04:52:51 -6.4442 129.7518 10.0 6.7 Banda Sea 106 29 51 2023-1101 21:04:48 -10.0561 123.7541 51.0 6.1 20 km NE of Kupang, Indonesi a 258 13 Note: The complete dataset includes 51 global earthquakes from November 2023 to March 2024, obtained from the USGS earthquake catalog. All events have magnitude ≥6.0 and were reviewed by the USGS. Additional parameters not shown in this table include: distance to nearest station (dmin), root mean square error (rms), horizontal error, depth error, magnitude error, and magnitude type details. Supplementary 2. Code R Analysis Scripts # ========================================== # REGRESSION ANALYSIS # ========================================== # Install and load required packages if (!require("pacman")) install.packages("pacman") pacman::p_load(tidyverse, broom, car, lmtest, moments, pROC, ggplot2, ggpubr, gridExtra) # ========================================== # 1. LOAD AND PREPARE DATA # ========================================== cat("1. DATA PREPARATION\n") cat("==========================================\n\n") data <- read.csv("/query (2).csv") # Convert and prepare data data <- data %>% mutate( depth = as.numeric(depth), mag = as.numeric(mag), nst = as.numeric(nst), gap = as.numeric(gap), dmin = as.numeric(dmin), rms = as.numeric(rms), horizontalError = as.numeric(horizontalError), magError = as.numeric(magError), # Create composite seismic score seismic_score = scale(mag) + scale(1/(depth + 1)) + scale(nst/100), # Normalized variables gap_normalized = gap / 100, horizontalError_normalized = horizontalError / 10, # For classification model major_earthquake = ifelse(mag >= 6.5, 1, 0) ) # Check for categorical variables with only one level cat("Checking data structure...\n") cat("Total observations:", nrow(data), "\n") cat("Major earthquakes (M ≥ 6.5):", sum(data$major_earthquake), "\n") cat("Minor earthquakes (M < 6.5):", sum(data$major_earthquake == 0), "\n\n") # ========================================== # 2. MODEL 1 # ========================================== cat("2. LINEAR REGRESSION MODEL\n") cat("======================================\n\n") # Fit the model model_significant <- lm(seismic_score ~ gap_normalized + log(dmin + 0.1) + rms + horizontalError_normalized + I(gap * rms / 1000), data = data) # ========================================== # 2.1 MODEL SUMMARY # ========================================== cat("MODEL SUMMARY\n") cat("-------------\n") model_summary <- summary(model_significant) print(model_summary) # Extract key statistics cat("\nKEY STATISTICS:\n") cat("F-statistic:", round(model_summary$fstatistic[1], 3), "\n") cat("p-value:", format.pval(pf(model_summary$fstatistic[1], model_summary$fstatistic[2], model_summary$fstatistic[3], lower.tail = FALSE), digits = 4), "\n") cat("R-squared:", round(model_summary$r.squared, 4), "\n") cat("Adjusted R-squared:", round(model_summary$adj.r.squared, 4), "\n") cat("Residual standard error:", round(model_summary$sigma, 4), "\n") cat("Observations:", nrow(data), "\n\n") # ========================================== # 2.2 ANOVA TABLE # ========================================== cat("ANOVA TABLE\n") cat("-----------\n") anova_table <- anova(model_significant) print(anova_table) cat("\n") # ========================================== # 2.3 REGRESSION COEFFICIENTS # ========================================== cat("REGRESSION COEFFICIENTS WITH 95% CONFIDENCE INTERVALS\n") cat("-----------------------------------------------------\n\n") coef_table <- broom::tidy(model_significant, conf.int = TRUE, conf.level = 0.95) coef_table_display <- coef_table %>% mutate( term = case_when( term == "(Intercept)" ~ "Intercept", term == "gap_normalized" ~ "Normalized Azimuthal Gap", term == "log(dmin + 0.1)" ~ "log(Distance to Nearest Station + 0.1)", term == "rms" ~ "Root Mean Square Error", term == "horizontalError_normalized" ~ "Normalized Horizontal Error", term == "I(gap * rms/1000)" ~ "Gap × RMS Interaction / 1000", TRUE ~ term ), estimate = round(estimate, 4), std.error = round(std.error, 4), statistic = round(statistic, 3), p.value = ifelse(p.value < 0.001, "<0.001", ifelse(p.value < 0.01, sprintf("%.3f", p.value), sprintf("%.3f", p.value))), conf.low = round(conf.low, 4), conf.high = round(conf.high, 4), significance = case_when( p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", p.value < 0.1 ~ ".", TRUE ~ "" ) ) %>% select(term, estimate, std.error, statistic, p.value, significance, conf.low, conf.high) print(coef_table_display) cat("\nSignificance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n") # ========================================== # 2.4 MODEL DIAGNOSTICS # ========================================== cat("MODEL DIAGNOSTICS\n") cat("================\n\n") residuals <- residuals(model_significant) fitted <- fitted(model_significant) # 2.4.1 Normality Test cat("1. NORMALITY OF RESIDUALS:\n") shapiro_test <- shapiro.test(residuals) cat(" Shapiro-Wilk test: W =", round(shapiro_test$statistic, 4), ", p =", format.pval(shapiro_test$p.value, digits = 4), "\n") cat(" Interpretation:", ifelse(shapiro_test$p.value > 0.05, "Residuals are normally distributed", "Residuals are not normally distributed"), "\n\n") # 2.4.2 Heteroscedasticity Test cat("2. HETEROSCEDASTICITY TEST:\n") bp_test <- lmtest::bptest(model_significant) cat(" Breusch-Pagan test: BP =", round(bp_test$statistic, 4), ", p =", format.pval(bp_test$p.value, digits = 4), "\n") cat(" Interpretation:", ifelse(bp_test$p.value > 0.05, "No evidence of heteroscedasticity", "Evidence of heteroscedasticity"), "\n\n") # 2.4.3 Multicollinearity cat("3. MULTICOLLINEARITY (VIF):\n") vif_values <- car::vif(model_significant) vif_table <- data.frame( Variable = names(vif_values), VIF = round(vif_values, 2) ) print(vif_table) cat("\n Interpretation: VIF < 5 indicates low multicollinearity\n\n") # 2.4.4 Autocorrelation cat("4. AUTOCORRELATION TEST:\n") dw_test <- lmtest::dwtest(model_significant) cat(" Durbin-Watson test: DW =", round(dw_test$statistic, 4), ", p =", format.pval(dw_test$p.value, digits = 4), "\n") cat(" Interpretation:", ifelse(dw_test$p.value > 0.05, "No evidence of autocorrelation", "Evidence of autocorrelation"), "\n\n") # ========================================== # 3. VISUALIZATIONS # ========================================== cat("3. VISUALIZATIONS\n") cat("================\n\n") # Create a directory for plots if (!dir.exists("figures")) dir.create("figures") # ========================================== # 3.1 Diagnostic Plots (4-in-1) # ========================================== cat("Creating diagnostic plots...\n") png("figures/diagnostic_plots.png", width = 10, height = 8, units = "in", res = 300) par(mfrow = c(2, 2), mar = c(4, 4, 3, 2), family = "sans") # Residuals vs Fitted plot(fitted, residuals, xlab = "Fitted Values", ylab = "Residuals", main = "(A) Residuals vs Fitted", pch = 19, col = "steelblue", cex = 0.8) abline(h = 0, col = "red", lty = 2, lwd = 2) lines(lowess(fitted, residuals), col = "darkgreen", lwd = 2) # Normal Q-Q Plot qqnorm(residuals, main = "(B) Normal Q-Q Plot", pch = 19, col = "steelblue", cex = 0.8) qqline(residuals, col = "red", lwd = 2) # Scale-Location Plot plot(fitted, sqrt(abs(residuals)), xlab = "Fitted Values", ylab = expression(sqrt("|Residuals|")), main = "(C) Scale-Location Plot", pch = 19, col = "steelblue", cex = 0.8) lines(lowess(fitted, sqrt(abs(residuals))), col = "darkgreen", lwd = 2) # Residuals vs Leverage plot(model_significant, which = 5, main = "(D) Residuals vs Leverage", pch = 19, col = "steelblue", cex = 0.8) dev.off() cat("✓ Diagnostic plots saved: figures/diagnostic_plots.png\n") # ========================================== # 3.2 Actual vs Predicted Plot # ========================================== cat("Creating actual vs predicted plot...\n") actual_vs_pred <- data.frame( Actual = data$seismic_score, Predicted = fitted, Magnitude = data$mag ) p1 <- ggplot(actual_vs_pred, aes(x = Predicted, y = Actual)) + geom_point(aes(color = Magnitude), size = 3, alpha = 0.7) + geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed", size = 1) + geom_smooth(method = "lm", se = TRUE, color = "blue", alpha = 0.2) + labs(title = "Actual vs Predicted Composite Seismic Score", subtitle = paste("R² =", round(model_summary$r.squared, 3), ", p =", format.pval(pf(model_summary$fstatistic[1], model_summary$fstatistic[2], model_summary$fstatistic[3], lower.tail = FALSE), digits = 3)), x = "Predicted Score", y = "Actual Score") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), plot.subtitle = element_text(size = 12), legend.position = "right") + scale_color_gradient(low = "blue", high = "red", name = "Magnitude") ggsave("figures/actual_vs_predicted.png", p1, width = 8, height = 6, dpi = 300) cat("✓ Actual vs predicted plot saved: figures/actual_vs_predicted.png\n") # ========================================== # 3.3 Coefficient Plot with Confidence Intervals # ========================================== cat("Creating coefficient plot...\n") coef_for_plot <- coef_table %>% filter(term != "(Intercept)") %>% mutate( term = case_when( term == "gap_normalized" ~ "Normalized Azimuthal Gap", term == "log(dmin + 0.1)" ~ "log(Distance to Station)", term == "rms" ~ "RMS Error", term == "horizontalError_normalized" ~ "Normalized Horizontal Error", term == "I(gap * rms/1000)" ~ "Gap × RMS Interaction", TRUE ~ term ), term = factor(term, levels = rev(unique(term))), significant = p.value < 0.05 ) p2 <- ggplot(coef_for_plot, aes(x = estimate, y = term)) + geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") + geom_errorbarh(aes(xmin = conf.low, xmax = conf.high, color = significant), height = 0.2, size = 1) + geom_point(aes(color = significant), size = 3) + labs(title = "Regression Coefficients with 95% CI", x = "Coefficient Estimate", y = "Predictor") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), axis.text = element_text(size = 10), legend.position = "none") + scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) ggsave("figures/coefficient_plot.png", p2, width = 8, height = 5, dpi = 300) cat("✓ Coefficient plot saved: figures/coefficient_plot.png\n") # ========================================== # 3.4 Residual Distribution Plot # ========================================== cat("Creating residual distribution plot...\n") resid_data <- data.frame(Residuals = residuals) p3 <- ggplot(resid_data, aes(x = Residuals)) + geom_histogram(aes(y = ..density..), bins = 15, fill = "steelblue", alpha = 0.7, color = "white") + geom_density(color = "red", size = 1) + stat_function(fun = dnorm, args = list(mean = mean(residuals), sd = sd(residuals)), color = "darkgreen", size = 1, linetype = "dashed") + labs(title = "Distribution of Residuals", subtitle = paste("Shapiro-Wilk p =", format.pval(shapiro_test$p.value, digits = 3)), x = "Residuals", y = "Density") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), plot.subtitle = element_text(size = 12)) ggsave("figures/residual_distribution.png", p3, width = 8, height = 5, dpi = 300) cat("✓ Residual distribution plot saved: figures/residual_distribution.png\n") cat("Minor earthquakes (M < 6.5):", sum(data$major_earthquake == 0), "\n\n") # ========================================== # 2. MODEL 1 # ========================================== cat("2. LINEAR REGRESSION MODEL\n") cat("======================================\n\n") # Fit the model model_significant <- lm(seismic_score ~ gap_normalized + log(dmin + 0.1) + rms + horizontalError_normalized + I(gap * rms / 1000), data = data) # ========================================== # 2.1 MODEL SUMMARY # ========================================== cat("MODEL SUMMARY\n") cat("-------------\n") model_summary <- summary(model_significant) print(model_summary) # Extract key statistics cat("\nKEY STATISTICS:\n") cat("F-statistic:", round(model_summary$fstatistic[1], 3), "\n") cat("p-value:", format.pval(pf(model_summary$fstatistic[1], model_summary$fstatistic[2], model_summary$fstatistic[3], lower.tail = FALSE), digits = 4), "\n") cat("R-squared:", round(model_summary$r.squared, 4), "\n") cat("Adjusted R-squared:", round(model_summary$adj.r.squared, 4), "\n") cat("Residual standard error:", round(model_summary$sigma, 4), "\n") cat("Observations:", nrow(data), "\n\n") # ========================================== # 2.2 ANOVA TABLE # ========================================== cat("ANOVA TABLE\n") cat("-----------\n") anova_table <- anova(model_significant) print(anova_table) cat("\n") # ========================================== # 2.3 REGRESSION COEFFICIENTS # ========================================== cat("REGRESSION COEFFICIENTS WITH 95% CONFIDENCE INTERVALS\n") cat("-----------------------------------------------------\n\n") coef_table <- broom::tidy(model_significant, conf.int = TRUE, conf.level = 0.95) coef_table_display <- coef_table %>% mutate( term = case_when( term == "(Intercept)" ~ "Intercept", term == "gap_normalized" ~ "Normalized Azimuthal Gap", term == "log(dmin + 0.1)" ~ "log(Distance to Nearest Station + 0.1)", term == "rms" ~ "Root Mean Square Error", term == "horizontalError_normalized" ~ "Normalized Horizontal Error", term == "I(gap * rms/1000)" ~ "Gap × RMS Interaction / 1000", TRUE ~ term ), estimate = round(estimate, 4), std.error = round(std.error, 4), statistic = round(statistic, 3), p.value = ifelse(p.value < 0.001, "<0.001", ifelse(p.value < 0.01, sprintf("%.3f", p.value), sprintf("%.3f", p.value))), conf.low = round(conf.low, 4), conf.high = round(conf.high, 4), significance = case_when( p.value < 0.001 ~ "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", p.value < 0.1 ~ ".", TRUE ~ "" ) ) %>% select(term, estimate, std.error, statistic, p.value, significance, conf.low, conf.high) print(coef_table_display) cat("\nSignificance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n") # ========================================== # 2.4 MODEL DIAGNOSTICS # ========================================== cat("MODEL DIAGNOSTICS\n") cat("================\n\n") residuals <- residuals(model_significant) fitted <- fitted(model_significant) # 2.4.1 Normality Test cat("1. NORMALITY OF RESIDUALS:\n") shapiro_test <- shapiro.test(residuals) cat(" Shapiro-Wilk test: W =", round(shapiro_test$statistic, 4), ", p =", format.pval(shapiro_test$p.value, digits = 4), "\n") cat(" Interpretation:", ifelse(shapiro_test$p.value > 0.05, "Residuals are normally distributed", "Residuals are not normally distributed"), "\n\n") # 2.4.2 Heteroscedasticity Test cat("2. HETEROSCEDASTICITY TEST:\n") bp_test <- lmtest::bptest(model_significant) cat(" Breusch-Pagan test: BP =", round(bp_test$statistic, 4), ", p =", format.pval(bp_test$p.value, digits = 4), "\n") cat(" Interpretation:", ifelse(bp_test$p.value > 0.05, "No evidence of heteroscedasticity", "Evidence of heteroscedasticity"), "\n\n") # 2.4.3 Multicollinearity cat("3. MULTICOLLINEARITY (VIF):\n") vif_values <- car::vif(model_significant) vif_table <- data.frame( Variable = names(vif_values), VIF = round(vif_values, 2) ) print(vif_table) cat("\n Interpretation: VIF < 5 indicates low multicollinearity\n\n") # 2.4.4 Autocorrelation cat("4. AUTOCORRELATION TEST:\n") dw_test <- lmtest::dwtest(model_significant) cat(" Durbin-Watson test: DW =", round(dw_test$statistic, 4), ", p =", format.pval(dw_test$p.value, digits = 4), "\n") cat(" Interpretation:", ifelse(dw_test$p.value > 0.05, "No evidence of autocorrelation", "Evidence of autocorrelation"), "\n\n") # ========================================== # 3. VISUALIZATIONS # ========================================== cat("3. VISUALIZATIONS\n") cat("================\n\n") # Create a directory for plots if (!dir.exists("figures")) dir.create("figures") # ========================================== # 3.1 Diagnostic Plots (4-in-1) # ========================================== cat("Creating diagnostic plots...\n") png("figures/diagnostic_plots.png", width = 10, height = 8, units = "in", res = 300) par(mfrow = c(2, 2), mar = c(4, 4, 3, 2), family = "sans") # Residuals vs Fitted plot(fitted, residuals, xlab = "Fitted Values", ylab = "Residuals", main = "(A) Residuals vs Fitted", pch = 19, col = "steelblue", cex = 0.8) abline(h = 0, col = "red", lty = 2, lwd = 2) lines(lowess(fitted, residuals), col = "darkgreen", lwd = 2) # Normal Q-Q Plot qqnorm(residuals, main = "(B) Normal Q-Q Plot", pch = 19, col = "steelblue", cex = 0.8) qqline(residuals, col = "red", lwd = 2) # Scale-Location Plot plot(fitted, sqrt(abs(residuals)), xlab = "Fitted Values", ylab = expression(sqrt("|Residuals|")), main = "(C) Scale-Location Plot", pch = 19, col = "steelblue", cex = 0.8) lines(lowess(fitted, sqrt(abs(residuals))), col = "darkgreen", lwd = 2) # Residuals vs Leverage plot(model_significant, which = 5, main = "(D) Residuals vs Leverage", pch = 19, col = "steelblue", cex = 0.8) dev.off() cat("✓ Diagnostic plots saved: figures/diagnostic_plots.png\n") # ========================================== # 3.2 Actual vs Predicted Plot # ========================================== cat("Creating actual vs predicted plot...\n") actual_vs_pred <- data.frame( Actual = data$seismic_score, Predicted = fitted, Magnitude = data$mag ) p1 <- ggplot(actual_vs_pred, aes(x = Predicted, y = Actual)) + geom_point(aes(color = Magnitude), size = 3, alpha = 0.7) + geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed", size = 1) + geom_smooth(method = "lm", se = TRUE, color = "blue", alpha = 0.2) + labs(title = "Actual vs Predicted Composite Seismic Score", subtitle = paste("R² =", round(model_summary$r.squared, 3), ", p =", format.pval(pf(model_summary$fstatistic[1], model_summary$fstatistic[2], model_summary$fstatistic[3], lower.tail = FALSE), digits = 3)), x = "Predicted Score", y = "Actual Score") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), plot.subtitle = element_text(size = 12), legend.position = "right") + scale_color_gradient(low = "blue", high = "red", name = "Magnitude") ggsave("figures/actual_vs_predicted.png", p1, width = 8, height = 6, dpi = 300) cat("✓ Actual vs predicted plot saved: figures/actual_vs_predicted.png\n") # ========================================== # 3.3 Coefficient Plot with Confidence Intervals # ========================================== cat("Creating coefficient plot...\n") coef_for_plot <- coef_table %>% filter(term != "(Intercept)") %>% mutate( term = case_when( term == "gap_normalized" ~ "Normalized Azimuthal Gap", term == "log(dmin + 0.1)" ~ "log(Distance to Station)", term == "rms" ~ "RMS Error", term == "horizontalError_normalized" ~ "Normalized Horizontal Error", term == "I(gap * rms/1000)" ~ "Gap × RMS Interaction", TRUE ~ term ), term = factor(term, levels = rev(unique(term))), significant = p.value < 0.05 ) p2 <- ggplot(coef_for_plot, aes(x = estimate, y = term)) + geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") + geom_errorbarh(aes(xmin = conf.low, xmax = conf.high, color = significant), height = 0.2, size = 1) + geom_point(aes(color = significant), size = 3) + labs(title = "Regression Coefficients with 95% CI", x = "Coefficient Estimate", y = "Predictor") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), axis.text = element_text(size = 10), legend.position = "none") + scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) ggsave("figures/coefficient_plot.png", p2, width = 8, height = 5, dpi = 300) cat("✓ Coefficient plot saved: figures/coefficient_plot.png\n") # ========================================== # 3.4 Residual Distribution Plot # ========================================== cat("Creating residual distribution plot...\n") resid_data <- data.frame(Residuals = residuals) p3 <- ggplot(resid_data, aes(x = Residuals)) + geom_histogram(aes(y = ..density..), bins = 15, fill = "steelblue", alpha = 0.7, color = "white") + geom_density(color = "red", size = 1) + stat_function(fun = dnorm, args = list(mean = mean(residuals), sd = sd(residuals)), color = "darkgreen", size = 1, linetype = "dashed") + labs(title = "Distribution of Residuals", subtitle = paste("Shapiro-Wilk p =", format.pval(shapiro_test$p.value, digits = 3)), x = "Residuals", y = "Density") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), plot.subtitle = element_text(size = 12)) ggsave("figures/residual_distribution.png", p3, width = 8, height = 5, dpi = 300) cat("✓ Residual distribution plot saved: figures/residual_distribution.png\n") # ========================================== # 4. SIMPLIFIED CLASSIFICATION MODEL # ========================================== cat("\n4. SIMPLIFIED CLASSIFICATION MODEL\n") cat("===================================\n\n") # Check if we have enough observations for classification if (sum(data$major_earthquake) >= 5 && sum(data$major_earthquake == 0) >= 5) { cat("Sufficient data for classification analysis.\n") cat("Major earthquakes (M ≥ 6.5):", sum(data$major_earthquake), "\n") cat("Minor earthquakes (M < 6.5):", sum(data$major_earthquake == 0), "\n\n") # Create simple numeric variables instead of factors data$depth_scaled <- scale(data$depth) data$gap_scaled <- scale(data$gap) data$nst_scaled <- scale(data$nst) # Fit logistic regression with numeric variables only logistic_model <- glm(major_earthquake ~ depth_scaled + gap_scaled + nst_scaled, data = data, family = binomial) cat("LOGISTIC REGRESSION SUMMARY\n") cat("---------------------------\n") logistic_summary <- summary(logistic_model) print(logistic_summary) # ROC Curve pred_probs <- predict(logistic_model, type = "response") roc_obj <- pROC::roc(data$major_earthquake, pred_probs) auc_value <- pROC::auc(roc_obj) cat("\nMODEL PERFORMANCE:\n") cat("AUC (Area Under ROC Curve):", round(auc_value, 3), "\n") # Confusion matrix (using optimal threshold) if (requireNamespace("pROC", quietly = TRUE)) { optimal_threshold <- pROC::coords(roc_obj, "best", ret = "threshold")$threshold[1] } else { optimal_threshold <- 0.5 } predictions <- ifelse(pred_probs > optimal_threshold, 1, 0) conf_matrix <- table(Predicted = predictions, Actual = data$major_earthquake) cat("\nConfusion Matrix (threshold =", round(optimal_threshold, 3), "):\n") print(conf_matrix) # Calculate metrics accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix) if (sum(conf_matrix[,2]) > 0) { sensitivity <- conf_matrix[2,2] / sum(conf_matrix[,2]) } else { sensitivity <- NA } if (sum(conf_matrix[,1]) > 0) { specificity <- conf_matrix[1,1] / sum(conf_matrix[,1]) } else { specificity <- NA } cat("\nPerformance Metrics:\n") cat("Accuracy:", round(accuracy, 3), "\n") if (!is.na(sensitivity)) cat("Sensitivity:", round(sensitivity, 3), "\n") if (!is.na(specificity)) cat("Specificity:", round(specificity, 3), "\n") # ========================================== # 4.1 ROC Curve Plot # ========================================== cat("\nCreating ROC curve plot...\n") # Create ROC curve data roc_data <- data.frame( fpr = 1 - roc_obj$specificities, tpr = roc_obj$sensitivities ) p4 <- ggplot(roc_data, aes(x = fpr, y = tpr)) + geom_line(color = "blue", size = 1) + geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray") + geom_polygon(aes(x = fpr), fill = "lightblue", alpha = 0.3) + labs(title = "ROC Curve for Earthquake Size Classification", subtitle = paste("AUC =", round(auc_value, 3), "(M ≥ 6.5 vs M < 6.5)"), x = "False Positive Rate (1 - Specificity)", y = "True Positive Rate (Sensitivity)") + theme_minimal() + theme(plot.title = element_text(face = "bold", size = 14), plot.subtitle = element_text(size = 12)) + coord_equal() ggsave("figures/roc_curve.png", p4, width = 6, height = 6, dpi = 300) cat("✓ ROC curve saved: figures/roc_curve.png\n") } else { cat("Insufficient data for reliable classification analysis.\n") cat("Major earthquakes (M ≥ 6.5):", sum(data$major_earthquake), "\n") cat("Minor earthquakes (M < 6.5):", sum(data$major_earthquake == 0), "\n") cat("Minimum 5 observations per class recommended.\n\n") } # ========================================== # 5. EXPORT TABLES # ========================================== cat("\n5. EXPORTING RESULTS\n") cat("===================\n\n") # Export regression coefficients write.csv(coef_table_display, "significant_regression_coefficients.csv", row.names = FALSE) cat("✓ Regression coefficients saved: significant_regression_coefficients.csv\n") # Export ANOVA table anova_df <- as.data.frame(anova_table) write.csv(anova_df, "significant_anova_table.csv", row.names = TRUE) cat("✓ ANOVA table saved: significant_anova_table.csv\n") # Export model summary model_stats <- data.frame( Statistic = c("R-squared", "Adjusted R-squared", "F-statistic", "p-value", "Residual SE", "Observations"), Value = c( round(model_summary$r.squared, 4), round(model_summary$adj.r.squared, 4), round(model_summary$fstatistic[1], 2), format.pval(pf(model_summary$fstatistic[1], model_summary$fstatistic[2], model_summary$fstatistic[3], lower.tail = FALSE), digits = 4), round(model_summary$sigma, 4), nrow(data) ) ) write.csv(model_stats, "significant_model_summary.csv", row.names = FALSE) cat("✓ Model summary saved: significant_model_summary.csv\n") # Export diagnostic test results diagnostics <- data.frame( Test = c("Shapiro-Wilk (Normality)", "Breusch-Pagan (Heteroscedasticity)", "Durbin-Watson (Autocorrelation)"), Statistic = c( round(shapiro_test$statistic, 4), round(bp_test$statistic, 4), round(dw_test$statistic, 4) ), p_value = c( format.pval(shapiro_test$p.value, digits = 4), format.pval(bp_test$p.value, digits = 4), format.pval(dw_test$p.value, digits = 4) ), Interpretation = c( ifelse(shapiro_test$p.value > 0.05, "Normal", "Not normal"), ifelse(bp_test$p.value > 0.05, "Homoscedastic", "Heteroscedastic"), ifelse(dw_test$p.value > 0.05, "No autocorrelation", "Autocorrelation") ) ) write.csv(diagnostics, "significant_model_diagnostics.csv", row.names = FALSE) cat("✓ Model diagnostics saved: significant_model_diagnostics.csv\n") # Export VIF values write.csv(vif_table, "significant_vif_values.csv", row.names = FALSE) cat("✓ VIF values saved: significant_vif_values.csv\n") # ========================================== # 6. FINAL SUMMARY # ========================================== cat("\n") cat(paste(rep("=", 60), collapse = ""), "\n") cat("SUMMARY\n") cat(paste(rep("=", 60), collapse = ""), "\n\n")