10/5 - Extended Data Analysis

Resources

• Handout

• Data

• Jen Spreadsheet & John Copy

• Colab Code

Last class you folks used data from 1 location to determine if precipitation and snow cover was trending from 1980 to 2026.  We used a simple approach to build your understanding from the ground up.  However, there are issues with using 1 location and using a linear regression for a trend analysis - namely:

A. Autocorrelation (Serial Correlation)

OLS regression assumes that all observational errors (residuals) are independent of one another. However, weather and hydrological data are sequentially ordered in time and often display autocorrelation—where the conditions in year t are correlated with year t-1 or t-2 due to persistent climate states (e.g., multi-year ocean warmth or groundwater storage memory). When autocorrelation is present, standard OLS formulas underestimate the true standard error of the slope. This artificially inflates t-statistics and yields deceptively small P-values, causing researchers to declare trends "statistically significant" when they are not.

B. Non-Stationarity

Standard linear regression assumes a stationary process where the mean, variance, and underlying data-generating distribution remain constant across the record. Under active global climate change, atmospheric mechanisms and baseline statistical distributions are actively shifting. A straight linear model cannot capture non-linear accelerations, sudden regime shifts, or tipping points.

C. Sensitivity to Extreme Outliers

Because OLS regression minimizes the squared vertical errors, a single extreme atmospheric river event (e.g., 1997) or an unprecedented drought year (e.g., 2001) near the start or end of a 40-year record can exert disproportionate leverage on the slope, heavily altering the fitted trendline.

D. Single Station issues & Climate Oscillations

Observational data from a single station like Big Red Mountain contain multiple overlaid signals operating across different temporal and spatial scales. What may look like random scatter on a 40-year chart is often driven by known natural climate oscillations.

Using Google Earthe Engine

To streamline this process we can use Google Earth Engine.  This way we have access to datasets while never having to download any data.  We will use python and Google Colab to access GEE.  Java script is natively supported by GEE.  You can also use R inside of Colab.  I have chosen to use Python because it is more generalized than R and allows the most efficient gathering and processing of data.

--- Linear Regression Results for WY_Prcp_mm ---

                            OLS Regression Results                            

==============================================================================

Dep. Variable:             WY_Prcp_mm   R-squared:                       0.003

Model:                            OLS   Adj. R-squared:                 -0.019

Method:                 Least Squares   F-statistic:                    0.1487

Date:                Sun, 04 Oct 2026   Prob (F-statistic):              0.702

Time:                        18:43:35   Log-Likelihood:                -279.29

No. Observations:                  46   AIC:                             562.6

Df Residuals:                      44   BIC:                             566.2

Df Model:                           1                                         

Covariance Type:            nonrobust                                         

==============================================================================

                 coef    std err          t      P>|t|      [0.025      0.975]

------------------------------------------------------------------------------

const       1508.4867   2384.381      0.633      0.530   -3296.917    6313.890

Water_Year    -0.4592      1.191     -0.386      0.702      -2.859       1.940

==============================================================================

Omnibus:                        2.696   Durbin-Watson:                   1.743

Prob(Omnibus):                  0.260   Jarque-Bera (JB):                1.528

Skew:                           0.125   Prob(JB):                        0.466

Kurtosis:                       2.143   Cond. No.                     3.02e+05

==============================================================================

Notes:

[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

[2] The condition number is large, 3.02e+05. This might indicate that there are

strong multicollinearity or other numerical problems.

--- Linear Regression Results for Apr1_SWE_kg_m2 ---

                            OLS Regression Results                            

==============================================================================

Dep. Variable:         Apr1_SWE_kg_m2   R-squared:                       0.318

Model:                            OLS   Adj. R-squared:                  0.302

Method:                 Least Squares   F-statistic:                     20.50

Date:                Sun, 04 Oct 2026   Prob (F-statistic):           4.51e-05

Time:                        18:43:35   Log-Likelihood:                -220.84

No. Observations:                  46   AIC:                             445.7

Df Residuals:                      44   BIC:                             449.3

Df Model:                           1                                         

Covariance Type:            nonrobust                                         

==============================================================================

                 coef    std err          t      P>|t|      [0.025      0.975]

------------------------------------------------------------------------------

const       3068.5214    669.251      4.585      0.000    1719.735    4417.308

Water_Year    -1.5131      0.334     -4.527      0.000      -2.187      -0.840

==============================================================================

Omnibus:                        9.719   Durbin-Watson:                   1.149

Prob(Omnibus):                  0.008   Jarque-Bera (JB):                9.028

Skew:                           1.009   Prob(JB):                       0.0110

Kurtosis:                       3.799   Cond. No.                     3.02e+05

==============================================================================

Notes:

[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

[2] The condition number is large, 3.02e+05. This might indicate that there are

strong multicollinearity or other numerical problems. 

=================== GAM MODEL SUMMARY ===================

Family: Tweedie(p=1.407) 

Link function: log 

Formula:

swe_apr1_kg_m2 ~ s(wy_prcp_mm) + s(elevation_m) + s(latitude) + 

    s(longitude) + s(winter_tmax_c)

Parametric coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept)  1.21754    0.09219   13.21   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:

                   edf Ref.df       F p-value    

s(wy_prcp_mm)    4.188  5.200   3.264 0.00446 ** 

s(elevation_m)   5.703  6.674  48.329 < 2e-16 ***

s(latitude)      8.423  8.894  21.760 < 2e-16 ***

s(longitude)     7.099  8.093   6.047 < 2e-16 ***

s(winter_tmax_c) 5.981  6.728 200.757 < 2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.745   Deviance explained = 73.1%

-REML =  11742  Scale est. = 12.248    n = 4600

=================== MODEL COMPARISON (AIC) ===================

Tweedie GAM AIC:  23410.21 

OLS Linear Model AIC:  54456.9 

AIC Difference:  31046.69  (Positive means GAM performs better) 

Plots of the smooth terms for each predictor from model_swe_gam have been generated. Each plot shows the estimated smooth function (the curve) and its 95% confidence interval (shaded area).

A flat curve close to zero suggests a weak or non-existent effect. A steeply sloping or highly curved line indicates a strong, non-linear relationship. The y-axis represents the contribution to the log(SWE) prediction.

Plots of the smooth terms on the actual SWE (response) scale have been generated above.

The y-axis now represents the partial effect on April 1st SWE in kg/m², making the practical significance more intuitive.

The Goal of a GAM

In a standard linear regression, you try to draw a straight line through your data. But in the real world, relationships (like how elevation affects snowpack, or how temperature affects plants) are rarely straight lines—they bend, curve, and plateau.

A Generalized Additive Model (GAM) fits a smooth, flexible curve to your data instead of a straight line. It does this by combining several simpler curves, called basis functions (or splines).

Left Plot: Step 1 – Building the Curve (Behind the Scenes)

Think of this plot as the "construction site" where the final curved trend is built.

• The Dashed Lines (Weighted Bases 1, 2, and 3): These are the raw building blocks. Individually, each represents a localized "bump" or shape at a specific part of your x-axis. In this example, the model has adjusted their heights and directions (multiplying them by a parameter, $\beta$$\beta$):

  • Basis 1 is pulled downward.

  • Basis 2 is pushed upward.

  • Basis 3 is pulled strongly downward.

• The Solid Red Line (Summed Smooth): This is what happens when you stack all three dashed shapes on top of each other. By adding up the values of the dashed lines at every point along the x-axis, you get a single, continuous, custom-curved line ($s(x)$$s(x)$). This is how GAMs capture complex, non-linear shapes mathematically.

Right Plot: Step 2 – Centered Partial Effect & Data Points (What we analyze)

This represents the final output of your model superimposed on your actual dataset.

• The Solid Dark Red Line (Centered Partial Effect): This is the finished, summed curve from the left plot, centered so its average is at $0$$0$. This line represents the "purified" effect of this specific variable. For example, as your covariate increases from $-3$$-3$ to $-1.5$$-1.5$, the effect drops; as it moves from $-1.5$$-1.5$ to $0$$0$, the effect rises sharply.

• The Red Shaded Area (95% Confidence Interval): This is the model's margin of error. Notice that the red band is narrowest in the center (where we have a lot of data/certainty) and gets wider at the far left and right edges (where data is sparse, so the model is less certain about the exact path of the curve).

• The Black Dots (Partial Residuals): These are your actual data points. They are called "partial" residuals because the influences of all other variables in the model have been subtracted out, allowing you to see the raw relationship between this specific variable and your response.

• The Blue Dots at the Bottom (Exact Zeros): This is a crucial feature of ecological data like Snow Water Equivalent (SWE). Many times, there is absolutely no snow on the ground ($0$$0$ mm of SWE). Because we analyze this data on a logarithmic scale, a value of $0$$0$ cannot be plotted directly (since $\log(0)$$\log(0)$ is mathematically undefined/negative infinity). To keep this valuable information from being thrown away, the model has grouped all the "true zero" data points together at the bottom of the graph (around $-15$$-15$). This helps you see where and how often the environmental conditions resulted in absolutely zero snow.

=================== CLIMATE ATTRIBUTION GAM SUMMARY ===================

Family: Tweedie(p=1.397) 

Link function: log 

Formula:

swe_apr1_kg_m2 ~ s(wy_prcp_mm, k = 5) + s(winter_tmax_c, k = 5) + 

    s(elevation_m, k = 5) + s(latitude, longitude, k = 15) + 

    s(Winter_PDO, k = 5) + s(Winter_ENSO_ONI, k = 5)

Parametric coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept)  1.15030    0.09189   12.52   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:

                         edf Ref.df       F p-value    

s(wy_prcp_mm)          2.869  3.374   2.179 0.06010 .  

s(winter_tmax_c)       3.740  3.941 299.437 < 2e-16 ***

s(elevation_m)         3.749  3.951 122.177 < 2e-16 ***

s(latitude,longitude) 13.488 13.953  24.536 < 2e-16 ***

s(Winter_PDO)          3.970  3.999  32.058 < 2e-16 ***

s(Winter_ENSO_ONI)     1.000  1.001  10.607 0.00113 ** 

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.714   Deviance explained = 74.9% 

Final Project Summary: Blue Mountains SWE Climate Attribution

All tasks and requested refinements have been successfully executed and organized within your notebook:

1. Data Caching & Integration: The monthly teleconnection datasets (PDO, ENSO ONI) are downloaded, aggregated to winter water years, merged with both regional and point-scale datasets, and stored in Google Drive. On all future runs, the files load instantaneously from cache, preventing unnecessary NOAA or Earth Engine queries.

2. Regression Diagnostic Validation: We validated the multiple regression model's assumptions. The VIF test returned factors near $1.0$$1.0$, proving that multicollinearity is not an issue (the standard statsmodels warnings were confirmed as a numerical artifact from the raw scale of the Water_Year term). The Breusch-Pagan test ($p \approx 0.66$$p \approx 0.66$) and visual analysis confirmed complete homoscedasticity.

3. Rigorous Point-Scale Modeling (GAM): We fit a highly advanced point-scale Tweedie Generalized Additive Model (GAM) using R's mgcv package. This model accounted for zero-inflated snow observations and revealed highly significant non-linear thresholds for winter temperatures, elevations, and winter climate anomalies (PDO and ENSO), explaining 74.9% of the deviance in SWE.