Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como resolver um problema de regressão linear com múltiplas variáveis explicativas utilizando a abordagem Bayesiana?

0 votos
48 visitas

1 Resposta

0 votos
respondida Jul 2 por Pedro Antero (26 pontos)  
editado Jul 6 por Pedro Antero

Essa resposta será construída tendo-se por base os conceitos apresentados nos Capítulos 2 e 3 do livro Bayesian Econometrics, de Gary Koop e o arcabouço teórico discutido na resposta da questão “Como utilizar o Python para ilustrar a aplicação da teoria Bayesiana sobre o caso mais simples de regressão linear?” (disponível em http://prorum.com/index.php/3243/utilizar-ilustrar-aplicacao-bayesiana-simples-regressao?show=3243#q3243). As equações que serão empregadas nesta explicação são idênticas àquelas discutidas na referida questão, que trata do caso de apenas uma variável explicativa. A única diferença é que aqui utilizaremos a notação matricial para representar as variáveis multidimensionais. Daqui pra frente, sempre que formos fazer menção às equações empregadas na resposta do caso unidimensional, utilizaremos a notação “UV”.

Tomemos um modelo de regressão linear com \( k \) variáveis explicativas e \( i = 1 ... N\) dado por:
\( y_i = \beta_1 + \beta_2 x_{i2} + ... + \beta_k x_{ik} + \epsilon_i \tag{1}\)

Para simplificar a notação, utilizaremos a notação matricial “empilhando” as observações da seguinte forma:

A imagem será apresentada aqui.

Assim, a Equação (1) pode ser reescrita como:
\( y = X \beta + \epsilon \tag{2}\)

Utilizando as mesmas hipóteses sobre a distribuição de \( \epsilon \) e sobre a independência entre \( \epsilon \) e \( X \), é possível mostrar que a função de máxima verossimilhança para a regressão linear da Equação (1) será dada por:

\( p(y| \beta, h) = \frac{h^{\frac{1}{2}}}{(2 \pi)^{\frac{N}{2}}} exp \left[ - \frac{h}{2} (\beta - \hat{ \beta})' X'X (\beta - \hat{ \beta}) \right] \left[ h^{\frac{\nu}{2}} exp \left( - \frac{h \nu}{ 2 s^{-2}} \right) \right] \tag{3} \)

Empregando a mesma sequência lógica adotada para o caso de uma variável, é fácil mostrar que as relações definidas na Equação (UV-6), para uma prior natural conjugada, e nas Equações (UV-7) a (UV-17), para a posterior, se mantêm. A única diferença é que, agora, o termo \( X \) contém mais de uma variável explicativa. Portanto, todas as operações descritas nestas Equações devem ser interpretadas sob a lógica matricial. Por exemplo, a operação \( \frac{1}{\sum{x_i^2}}\) deve ser entendida como \( (X'X)^{-1} \).

Seguindo a Seção 3.9 do livro Bayesian Econometrics, de Gary Koop, vamos ilustrar empiricamente a utilização da abordagem Bayesiana sobre uma regressão linear multivariada por meio do banco de dados disponibilizado por Anglin e Gencay (1996). A referida base de dados ilustra as seguintes características de \( N = 546 \) casas no Canadá:

  • \( y = \) preço de venda da i-ésima casa, em dólares Canadenses
  • \( x_{i2} = \) tamanho do terreno da i-ésima casa, em pés quadrados
  • \( x_{i3} = \) número de quartos da i-ésima casa
  • \( x_{i4} = \) número de banheiros da i-ésima casa
  • \( x_{i5} = \) número de andares da i-ésima casa

Observe que utilizaremos o coeficiente \( \beta_{i1} \) para representarmos o intercepto. Isto é, \( x_{i1} =1 \) para \( i = 1...N\).

O código abaixo realiza a importação da base de dados e de algumas bibliotecas que serão utilizadas adiante.

import xlrd
import numpy as np  # imports numpy library
from scipy import stats
from scipy.stats import t

file_location = "C:/Users/pantero/PycharmProjects/metodos_computacionais/src/ag_data.xlsx"
workbook = xlrd.open_workbook(file_location)
sheet = workbook.sheet_by_index(0)

data = [[sheet.cell_value(r, c) for c in range(sheet.ncols)] for r in range(sheet.nrows)]
data = np.asarray(data, float)
N = sheet.nrows
k = 5
y = data[:, 0]
X = data[:, 1:k]
X = np.c_[np.ones((N, 1)), X]

Passemos agora a definição de uma prior para nossa análise. Para tanto, conforme discutido no caso univariado, temos que estabelecer os chamados hiperparâmetros da prior: \( \underline{\beta}, \underline{V}, \underline{s}^{-2}, \underline{\nu} \).

Vamos supor que um corretor imobiliário na região nos indique que o preço das casas de Windsor variam, em média, entre $ 50.000 e $ 150.000. Portanto, um bom modelo de ajuste não deveria ter erros superiores a $ 10.000. Ou seja, devemos esperar que \( \sigma \approx 5000) \). Como mantivemos a hipótese de que \( \epsilon_i \sim N(0, 1) \), então os erros serão menores que \( 1.96 \times 5000 = $ 9800\), em valores absolutos, para \( \sigma = 5000 \). Ou seja, \( \frac{1}{\sigma^2} = \frac{1}{5000^2} = 4,0 \times 10^{-8} \) seria um bom chute para o parâmetro \( \underline{s}^{-2} \) da prior. Como esse chute é um pouco grosseiro, não seria prudente adotar um peso muito grande para esse parâmetro. Assim, escolheremos \( \underline{\nu} = 5 << N = 546\). De modo aproximado, isso quer dizer que o peso de <em>prior deve ser de cerca de 1% daquele atribuído aos dados. Vamos supor ainda, que uma escolha razoável de \( \underline{\beta} \) possa ser dada por:

A imagem será apresentada aqui.

Como o "chute" utilizado para definição dos parâmetros acima também foi um tanto grosseiro, é prudente que se associe uma variância alta a esses parâmetros. Por exemplo, podemos adotar:
A imagem será apresentada aqui.

Observe que todas as covariâncias foram admitidas como nulas. Isso normalmente é feito no desenho de priors, já que é difícil assumirmos previamente valores para estes parâmetros.

A partir das definições anteriores de \(\underline{s}\), \(\underline{\nu}\), \( var(\beta) \) e da Equação (UV-15) chegamos, então, ao seguinte resultado para \( \underline{V} \):

A imagem será apresentada aqui.

O código abaixo nos permite chegar aos resultados apresentados.

# Hyper-parameters of the informative prior
beta_uline = np.array([0, 10, 5000, 10000, 10000]).reshape((k, 1))
var_beta = [10000 ** 2, 25, 2500 ** 2, 5000 ** 2, 5000 ** 2]
var_beta_m = np.zeros((k, k))
for r in range(len(var_beta_m[0][:])):
    var_beta_m[r][r] = var_beta[r]
nu_uline = k
s_minus_sqd_uline = 4 * (10 ** (-8))
V_uline = np.asarray((nu_uline - 2) / (nu_uline * (1/s_minus_sqd_uline)) * var_beta_m)

Com isso, completamos a especificação da nossa prior natural conjugada.

As Tabelas 1 e 2 abaixo ilustram os resultados das características dos parâmetros, calculados conforme a teoria já discutida no caso univariado.

A imagem será apresentada aqui.
Tabela 1 - Características de \( \beta \)

A imagem será apresentada aqui.
Tabela 2 - Características de \( h \)

O código abaixo foi utilizado para o cálculo dos parâmetros acima.

nu = N - k
beta_hat = np.dot(np.dot(np.linalg.inv(np.dot(X.transpose(), X)), X.transpose()), y)
s_sqd = np.dot(np.transpose(y - np.dot(X, beta_hat)), y - np.dot(X, beta_hat)) / nu

V_bar_nip = np.linalg.inv(np.dot(X.transpose(), X))
beta_bar_nip = beta_hat
nu_bar_nip = N
s_bar_sqd_nip = (nu * s_sqd) / nu_bar_nip

e_beta_post_nip = beta_bar_nip
var_beta_post_nip = ((nu_bar_nip * s_bar_sqd_nip) / (nu_bar_nip - 2)) * V_bar_nip  # Eq. (15)
std_dev_beta_post_nip = np.sqrt(np.diagonal(var_beta_post_nip))

# print('e_beta_post_nip')
# print(e_beta_post_nip)
# print('std_dev_beta_post_nip')
# print(std_dev_beta_post_nip)

e_h_post_nip = 1 / s_bar_sqd_nip  # Eq. (16)
var_h_post_nip = (2 * (1 / s_bar_sqd_nip)) / nu_bar_nip  # Eq. (17)
std_dev_h_post_nip = np.sqrt(var_h_post_nip)

# print('e_h_post_nip')
# print(e_h_post_nip)
# print('std_dev_h_post_nip')
# print(std_dev_h_post_nip)

e_beta_prior = beta_uline
var_beta_prior = ((nu_uline * (1 / s_minus_sqd_uline)) / (nu_uline - 2)) * V_uline
std_dev_beta_prior = np.sqrt(np.diagonal(var_beta_prior))

# print('e_beta_prior')
# print(e_beta_prior)
# print('std_dev_beta_prior')
# print(std_dev_beta_prior)

e_h_prior = s_minus_sqd_uline
var_h_prior = (2 * s_minus_sqd_uline) / nu_uline
std_dev_h_prior = np.sqrt(var_h_prior)  # há um erro no valor reportado pelo livro, que fornece a variância ao invés do desvio

# print('e_h_prior')
# print(e_h_prior)
# print('std_dev_h_prior')
# print(std_dev_h_prior)

V_bar_ip = np.linalg.inv((np.linalg.inv(V_uline) + np.dot(X.transpose(), X)))
beta_hat = beta_hat.reshape((k, 1))
beta_bar_ip = np.dot(V_bar_ip, np.dot(np.linalg.inv(V_uline), beta_uline) + np.dot(np.dot(X.transpose(), X), beta_hat))  # Eq. (9)
nu_bar_ip = nu_uline + N  # Eq. (10)
s_bar_sqd_ip = (1 / nu_bar_ip) * ((nu_uline * (1 / s_minus_sqd_uline)) + (nu * s_sqd) + np.dot(np.dot(np.transpose(beta_hat - beta_uline), np.linalg.inv(V_uline + np.linalg.inv(np.dot(X.transpose(), X)))), beta_hat - beta_uline))  # Eq. (11)

e_beta_post_ip = beta_bar_ip
var_beta_post_ip = ((nu_bar_ip * s_bar_sqd_ip) / (nu_bar_ip - 2)) * V_bar_ip
std_dev_beta_post_ip = np.sqrt(np.diagonal(var_beta_post_ip))

# print('e_beta_post_ip')
# print(e_beta_post_ip)
# print('std_dev_beta_post_ip')
# print(std_dev_beta_post_ip)

e_h_post_ip = 1 / s_bar_sqd_ip
var_h_post_ip = (2 * (1 / s_bar_sqd_ip)) / nu_bar_ip
std_dev_h_post_ip = np.sqrt(var_h_post_ip)

# print('e_h_post_ip')
# print(e_h_post_ip)
# print('std_dev_h_post_ip')
# print(std_dev_h_post_ip)

Como se pode perceber, os resultados são idênticos àqueles apresentados na Seção 3.9 o livro Bayesian Econometrics, de Gary Koop. (A única exceção é o parâmetro de desvio padrão da prior. Para esse caso, há um pequeno erro no livro, que apresenta a variância, e não o desvio do parâmetro.)

Passemos agora ao estudo da comparação entre modelos, já discutida para o caso univariado.

Inicialmente, tomemos 2 modelos concorrentes \(M_1\) restrito e \(M_2\), irrestrito. Vamos supor que a restrição de \(M_1\) envolve uma desigualdade da forma:
\( M_1: R \beta \geq r \tag{MC-1}\)

Em que \( R_{J \times k} \) é uma matriz dada de restrições e \( r \) um vetor de dimensão \(J \).

Para modelos deste tipo a razão de probabilidades da posterior, definida pela Equação (UV-MC-5) pode ser escrita por:

\( PO_{12} = \frac{p(M_1 | y)}{p(M_2 | y)} = \frac{p(R \beta \geq r | y)}{p(R \beta \ngeq r | y)} \tag{MC-2} \)

Como \( \beta \) segue uma distribuição t de Student multivariada (nos moldes na Equação (UV-12)), pode-se mostrar que \( p(R \beta | y ) \) também seguirá uma distribuição deste tipo. Dessa forma, o cálculo de \( p(R \beta > r | y ) \) é bastante trivial.

Consideremos agora um outro tipo de controle, dado por uma restrição de "igualdade". A Equação abaixo sintetiza esse tipo de restrição para \( j = 1,2\) modelos concorrentes.

\( M_j: y_j = X_j \beta_{(j)} + \epsilon_j \tag{MC-3} \)

Novamente respaldados pela discussão apresentada no caso univariado, podemos reescrever a Equação (UV-MC-5) como:

\( PO_{12} = \frac{c_1 \left( \frac{|\bar{V}_1|}{|\underline{V}_1|} \right)^{\frac{1}{2}} (\bar{\nu}_1 \bar{s}_1^2)^{-{\frac{\bar{\nu}_1}{2}}} p(M_1)}{c_1 \left( \frac{|\bar{V}_2|}{|\underline{V}_2|} \right)^{\frac{1}{2}} (\bar{\nu}_2 \bar{s}_2^2)^{-{\frac{\bar{\nu}_2}{2}}}p(M_2)} \tag{MC-4} \)

Com:

\( c_j = \frac{\Gamma \left( \frac{\bar{\nu}_j}{2} \right) (\underline{\nu}_j \underline{s}_j^2)^{-\frac{\underline{\nu}_j}{2}} }{\Gamma \left( \frac{\underline{\nu}_j}{2} \right)\pi^{\frac{N}{2}} } \tag{MC-5} \)

Para \( j = 1,2 \)

Na comparação entre modelos diferentes, também é comum se utilizar a noção de Highest Posterior Density Intervals (HPDI). Tomemos novamente dois modelos concorrentes, em que \( M_2 \) é irrestrito e \(M_1\) apresenta uma restrição de igualdade do tipo:
\( M_1: \beta_j = 0 \tag{MC-6}\)

Usando novamente o resultado de que \( p( R \beta | y) \) segue uma distribuição t de Student multivariada, é muito fácil calcularmos intervalos de confiança para o modelo restrito.

Para exemplificarmos a discussão acima, tomemos, por exemplo as restrições \( p(\beta_j >0 | y) \) e \( p(\beta_j =0 | y) \), com \( j = 1...k \).

A Tabela 3 ilustra as características de \( \beta \) para cada uma destas restrições utilizando as noções discutidas acima para o cálculo da razão de probabilidades da posterior (\( PO_{12}\)) e do HPDI.

A imagem será apresentada aqui.
Tabela 3 - Comparação de modelos de restrição à \( \beta \)

O código abaixo foi utilizado para obtenção dos resultados.

y = y.reshape((N, 1))
for i in range(0, k):

# Column 1
    print('p(beta_', i, '> 0 |y')
    print(1 - stats.t.cdf(0, nu_bar_ip, e_beta_post_ip[i], np.sqrt(var_beta_post_ip[i][i])))

# Column 2
    print('95% HPDI beta_', i, '> 0 |y')
    print(stats.t.interval(0.95, nu_bar_ip, e_beta_post_ip[i], np.sqrt(var_beta_post_ip[i][i])))

# Column 3
    print('99% HPDI beta_', i, '> 0 |y')
    print(stats.t.interval(0.99, nu_bar_ip, e_beta_post_ip[i], np.sqrt(var_beta_post_ip[i][i])))

# Column 4
    X_rest = np.delete(X, i, 1)
    beta_hat_rest = np.dot(np.dot(np.linalg.inv(np.dot(X_rest.transpose(), X_rest)), X_rest.transpose()), y)
    beta_uline_rest = np.delete(beta_uline, i, 0)
    V_uline_rest = np.delete(V_uline, i, 0)
    V_uline_rest = np.delete(V_uline_rest, i, 1)
    V_bar_ip_rest = np.linalg.inv((np.linalg.inv(V_uline_rest) + np.dot(X_rest.transpose(), X_rest)))
    s_sqd_rest = np.dot(np.transpose(y - np.dot(X_rest, beta_hat_rest)), y - np.dot(X_rest, beta_hat_rest)) / nu
    s_bar_sqd_ip_rest = (1 / nu_bar_ip) * ((nu_uline * (1 / s_minus_sqd_uline)) + (nu * s_sqd_rest) + np.dot(np.dot(np.transpose(beta_hat_rest - beta_uline_rest), np.linalg.inv(V_uline_rest + np.linalg.inv(np.dot(X_rest.transpose(), X_rest)))), beta_hat_rest - beta_uline_rest))  # Eq. (11)
    PO12_num = (( np.linalg.det(V_bar_ip_rest) / np.linalg.det(V_uline_rest)) ** (1 / 2)) * ((s_bar_sqd_ip_rest * 10 ** (-8)) ** (- nu_bar_ip / 2))
    PO12_den = (( np.linalg.det(V_bar_ip) / np.linalg.det(V_uline) ) ** (1 / 2)) * ((s_bar_sqd_ip * 10 ** (-8)) ** (- nu_bar_ip / 2))

    print(PO12_num / PO12_den)

Conforme pode ser observado, os resultados são idênticos aos apresentados na Seção 3.9 do livro Bayesian Econometrics, de Gary Koop. As pequenas diferenças observadas nos intervalos do HDPI se devem ao algoritmo utilizado para o cálculo da CDF da função t de Student. Ao que parece, Gary Koop utilizou o MATLAB, enquanto nossas computações foram realizadas em Python.

Por fim, iremos calcular as características da posterior para \( \beta_2 \) por meio do método de Monte Carlo. Iremos fazer isso para \( S = {10, 100, 1.000, 10.000 e 100.000} \) replicações.

A imagem será apresentada aqui.

O seguinte código foi utilizado para obtenção dos resultados:

for i in range(0, 5):
    mean_acc = 0
    S = 10 * (10 ** i)
    for j in range(S):
        mean_acc = mean_acc + t.rvs(nu_bar_ip , e_beta_post_ip[1], np.sqrt(var_beta_post_ip[1][1]))
    mean_mc = mean_acc / S
    print('mean_mc')
    print(mean_mc)

for i in range(0, 5):
    var_acc = 0
    S = 10 * (10 ** i)
    for k in range(S):
        var_acc = var_acc + ((t.rvs(nu_bar_ip, e_beta_post_ip[1], np.sqrt(var_beta_post_ip[1][1])) - (mean_mc)) ** 2)
    std_dev_mc = np.sqrt(var_acc) / np.sqrt(S)
    print('std_dev_mc')
    print(std_dev_mc)
    print('num_std_dev_mc')
    print(std_dev_mc / np.sqrt(S))

Analiticamente, pode-se mostrar que a média e o desvio da posterior para \(\beta_2 \) são 5,4316 e 0,3662 (vide Seção 3.9 do livro Bayesian Econometrics, de Gray Koop). Como podemos perceber, o Método de Monte Carlo consegue aproximar bem estes momentos.

comentou Jul 10 por Carlos Eduardo Véras (11 pontos)  
Lendo a sua explicação em conjunto com o código consegui entender bem melhor que lendo o livro. Parabéns!
...