Spatial Analysis of Vegetation Indices Across Kakamega County

Important: scope and interpretation

This presentation implements the spatial-analysis stages in the supplied workflow using the uploaded zonal-statistics CSV.

Important data constraint: the CSV contains 12 sub-county summaries for NDVI, EVI and NDWI, plus polygon geometry and centroids. It does not contain rain_mean, repeated dates, independent field measurements, or documented index-processing metadata. Accordingly:

• For exploratory spatial association/regression, the available pair is NDVI mean and EVI mean; this is not a rainfall-effect model.
• Spatial statistics can be calculated only if the geometry parses correctly and the neighbor graph is usable.
• Ground-truth validation and a genuine spatial predictive-accuracy assessment cannot be completed without independent observations and a prespecified prediction target.
• 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.

1. Setup and load data

2. Data audit and descriptive summary

The mean of the sub-county means is unweighted: each sub-county contributes equally regardless of area or valid pixel count. It is not automatically the area-weighted county-wide mean. A zonal standard deviation describes dispersion in the underlying pixel values as calculated by the original workflow; it is not the standard error of a sub-county mean.

3. Map and rank NDVI

4. Construct spatial weights

Interpretation: Queen contiguity is a defensible starting point for administrative polygons, but the results depend on the neighbor definition. Review the neighbor graph and any isolated polygons before interpreting statistics. zero.policy = TRUE permits isolates to be handled computationally; it does not make their spatial information equivalent to connected zones.

5. Global spatial autocorrelation — Moran’s I

Interpretation rule: Moran’s I above its null expectation, together with a small p-value, supports positive global spatial autocorrelation (similar NDVI values tend to be near one another). A significantly negative value supports neighboring dissimilarity. A non-significant result means the analysis does not provide sufficient evidence against the null; it does not prove spatial randomness. Report the statistic, p-value, neighbor definition, and the small sample size together.

Caution: the conventional test’s inference depends on its assumptions and the spatial weights. With 12 zones, treat the p-value cautiously.

6. Local clusters and outliers — LISA / Local Moran’s I

7. Hotspot analysis — Getis–Ord Gi*

8. Spatial association: NDVI and EVI

The supplied CSV has no rain_mean. Therefore, the following association analysis uses NDVI mean and EVI mean, both available for the same zones. This is a relationship between two derived spectral indices, not evidence of a causal relationship or a rainfall effect.

8.1 Ordinary Pearson correlation (baseline)

Interpretation: Pearson’s r describes the strength and direction of a linear relationship among the 12 zonal summaries. The confidence interval and p-value rely on assumptions that may be questionable for spatially dependent zones. Report it as an unadjusted baseline, not the only inferential result.

8.2 Spatially corrected bivariate association — Lee’s L

    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 combines information about the two variables and their spatial lags. Interpret the reported estimate and p-value using the exact output from spdep; do not treat Lee’s L as interchangeable with Pearson’s r. If the test fails, do not substitute an invented result: check package version, missing values, weight-list alignment and documentation. With only 12 zones, any spatially corrected inference is fragile.

8.3 Categorical spatial association — join counts

    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: join counts test whether same-category neighboring joins (such as High-High or Low-Low) depart from what the null model expects. The median split is a data-driven categorisation for exploration and discards information from the continuous NDVI values. State the category rule and inspect the returned tests rather than reporting only a map. The exact output depends on the spdep version and weight graph.

8.4 Visualize the bivariate relationship

9. Spatial regression modeling

9.1 OLS as a diagnostic baseline

Because rainfall is absent, the example model is NDVI_mean ~ EVI_mean. This is an exploratory index-to-index regression, not a model of rainfall effects. Do not choose covariates simply because they are available; the scientific question should drive the model.

The EVI coefficient is the estimated average change in zonal mean NDVI associated with a one-unit difference in zonal mean EVI under the linear model. It is not causal. Conventional OLS standard errors assume an appropriate error structure and can be unreliable if residuals are spatially autocorrelated.

9.2 Test OLS residuals for spatial autocorrelation

    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 

If residual Moran’s I is significant, OLS leaves spatial structure unexplained and its usual inferential results may be unreliable. A non-significant result does not prove all assumptions are satisfied, particularly with only 12 zones.

9.3 LM diagnostics: spatial lag vs spatial error

    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

Decision guide: compare the robust lag and robust error tests, not only the unadjusted tests. Robust lag significance supports considering a spatial-lag model; robust error significance supports considering a spatial-error model. If both or neither are significant, theory, diagnostics, sensitivity to weights and sample limitations matter. Do not select a model solely by a mechanical rule.

9.4 Fit spatial lag and spatial error models

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

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: the spatial-lag model includes a spatially lagged outcome term, often reported as rho (ρ), to represent dependence on neighboring outcomes under the model. The spatial-error model represents spatial dependence in the unobserved error process, often summarized by lambda (λ). These terms are not proof of a physical causal spillover or a named omitted variable. With 12 observations, coefficient estimates can be unstable.

9.5 Compare AIC and select cautiously

Lower AIC indicates a better trade-off between fit and model complexity among comparable models fitted to the same response and observations. AIC alone does not prove a model is correct. Use residual diagnostics, the scientific question, sample size and sensitivity to spatial weights alongside AIC. Do not call the first-ranked model “final” without those checks.

9.6 Re-test spatial structure in the selected spatial-model residuals

 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 

A non-significant residual Moran’s I is a desirable diagnostic, but not a guarantee of a correct model. A significant residual result suggests remaining spatial structure; review the weights, covariates, functional form and data quality.

10. Geographically Weighted Regression (GWR)

GWR estimates locally varying coefficients. It is data-hungry and may be unstable with only 12 polygons. The following is included as a diagnostic/exploratory workflow only; do not publish local coefficients as reliable evidence without adequate spatial units, bandwidth checks, diagnostics and validation.

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 18:57:36.498946 
   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 18:57:36.503305 

If GWR runs, examine the local EVI coefficients, local standard errors, local collinearity diagnostics, effective number of parameters and corrected AIC. Spatially varying coefficients are not automatically meaningful differences; local multiple testing and bandwidth choice can produce unstable patterns. Compare GWR to the global model and validate out-of-sample before making local claims.

11. Validation and accuracy assessment

11.1 What can be validated with the supplied file?

The supplied CSV has 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 — requirements and limitation

Spatial cross-validation should hold out spatial blocks or groups of neighboring zones so nearby observations do not leak across training and test sets. Block size should be justified by the spatial scale of dependence and the study geography, not chosen arbitrarily.

With only 12 sub-counties, five spatial folds may leave too few zones for stable model fitting. Do not report cross-validated performance until a valid fold design and an explicit prediction workflow have been specified. The example below creates spatial blocks as a design aid; it does not train or validate a model.

blockCV is installed. Design and inspect folds with the installed blockCV version before running cross-validation.
No CV score is calculated in this presentation because the dataset has only 12 zones and no pre-specified prediction workflow.

11.3 Accuracy metrics for a continuous outcome

For a continuous outcome such as NDVI, common metrics include RMSE and MAE, calculated from held-out predictions and observed outcomes. In-sample values are optimistic and should not be labelled predictive accuracy.

These values describe the OLS fit to the same 12 observations used to estimate the model. They are not independent accuracy estimates and should not be used to claim predictive performance. For a final model, calculate RMSE/MAE from predictions for held-out spatial blocks and compare them to a transparent baseline.

11.4 Ground-truth validation using independent field data

A defensible validation dataset should contain independent field observations with coordinates, observation dates, a clearly defined measured vegetation variable, and a matching spatial/temporal scale. Extract or match the corresponding satellite-derived index at those locations, then report sample size, bias, MAE/RMSE, correlation and uncertainty as appropriate.

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

Do not use the zonal NDVI summaries as their own ground truth. Correlating a model-derived value with itself or with another product that shares the same inputs does not provide independent validation.

12. Conclusions and reporting checklist

• Report Moran’s I with its value, p-value, Queen-contiguity definition and sample size.
• Report LISA class and local p-value; account for multiple testing.
• Report Gi* z-scores and the chosen significance threshold, while acknowledging approximation and multiple testing.
• Report Pearson association as a baseline and Lee’s L as a distinct spatial statistic; state that the available pair is NDVI–EVI, not NDVI–rainfall.
• Use join counts only with the stated categorical split and a clear description of which joins are tested.
• Treat OLS as a diagnostic baseline; use residual Moran’s I and robust LM diagnostics before considering spatial regression.
• Compare spatial lag/error models using diagnostics and AIC together, not AIC alone.
• Treat GWR as exploratory with only 12 zones.
• Do not claim spatial cross-validation or ground-truth accuracy until independent held-out/field observations are available.

13. Reproducibility and limitations

The analysis depends on the geometry parsing correctly and on the chosen Queen-contiguity neighbor graph. Confirm that the polygons are the intended administrative boundaries, that no unexpected overlaps/gaps affect contiguity, and that all index values were generated from comparable imagery and processing rules.

Before publication, document the satellite product, acquisition/composite dates, band formula for NDWI, cloud/shadow mask, resolution, valid-pixel counts, and zonal aggregation method. The present dataset is cross-sectional and has only 12 sub-county summaries; inference and especially complex spatial models should therefore be described as exploratory.