<jupyter_start><jupyter_text>Data analysis in the Amazon auditMost of the materials in this notebook were kindly provided by <NAME> (www.sapiezynski.com)  The motivation behind the analysisRead it herehttps://themarkup.org/amazons-advantage/2021/10/14/amazon-puts-its-own-brands-first-above-better-rated-products<jupyter_code># installing the packages for today (if needed)
!pip install -U pandas
!pip install -U scikit-learn
!pip install -U statsmodels
import pandas as pd
from sklearn.model_selection import train_test_split
# this is the data from the markup audit, downloaded from their github: 
# https://github.com/the-markup/investigation-amazon-brands 

data = pd.read_csv('https://www.sapiezynski.com/cs4910/markup/data.csv')
data<jupyter_output><empty_output><jupyter_text>![features](https://mrkp-static-production.themarkup.org/graphics/amazon-methodology-product-comparison/1634143723511/assets/product-pairs-subtraction-desktop.png)<jupyter_code>columns = ['stars_delta', 'reviews_delta', 'is_amazon',
        'is_shipped_by_amazon', 'is_sold_by_amazon',
        'is_top_clicked', 'placed_higher']
dataset = data[columns]
train, test = train_test_split(dataset, test_size=0.2, stratify=dataset['placed_higher'])
features = ['is_amazon',
        'is_shipped_by_amazon', 'is_sold_by_amazon',
        'is_top_clicked', 'stars_delta', 'reviews_delta', ]
from sklearn.ensemble import RandomForestClassifier
clf = RandomForestClassifier(n_estimators=500, max_depth=3, max_features = 6)
clf.fit(train[features], train['placed_higher'])

# let's see how well it does in predicting the test set
# 1. predict the test data
y_pred = clf.predict(test[features])

# 2. calculate accuracy, i.e. the fraction of time that the prediction was correct out of all predictions
(y_pred == test['placed_higher']).sum()/len(y_pred)
for feature, importance in zip(features, clf.feature_importances_):
    print(feature, importance)<jupyter_output>is_amazon 0.9032928487280287
is_shipped_by_amazon 0.007487632851518247
is_sold_by_amazon 0.002829947430288105
is_top_clicked 0.002445399880340501
stars_delta 0.019225107462222012
reviews_delta 0.06471906364760226<jupyter_text>![features](https://www.sapiezynski.com/cs4910/markup/importances.png) So clearly, `is_amazon` is the most important, but what exactly does it mean, how does it translate the the probability of being the top result? RandomForests won't tell us... But logistic regression can.First, quick reminder on regression. **Linear** regression has the following formula:$$ y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... $$This means that:* when all variables are equal to 0, the outcome variable $y$ is equal to $\beta_0$, the intercept.* a **unit** change in $x_1$ corresponds to $y$ changing by $\beta_1$ It might not be a powerful model, but it offers a clear interpretation/explanation in regression problems.**Logistic** regression is mostly used for binary classification and it looks similar to the linear regression, but now the outcome variable is the logarithm of odds of success.$$ ln(\frac{\pi}{1-\pi}) = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... $$where $pi$ is the probability of success, so $1-\pi$ is the probability of failure.Then, a unit change in $x_1$ corresponds to $\beta_1$ change in log odd ratios... but that is not intuitive at all. Let's rewrite it to calculate the probability directly:$$\begin{align}ln(\frac{\pi}{1-\pi}) &= \beta_0 + \beta_1x_1 + \beta_2x_2 + ... \\\frac{\pi}{1-\pi} &= e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...} \\{\pi} &= ({1-\pi})e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...} \\{\pi} &= e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...}-\pi e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...}\\{\pi}(1+e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...}) &= e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...} \\{\pi} &= \frac{e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...}}{1+e^{\beta_0 + \beta_1x_1 + \beta_2x_2 + ...}}\end{align}$$So, while we cannot directly say from that a unit change of $x_1$ corresponds to a certain change in the probability of the outcome, we can just plug zeros and ones into the equation above and compute the changes. Let's load the data from scratch and do some cleaning (I'm dividing by two, because The Markup left the values of -2, 0, 2 so they never have a "unit" change. If we divide them by two, the unit change actually is meaningful).<jupyter_code>data = pd.read_csv('https://www.sapiezynski.com/cs4910/markup/data.csv')
data['is_shipped_by_amazon'] /= 2
data['is_sold_by_amazon'] /= 2
data['is_amazon'] /= 2
data['is_top_clicked'] /= 2
dataset = data[columns]

train, test = train_test_split(dataset, test_size=0.2, stratify=dataset['placed_higher'])<jupyter_output><empty_output><jupyter_text>Let's start with a model that only has the intercept (constant) and `is_amazon` as features.<jupyter_code>import statsmodels.api as sm
from statsmodels.tools.tools import add_constant

# adding the intercept as a variable
train = add_constant(train)

# selecting which features to train on
features = ['const', 'is_amazon']

# training a logistic regression model. First the dependent (outcome variable), then the features
log_reg = sm.Logit(train['placed_higher'], train[features]).fit()
log_reg.summary()<jupyter_output><empty_output><jupyter_text>We will mostly be looking at the bottom table.`const` (the intercept) translates to the probability if all other variables are equal to 0$$ p = \frac{e^{\beta_{const}}}{1+e^{\beta_{const}}} $$<jupyter_code>import numpy as np

# the beta coefficients are stored in the .params list, with the constant being first
# get the constant out of the fit model:
const = log_reg.params[0]

np.e**const/(1+np.e**const)<jupyter_output><empty_output><jupyter_text>This should be very close to 50%. In our case if all variables are zero, it means that either:1. The products are both from amazon2. The products are both not from amazonWe don't have a way to differentiate them, so it's a coin toss (50%) The other coefficients describe the change in odds that corresponds to a unit change in the variable.* if the coefficient is positive, it means that an **increase** of variable corresponds to an **increase** in the odds of the positive outcome* if the coefficient is positive, it means that an **increase** of variable corresponds to an **decrease** in the odds of the positive outcomeSo what happens when `is_amazon` is equal to 1, i.e. one product is an amazon product and the other isn't?By looking at the `coef` we already know it's going to go up (because the coefficient is positive). But by How much? Let's calculate:$$ p = \frac{e^{\beta_{const} + \beta_{is\_amazon}}}{1+e^{\beta_{const} + \beta_{is\_amazon}}} $$<jupyter_code>numerator = np.e**(log_reg.params[0] + log_reg.params[1])
print(numerator/(1+numerator))<jupyter_output>0.9293369710808036<jupyter_text>That means that in this dataset 93\% of the time we have two products where one is from amazon and the other isn't, it's the amazon product that will be on the first place.Solved? Well no, because maybe this is just because amazon products are better. Let's include the star difference in the analysis to **control** for this:<jupyter_code>features = ['const', 'is_amazon','stars_delta', ]
log_reg = sm.Logit(train['placed_higher'], train[features]).fit()
log_reg.summary()<jupyter_output>Optimization terminated successfully.
         Current function value: 0.494329
         Iterations 6<jupyter_text>Now notice the last two columns - that's the 95% confidence interval. If one value is positive, and the other negative, the interval includes 0, which means we can't really tell if the effect is positive, or negative - it's not significant. Turns out, controling for the star rating doesn't change anything.Let's add the review delta:<jupyter_code>features = ['const', 'is_amazon','stars_delta', 'reviews_delta']
log_reg = sm.Logit(train['placed_higher'], train[features]).fit()
log_reg.summary()<jupyter_output>Optimization terminated successfully.
         Current function value: 0.493737
         Iterations 6<jupyter_text>Same story, so let's include the rest of our variables.<jupyter_code>features = ['const', 'is_amazon','stars_delta', 'reviews_delta',
        'is_shipped_by_amazon', 'is_sold_by_amazon',
        'is_top_clicked']
log_reg = sm.Logit(train['placed_higher'], train[features]).fit()
log_reg.summary()<jupyter_output>Optimization terminated successfully.
         Current function value: 0.485827
         Iterations 6<jupyter_text>The only other significant variable is `is_sold_by_amazon`. How do we interpret its value?For exaple, both products are from amazon or are not from amazon, but one is sold by amazon, so `is_amazon = 0` and `is_sold_by_amazon = 1`<jupyter_code>numerator = np.e**(log_reg.params[0] + log_reg.params[1] * 0 + log_reg.params[5] * 1)
print(numerator/(1+numerator))<jupyter_output>0.6933580136415483<jupyter_text>That's the probability the product sold by amazon will be first.For another example, first product is from amazon and sold by amazon, the other is neither:<jupyter_code>numerator = np.e**(log_reg.params[0] + log_reg.params[1] * 1 + log_reg.params[5] * 1)
print(numerator/(1+numerator))<jupyter_output>0.9549840944621973<jupyter_text>That's the probability the product sold by amazon will be first - it's slightly higher than when the first product is from amazon and the other isn't but we don't know who sells it. What about the predictive performance? The Markup uses RandomForests with hyperparameter tuning to get the best prediction and ends up with 73.2% accuracy<jupyter_code>test = add_constant(test)
y_pred = log_reg.predict(test[features])
((y_pred > 0.5) == test['placed_higher']).mean()
train.shape
train
test<jupyter_output><empty_output><jupyter_text>That's not any worse than random forest - we can an equally good fit but we can actually interpret just how important the different features are! The problem that we looked at so far is defined as follows: given the top two items, which one is going to be be placed in the first position?Now, let's solve a slightly different problem formulation (not covered in the Markup write up): given characteristics of the product, how likely is it the top product on its search page?<jupyter_code># let's get the data first:
dataset = pd.read_csv('https://www.sapiezynski.com/cs4910/markup/dataset_full.csv')
dataset<jupyter_output><empty_output><jupyter_text>This dataset is an aggregation of all first search result pages from the Markup audit that had any amazon products on them. The sponsored results are removed.In this dataset each row is a product with its star rating, the number of reviews, the binary indicators of whether it's an amazon product, it's sold or shipped by amazon, whether it was among the most clicked products, and the outcome variable - whether it was the top product in its search results page.We're looking at at most 24 first results, so on average only one in 24 rows is the top placed product:<jupyter_code>dataset['is_top'].mean()
# the lowest star rating is 1, not 0, let's adjust the rating such that it starts at 0:
dataset['stars'] -= 1<jupyter_output><empty_output><jupyter_text>Excercise 1: Interpreting the meaning of logistic regression coefficientsIn this exercise you will interpret the values of the coefficients of a logistic regression model that tries to predict whether a product was the top product on its search page.In particular, please answer these questions:1. What is the probability that a hypothetical product with the lowest rating, no reviews, not sold/shipped/produced by amazon, and not among the top clicked products will be placed at the top?1. Increase in which variables corresponds to a higher probability of being placed at the top?1. How do you interpret the `stars` coefficient and its statistical significance?1. If a product receives one more review without changing any other of its characteristics, how does it affect its chances to placed at the top?     * How about 1000000 more reviews?1. What is the probability that a five-star product with 1000 reviews, produced, sold, and shipped by amazon that was among the top clicked results will be placed on the top of the result list?     * Is this number higher, or lower than you expected? Does it change how you interpret the results of the model for only the top two positions? Why?<jupyter_code>dataset = add_constant(dataset)
cols =['const', 'stars', 'reviews', 'is_amazon', 'is_sold_by_amazon', 'is_shipped_by_amazon', 'top_clicked']
log_reg = sm.Logit(dataset['is_top'],
                   dataset[cols]).fit()
log_reg.summary()
# What is the probability that a hypothetical product with the lowest rating, no reviews, 
# not sold/shipped/produced by amazon, and not among the top clicked products will be placed at the top?
np.e**(log_reg.params[0])/(1+ np.e**(log_reg.params[0]))<jupyter_output><empty_output><jupyter_text>- Increase in which variables corresponds to a higher probability of being placed at the top?>reviews, is_amazon, is_shipped_by_amazon, top_clicked (all ones that have positive coefs and are significant) - How do you interpret the stars coefficient and its statistical significance?> It is not statistically significant, so it appears not to have an effect on the placement.> its value is negative, so if it was significant it would indicate that higher rated items are less likely to be on top - If a product with no reviews receives its first one without changing any other of its characteristics, how does it affect its chances to placed at the top? How about when it receives 1000000 reviews instead?<jupyter_code>p0 = np.e**(log_reg.params[0])/(1+ np.e**(log_reg.params[0]))
p1 = np.e**(log_reg.params[0] + log_reg.params[2])/(1+ np.e**(log_reg.params[0]+log_reg.params[2]))
print(p1-p0)
p2 = np.e**(log_reg.params[0] + log_reg.params[2] * 1000000)/(1+ np.e**(log_reg.params[0]+log_reg.params[2]*1000000))
print(p2-p0)
p0
p2<jupyter_output><empty_output><jupyter_text>> The probability of being placed at the top grows by 7 percentage points from 3% to 10% - What is the probability that a five-star product with 100000 reviews, produced, sold, and shipped by amazon that was among the top clicked results will be placed on the top of the result list?<jupyter_code>numerator = np.e**(log_reg.params[0]\
                  + log_reg.params[1]*4\
                  + log_reg.params[2]*100000\
                  + log_reg.params[3]\
                  + log_reg.params[4]\
                  + log_reg.params[5]\
                  + log_reg.params[6])
numerator/(1+ numerator)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1015009252278802, 'fraction_numerical': 0.0304982523473374, 'mean_word_length': 5.014839241549876, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 68, 'https://': 7, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '185681', 'license': 'MIT', 'n_tokens_mistral': 4740, 'n_tokens_neox': 4341, 'n_words': 1656, 'path': 'Week 10/solution.ipynb', 'repo_name': 'carlomarxdk/algorithmic-fairness', 'size': '14602'}