Spatial Analysis of Vegetation Indices Across Kakamega County

Important: scope and interpretation

Important data constraint: the CSV contains 12 sub-county summaries for NDVI, EVI and NDWI, plus polygon geometry and centroids. Accordingly:

• For exploratory spatial association/regression, the available pair is NDVI mean and EVI mean; this is not a rainfall-effect model.
• With only 12 areal units, local statistics, spatial regression and GWR have low power and can be unstable. Results should be treated as exploratory, not definitive evidence.

Research questions

1. Is NDVI spatially autocorrelated across Kakamega’s sub-counties? Global Moran’s I.
2. Where are the local clusters and outliers? LISA / Local Moran’s I and Getis–Ord Gi*.
3. Is NDVI associated with EVI across zones? Pearson correlation, Lee’s L, and join counts.
4. Does accounting for spatial dependence change the NDVI–EVI relationship? OLS compared with spatial lag and spatial error models.
5. Does the relationship vary geographically? Geographically Weighted Regression (GWR).
6. Can any of this be validated? Spatial cross-validation and independent ground-truth validation.

These questions are exploratory. The available dataset contains only 12 sub-county summaries and does not include independent field measurements.

1. Setup and load data

library(sf)
library(spdep)
library(spatialreg)

has_gwmodel <- requireNamespace("GWmodel", quietly = TRUE)
if (has_gwmodel) library(GWmodel)

library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)
library(jsonlite)
library(knitr)

data_file <- "NDVI_zonal_stats_2024 (1).csv"
raw <- read.csv(data_file, stringsAsFactors = FALSE, check.names = FALSE)

required <- c(
  "adm2_name", "NDVI_mean", "EVI_mean", "NDWI_mean",
  "NDVI_stdDev", "EVI_stdDev", "NDWI_stdDev", ".geo"
)

missing_cols <- setdiff(required, names(raw))
if (length(missing_cols)) {
  stop(
    "Missing required columns: ",
    paste(missing_cols, collapse = ", ")
  )
}

zones_tbl <- raw |>
  mutate(
    across(
      c(
        NDVI_mean, EVI_mean, NDWI_mean,
        NDVI_stdDev, EVI_stdDev, NDWI_stdDev,
        area_sqkm, center_lat, center_lon
      ),
      as.numeric
    )
  )

geo_text <- zones_tbl[[".geo"]]

geom <- tryCatch(
  sf::st_as_sfc(geo_text, GeoJSON = TRUE, crs = 4326),
  error = function(e) NULL
)

if (is.null(geom)) {
  if (requireNamespace("geojsonsf", quietly = TRUE)) {
    geom <- geojsonsf::geojson_sfc(geo_text, crs = 4326)
  } else {
    stop(
      "Could not parse .geo as GeoJSON. ",
      "Install geojsonsf and rerun."
    )
  }
}

zones <- st_sf(
  zones_tbl |> select(-all_of(".geo")),
  geometry = geom
)

zones <- st_make_valid(zones)
zones <- zones |> filter(!is.na(NDVI_mean), !is.na(EVI_mean))

if (nrow(zones) < 4) {
  stop("Too few complete zones for the planned exploratory analyses.")
}

zones_m <- st_transform(zones, 32736)

2. Data audit and descriptive summary

audit <- data.frame(
  Item = c(
    "Rows in source",
    "Zones retained",
    "Unique sub-county names",
    "Missing NDVI",
    "Missing EVI",
    "Missing NDWI",
    "Geometry validity after st_make_valid"
  ),
  Result = c(
    nrow(raw),
    nrow(zones),
    dplyr::n_distinct(zones$adm2_name),
    sum(is.na(zones$NDVI_mean)),
    sum(is.na(zones$EVI_mean)),
    sum(is.na(zones$NDWI_mean)),
    sum(st_is_valid(zones))
  )
)

kable(audit)

summary_indices <- zones |>
  st_drop_geometry() |>
  summarise(
    n_zones = n(),
    NDVI_mean_of_zones = mean(NDVI_mean, na.rm = TRUE),
    NDVI_median = median(NDVI_mean, na.rm = TRUE),
    SD_between_zone_NDVI_means = sd(NDVI_mean, na.rm = TRUE),
    EVI_mean_of_zones = mean(EVI_mean, na.rm = TRUE),
    NDWI_mean_of_zones = mean(NDWI_mean, na.rm = TRUE),
    NDVI_min = min(NDVI_mean, na.rm = TRUE),
    NDVI_max = max(NDVI_mean, na.rm = TRUE)
  )

kable(round(summary_indices, 4))

Interpretation: Across the 12 sub-counties, mean NDVI was 0.6613 (median 0.6575; between-zone SD 0.0375). Values ranged from 0.5891 in Likuyani to 0.7519 in Shinyalu.

The county-level mean reported here is unweighted: each sub-county contributes equally regardless of its area or valid pixel count. A zonal standard deviation describes dispersion among underlying pixel values; it is not the standard error of a sub-county mean.

3. Map and rank NDVI

ggplot(zones) +
  geom_sf(
    aes(fill = NDVI_mean),
    colour = "white",
    linewidth = 0.25
  ) +
  scale_fill_viridis_c(
    option = "C",
    name = "Mean NDVI"
  ) +
  labs(
    title = "Mean NDVI by sub-county",
    caption = "Unweighted sub-county means."
  ) +
  theme_minimal(base_size = 13)

ndvi_rank <- zones |>
  st_drop_geometry() |>
  arrange(NDVI_mean) |>
  select(adm2_name, NDVI_mean, NDVI_stdDev, area_sqkm)

kable(
  mutate(
    ndvi_rank,
    across(where(is.numeric), ~round(.x, 4))
  )
)

4. Construct spatial weights

nb <- poly2nb(zones, queen = TRUE)
n_neighbours <- card(nb)
isolates <- which(n_neighbours == 0)

weights <- nb2listw(
  nb,
  style = "W",
  zero.policy = TRUE
)

weight_summary <- data.frame(
  Measure = c(
    "Zones",
    "Zones with no Queen neighbors",
    "Minimum neighbors",
    "Median neighbors",
    "Maximum neighbors"
  ),
  Result = c(
    length(nb),
    length(isolates),
    min(n_neighbours),
    median(n_neighbours),
    max(n_neighbours)
  )
)

kable(weight_summary)

Interpretation: All 12 sub-counties had at least one Queen neighbour; there were no isolates. The median number of neighbours was 3, with a minimum of 1 and maximum of 7.

Queen contiguity defines neighbours as polygons sharing either an edge or a vertex. Results may change under alternative neighbour definitions.

5. Global spatial autocorrelation: Moran’s I

moran_ndvi <- moran.test(
  zones$NDVI_mean,
  weights,
  zero.policy = TRUE,
  alternative = "two.sided"
)

moran_out <- data.frame(
  Statistic = c(
    "Moran's I",
    "Expected I under randomisation",
    "p-value"
  ),
  Value = c(
    unname(moran_ndvi$estimate[["Moran I statistic"]]),
    unname(moran_ndvi$estimate[["Expectation"]]),
    moran_ndvi$p.value
  )
)

kable(moran_out, digits = 5)

Interpretation: Global Moran’s I for mean NDVI was 0.21986, compared with an expected value of -0.09091 under spatial randomisation. The two-sided p-value was 0.05011.

This indicates weak-to-moderate positive spatial autocorrelation: sub-counties with similar NDVI values tended to be geographically close. However, because (p = 0.05011), the result was not statistically significant at the conventional 5% threshold. It should therefore be reported as borderline or suggestive evidence, not conclusive evidence of global spatial clustering.

6. Local clusters and outliers: LISA

local_res <- localmoran(
  zones$NDVI_mean,
  weights,
  zero.policy = TRUE
)

zones$local_I <- as.numeric(local_res[, "Ii"])
zones$local_p <- as.numeric(local_res[, "Pr(z != E(Ii))"])

x <- zones$NDVI_mean - mean(zones$NDVI_mean, na.rm = TRUE)
lag_x <- lag.listw(weights, x, zero.policy = TRUE)

zones$lisa_type <- case_when(
  x >= 0 & lag_x >= 0 ~ "High-High",
  x < 0 & lag_x < 0 ~ "Low-Low",
  x >= 0 & lag_x < 0 ~ "High-Low outlier",
  x < 0 & lag_x >= 0 ~ "Low-High outlier",
  TRUE ~ "Unclassified"
)

zones$lisa_significant <- zones$local_p < 0.05
zones$lisa_label <- ifelse(
  zones$lisa_significant,
  zones$lisa_type,
  "Not significant"
)

ggplot(zones) +
  geom_sf(
    aes(fill = lisa_label),
    colour = "white",
    linewidth = 0.25
  ) +
  scale_fill_manual(
    values = c(
      "High-High" = "#d73027",
      "Low-Low" = "#4575b4",
      "High-Low outlier" = "#fc8d59",
      "Low-High outlier" = "#91bfdb",
      "Not significant" = "grey85"
    ),
    name = "LISA class"
  ) +
  labs(
    title = "LISA cluster/outlier map for NDVI",
    subtitle = "Nominal p < 0.05; no multiple-testing adjustment applied"
  ) +
  theme_minimal(base_size = 12)

lisa_table <- zones |>
  st_drop_geometry() |>
  select(
    adm2_name,
    local_I,
    local_p,
    lisa_type,
    lisa_significant
  ) |>
  arrange(local_p)

kable(
  mutate(
    lisa_table,
    across(where(is.numeric), ~round(.x, 4))
  )
)

7. Hotspot analysis: Getis–Ord Gi*

gi <- localG(
  zones$NDVI_mean,
  weights,
  zero.policy = TRUE
)

zones$gi_star_z <- as.numeric(gi)

zones$gi_class <- case_when(
  zones$gi_star_z >= 1.96 ~ "Hot spot (z >= 1.96)",
  zones$gi_star_z <= -1.96 ~ "Cold spot (z <= -1.96)",
  TRUE ~ "Not significant at 5%"
)

ggplot(zones) +
  geom_sf(
    aes(fill = gi_class),
    colour = "white",
    linewidth = 0.25
  ) +
  scale_fill_manual(
    values = c(
      "Hot spot (z >= 1.96)" = "#d73027",
      "Cold spot (z <= -1.96)" = "#4575b4",
      "Not significant at 5%" = "grey85"
    ),
    name = "Gi* class"
  ) +
  labs(
    title = "Getis–Ord local G hotspot/coldspot map",
    subtitle = "Approximate two-sided 5% z-score cut-off: ±1.96"
  ) +
  theme_minimal(base_size = 12)

8. Spatial association: NDVI and EVI

The supplied CSV has no rain_mean. Therefore, the association analysis uses NDVI mean and EVI mean. This is a relationship between two remotely sensed vegetation indices, not evidence of causation or a rainfall effect.

8.1 Pearson correlation

pair <- zones |>
  st_drop_geometry() |>
  filter(is.finite(NDVI_mean), is.finite(EVI_mean))

pearson <- cor.test(
  pair$NDVI_mean,
  pair$EVI_mean,
  method = "pearson"
)

pearson_ci <- as.numeric(pearson$conf.int)
pearson_ci_lower <- head(pearson_ci, 1)
pearson_ci_upper <- tail(pearson_ci, 1)

pearson_table <- data.frame(
  Quantity = c(
    "Pearson r",
    "95% CI lower",
    "95% CI upper",
    "p-value",
    "n zones"
  ),
  Value = c(
    unname(pearson$estimate),
    pearson_ci_lower,
    pearson_ci_upper,
    pearson$p.value,
    nrow(pair)
  )
)

kable(
  mutate(
    pearson_table,
    Value = round(Value, 5)
  )
)

Interpretation: There was a very strong positive linear association between sub-county mean NDVI and mean EVI: Pearson (r = 0.95150), 95% CI 0.83182–0.98663, (p < 0.00001), (n = 12).

Sub-counties with higher EVI generally had higher NDVI. This correlation does not establish causality, and the two indices may share spectral inputs or respond to similar vegetation and land-cover conditions. The conventional confidence interval and p-value also assume independent observations, which may not hold for spatially arranged zones.

8.2 Lee’s L spatial bivariate association

lee_result <- tryCatch(
  lee.test(
    pair$NDVI_mean,
    pair$EVI_mean,
    listw = weights,
    zero.policy = TRUE
  ),
  error = function(e) e
)

if (inherits(lee_result, "error")) {
  cat("Lee's L could not be computed:", lee_result$message)
} else {
  print(lee_result)
}

    Lee's L statistic randomisation

data:  pair$NDVI_mean ,  pair$EVI_mean 
weights: weights  

Lee's L statistic standard deviate = -0.34143, p-value = 0.6336
alternative hypothesis: greater
sample estimates:
Lee's L statistic       Expectation          Variance 
       0.21241588        0.26320591        0.02212874 

Interpretation: Lee’s L was 0.21242, with an expectation of 0.26321 under the randomisation null. The one-sided test gave (p = 0.6336).

Therefore, after accounting for the spatial arrangement and autocorrelation of the two indices, there was no statistically significant bivariate spatial association between NDVI and EVI.

This does not contradict the very high Pearson correlation. Pearson’s r measures the overall linear relationship between the 12 zonal values, whereas Lee’s L assesses whether the two variables show similar spatial co-patterning after incorporating spatial dependence.

8.3 Join-count analysis

zones$ndvi_class <- factor(
  ifelse(
    zones$NDVI_mean > median(zones$NDVI_mean, na.rm = TRUE),
    "High",
    "Low"
  ),
  levels = c("High", "Low")
)

join_result <- tryCatch(
  joincount.test(
    zones$ndvi_class,
    weights,
    zero.policy = TRUE
  ),
  error = function(e) e
)

if (inherits(join_result, "error")) {
  cat(
    "Join-count test could not be computed:",
    join_result$message
  )
} else {
  print(join_result)
}

    Join count test under nonfree sampling

data:  zones$ndvi_class 
weights: weights 

Std. deviate for High = 2.7819, p-value = 0.002702
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            2.3571429             1.3636364             0.1275401 


    Join count test under nonfree sampling

data:  zones$ndvi_class 
weights: weights 

Std. deviate for Low = 2.5286, p-value = 0.005726
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
            2.2666667             1.3636364             0.1275401 

Interpretation: Using a median split, both High–High and Low–Low neighbouring joins occurred significantly more often than expected under the null hypothesis.

• High–High joins: observed 2.357 versus 1.364 expected; standardised deviate 2.7819; (p = 0.002702).
• Low–Low joins: observed 2.267 versus 1.364 expected; standardised deviate 2.5286; (p = 0.005726).

Thus, high-NDVI sub-counties tended to border other high-NDVI sub-counties, and low-NDVI sub-counties tended to border other low-NDVI sub-counties. This provides categorical evidence of spatial clustering, even though the continuous global Moran’s I test was only borderline.

The median split is data-driven and discards information from the continuous NDVI values, so this result is exploratory.

8.4 NDVI–EVI scatterplot

ggplot(
  pair,
  aes(
    x = NDVI_mean,
    y = EVI_mean,
    label = adm2_name
  )
) +
  geom_point(
    size = 3,
    colour = "#2166ac"
  ) +
  geom_text(
    nudge_y = 0.003,
    check_overlap = TRUE,
    size = 3.2
  ) +
  geom_smooth(
    method = "lm",
    se = TRUE,
    colour = "#b2182b"
  ) +
  labs(
    title = "Association between sub-county mean NDVI and EVI",
    subtitle = "Pearson r = 0.95150; 95% CI 0.83182–0.98663",
    x = "Mean NDVI",
    y = "Mean EVI"
  ) +
  theme_minimal(base_size = 13)

9. Spatial regression modelling

9.1 OLS baseline

ols <- lm(
  NDVI_mean ~ EVI_mean,
  data = zones
)

kable(
  as.data.frame(coef(summary(ols))),
  digits = 5
)

Interpretation: The OLS model was:

[ = 0.12021 + 1.23847 . ]

A one-unit increase in sub-county mean EVI was associated with a 1.23847-unit increase in mean NDVI ((p < 0.00001)). The intercept was 0.12021 ((p = 0.05532)).

The global OLS model explained 90.53% of the between-sub-county variation in NDVI ((R^2 = 0.9053); adjusted (R^2 = 0.8959)). This is an associational, exploratory index-to-index model, not a causal or rainfall-effect model.

9.2 Residual spatial autocorrelation

zones$ols_resid <- residuals(ols)

moran_resid <- moran.test(
  zones$ols_resid,
  weights,
  zero.policy = TRUE
)

print(moran_resid)

    Moran I test under randomisation

data:  zones$ols_resid  
weights: weights    

Moran I statistic standard deviate = 0.87172, p-value = 0.1917
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.06828930       -0.09090909        0.03335210 

Interpretation: Moran’s I for OLS residuals was 0.06829 ((p = 0.1917)). This was not statistically significant, indicating that the OLS residuals did not show detectable spatial autocorrelation under the Queen-contiguity weights.

Therefore, there was no strong statistical evidence that OLS left substantial spatial structure in the residuals.

9.3 LM diagnostics

lm_diag <- tryCatch(
  lm.RStests(
    ols,
    weights,
    test = c(
      "RSlag",
      "RSerr",
      "adjRSlag",
      "adjRSerr"
    ),
    zero.policy = TRUE
  ),
  error = function(e) e
)

if (inherits(lm_diag, "error")) {
  cat(
    "LM diagnostics failed:",
    lm_diag$message
  )
} else {
  print(lm_diag)
}

    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = NDVI_mean ~ EVI_mean, data = zones)
test weights: weights

RSlag = 0.10105, df = 1, p-value = 0.7506


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = NDVI_mean ~ EVI_mean, data = zones)
test weights: weights

RSerr = 0.090759, df = 1, p-value = 0.7632


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = NDVI_mean ~ EVI_mean, data = zones)
test weights: weights

adjRSlag = 0.37444, df = 1, p-value = 0.5406


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = NDVI_mean ~ EVI_mean, data = zones)
test weights: weights

adjRSerr = 0.36415, df = 1, p-value = 0.5462

Interpretation: Neither robust diagnostic was significant:

• Robust spatial-lag test: (p = 0.5406).
• Robust spatial-error test: (p = 0.5462).

Therefore, the diagnostics did not support a spatial-lag or spatial-error specification. The simpler OLS model was statistically adequate based on these residual diagnostics.

9.4 Spatial lag and error models

lag_model <- tryCatch(
  lagsarlm(
    NDVI_mean ~ EVI_mean,
    data = zones,
    listw = weights,
    zero.policy = TRUE
  ),
  error = function(e) e
)

error_model <- tryCatch(
  errorsarlm(
    NDVI_mean ~ EVI_mean,
    data = zones,
    listw = weights,
    zero.policy = TRUE
  ),
  error = function(e) e
)

if (!inherits(lag_model, "error")) {
  cat("SPATIAL LAG MODEL\n")
  print(summary(lag_model))
} else {
  cat(
    "Spatial lag model failed:",
    lag_model$message,
    "\n"
  )
}

SPATIAL LAG MODEL

Call:lagsarlm(formula = NDVI_mean ~ EVI_mean, data = zones, listw = weights, 
    zero.policy = TRUE)

Residuals:
        Min          1Q      Median          3Q         Max 
-0.01465650 -0.00769188 -0.00096299  0.00472519  0.02488147 

Type: lag 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.15326    0.11874  1.2907   0.1968
EVI_mean     1.25841    0.14082  8.9362   <2e-16

Rho: -0.062964, LR test value: 0.092898, p-value: 0.76052
Asymptotic standard error: 0.21517
    z-value: -0.29262, p-value: 0.76981
Wald statistic: 0.085629, p-value: 0.76981

Log likelihood: 37.10212 for lag model
ML residual variance (sigma squared): 0.00012064, (sigma: 0.010984)
Number of observations: 12 
Number of parameters estimated: 4 
AIC: -66.204, (AIC for lm: -68.111)
LM test for residual autocorrelation
test value: 0.34049, p-value: 0.55955

if (!inherits(error_model, "error")) {
  cat("SPATIAL ERROR MODEL\n")
  print(summary(error_model))
} else {
  cat(
    "Spatial error model failed:",
    error_model$message,
    "\n"
  )
}

SPATIAL ERROR MODEL

Call:errorsarlm(formula = NDVI_mean ~ EVI_mean, data = zones, listw = weights, 
    zero.policy = TRUE)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0136622 -0.0091167 -0.0017879  0.0056070  0.0252391 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.113455   0.052484  2.1617  0.03064
EVI_mean    1.254428   0.119944 10.4585  < 2e-16

Lambda: 0.13907, LR test value: 0.11491, p-value: 0.73462
Asymptotic standard error: 0.35056
    z-value: 0.39672, p-value: 0.69157
Wald statistic: 0.15739, p-value: 0.69157

Log likelihood: 37.11312 for error model
ML residual variance (sigma squared): 0.00011989, (sigma: 0.010949)
Number of observations: 12 
Number of parameters estimated: 4 
AIC: -66.226, (AIC for lm: -68.111)

Interpretation: In the spatial-lag model, the autoregressive coefficient was (= -0.06296) ((p = 0.76981)), indicating no evidence that a sub-county’s NDVI depended on neighbouring NDVI after adjusting for EVI.

In the spatial-error model, (= 0.13907) ((p = 0.69157)), indicating no evidence of spatial dependence in the regression errors after adjusting for EVI.

EVI remained strongly and positively associated with NDVI in both spatial models:

• Spatial lag model: EVI coefficient 1.25841 ((p < 0.00001)).
• Spatial error model: EVI coefficient 1.25443 ((p < 0.00001)).

These spatial terms should not be interpreted as physical spillovers or as evidence of a specific omitted variable, particularly with only 12 observations.

9.5 Model comparison

Interpretation: OLS had the lowest AIC (-68.111), followed by the spatial-error model (-66.226) and spatial-lag model (-66.204).

The spatial models had essentially equivalent AIC values and were about 1.9 AIC units higher than OLS. Together with the non-significant robust LM tests, non-significant () and (), and non-significant residual Moran’s I, OLS is the most defensible global model for these data.

9.6 Residual checks for spatial models

 Spatial lag residuals 

    Moran I test under randomisation

data:  as.numeric(residuals(lag_model))  
weights: weights    

Moran I statistic standard deviate = 1.086, p-value = 0.1387
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.10800096       -0.09090909        0.03354823 


 Spatial error residuals 

    Moran I test under randomisation

data:  as.numeric(residuals(error_model))  
weights: weights    

Moran I statistic standard deviate = 0.50133, p-value = 0.3081
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.0007484675     -0.0909090909      0.0334263209 

Interpretation: Residual Moran’s I was non-significant for both the spatial-lag model ((I = 0.10800), (p = 0.1387)) and the spatial-error model ((I = 0.00075), (p = 0.3081)).

Thus, neither spatial model left significant residual spatial autocorrelation. However, this does not make either spatial model preferable to OLS, because the spatial dependence terms themselves were non-significant and AIC favoured OLS.

10. Geographically Weighted Regression

Adaptive bandwidth (number of nearest neighbours): 15 AICc value: -67.65311 
Adaptive bandwidth (number of nearest neighbours): 17 AICc value: -67.34994 
Adaptive bandwidth (number of nearest neighbours): 14 AICc value: -67.62102 
Adaptive bandwidth (number of nearest neighbours): 15 AICc value: -67.65311 

Selected adaptive bandwidth: 15 
   ***********************************************************************
   *                       Package   GWmodel                             *
   ***********************************************************************
   Program starts at: 2026-10-15 20:24:23.120889 
   Call:
   gwr.basic(formula = NDVI_mean ~ EVI_mean, data = zones_sp, bw = gwr_bw, 
    kernel = "bisquare", adaptive = TRUE)

   Dependent (y) variable:  NDVI_mean
   Independent variables:  EVI_mean
   Number of data points: 12
   ***********************************************************************
   *                    Results of Global Regression                     *
   ***********************************************************************

   Call:
    lm(formula = formula, data = data)

   Residuals:
      Min        1Q    Median        3Q       Max 
-0.013996 -0.008471 -0.001614  0.005401  0.025398 

   Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
   (Intercept)  0.12021    0.05544   2.168   0.0553 .  
   EVI_mean     1.23847    0.12664   9.780 1.95e-06 ***

   ---Significance stars
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
   Residual standard error: 0.01209 on 10 degrees of freedom
   Multiple R-squared: 0.9053
   Adjusted R-squared: 0.8959 
   F-statistic: 95.65 on 1 and 10 DF,  p-value: 1.949e-06 
   ***Extra Diagnostic information
   Residual sum of squares: 0.00146059
   Sigma(hat): 0.01208549
   AIC:  -68.11133
   AICc:  -65.11133
   BIC:  -71.20189
   ***********************************************************************
   *          Results of Geographically Weighted Regression              *
   ***********************************************************************

   *********************Model calibration information*********************
   Kernel function: bisquare 
   Adaptive bandwidth: 15 (number of nearest neighbours)
   Regression points: the same locations as observations are used.
   Distance metric: Euclidean distance metric is used.

   ****************Summary of GWR coefficient estimates:******************
                 Min.  1st Qu.   Median  3rd Qu.   Max.
   Intercept 0.017142 0.029474 0.037755 0.052166 0.1441
   EVI_mean  1.192015 1.391906 1.421371 1.439572 1.4705
   ************************Diagnostic information*************************
   Number of data points: 12 
   Effective number of parameters (2trace(S) - trace(S'S)): 3.402485 
   Effective degrees of freedom (n-2trace(S) + trace(S'S)): 8.597515 
   AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): -67.65311 
   AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): -77.98694 
   BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): -85.66297 
   Residual sum of squares: 0.0008296191 
   R-square value:  0.946235 
   Adjusted R-square value:  0.9221568 

   ***********************************************************************
   Program stops at: 2026-10-15 20:24:23.127988 

Interpretation: The selected adaptive bandwidth was 15 nearest neighbours. Local EVI coefficients ranged from 1.19202 to 1.47050, with a median of 1.42137. This narrow range suggests that the NDVI–EVI relationship was broadly similar across Kakamega’s sub-counties, rather than showing strong geographical variation.

GWR produced a higher in-sample (R^2) (0.9462) than global OLS (0.9053) and a lower AICc (-67.65311 versus -65.11133). However, these are in-sample measures.

With only 12 spatial units, an adaptive bandwidth of 15 neighbours, and no out-of-sample validation, the GWR results should be treated as exploratory and potentially unstable. They do not justify strong claims about locally varying NDVI–EVI relationships.

11. Validation and accuracy assessment

11.1 What can be validated?

The supplied CSV contains no independent ground-truth vegetation measurements and no separate test dataset. Therefore, it cannot support a genuine external validation claim.

Fitted values compared with the same observations used to fit a model are in-sample fit, not independent predictive validation.

11.2 Spatial cross-validation

Spatial cross-validation should hold out spatial blocks or groups of neighbouring zones to reduce leakage from spatially autocorrelated observations.

With only 12 sub-counties, five spatial folds may leave too few zones for stable model fitting. No cross-validated performance estimate should be reported from the current data.

11.3 In-sample fit metrics

Interpretation: The OLS model had an in-sample RMSE of 0.01103 NDVI units and an in-sample MAE of 0.00877 NDVI units.

These values describe how closely the fitted line reproduced the 12 observed sub-county means. They are not predictive accuracy estimates and must not be presented as validation accuracy.

11.4 Ground-truth validation

No field_validation_points.shp was supplied.  Ground-truth validation cannot be performed from the current files.

Interpretation: No field_validation_points.shp was supplied. Therefore, ground-truth validation could not be performed.

A defensible validation dataset would require independent field observations with coordinates, dates, a clearly defined vegetation measurement, and matching spatial and temporal scales.

12. Overall conclusions

• Mean NDVI varied moderately across Kakamega’s 12 sub-counties, from 0.5891 in Likuyani to 0.7519 in Shinyalu.
• Global Moran’s I was 0.21986 ((p = 0.05011)), providing borderline, non-significant evidence of positive global spatial autocorrelation.
• Ikolomani was the only significant local NDVI cluster: it was a significant High–High LISA cluster and a significant Getis–Ord Gi* hotspot.
• Join-count analysis found significant clustering of both high- and low-NDVI sub-counties.
• NDVI and EVI were very strongly positively correlated ((r = 0.95150), (p < 0.00001)), but Lee’s L did not show significant bivariate spatial co-patterning ((L = 0.21242), (p = 0.6336)).
• OLS indicated that EVI was strongly associated with NDVI ((= 1.23847), (p < 0.00001); (R^2 = 0.9053)).
• Residual Moran’s I and robust LM diagnostics did not support spatial lag or spatial error terms; OLS had the lowest AIC and was the preferred global model.
• GWR suggested limited geographical variation in the EVI coefficient, but the result is exploratory and unreliable with only 12 zones.
• No independent predictive validation or ground-truth validation was possible with the supplied data.

13. Reproducibility and limitations

The analysis depends on correct geometry parsing and the Queen-contiguity neighbour graph. Before publication, document the satellite product, acquisition or composite dates, NDWI band formula, cloud/shadow masking, spatial resolution, valid-pixel counts, and zonal aggregation method.

The present analysis is cross-sectional and based on only 12 sub-county summaries. It has no rainfall or other environmental covariates, no independent field data, and limited statistical power for local statistics, spatial regression, GWR, and cross-validation. All inferential results should therefore be described as exploratory.