Double/Debiased Machine Learning

Conteúdo

1. Introdução
2. Processo Gerador de Dados
3. Regressão Linear
4. Aprendizado De Máquina
5. Double/Debiased Machine Learning
6. Referências

Introdução

Digamos que você trabalhe numa empresa de cobranças e construiu um modelo de Aprendizado de Máquina mega deep poderoso que prevê quase perfeitamente que vai pagar ou não a dívida - no set de teste, claro. Ou que você trabalha numa empresa de logística e consegue prever todos os atrasos que irão acontecer. Você então vai todo contente mostrar seu sucesso pro seu chefe ou para a gerência. E, quem diria!, descobre que eles não estão muito interessado em prever essas coisas. Se você pensar bem até que faz sentido. A gerência está muito mais interessada em descobrir como aumentar o pagamento das dívidas ou como diminuir os atrasos. Nessa caso, previsão pode até ser interesse, mas não é o objetivo principal.

O problema de estimação de impacto está dentro de toda uma literatura sobre Inferência Causal. Há anos que cientistas, principalmente da econometria e bioestatística, estudam inferência causal, mas só recentemente o campo vem ganhando relevância na comunidade de aprendizado de máquina.

Hoje, a maioria dos modelos “de prateleira” de aprendizado de máquina não são pensados para otimizar estimação de impacto ou inferência causal. Assim, quem quer responder perguntas de impacto (“como X impacta y”) ou contrafactuais (“como y seria se X fosse diferente do que é”) precisa ainda implementar modelos de inferência causa na mão.

Aqui, vamos ver uma das técnicas que misturam modelos complexos de aprendizado de máquina com fundamentos teóricos de inferência causal. Como exemplo, vamos tentar descobrir como a dosagem de um medicamento impacta no tempo de recuperação de pacientes doentes.

Processo Gerador de Dados

Imagine que você foi contratado para descobrir o impacto da dosagem de um remédio no tempo de recuperação de pacientes doentes. Parece um problema fácil. Bastaria ver se pessoa que receberam uma alta dosagem recuperam mais ou menos rápido do que pacientes que receberam uma dosagem menos. Seria fácil… Mas os médicos tendem a dar uma dosagem maior para pacientes em estados mais críticos. Por isso, se a gente simplesmente olhar correlações simples iremos obter algo bem contra intuitivo: mais remédio significa mais tempo no hospital.

Para o propósito desse tutorial, é melhor usar dados simulados para os quais sabemos exatamente o efeito causal. I ideia será ver como nosso modelo consegue identificar esse efeito. Os dados são simulados de forma a gerar estruturas de confusão, em que uma variável causa tanto o tratamento como a variável resposta. No caso, temos que tanto Sexo do paciente quanto Severidade da doença causam tanto a recuperação quanto o tratamento. Vamos também colocar algumas relações não lineares que são pouco triviais.

O que é importante perceber no processo acima é que o efeito causal do medicamento é -1, isto é, aumentar a dose do medicamento 1mg diminui o tempo no hospital em um dia.

Regressão Linear

Se você tem contato com a área de econometria ou estatística seu pais devem ter te falado que sempre que você quer controlar para algum fator em dados observacionais você precisa rodar uma regressão linear. Essa é uma boa intuição. Regressão linear tem uma interpretação de controlar para os fatores incluídos no modelo que simula dados experimentais. Assim, é de se esperar que se a gente controlar as variáveis de confusão, Sexo e Severidade, deveríamos ser capazes de achar o efeito causal do medicamento. Em outras palavras, regressão linear nos daria o efeito do medicamento mantidos Sexo e Severidade constante. Certo? Vamos ver…

import statsmodels.formula.api as smf
forumla = "recovery ~ sex + age + severity + medication"
model = smf.ols(forumla, data=df_obs).fit()
model.summary()

Ops… Aparentemente não. Nosso modelo de regressão está estimando um efeito causal de -2.0233 em vez de -1, ou seja, ele está superestimando o efeito causal. Isso acontece porque nosso modelo não está capturando todas as relações não lineares no nosso processo gerador de dados. Se nós fizéssemos um modelo que corresponde com o PGD isso provavelmente não aconteceria e conseguiríamos estimar o impacto causal corretamente. Só que, na vida real, nós raramente conhecemos o PGD.

Por sinal, não linearidade complexas é um dos grandes motivos para usarmos modelos de aprendizado de máquina.

Aprendizado De Máquina

Se não linearidades são o problema, vamos usar um modelo de aprendizado de máquina, certo? Bom, vamos ver…

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score
np.random.seed(42)

X = df_obs.drop("recovery", axis=1)
y = df_obs["recovery"]

model = RandomForestRegressor(max_depth=5, n_estimators=100, n_jobs=-1)

print(cross_val_score(model, X, y, cv=5))

[0.9709928  0.96984636 0.9733455  0.96893817 0.96942233]

Acima, nós treinamos um modelo de floresta aleatória, com 100 árvores e profundidade máxima de 5. Nossa performance num set de validação é de aproximadamente 0,98, o que é um bom indicativo que nosso modelo está bem ajustado.

Para estimar o efeito causal nós precisaremos ser um pouco espertos. Como estamos lidando com um modelo de aprendizado de máquina, o efeito não precisa ser linear nem constante. Também não é possível simplesmente “ler” o modelo para ver como ele está entendendo o funcionamento do grau de medicamento.

O que vamos fazer então é gerar artificialmente alguns níveis de medicamento, de 0 a 4 em intervalos de 0.1mg. Depois, para cada nível, vamos ver a previsão média do modelo na base toda. Por fim, vamos plotar a previsão média contra o nível de medicamento

np.random.seed(42)
model.fit(X, y)
med_level = np.arange(0, 4, 0.1)
y_hats = [model.predict(X.assign(medication = med)).mean() for med in med_level]

plt.plot(med_level, y_hats)
plt.xlabel("Medication Level")
plt.ylabel("Average Prediction of Recovery")
plt.show()

E mais uma vez, parece que nossa ideia não deu certo. O modelo está prevendo que o medicamento diminui o tempo de recuperação, muito pouco. A queda prevista no tempo de recuperação é apenas de 1,99 dias para 1,96 dias e só acontece após o paciente receber 2mg de medicamento.

Double/Debiased Machine Learning

Nossa abordagem final envolve um pouco de teoria de regressão linear e modelos de aprendizado de máquina. Em um artigo de 2017, Chernozhukov et al propõe o uso de ortogonalização usando modelos de AM como uma forma de estimar o efeito causal. Começamos definindo \(X\) o conjunto de variáveis de confusão, no nosso caso idade, sexo de severidade. Definindo \(D\) como a nossa variável de tratamento, no nosso caso o nível de medicamento dado ao paciente. E, finalmente, definimos \(Y\) como a nossa variável resposta, dias de recuperação, no nosso caso.

Nós podemos retirar o efeito de \(X\) em \(D\) construindo um regressor \(\hat{V} = D - m(X)\), em que \(m\) é um modelo de regressão estimador de \(D\). \(\hat{V}\) nesse contexto são os resíduos da regressão de \(D\) em \(X\), de forma que \(\hat{V}\) será ortogonal a \(X\).

Analogamente, nós podemos retirar o efeito de \(X\) em \(Y\) construindo um regressor \(\hat{Y} = Y - l(X)\), em que \(l\) é um modelo de regressão estimador de \(Y\). \(\hat{U}\) nesse contexto são os resíduos da regressão de \(Y\) em \(X\), de forma que \(\hat{U}\) será ortogonal a \(X\).

Com isso, nós podemos construir o estimador do efeito causal da seguinte forma:

Vamos tentar…

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score

X = df_obs.drop(columns=["recovery", "medication"]).values
Y = df_obs["recovery"].values
D = df_obs["medication"].values

np.random.seed(42)
Lhat = RandomForestRegressor(max_depth=5, n_estimators=100, n_jobs=-1).fit(X, Y).predict(X)
Mhat = RandomForestRegressor(max_depth=5, n_estimators=100, n_jobs=-1).fit(X, D).predict(X)

theta = np.dot((D - Mhat), (Y - Lhat)) / np.dot(D - Mhat, D - Mhat)

print("Causal Effect:", theta)

-0.9601128788228591

Agora sim! O efeito estimado do medicamento é de -0.96. Isso é bem próximo do efeito real, que era -1. Mas por que isso funciona? Bom, para isso precisamos ir um pouco mais fundo em teoria econométrica.

Por Que Funciona?

A primeira coisa que precisamos reconhecer é que o estimador é simplesmente o coeficiente \(\hat{\beta_1}\) da seguinte regressão linear univariada

Se você não acredita em mim, pode verificar isso rodando essa regressão no Python e vendo que o coeficiente é exatamente o mesmo que estimamos acima.

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit((D - Mhat).reshape(-1,1), Y - Lhat)

print(model.coef_)

[-0.96011223]

A demonstração matemática é um pouco mais longa, mas você pode conferi-la neste tutorial que eu fiz. De maneira resumida, o que isso nos mostra que que \(\beta_1\) pode ser entendido como a regressão de y em \(D\) após os termos considerados (projetado) a informação em \(X\). Sendo um pouco mais técnico, podemos usar o teorema de Frisch-Waugh-Lovell para mostrar que parâmetro \(\beta_1\) pode ser obtido rodando duas regressões, isto é, regressão por partes. Esse parâmetro é matematicamente igual a regredir os resíduos da regressão de \(X\) em \(D\) (ou seja, \(\hat{U}\)) nos resíduos da regressão de \(X\) em \(Y\) (ou seja, \(\hat{V}\)).

Particionando os Dados

Podemos ainda ter uma última melhora em cima do estimador acima particionando os dados. Isso envolve separar os dados em treino e teste, treinar os modelos de aprendizado de máquina no dataset de treino, mas computar os resíduos usando o dataset de teste.

No código abaixo, nos fazemos esse particionamento no estilo k-fold, onde criamos 10 datasets de treino e 10 datasets de teste e realizamos 10 iterações de treino, teste e estimação do parâmetro \(\theta\). Por fim, tiramos a média de todos os parâmetros estimados para obter uma estimativa final.

import numpy as np
from sklearn.model_selection import KFold

kf = KFold(n_splits=10)

thetas = []

np.random.seed(42)
for train, test in kf.split(X):
    
    Ghat = RandomForestRegressor(max_depth=5, n_estimators=100, n_jobs=-1).fit(X[train],Y[train]).predict(X[test])
    Mhat = RandomForestRegressor(max_depth=5, n_estimators=100, n_jobs=-1).fit(X[train],D[train]).predict(X[test])
    
    Vhat = D[test] - Mhat
    Uhat = Y[test] - Ghat
    
    thetas += [np.mean(np.dot(Vhat,Uhat))/np.mean(np.dot(Vhat,Vhat))]


print("Causal Effect:", np.mean(thetas))    

Causal Effect: -0.977882592667485

O efeito causal estimado é ainda mais próximo do efeito causal real de -1. Isso funciona pois previne problemas de sobre-ajustamento de modelos. Voltando a fórmula do estimador, é intuitivo ver porque isso seria um problema.

Se alguns dos modelos de aprendizado de máquina sobre-ajustarem, ou \(\hat{V}\) ou \(\hat{U}\) serão artificialmente próximos de zero. Se \(\hat{V} < \hat{U}\) porque \(m\) está sobre ajustado, o denominador da fórmula acima seria muito pequeno e o estimador explodiria para um número muito grande. Se \(\hat{U} < \hat{V}\) porque \(g\) está sobre ajustado, o numerador da fórmula acima seria muito pequeno e o estimador seria jogado para zero.

Separas a amostra de estimação da de computação de resíduo previsse esses problemas.

Referências

Este post implementa o modelo descrito em Double/Debiased Machine Learning for Treatment and Causal Parameters, de Chernozhukov et al, 2016. Uma boa referência para quem quiser aprender mais sobre esse modelo é a palestra dada por Victor Chernozhukov. Como de costume, o código deste post está no meu GitHub.

• ← Previous Post
• Next Post →