Regressão Múltipla

Autor

Douglas Braga

Nota

Esta seção corresponde à seção 3.2 de James et al. (2023).

Advertising traz três mídias, e a seção 8.1 usou só uma. A pergunta mais direta é repetir a mesma reta duas vezes mais — uma para radio, outra para jornal — e ler os três coeficientes lado a lado. Esta seção começa por essa pergunta, e termina desfazendo a resposta que ela parece dar.

Jornal, sozinho, parece importar

X_jornal = propaganda[["jornal"]]
modelo_jornal = LinearRegression().fit(X_jornal, y)
beta1_jornal_sozinho = float(modelo_jornal.coef_[0])

round(beta1_jornal_sozinho, 4), round(beta1_jornal_sozinho * 1000, 1)

(0.0547, 54.7)

Ajustada sozinha contra jornal, a reta de mínimos quadrados dá um coeficiente de 0,0547: cada mil dólares a mais gastos em jornal está associado, em média, a 54,7 unidades a mais vendidas. É um coeficiente positivo, do mesmo tipo que a seção 8.1 encontrou para tv — nada nesta regressão isolada distingue jornal das outras duas mídias.

Só que ajustar jornal sozinho ignora tv e radio por completo. Se as três mídias variam juntas de mercado para mercado, o coeficiente de uma pode estar levando crédito por vendas que outra produziu — e a única forma de saber é colocar as três na mesma equação.

As três mídias na mesma equação

O modelo múltiplo estende a mesma ideia de mínimos quadrados a mais de um preditor: em vez de uma reta, um hiperplano que minimiza a soma dos quadrados dos resíduos sobre tv, radio e jornal ao mesmo tempo.

X_multiplo = propaganda[["tv", "radio", "jornal"]]
modelo_multiplo = LinearRegression().fit(X_multiplo, y)

coeficientes = pd.Series(modelo_multiplo.coef_, index=modelo_multiplo.feature_names_in_)
coeficientes_mil_dolares = (coeficientes * 1000).round(1)

coeficientes.round(4), coeficientes_mil_dolares

(tv        0.0458
 radio     0.1885
 jornal   -0.0010
 dtype: float64,
 tv         45.8
 radio     188.5
 jornal     -1.0
 dtype: float64)

feature_names_in_ guarda os três nomes de coluna que o DataFrame carregava, na mesma ordem dos coeficientes — por isso dá para juntar os dois numa única Series em vez de decorar qual número é qual. tv sai com 0,0458 (45,8 por mil dólares) e radio com 0,1885 (188,5 por mil dólares), os dois positivos. jornal sai com -0,0010 (-1,0 por mil dólares): praticamente zero, e de sinal oposto ao 0,0547 que a regressão sozinha tinha encontrado para ele.

Num modelo múltiplo, cada coeficiente \(\hat\beta_j\) se lê como o efeito médio sobre a resposta de aumentar o preditor \(X_j\) em uma unidade, mantendo todos os outros preditores fixos. É essa cláusula — mantendo os demais fixos — que separa a leitura de um coeficiente múltiplo da de uma regressão simples, e é ela que muda o que se pode dizer sobre jornal.

Lido dessa forma: gastar mais mil dólares em tv, mantendo radio e jornal fixos, está associado a 45,8 unidades a mais de venda; mais mil dólares em radio, mantendo tv e jornal fixos, a 188,5 unidades a mais. Mais mil dólares em jornal, mantendo tv e radio fixos, está associado a uma variação de -1,0 unidade — perto o bastante de zero para não sobrar efeito de jornal depois que as outras duas mídias já estão na conta.

Por que jornal some: quem anda com quem

A explicação não está em mais uma regressão — está em como as três mídias se relacionam entre si, antes de qualquer venda entrar na conta.

correlacoes = propaganda[["tv", "radio", "jornal", "vendas"]].corr()
correlacoes.round(4)

	tv	radio	jornal	vendas
tv	1.0000	0.0548	0.0566	0.7822
radio	0.0548	1.0000	0.3541	0.5762
jornal	0.0566	0.3541	1.0000	0.2283
vendas	0.7822	0.5762	0.2283	1.0000

corr_jornal_radio = float(correlacoes.loc["jornal", "radio"])
corr_jornal_tv = float(correlacoes.loc["jornal", "tv"])
jornal_mais_correlacionado_com_radio = bool(corr_jornal_radio > corr_jornal_tv)

round(corr_jornal_radio, 4), round(corr_jornal_tv, 4), jornal_mais_correlacionado_com_radio

(0.3541, 0.0566, True)

jornal correlaciona com radio a 0,3541 e com tv a só 0,0566 — jornal_mais_correlacionado_com_radio confirma que a primeira supera a segunda. Mercados que gastam mais em rádio tendem a gastar mais em jornal também: os dois investimentos sobem e descem juntos, sem que isso implique nada causal entre eles.

É essa correlação que explica a reviravolta do coeficiente de jornal, de 0,0547 sozinho para -0,0010 na companhia das outras duas mídias. Se é o rádio — não o jornal — que de fato move vendas, então nos mercados onde se gasta mais em rádio as vendas tendem a ser maiores, e a tabela de correlação mostra que esses são os mesmos mercados que gastam mais em jornal. Uma regressão que olha só para jornal, sem radio por perto, não tem como separar as duas coisas: ela atribui a jornal parte do crédito que é do rádio. Colocar as duas mídias na mesma equação é o que permite distinguir — e é aí que o coeficiente de jornal cai a praticamente zero.

R² e RSE: o quanto o modelo múltiplo melhora

R² não muda de definição com mais preditores; o RSE muda de denominador. Com \(p\) preditores, a raiz é de \(\text{RSS}/(n - p - 1)\) — um grau gasto pelo intercepto e um por coeficiente —, e o \(n-2\) da seção anterior é o caso \(p = 1\) dessa mesma forma.

X_tv = propaganda[["tv"]]
modelo_tv = LinearRegression().fit(X_tv, y)
yhat_tv = modelo_tv.predict(X_tv)

n = len(y)
rss_tv = float(((y - yhat_tv) ** 2).sum())
rse_tv = float(np.sqrt(rss_tv / (n - 2)))
r2_tv = float(r2_score(y, yhat_tv))

round(rse_tv, 4), round(r2_tv, 4)

(3.2587, 0.6119)

yhat_multiplo = modelo_multiplo.predict(X_multiplo)

p = X_multiplo.shape[1]
rss_multiplo = float(((y - yhat_multiplo) ** 2).sum())
rse_multiplo = float(np.sqrt(rss_multiplo / (n - p - 1)))
r2_multiplo = float(r2_score(y, yhat_multiplo))

round(rse_multiplo, 4), round(r2_multiplo, 4)

(1.6855, 0.8972)

comparacao = pd.DataFrame(
    {"R²": [r2_tv, r2_multiplo], "RSE": [rse_tv, rse_multiplo]},
    index=["tv sozinho", "tv + radio + jornal"],
).round(4)
comparacao

	R²	RSE
tv sozinho	0.6119	3.2587
tv + radio + jornal	0.8972	1.6855

O modelo de tv sozinho explica 0,6119 da variância de vendas, com erro típico de 3,2587 mil unidades. Somar radio e jornal sobe o R² para 0,8972 e derruba o RSE para 1,6855: as três mídias juntas explicam mais da variância de vendas do que tv sozinho, e erram menos a cada previsão. A régua não muda — R² e RSE continuam sendo as mesmas duas medidas da seção anterior —, só o modelo que elas avaliam.

A superfície ajustada com tv e radio

Como jornal quase não muda a previsão, a superfície que os coeficientes de tv e radio desenham já carrega quase toda a informação do modelo múltiplo — e, com só dois preditores, dá para desenhar essa superfície inteira.

X_tv_radio = propaganda[["tv", "radio"]]
modelo_tv_radio = LinearRegression().fit(X_tv_radio, y)
yhat_tv_radio = modelo_tv_radio.predict(X_tv_radio)

r2_tv_radio = float(r2_score(y, yhat_tv_radio))
round(r2_tv_radio, 4)

0.8972

grade_tv = np.linspace(propaganda["tv"].min(), propaganda["tv"].max(), 60)
grade_radio = np.linspace(propaganda["radio"].min(), propaganda["radio"].max(), 60)
malha_tv, malha_radio = np.meshgrid(grade_tv, grade_radio)
grade = pd.DataFrame({"tv": malha_tv.ravel(), "radio": malha_radio.ravel()})
superficie = modelo_tv_radio.predict(grade).reshape(malha_tv.shape)

azuis = LinearSegmentedColormap.from_list(
    "azuis", plt.get_cmap("Blues")(np.linspace(0.60, 0.74, 256))
)

fig, ax = plt.subplots()
mapa = ax.contourf(malha_tv, malha_radio, superficie, levels=14, cmap=azuis)
ax.scatter(propaganda["tv"], propaganda["radio"], color="C1", s=16, edgecolor="white", linewidth=0.5)
ax.set_xlabel("tv (milhares de dólares)")
ax.set_ylabel("radio (milhares de dólares)")
barra = fig.colorbar(mapa, ax=ax, shrink=0.9, pad=0.02)
barra.set_label("vendas prevista (mil unidades)")
plt.tight_layout()
plt.show()

Figura 45.1: Vendas previstas para toda combinação de investimento em tv e radio, pelo modelo ajustado com as duas mídias. Os pontos são os duzentos mercados observados.

Cada curva de nível reúne as combinações de tv e radio que o modelo prevê com a mesma venda; a superfície cresce para a direita e para cima, do mesmo jeito que os dois coeficientes positivos calculados acima já anunciavam.

O que a superfície plana ainda erra

O modelo com tv e radio explica 0,8972 da variância de vendas — o mesmo valor, até a quarta casa, do modelo de três mídias da tabela acima, já que jornal quase não muda a previsão —, mas sobra no resíduo um padrão que essa superfície plana não captura.

df_resid = pd.DataFrame({
    "previsto": yhat_tv_radio,
    "residuo": y.to_numpy() - yhat_tv_radio,
})
df_resid["terco"] = pd.qcut(df_resid["previsto"], 3, labels=["baixo", "medio", "alto"])

medias_por_terco = df_resid.groupby("terco", observed=True)[["previsto", "residuo"]].mean()
medias_por_terco.round(4)

	previsto	residuo
terco
baixo	8.5201	0.4307
medio	13.9244	-0.8077
alto	19.6216	0.3650

Dividindo os duzentos mercados em três grupos pelo valor previsto — o terço com a previsão mais baixa, o do meio, o mais alto —, o resíduo médio não fica perto de zero nos três grupos: 0,4307 no terço mais baixo, -0,8077 no do meio, e de volta a 0,3650 no mais alto. O sinal do erro típico muda com a faixa de previsão, em vez de se espalhar ao acaso ao redor de zero — um padrão que uma superfície plana, por definição, não reproduz.

fig, ax = plt.subplots()
ax.axhline(0, color="C2", linewidth=1, linestyle="--")
ax.scatter(df_resid["previsto"], df_resid["residuo"], color="C0", s=14, alpha=0.6)
ax.plot(
    medias_por_terco["previsto"],
    medias_por_terco["residuo"],
    color="C1",
    linewidth=2,
    marker="o",
    markersize=7,
)
ax.set_xlabel("vendas prevista (mil unidades)")
ax.set_ylabel("resíduo (vendas observada − prevista)")
plt.tight_layout()
plt.show()

Figura 45.2: Resíduo contra o valor previsto, no modelo de tv e radio. A linha tracejada marca resíduo zero; os três marcadores maiores são a média de cada terço de previsão — abaixo de zero no meio, acima nas duas pontas.

Essa curvatura — negativa no meio da faixa de previsão, positiva nas duas pontas — é o sinal de que tv e radio não agem de forma independente sobre vendas. A seção 8.5 retoma este mesmo modelo para mostrar o que muda ao deixar as duas mídias interagirem.

James, Gareth, Daniela Witten, Trevor Hastie, Robert Tibshirani, e Jonathan Taylor. 2023. An Introduction to Statistical Learning with Applications in Python. Springer.