# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Extinction risks for cartilaginous fish
#
# An exploration of some of the results in [Extinction risk is most acute for the world’s largest and smallest vertebrates](https://www.pnas.org/content/114/40/10678), Ripple et al., PNAS October 3, 2017 114 (40) 10678-10683
#
# Specifically, we'll investigate how extinction risks vary by weight for cartilaginous fish. This provides some nice practice with simple linear and logistic regression, with the overall goal of explaining basic diagnostics for both methods. 
#
# All of this (and more!) is in Chapters 2-5 of my Manning book, [Regression: A friendly guide](https://www.manning.com/books/regression-a-friendly-guide). 
#
# This notebook and the relevant CSVs are available in my [regression repo on github](https://github.com/mbrudd/regression), along with other code and data for the book. Clone and fork at will!

# ### Imports and settings

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import chi2
from sklearn import metrics

import warnings
warnings.filterwarnings('ignore')

sns.set_theme()
plt.rcParams['figure.figsize'] = [8,8]

# ### Data
#
# First things first: we need data! Thanks to the good people at [ROpenSci](https://ropensci.org/), the data from [Fishbase.org](https://www.fishbase.se/search.php) is easily available in the [rfishbase](https://docs.ropensci.org/rfishbase/index.html) package.
#

fish = pd.read_csv("fish.csv")
fish.shape

fish.columns

fish = fish.filter(["Species","Length","Weight"])
fish

train = fish.dropna(axis='index')
train

sns.scatterplot(data=train, x="Length", y="Weight")
plt.title("Fish weights versus fish lengths")

train = train.assign(LogLength = np.log(train.Length), LogWeight = np.log(train.Weight))

sns.scatterplot(data=train, x="LogLength", y="LogWeight")
plt.title("Log(Weight) versus Log(Length)")
plt.axvline( np.mean( train["LogLength"] ), linestyle='--')
plt.axhline( np.mean( train["LogWeight"] ), linestyle='--')

# ### Linear regression
#
# The basic metric for the strength of a _linear_ relationship is the _correlation coefficient_:
#

train.LogLength.corr( train.LogWeight )

# This is a very strong correlation! In real life, especially in the social sciences, correlations between .3 and .7 in magnitude are much more common. Having checked the linear relationship, let's fit the regression line:
#

train_model = smf.ols( "LogWeight ~ LogLength", data=train)
train_fit = train_model.fit()
train_fit.params

#
# This model says that 
#
# $$ \log{ \left( \text{Weight} \right) } ~ = ~ -3.322617 + 2.681095 \times \log{ \left( \text{Length} \right) } \ ,$$
#
# which is easier to digest after exponentiating:
#
# $$ \text{Weight} ~ = ~ e^{-3.322617} \times \text{Length}^{2.681095} ~ = ~ 0.036 \times \text{Length}^{2.681095} \ .$$
#
# This _power law relationship_ says that weight is roughly proportional to the cube of the length!
#
# ***
#
# The _null model_ predicts that _every_ needed/unseen weight equals the average of the known weights:
#

np.mean( fish["Weight"] )

np.log( np.mean( fish["Weight"] ) )

# Is the regression model better than this at predicting weights? Answering this specific question is the job of the _coefficient of determination_, denoted $R^2$. 
#
# $$R^2 ~ = ~ \frac{ \text{TSS} - \text{SSR} }{ \text{TSS} } ~ = ~ 1 - \frac{ \text{SSR} }{ \text{TSS} }$$
#
# You could compute it this way...

( train_fit.centered_tss - train_fit.ssr) / train_fit.centered_tss

#
# but don't! It's already provided:
#

train_fit.rsquared

# ### Sharks!
#
# The information we need about [cartilaginous fish](https://en.wikipedia.org/wiki/Chondrichthyes) (sharks, rays, skates, sawfish, ghost sharks) comes from the [IUCN Red List](https://www.iucnredlist.org/):

sharks = pd.read_csv("chondrichthyes.csv")
sharks

sharks = sharks.join( fish.set_index("Species"), on="Species")
sharks

sharks = sharks[ sharks.Length.notna() ]
sharks = sharks[ sharks.Category != "Data Deficient" ]
sharks

# ### Data imputation
#
# Use the power law relationship to _impute_ the missing weights:
#

imp = np.exp( train_fit.params.Intercept )*np.power( sharks.Length, train_fit.params.LogLength )
sharks.Weight = sharks.Weight.where( sharks.Weight.notna(), imp )
sharks

sharks = sharks.assign(LogLength = np.log(sharks.Length), 
                       LogWeight = np.log(sharks.Weight))

sns.scatterplot( data=sharks, x="LogLength", y="LogWeight")
plt.title("Log(Weight) versus Log(Length) for sharks")

sharks

threatened = ["Critically Endangered","Endangered","Vulnerable"]

sharks["Threatened"] = sharks["Category"].isin( threatened ).astype('int')
sharks = sharks.drop(columns = "Category")
sharks

null_prob = np.mean(sharks["Threatened"])
null_prob

sharks_model = smf.glm("Threatened ~ LogWeight", data=sharks, family=sm.families.Binomial())
sharks_fit = sharks_model.fit()
sharks_fit.params

# This model says that
#
# $$\log{ \left( \text{Odds of being threatened} \right) } ~ = ~ -3.173571 + 0.293120 \times \log{\left( \text{Weight} \right) } \ ,$$
#
# which is equivalent to a power law:
#
# $$\text{Odds of being threatened} ~ = ~ .042 \times \text{Weight}^{.293120} \ .$$
#
# In other words, bigger fish are more likely to be threatened.
#
# ***
#
# - For logistic models, the _deviance_ is analogous to the sum of squared residuals in linear regression analysis; logistic model coefficients minimize the deviance. 
#
# - The _likelihood ratio statistic_ compares the deviances of the simple logistic model and the null model; it's analogous to the coefficient of determination.  
#
# - Unlike the coefficient of determination, the likelihood ratio statistic defies easy interpretation. It's easy to gauge its size, though: it's a $\chi^2$ statistic with $df=1$ (_why_ this is true is another story...).
#

sharks_fit.null_deviance - sharks_fit.deviance

1 - chi2.cdf(sharks_fit.null_deviance - sharks_fit.deviance, df=1)

#
# This is astronomically small -- the logistic model that includes `LogLength` is better than the null model that ignores it!
#
# And if we plot the results, things look pretty good:
#

sns.regplot(data=sharks, x="LogWeight", y="Threatened", logistic=True, ci=None)
plt.savefig("sharks_fit.png")

# ### Model assessment 
#
# Ripple et al. stop here with this particular model, but they should have assessed it carefully! We'll look at two options for what to do next.
#

# #### Binary classification and the ROC curve

#
# The naive thing to do is to compare the model's fitted probabilities to a threshold of 50% : classify the fish as `Threatened` if the fitted probability is higher than 50%, as `Not threatened` otherwise.
#

sharks["Class"] = (sharks_fit.fittedvalues > 0.50).astype(int)
sharks

pd.crosstab(sharks["Threatened"], sharks["Class"])

np.mean( sharks["Threatened"] == sharks["Class"] )

fpr, tpr, thresholds = metrics.roc_curve(sharks["Threatened"], sharks_fit.fittedvalues)
chronic_auc = metrics.auc(fpr, tpr)
chronic_auc

plt.figure()
plt.plot(fpr, tpr, label='ROC curve    AUC: %0.2f' % chronic_auc)
plt.plot([0,1], [0,1], 'r--', label='Random classification')
# plt.xlim([0, 1]) 
# plt.ylim([0, 1.05])
plt.xlabel('False Positive Rate (1-Specificity)')
plt.ylabel('True Positive Rate (Sensitivity)')
plt.title('ROC curve for shark extinction risk classifier')
plt.legend(loc="lower right")

# #### Logistic analogues of $R^2$

# _McFadden's pseudo_-$R^2$ : replace sums of squares with deviances to measure the proportional reduction in the deviance

R2_M = 1 - (sharks_fit.deviance / sharks_fit.null_deviance)
R2_M

#
# Or use the native sums of squares in this context:
#

sharks["Null_residual"] = sharks["Threatened"] - null_prob
sharks["Residual"] = sharks["Threatened"] - sharks_fit.fittedvalues
sharks["Difference"] = sharks_fit.fittedvalues - null_prob
R2_S = np.sum(sharks["Difference"]**2) / np.sum(sharks["Null_residual"]**2)
R2_S

#
# Or compute _Tjur's coefficient of discrimination_: a good model should, on average, assign high probabilities to observed successes (1's) and low probabilities to observed failures (0's)
#

sharks["Fit_prob"] = sharks_fit.fittedvalues
sns.displot( data=sharks, x="Fit_prob", col="Threatened", binwidth=0.2)

fit_avgs = sharks.groupby("Threatened").agg(Fit_average=('Fit_prob','mean'))
fit_avgs

R2_D = fit_avgs["Fit_average"][1] - fit_avgs["Fit_average"][0]
R2_D

# Yikes! Not a very good model after all. :(



StarcoderdataJupyterScriptsDedupFiltered