This article is the first one in a series of publications dedicated to explaining various aspects of Logistic Regression as a substitute for the “classic” A/B test. Unlike t-tests, logistic regression allows performing the same hypothesis testing plus gives you tools for estimating interaction effects and primacy/novelty effects. Unlike Linear Regression, which deals with continuous outcomes (revenue per user, loading time etc.), Logistic Regression is applicable to binary outcomes (conversion rates, retention etc.).
Here’s a typical logistic regression output (as given by statsmodels.api.Logit):
Optimization terminated successfully.
Current function value: 0.543538
Iterations 5
Logit Regression Results
==============================================================================
Dep. Variable: Y No. Observations: 12
Model: Logit Df Residuals: 10
Method: MLE Df Model: 1
Date: Wed, 13 Mar 2024 Pseudo R-squ.: 0.03343
Time: 11:12:10 Log-Likelihood: -6.5225
converged: True LL-Null: -6.7480
Covariance Type: nonrobust LLR p-value: 0.5018
================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
const -1.6094 1.095 -1.469 0.142 -3.756 0.538
IS_TREATMENT 0.9163 1.396 0.656 0.512 -1.821 3.653
================================================================================There are plenty of things to digest, but let’s concentrate on these three for the time being: Log-Likelihood, LL-Null, LLR p-value. After all, these three guys are what we are interested in this article. These three will be enough for us to either reject or fail to reject the null hypothesis, which is — there’s no difference between control and treatment.
But, before getting into those scary numbers, let’s take a step back and look at this code.
import pandas as pd
import numpy as np
# Set a seed for reproducibility
np.random.seed(36) # You can choose any number as a seed
sample_size = 6
prob_control = 0.10
prob_treatment = 0.20
data_control = np.random.binomial(1, prob_control, sample_size)
data_treatment = np.random.binomial(1, prob_treatment, sample_size)
df_control = pd.DataFrame({
"IS_TREATMENT": [0] * sample_size,
"Y": data_control
})
df_treatment = pd.DataFrame({
"IS_TREATMENT": [1] * sample_size,
"Y": data_treatment
})
df = pd.concat([
df_control, df_treatment
], ignore_index=True)
mean_p = df['Y'].mean()
def assign_sample_probability(row, mean_prob):
return mean_prob if row['Y'] == 1 else 1 - mean_prob
df['p_sample'] = df.apply(assign_sample_probability, mean_prob=mean_p, axis=1)
df['p_sample_log'] = np.log(df['p_sample'])
df.head()This code:
- Generates a small sample by drawing 12 observations from Bernoulli distribution, 6 for control and 6 for treatment with the probability of success being 0.10 and 0.20
- Calculates the average probability of successes across the entire dataset
- Calculates the probability of an individual values occurring given the average probability of success across the dataset as seen in p_sample
- Takes the log of each individual probability (as seen in p_sample_log)
Code output
| IS_TREATMENT | Y | p_sample | p_sample_log | |
|---|---|---|---|---|
| 0 | 0 | 0 | 0.75 | −0.287682 |
| 1 | 0 | 0 | 0.75 | −0.287682 |
| 2 | 0 | 1 | 0.25 | −1.386294 |
| 3 | 0 | 0 | 0.75 | −0.287682 |
| 4 | 0 | 0 | 0.75 | −0.287682 |
| 5 | 0 | 0 | 0.75 | −0.287682 |
| 6 | 1 | 0 | 0.75 | −0.287682 |
| 7 | 1 | 0 | 0.75 | −0.287682 |
| 8 | 1 | 0 | 0.75 | −0.287682 |
| 9 | 1 | 0 | 0.75 | −0.287682 |
| 10 | 1 | 1 | 0.25 | −1.386294 |
| 11 | 1 | 1 | 0.25 | −1.386294 |
What this dataset mimics is results of a tiny experiment with only 12 observations (this tiny sample is chosen for teaching purposes so that I can better demonstrate the point).
Imagine we performed an experiment (improved something on our website). But instead of a simple A/B test, we chose to apply Logistic Regression to the data (for the reasons described in the introduction for this publication), but we are still interested only in one thing — is there actually any difference between visitors, who received the treatment and those who didn’t? For our case, simple A/B testing is more than suitable, but this simplification of application of logistic regression is necessary to better get the point across and get us started.
Log-Likelihood Null (LL-Null) and its “father” — Likelihood Null
So, where does LL-Null come from, or Log-Likelihood-Null for that matter? If you thought “null hypothesis” by looking at the name, you’d be absolutely right! This has everything to do with null hypothesis. And you have probably heard the word “likelihood” in statistics before. Let’s make things easier and drop “Log-” and only think about Likelihood-Null for a moment. For that, we need to exponentiate the value of −6.7480. Exponentiation simply converts log-value back to the conventional, intuitive value, which is 0.0011732237233539518 in our case.
So, 0.0011732237233539518 is our null-likelihood.
What does this value represent and what does it tell us? Likelihood of what are we dealing with exactly? It tells us that, given the assumed probability of success, we have roughly one in a thousand chance to get this combination of outcomes (these are the values of Y in our dataset):
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1]
What’s the assumed probability of success in our case? Since we are dealing with null-hypothesis:
- we assume that there’s no difference between treatment and control and the probability of success between the two groups is the same, any discrepancy is due to sampling error
- our best shot at estimating the population probability is taking average across all groups (just as we do in the classic A/B testing)
- this average becomes our assumed probability of success under null hypothesis
Like many things in statistics, the only way we can guess the real population probability of success is by taking the average of the entire dataset, which is 12 observations. Sometimes our guess will be far off from the real population probability, but that’s the harsh reality. The average of the dataset in our case is 0.25. In our situation we of course have the population parameters, but that’s not something we usually have in real life, so our best estimate is derived from the available data.
So, getting back to our null-likelihood. What it tells us is that:
If control and treatment actually have no difference and hence, belong to the same population with the probability of success is 0.25, there’s only 0.0011732237233539518 likelihood that we will get the given set of values from this population. Or roughly 12 times out of 10000, if we round things up.
How do we actually arrive at this value? Remember, the probability of drawing 1 is 0.25 (that’s what we established from the dataset), so, the probability of 0 is 0.75.
If you are a bit familiar with probabilities, you may know, that the probability of getting all the values together is their individual probabilities multiplied. So:
0.75 × 0.75 × 0.75 × 0.75 × 0.25 × 0.75 × 0.25 × 0.75 × 0.75 × 0.75 × 0.25 × 0.25 = 0.0011731982231140137
which is pretty close to 0.0011732237233539518, the discrepancy is due to rounding error.
Let’s not take my word for it and do a little simulation:
import numpy as np
# Parameters for the Bernoulli distribution
p_success = 0.25 # Probability of success
n_draws = 12 # Number of values to draw
identical_draws = 0
iterations = 10000000
for i in range(iterations) :
# Draw 12 values from a Bernoulli distribution
bernoulli_draws = np.random.binomial(n=1, p=p_success, size=n_draws)
if(bernoulli_draws.tolist() == [0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0]) :
identical_draws += 1
print(identical_draws/iterations)What this code does is:
- makes 12 draws from Bernoulli distribution with the probability of success at 0.25;
- compares the outcome to the outcome we have in our dataset;
- if the outcome is identical to the one we have in our dataset, our counter identical_draws is increased by one;
- we perform this procedure 10 000 000 times for precision.
- we calculate the proportion of iterations which resulted in the desired outcome.
Hopefully, when you run this code, you’ll get something close to 0.0011731982231140137 — I ended up with 0.00117304.
So, to reiterate, if we assume that:
- every 0 across two groups has the probability to be drawn at 0.75
- every 1 across two groups has the probability to be drawn at 0.25
- the likelihood of drawing this particular set of values is 0.0011731982231140137
This is pretty small, isn’t it?
Log-Likelihood, or Likelihood or Likelihood-Alternative
Our regression output also has something called “Log-Likelihood”, also known as Log-Likelihood Model, but since we have Log-Likelihood Null, I prefer the name Log-Likelihood Alternative, since it better aligns with my background in A/B testing and better maps to concepts of hypothesis testing.
When we perform logistic regression with categorical predictors, we ask ourselves:
What if the probability of success is not the same among the groups (“IS_TREATMENT = 0”, “IS_TREATMENT = 1” in our case)?
And if it’s not the same, how do we find out which probability of success the control has, and which probability of success the treatment has? We assume that the probability of drawing 1, and, by extension, 0 is different for control and treatment.
Is there a specific combination of probabilities for control and treatment, that would make our given values more likely to occur? Moreover, what we want to find out is not only a combination of two probabilities of success that would make our combination more likely to occur, but such a combination, for which the likelihood of our outcome is at max. In other words, there is no such set of probability of successes for control and probability of success for treatment that would make the given combination of values more likely. That’s what maximizing the likelihood actually means.
The goal of logistic regression for hypothesis testing is to find such a combination of probabilities, for which the likelihood of given values occurring is maximal.
If we were simpler creatures, and still had a computer, we would try to go through all possible combinations of probabilities and see which combination of the probabilities has the highest likelihood of giving us the desired outcome. We might even go about the problem with a code like this:
import pandas as pd
import numpy as np
from functools import reduce
# Generate a range of probability values
values = np.arange(0.01, 0.99, 0.01)
# Convert numpy array to DataFrame
df1 = pd.DataFrame(values, columns=['Control'])
df2 = pd.DataFrame(values, columns=['Treatment'])
# Perform a cross join by creating a key to join on
df1['key'] = 1
df2['key'] = 1
# Merge the two DataFrames on the key column to perform the cross join
probability_combinations = pd.merge(df1, df2, on='key').drop('key', axis=1)
# cross_joined_df now contains the cross-joined values
control_values = [0, 1, 0, 0, 0, 0]
treatment_values = [0, 1, 0, 0, 0, 1]
likelihood = []
results = []
# Estimating Likelihood for all combinations of probabilities
for index, row in probability_combinations.iterrows():
p_success_control = row['Control']
p_success_treatment = row['Treatment']
probs = []
for value in control_values:
probs.append(p_success_control if value == 1 else 1 - p_success_control)
for value in treatment_values:
probs.append(p_success_treatment if value == 1 else 1 - p_success_treatment)
likelihood_value = reduce(lambda x, y: x * y, probs)
likelihood.append(likelihood_value)
# Append the combination along with its likelihood to the results list
results.append({
'Control Probability': p_success_control,
'Treatment Probability': p_success_treatment,
'Likelihood': likelihood_value
})
# Convert the results list to a DataFrame
results_df = pd.DataFrame(results)
# Find the index of the row with the maximum likelihood value
max_likelihood_index = results_df['Likelihood'].idxmax()
# Retrieve the row with the maximum likelihood value
max_likelihood_row = results_df.loc[max_likelihood_index]
# Print the row with the maximum likelihood value
print(max_likelihood_row)What this code does is:
- Creates all possible combinations of population probabilities for Control and Treatment (rounded to 2 decimal places), 9604 combinations in total
- Estimates the likelihood of observing the given combination of values under the specified probabilities
- Finds the row with the combination that resulted in the highest likelihood
In our case, the code is going to produce:
Control Probability 0.170000
Treatment Probability 0.330000
Likelihood 0.001469Which reads as “The highest likelihood of observing the given arrangement of values: [0, 1, 0, 0, 0, 0] for control and [0, 1, 0, 0, 0, 1] for treatment is 0.001469 (or roughly 15 times out of 10000). This likelihood corresponds to the probability of success for control being 0.17 and probability of success for treatment being 0.33.” Or:
0.83 × 0.17 × 0.83 × 0.83 × 0.83 × 0.83 × 0.67 × 0.33 × 0.67 × 0.67 × 0.67 × 0.33 = 0.001469489454588773
Now, let’s compare this output to our regression results presented earlier.
Logit Regression Results
==============================================================================
Dep. Variable: Y No. Observations: 12
Model: Logit Df Residuals: 10
Method: MLE Df Model: 1
Date: Wed, 13 Mar 2024 Pseudo R-squ.: 0.03343
Time: 11:12:10 Log-Likelihood: -6.5225
converged: True LL-Null: -6.7480
Covariance Type: nonrobust LLR p-value: 0.5018
==============================================================================Oh well, 0.001469 is nowhere to be found. But remember, we are dealing with logs in our logistic regression. So, if we take the natural logarithm of 0.001469489454588773, what we’ll end up with is… −6.522840248321066! (This is pretty darn close to the Log-Likelihood in the regression results).
import numpy as np
lla = np.log(0.83 * 0.17 * 0.83 * 0.83 * 0.83 * 0.83 * 0.67 * 0.33 * 0.67 * 0.67 * 0.67 * 0.33)
print(lla)So, our little simple code and the regression basically produce the same results, but expressed differently.
Why logs?
“But wait”, you may say, “why do we have to deal with log-likelihoods, log-odds and all that other fancy stuff instead of finding the probabilities directly like we did above?”
The reasons are multifold, but here are some of them.
As our sample size increases, (instead of 12 we may even have to work with something like a thousand or hundreds of thousands of observations), the likelihoods of observing this particular combination of values under any combination of probabilities gets extremely small. In fact, it becomes so small that even modern computers can’t handle them and simply give up, converting the values to their logs makes them more “compact” and more stable.
Here’s an updated likelihood-estimation code, where we increase the sample sizes to 1000 observations per group and try to run the code:
import pandas as pd
import numpy as np
from functools import reduce
# Generate a range of probability values
values = np.arange(0.01, 0.99, 0.01)
# Convert numpy array to DataFrame
df1 = pd.DataFrame(values, columns=['Control'])
df2 = pd.DataFrame(values, columns=['Treatment'])
# Perform a cross join by creating a key to join on
df1['key'] = 1
df2['key'] = 1
# Merge the two DataFrames on the key column to perform the cross join
probability_combinations = pd.merge(df1, df2, on='key').drop('key', axis=1)
# cross_joined_df now contains the cross-joined values
np.random.seed(36)
control_values = np.random.binomial(n=1, p=0.17, size=1000)
treatment_values = np.random.binomial(n=1, p=0.33, size=1000)
likelihood = []
results = []
# Estimating Likelihood for all combinations of probabilities
for index, row in probability_combinations.iterrows():
p_success_control = row['Control']
p_success_treatment = row['Treatment']
probs = []
for value in control_values:
probs.append(p_success_control if value == 1 else 1 - p_success_control)
for value in treatment_values:
probs.append(p_success_treatment if value == 1 else 1 - p_success_treatment)
likelihood_value = reduce(lambda x, y: x * y, probs)
likelihood.append(likelihood_value)
# Append the combination along with its likelihood to the results list
results.append({
'Control Probability': p_success_control,
'Treatment Probability': p_success_treatment,
'Likelihood': likelihood_value
})
# Convert the results list to a DataFrame
results_df = pd.DataFrame(results)
# Find the index of the row with the maximum likelihood value
max_likelihood_index = results_df['Likelihood'].idxmax()
# Retrieve the row with the maximum likelihood value
max_likelihood_row = results_df.loc[max_likelihood_index]
# Print the row with the maximum likelihood value
print(max_likelihood_row)We will get something like this:
Control Probability 0.01
Treatment Probability 0.01
Likelihood 0.00Likelihoods are never zero (same as probabilities) and this output is just a sign that our computer just couldn’t handle likelihoods and probabilities that small.
Here’s the output dataframe:
Results of the brute-force likelihood estimation
| Control probability | Treatment probability | Likelihood | |
|---|---|---|---|
| 0 | 0.01 | 0.01 | 0.0 |
| 1 | 0.01 | 0.02 | 0.0 |
| 2 | 0.01 | 0.03 | 0.0 |
| 3 | 0.01 | 0.04 | 0.0 |
| 4 | 0.01 | 0.05 | 0.0 |
| … | … | … | … |
| 9599 | 0.98 | 0.94 | 0.0 |
| 9600 | 0.98 | 0.95 | 0.0 |
| 9601 | 0.98 | 0.96 | 0.0 |
| 9602 | 0.98 | 0.97 | 0.0 |
| 9603 | 0.98 | 0.98 | 0.0 |
So remember, using our simple algorithm isn’t feasible with large sample sizes and with more than two groups — that’s why I chose a small sample size to demonstrate the idea behind the maximum likelihood estimation.
To make their work more efficient and less resource-heavy, humans invented sophisticated algorithms for finding this exact combination of probability values without going through all of them iteratively and these algorithms require converting our probabilities to something the algorithm can work with, but they basically do the same thing, it’s just that for us, humans, it doesn’t look as transparent and intuitive as working with the probabilities directly.
Log-Likelihood Ratio
Now, to the Log-Likelihood Ratio. As with any statistical test we can’t simply take values and reject the null hypothesis just because our likelihoods aren’t identical. That’s because we know already what sampling error and variability are and that, even if there’s no real difference between two groups, we sometimes may observe a considerable difference in the samples. That’s why we also need a method here to estimate how much our difference deviates from the one we might see under the null hypothesis, and if the deviation is quite extreme we go ahead and reject the null.
Log-likelihood ratio is a way to estimate how much our alternative hypothesis deviates from the null hypothesis. We estimate this deviation with p-value, just like we do with the classic A/B test.
First, what in the world is Log-Likelihood Ratio?
Let’s start with the formula in its traditional form. You may wonder what’s going on here and why we call it ratio if we subtract one log-likelihood from the other. This question may arise if you are not familiar with how logs work.
It may be confusing at first (it definitely is for me), but consider this expression:
Those are equal expressions, due to the properties of logarithm, and the second one is indeed a ratio of squared likelihoods, which kind of justifies its name. But remember, the one that contains the ratio is computationally unstable and that’s why we have to completely eliminate pure likelihoods from our calculations and work with their logs instead.
Incidentally, this expression will also produce the same result:
Turns out, squaring the likelihoods and then finding the difference of their logs gives us a statistic which, under null hypothesis, follows a chi2 distribution (like with a usual Student’s T-test, we have a distribution which we use to compare our given statistic to).
likelihood_null = 0.0011732237233539518
likelihood_alt = 0.0014699895309787728
print("Intuitive LLR", np.log(likelihood_alt**2 / likelihood_null**2))
print("Traditional LLR",-2 * (np.log(likelihood_null)-np.log(likelihood_alt)))
print("Another form", np.log(likelihood_alt**2)-np.log(likelihood_null **2))Null-Distribution
Let’s approximate our null-distribution by drawing values for treatment and control from the same distribution this time and do a logistic regression on the resulting data. The regression will assume there’s a difference between treatment and control. We will repeat this procedure 1000 times, each time storing the value of our log-likelihood ratio to a list. For regression analysis we will use statsmodels.api.Logit
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy.stats import chi2
from scipy.stats import kstest
llr_all = []
for i in range (1000) :
sample_size = 1000
prob_old_control, prob_new_control = 0.10, 0.10
prob_old_treatment, prob_new_treatment = 0.10, 0.10
data_old_users_control = np.random.binomial(1, prob_old_control, sample_size)
data_new_users_control = np.random.binomial(1, prob_new_control, sample_size)
data_old_users_treatment = np.random.binomial(1, prob_old_treatment, sample_size)
data_new_users_treatment = np.random.binomial(1, prob_new_treatment, sample_size)
df_old_control = pd.DataFrame({
"IS_NEW": [0] * sample_size,
"IS_TREATMENT": [0] * sample_size,
"Y": data_old_users_control
})
df_new_control = pd.DataFrame({
"IS_NEW": [1] * sample_size,
"IS_TREATMENT": [0] * sample_size,
"Y": data_new_users_control
})
df_old_treatment = pd.DataFrame({
"IS_NEW": [0] * sample_size,
"IS_TREATMENT": [1] * sample_size,
"Y": data_old_users_treatment
})
df_new_treatment = pd.DataFrame({
"IS_NEW": [1] * sample_size,
"IS_TREATMENT": [1] * sample_size,
"Y": data_new_users_treatment
})
df = pd.concat([
df_old_control, df_new_control,
df_old_treatment, df_new_treatment
], ignore_index=True)
mean_p = df['Y'].mean()
def assign_sample_probability(row, mean_prob):
return mean_prob if row['Y'] == 1 else 1 - mean_prob
df['p_sample'] = df.apply(assign_sample_probability, mean_prob=mean_p, axis=1)
df['p_sample_log'] = np.log(df['p_sample'])
# Create a design matrix
X = sm.add_constant(df[['IS_TREATMENT']])
# Fit the logistic regression model
model = sm.Logit(df['Y'], X)
result = model.fit(maxiter=100, disp=0)
# Get confidence intervals for the coefficients
confidence_intervals = result.conf_int()
lln = result.llnull
llf = result.llf
llr = 2 * (result.llf - result.llnull)
llr_all.append(llr)
plt.hist(llr_all, bins=30, density=True, alpha=0.5, label='LLR Histogram')
x_range = np.linspace(0, 10, 1000) # More typical range for chi2 with df=1
plt.plot(x_range, chi2.pdf(x_range, df=1), 'r-', lw=2, label='Chi2 PDF (df=1)')
plt.xlabel('Log-Likelihood Ratio')
plt.ylabel('Density')
plt.title('LLR Histogram vs. Chi2 PDF')
plt.legend()
plt.show()
# Perform the one-sample KS test
ks_statistic, p_value = kstest(llr_all, 'chi2', args=(1,))
print(f"KS Statistic: {ks_statistic}")
print(f"P-value: {p_value}")This time around we draw a sample from the same distribution, which means, the treatment and control have no difference in the probability of success:
prob_old_control, prob_new_control = 0.10, 0.10
prob_old_treatment, prob_new_treatment = 0.10, 0.10- We draw 1000 observations from each of four groups and run a logistic regression on them with the predictor “is_treatment”.
- We repeat this procedure 1000 times, each time storing the log-likelihood ratios to a list for further analysis.
- We plot all LLRs as a histogram against pdf-curve of chi2 distribution with 1 degree of freedom
- We also perform a one-sample Kolmogorov-Smirnov test, comparing our values to a theoretical chi2 distribution with 1 degree of freedom.
When you run the code, the output is going to look something like this:
LLR histogram vs. χ² density (df = 1)
- KS statistic
- 0.0325
- KS p-value
- 0.2360
This simple trick hopefully demonstrates that we are indeed dealing with chi2 distribution with 1 degree of freedom. Degrees of freedom change based on the number of parameters we estimate in our logistic regression. In our case it’s one (is the user in treatment or not).
And this is it as far as LL-null, LL-alternative and LLR are concerned, but there are still some things to cover. Namely, how odds and their logs fit into the picture and more importantly, why we need them in the first place.
First published on Medium on Mar 29, 2024.