<jupyter_start><jupyter_text>Tutorial: Logistic Regression for Factor Mapping This tutorial replicates a scenario discovery analysis performed in [Hadjimichael et al. (2020)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020EF001503). BackgroundPlanners in the the Upper Colorado River Basin (UCRB, shown in the figure below) are seeking to understand the vulnerability of water users to uncertainties stemming from climate change, population growth and water policy changes. The UCRB spans 25,682 km2 in western Colorado and is home to approximately 300,000 residents and 1,012 km2 of irrigated land. Several thousand irrigation ditches divert water from the main river and its tributaties for irrigation (shown as small black dots in the figure). Transmountain diversions of approximately 567,4000,000 m3 per year are exported for irrigation, industrial and municipal uses in northern and eastern Colorado, serving the major population centers of Denver and Colorado Springs. These diversions are carried through tunnels, shown as large black dots in the figure.An important planning consideration is the water rights of each user, defined by seniority across all water uses (irrigation diversions, transboundary diversions, power plants etc.) in the basin. To assess the vulnerability of users with varying degrees of water rights seniority, planners simulate the system across an ensemble of scenarios using the state of Colorado's StateMod platform. The model simulates streamflow, diversions, instream demands, and reservoir operations.  Focusing on decision-relevant metrics, the scenario discovery is applied to the water shortages experienced by each individual user (i.e., not on a single basin-wide or sector-wide metric). For this training example, we'll be performing scenario discovery for three different water users, two irrigation users and one municipal user.Before we start our analysis, we'll load the relevant Python libraries and create lists storing the names of the users of interest.<jupyter_code>import pandas as pd
from utils import *
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

all_IDs = ['7000550','7200799','3704614'] # IDs for three that we will perform the analysis for
usernames = ['Medium seniority irrigation',
             'Low seniority irrigation',
             'Transbasin municipal diversion']
nStructures = len(all_IDs)<jupyter_output><empty_output><jupyter_text>Step 1:  Load Latin Hypercube Sample and set up problemTo examine regional vulnerability, we generate an ensemble of plausible future states of the worlds (SOWs) using Latin Hypercube Sampling. For this tutorial, we'll load a file containing 1,000 parameter samples.<jupyter_code>LHsamples = np.loadtxt('./data/LHsamples_original_1000.txt')
param_bounds=np.loadtxt('./data/uncertain_params_bounds.txt', usecols=(1,2))

param_names=['Irrigation demand multiplier','Reservoir loss','Transbasin demand multiplier',
             'Municipal & industrial multiplier', 'Shoshone','Environmental flows',
             'Evaporation change','Mean dry flow','Dry flow variance',
             'Mean wet flow','Wet flow variance','Dry-dry probability',
             'Wet-wet probability', 'Snowmelt shift']<jupyter_output><empty_output><jupyter_text>Step 2: Define decision-relevant metrics for illustrationScenario discovery attempts to identify parametric regions that lead to 'success' and 'failure'. For this demonstration we'll be defining 'success' as states of the world where a shortage level doesn't exceed its historical frequency.  Step 3: Run the logistic regressionLogistic regression estimates the probability that a future SOW will be classified as a success or failure given a set of performance criteria. A logistic regression model is defined by:$$ln \bigg (\frac{p_i}{1-p_i} \bigg ) = X^T_i \beta$$where $p_i$ is the probability the performance in the $i^{th}$ SOW will be classified as a success, $X_i$ is the vector of covariates describing the $i^{th}$ SOW, and $\beta$ is the vector of coefficients describing the relationship between the covariates and the response, which here will be estimated using maximum likelihood estimation. A logistic regression model was fit to the ensemble of SOWs using the performance criteria defined in step 2. Logistic regression modeling was conducted using the [Statsmodel Python](https://www.statsmodels.org/stable/index.html) package. The data required for the full analysis is too large to include in this tutorial, but results can be found in the data file loaded below.The results files contain the occurence of different frequency and magnitude combinations under the experiment, in increments of 10, between 0 and 100. These combinations (100 for each user) are alternative decision-relevant metrics that can be used for scenario discovery.<jupyter_code># Set arrays for shortage frequencies and magnitudes
frequencies = np.arange(10, 110, 10)
magnitudes = np.arange(10, 110, 10)
realizations = 10


# Load performance and pseudo r scores for each of the users
results = [load_performance(all_IDs[i])/100 for i in range(len(all_IDs))]<jupyter_output><empty_output><jupyter_text>Step 3: Factor rankingTo rank the importance of each uncertain factor, we utilize McFadden's psuedo-R2, a measure that quantifies the improvement of the model when utilizing each given predictor as compared to prediction using the mean of the data set:$$R^2_{McFadden}=1-\frac{ln \hat{L}(M_{full})}{ln \hat{L}(M_{intercept})}$$Where $ln \hat{L}(M_{full})$ is the log likelihood of the full model (including the predictor) and $ln \hat{L}(M_{intercept})$ is the log likelihood of the intercept model (which predicts the mean probability of success across all SOWs). Higher values of McFadden's psuedo-R2 indicate higher factor importance (when the likelihood of the full model approaches one, the ratio of the likelihood of the full model compared to the intercept model will get very small).<jupyter_code>scores = [pd.read_csv('./data/' + all_IDs[i] + '_pseudo_r_scores.csv', sep=",") for i in range(len(all_IDs))]

freq = [1,0,0]
mag = [7,3,7]<jupyter_output><empty_output><jupyter_text>Step 4: Draw factor mapsThe McFadden's psuedo-R2 scores files contain preliminary logistic regression results on parameter importance for each of these combinations. Using these psuedo-R2 scores we will identify the two most important factors for each metric which we'll use to generate the final scenario discovery maps (note: there may be more than two important metrics for each user, but here we will demonstrate by mapping two).<jupyter_code>fig, axes = plt.subplots(3,1, figsize=(6,18), tight_layout=True)
fig.patch.set_facecolor('white')
for i in range(len(axes.flat)):
    ax = axes.flat[i]
    allSOWsperformance = results[i]
    all_pseudo_r_scores = scores[i]
    dta = pd.DataFrame(data = np.repeat(LHsamples, realizations, axis = 0), columns=param_names)
    dta['Success']=allSOWsperformance[freq[i],mag[i],:]
    pseudo_r_scores=all_pseudo_r_scores[str(frequencies[freq[i]])+'yrs_'+str(magnitudes[mag[i]])+'prc'].values
    top_predictors = np.argsort(pseudo_r_scores)[::-1][:2] #Sort scores and pick top 2 predictors
    
    # define color map for dots representing SOWs in which the policy
    # succeeds (light blue) and fails (dark red)
    dot_cmap = mpl.colors.ListedColormap(np.array([[227,26,28],[166,206,227]])/255.0)
    
    # define color map for probability contours
    contour_cmap = mpl.cm.get_cmap('RdBu')
    
    # define probability contours
    contour_levels = np.arange(0.0, 1.05,0.1)
    
    # define base values of the predictors
    SOW_values = np.array([1,1,1,1,0,0,1,1,1,1,1,0,0,0]) # default parameter values for base SOW
    base = SOW_values[top_predictors]
    ranges = param_bounds[top_predictors]
    
    # define grid of x (1st predictor), and y (2nd predictor) dimensions
    # to plot contour map over
    xgrid = np.arange(param_bounds[top_predictors[0]][0],
                      param_bounds[top_predictors[0]][1], np.around((ranges[0][1]-ranges[0][0])/500,decimals=4))
    ygrid = np.arange(param_bounds[top_predictors[1]][0],
                      param_bounds[top_predictors[1]][1], np.around((ranges[1][1]-ranges[1][0])/500,decimals=4))
    all_predictors = [ dta.columns.tolist()[i] for i in top_predictors]
    dta['Interaction'] = dta[all_predictors[0]]*dta[all_predictors[1]]
    
    # logistic regression here
    result = fitLogit(dta, [all_predictors[i] for i in [0,1]])
    
    
    contourset = plotContourMap(ax, result, dta, contour_cmap,
                                dot_cmap, contour_levels, xgrid,
                                ygrid, all_predictors[0], all_predictors[1], base)
    ax.set_title(usernames[i])
cbar_ax = fig.add_axes([0.98, 0.15, 0.05, 0.7])
cbar = fig.colorbar(contourset, cax=cbar_ax)
cbar_ax.set_ylabel('Probability of Success', fontsize=16)
cbar_ax.tick_params(axis='y', which='major', labelsize=12)<jupyter_output>Optimization terminated successfully.
         Current function value: 0.378619
         Iterations 8
Optimization terminated successfully.
         Current function value: 0.397285
         Iterations 8
Optimization terminated successfully.
         Current function value: 0.377323
         Iterations 8
{'fraction_non_alphanumeric': 0.08439181916038752, 'fraction_numerical': 0.02658772874058127, 'mean_word_length': 4.788785046728972, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '312964', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 2822, 'n_tokens_neox': 2557, 'n_words': 1017, 'path': 'notebooks/basin_users_logistic_regression.ipynb', 'repo_name': 'IMMM-SFA/nanites', 'size': '9290'}