Predicting a fair coin flip outcome with logistic regression

Question

I am trying to better understand the results of logistic regression models and I wanted to apply a logistic regression model on a trivial "fair" coin flip simulation example. I want to simulate a number of coin tosses and fit a logistic regression model with whether we had at least one heads as the outcome and the number of coin tosses as the only predictor.

In this scenario, I know that a logistic regression model is not very useful, but I am just trying to understand how to interpret things and apply the model to real-life examples better.

Here is a simple explanation of what the simulation does:

• There are 10000 examples (rows)
• Each example has a random number of coin tosses (between 0 and 30, uniform distribution)
• Each coin toss has a 0.5 chance of getting a heads
• For each example, result = 1 if at least one toss resulted in a heads, result = 0 otherwise
• We end up with a 10000 by 2 dataframe with two columns: nb_toss and result

My simulation uses the following code:

import numpy as np
import pandas as pd

# Set random seed to always get the same results
np.random.seed = 2

examples = [str(x) for x in range(10000)]
results = []
nb_tosses = []

# For each of my 10000 example rows...
for e in examples:

  # Assign a random (uniform distribution) number of tosses to each example, from 0 to 30
  nb_toss = int(np.random.uniform(0,30, size=1))
  nb_tosses.append(nb_toss)

  result = 0
  for toss in range(nb_toss):

    # Each coin toss has a 0.5 chance of success (binomial random variable)
    # If at least one coin results in a success, assign 1 to result, otherwise assign 0
    result = np.random.binomial(1, 0.5, 1)[0]

    if result == 1:
      break
  results.append(result)

# Build a dataframe with the number of tosses and the result (1 or 0) for each row
toss_df = pd.DataFrame({"nb_toss" : nb_tosses, "result" : results})

So as you can see the coin is completely fair (50% probability of heads).

Now I fit my logistic regression model:

from patsy import dmatrices
import statsmodels.api as sm

# Build my X (10000 rows, cols: intercept = 1, nb_toss)
# and y (result)
y, X = dmatrices("result ~ nb_toss", data=toss_df)

# Fit a logistic regression model that predicts toss result given number of coin tosses
logit_mod = sm.Logit(y, X)
logit_res = logit_mod.fit(disp=0)

# Print out the summary
logit_res.summary()

Logit Regression Results
Dep. Variable:  result  No. Observations:   10000
Model:  Logit   Df Residuals:   9998
Method: MLE Df Model:   1
Date:   Mon, 15 May 2017    Pseudo R-squ.:  0.6503
Time:   13:47:43    Log-Likelihood: -849.40
converged:  True    LL-Null:    -2429.0
LLR p-value:    0.000
coef    std err z   P>|z|   [95.0% Conf. Int.]
Intercept   -1.5767 0.102   -15.490 0.000   -1.776 -1.377
nb_toss 1.1099  0.046   24.204  0.000   1.020 1.200

I am particularly interested in the interpretation of the pseudo R-squared.

Also, I find a predicted probability of 17% for 0 coin tosses, which is obviously impossible. Is there a way to adjust the model to take this into account?

I am finally trying to understand the following: say I have a dataset where the result of the coin toss has an influence on another binary outcome y. If my example gets at least one heads, then the probability of y = 1 is 80%. Otherwise, the probability of y = 1 is 40%. If I am unaware of these true probabilities, how well can I fit a logistic regression model that estimates the effect of getting at least a heads or tails, or the effect of the number of coin tosses on this binary outcome?

Answer 1

Below is part of an answer for how to use a different link function to capture the nonlinearity.

As discussed in the comments, the relationship between the number of tosses, nb_toss, and the probability of observing a success is nonlinear.

If $p_0$ is the probability of observing a failure in a single toss, then the probability of observing only failures in n tosses is $p_0^n$, and of observing a success overall is $1-p_0^n$

For simplicity I switch the definition of the outcome so that failure is 1 and success is zero. Then the probability of observing a failure is $Pr(failure) = p_0^n$ which we can rewrite as.

$Pr(failure) = exp(n * log(p_0)$

This is just a log-link with n, i.e. nb_toss in the linear prediction part. The estimated coefficient is the log probability of a toss, so we need to take exp to recover $p_0$.

Below I use statsmodels with the data and imports from the question. The estimated $p_0$ is 0.508, close to 0.5, and the prediction also match closely the true probabilities.

res_glm = sm.GLM(1 - y, X[:, 1], 
                 family=sm.families.Binomial(link=sm.families.links.log())).fit()

print(res_glm.summary())
print(np.exp(res_glm.params))
nbt = X[:, 1]
ii = np.arange(20)

table = np.column_stack((0.5**ii, res_glm.predict(ii),
                         [1 - y[nbt == i].mean() for i in ii]))
print(pd.DataFrame(table, columns=['true', 'predicted', 'sample']))

this prints

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                10000
Model:                            GLM   Df Residuals:                     9999
Model Family:                Binomial   Df Model:                            0
Link Function:                    log   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -765.47
Date:                Mon, 15 May 2017   Deviance:                       1530.9
Time:                        19:57:36   Pearson chi2:                 1.88e+04
No. Iterations:                    10                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
var_0         -0.6762      0.020    -34.626      0.000      -0.714      -0.638
==============================================================================
[ 0.50856518]

        true  predicted    sample
0   1.000000   1.000000       NaN
1   0.500000   0.508565  0.521368
2   0.250000   0.258639  0.290141
3   0.125000   0.131535  0.120944
4   0.062500   0.066894  0.055385
5   0.031250   0.034020  0.027778
6   0.015625   0.017301  0.008380
7   0.007812   0.008799  0.014663
8   0.003906   0.004475  0.002770
9   0.001953   0.002276  0.000000
10  0.000977   0.001157  0.000000
11  0.000488   0.000589  0.000000
12  0.000244   0.000299  0.002924
13  0.000122   0.000152  0.000000
14  0.000061   0.000077  0.002915
15  0.000031   0.000039  0.000000
16  0.000015   0.000020  0.000000
17  0.000008   0.000010  0.000000
18  0.000004   0.000005  0.000000
19  0.000002   0.000003  0.000000

aside: statsmodels prints a DomainWarning
DomainWarning: The log link function does not respect the domain of the Binomial family.

In general, there can be problems when using log-link with Binomial, i.e. log-Binomial, because the log-link does not force the predicted values to be in the range [0, 1]. However, the way this is set up in this example, the prediction is limited by 1 because the explanatory variable is nonnegative.