# Estimate the effect of dummies with unbalanced data

Imagine having a factory with jobs of different complexity, and a few machines that can execute these jobs; a job can succeed or fail. However, jobs are not assigned completely randomly, there is some unknown process that distributes the jobs over the machines. I am interested in the effect of the machines on the probability of solving a job.

To illustrate, let's look at a sample dataset. I constructed one in Python, the code is given at the bottom. The fraction of solved jobs per machine is given below.

+---------+-------------+
| machine | mean_solved |
+---------+-------------+
| A       |        0.30 |
| B       |        0.34 |
| C       |        0.31 |
+---------+-------------+


If we look at this table, we might be tempted to say that machine A is the worst; it has the lowest success rate. However, if we dive deeper into this data, we see that it performs better or equal than the other machines in all categories, it was just given more difficult tasks:

+---------+-----------+-------------+-----+
| machine | difficult | mean_solved |  n  |
+---------+-----------+-------------+-----+
| A       |         0 |        0.50 |  20 |
|         |         1 |        0.25 |  80 |
| B       |         0 |        0.50 | 100 |
|         |         1 |        0.18 | 100 |
| C       |         0 |        0.50 |  50 |
|         |         1 |        0.12 |  50 |
+---------+-----------+-------------+-----+


Therefore, I want to measure the performance of the machines, assuming they would obtain a similar set of jobs. In this case, I could imagine sampling the data based on the difficult variable. However, in my actual problem, I encounter the problem that the variable difficult is unknown. In practice, I have up to ten different binary and continuous variables, that may influence the difficulty of a problem.

My initial idea was to create dummy variables for the machines, and apply a logistic regression to the data to find the coefficients of the created dummies, but I find that it does not accurately capture the effect that I am looking for.

Does anyone have an idea on how to approach this kind of problem? I'm quite sure there should be plenty of literature or information on this, but I feel I lack the knowledge to know exactly which terms to search for.

import pandas as pd
import numpy as np

np.random.seed(seed=1234)

df_A = pd.DataFrame({'machine' : ['A'] * 100,
'difficult' : [1] * 80 + [0] * 20,
'solved' : list(np.random.choice([0, 1], size=(80,), p=[3./4, 1./4])) + [1]*10 + [0]*10 })
df_B = pd.DataFrame({'machine' : ['B'] * 200,
'difficult' : [1] * 100 + [0] * 100,
'solved' : list(np.random.choice([0, 1], size=(100,), p=[7./8, 1./8])) + [1]*50 + [0]*50 })
df_C = pd.DataFrame({'machine' : ['C'] * 100,
'difficult' : [1] * 50 + [0] * 50,
'solved' : list(np.random.choice([0, 1], size=(50,), p=[7./8, 1./8])) + [1]*25 + [0]*25 })
df = pd.concat([df_A,df_B, df_C])

# table 1
display(df.groupby('machine')['solved'].mean())

# table 2
display(df.groupby(['machine','difficult']).agg(mean_solved=('solved', np.mean), n = ('solved',len)))

# create dummies
df_dummies = pd.get_dummies(df['machine'], prefix ='machine')
df = pd.concat([df, df_dummies], axis=1)
df = df.drop(['machine'],axis=1)


I realized I fell for the dummy variable trap, creating variables with perfect multicollinearity. So initially, I tried

X = df.drop('solved',axis=1)
y = pd.DataFrame(df['solved'])

model = sm.Logit(y,X)
result=model.fit()
print(result.summary2())


which results in

------------------------------------------------------------------------
Coef.    Std.Err.      z    P>|z|      [0.025       0.975]
------------------------------------------------------------------------
const      0.0474 2301815.0453  0.0000 1.0000 -4511474.5405 4511474.6353
difficult -1.5380       0.2414 -6.3719 0.0000       -2.0111      -1.0649
machine_A  0.2805 2301815.0453  0.0000 1.0000 -4511474.3073 4511474.8684
machine_B -0.0394 2301815.0453 -0.0000 1.0000 -4511474.6273 4511474.5485
machine_C -0.1937 2301815.0453 -0.0000 1.0000 -4511474.7816 4511474.3941
========================================================================


I solved it by dropping one of the dummies, taking one of the machines as baseline:

X = X.drop('machine_C',axis = 1)

model = sm.Logit(y,X)
result=model.fit()
print(result.summary2())


which results in

------------------------------------------------------------------
Coef.   Std.Err.     z     P>|z|    [0.025   0.975]
------------------------------------------------------------------
const         -0.1463    0.2440  -0.5997  0.5487  -0.6246   0.3319
difficult     -1.5380    0.2414  -6.3719  0.0000  -2.0111  -1.0649
machine_A      0.4743    0.3381   1.4028  0.1607  -0.1884   1.1369
machine_B      0.1544    0.2789   0.5534  0.5800  -0.3923   0.7010
=================================================================


This approach also improved results on my larger dataset. If anyone else has other approaches to better model this, I am still interested in hearing about them.