2
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ Unless you explain in words what your code is doing, this question will be inaccessible to anyone without a knowledge of Python and the time and interest to review the code, run it, and study its output. That would exclude many of the best candidates to answer it. $\endgroup$ – whuber May 15 '17 at 14:25
  • 1
    $\begingroup$ One reason why an explanation in words is so important is to check that the code is doing what you think it is. The explanation you added is very helpful (+1). Why do you suppose, though, that the log odds of at least one head ought to be a linear function of the number of tosses? You might want to carefully consider the cases $n=0$ and $n=1$ especially. Ordinarily when one is learning a technique to deal with a particular model, it is best to begin by generating data that conform to that model--otherwise the results are difficult or impossible to interpret. $\endgroup$ – whuber May 15 '17 at 16:34
  • 1
    $\begingroup$ I ran your example and the prediction for nb_toss=1, logit_res.predict([1,1]), is around 0.5 which corresponds to your simulation. If you don't want to assume a specific nonlinear function for the effect of nb_toss, then one possibility is to add log(nb_toss), nb_toss**2 or similar polynomial nonlinear terms. (note: if nb_toss=0, then we can perfectly predict that y has to be zero.) $\endgroup$ – Josef May 15 '17 at 17:56
  • 2
    $\begingroup$ As whuber pointed out, the true function for your model is specific nonlinear function of nb_toss and the single event probability. The linear function (plus logit link) will be a reasonable approximation over some range. However, not for nb_toss=0. The probability jumps from 0 to 0.5 going from nb_toss=0 to 1 and then increases more slowly. (A symptom of the nonlinearity is that we can perfectly predict the outcome if nb_toss=0 and when nb_toss gets large, the probability is essentially 1. statsmodels gives a perfect separation warning because a large number of predictions are close to 0 or 1 $\endgroup$ – Josef May 15 '17 at 20:52
  • 1
    $\begingroup$ AFAICS, log-link is the correct nonlinear function between nb_toss and the probability for the outcome. e.g. GLM family=Binomial link=log $\endgroup$ – Josef May 15 '17 at 23:33
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ +1 This is an excellent analysis, showing clearly how understanding the underlying situation can lead to an appropriate and tractable model. $\endgroup$ – whuber May 16 '17 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.