1 Introdução

A base proposta, Advertising, apresenta variáveis referentes a vendas associadas a veiculações de marteking. Tem como característica variáveis contínuas estritamente positivas. É de interesse avaliar se é possível ajustar um modelo bayesiano gamma através da covariável jornal. Por fim, comparar os resultados com uma árvore de decisão.

2 Modelo proposto

O modelo ajustado será um Modelo Linear Generalizado através da abordagem bayesiana, com a seguinte equação \(Y_i = \beta_0+\beta_1x_i+\beta_2x_i+\beta_3x_i + e_i\) ; \(e_i \sim N(0,\sigma^2)\)

3 Verossimilhança

Posto que a distribuiçãoo da variável resposta não pode ser assumida como normal, o conceito de subjetividade da filosofia Bayesiana nos permite atribuir uma distribuição para verossimilhan?a \(Y|\phi, \beta_0,\beta_1,\beta_2, \sigma^2\), onde todos os parâmetros mencionados são desconhecidos.

Logo, será assumida para a mesma uma distribuiçao Gamma reparametrizada, com parâmetro de escala \(\phi\) e parâmetro de dispersão \(\cfrac{\phi}{\mu}\)

De modo que \(\mu = e^{\beta_0+\beta_1x_i+\beta_2x_i+\beta_3x_i}\), isto é, o preditor linear de uma regressão linear gama usual supondo a inclusão na família exponencial para garantia das propriedades de convergência.

4 Prioris hierárquicas

Em inferência bayesiana, quantidades desconhecidas são atribuídas com distribuição em probabilidade, chamadas de priori.

Para os casos estudados:

\(\phi \sim Gamma(1;0,1)\)

\(\beta_0 \sim Gamma(1;0,1)\)

\(\beta_1 \sim Gamma(1;0,1)\)

\(\beta_2 \sim Gamma(1;0,1)\)

\(\beta_3 \sim Gamma(1;0,1)\)

\(\sigma^2 \sim Gamma(1;0,1)\)

Atribui-se os seguintes parâmetros de forma a obter-se média 10 e variãncia 100, garantindo o que é chamado de priori não informativa, utilizada quando não há informação sob a base trabalhada, o que é o caso.

5 Aproximação de Laplace

Segue do Teorema de Bayes que

\(\pi(\theta|x) = \cfrac{\pi(x|\theta)\pi(\theta)}{\pi(x)}\)

\(\propto\)

\(\pi(x|\theta)\pi(\theta)\)

Onde

\(\pi(\theta|x)\) é um vetor paramétrico de distribuição a posteriori, desconhecida

\(\pi(x|\theta)\) é a verossimilhançaa dos dados, observado

\(\pi(x)\) é um vetor paramétrico de distribuições à priori

Em geral, não é possível a obtenção direta da distribuição a posteriori, que é o grande interesse. Para estas situações, a metodologia de Aproximação de Laplace oferece uma convergência garantida via ajuste uma normal multivariada que assemelha-se a posteriori estudada, de tal forma que:

\(\pi(\theta_p|x) \sim N_p(\theta^*,\sum)\)

Onde a média desta normal multivariada é obtida através da moda de cada parâmetro e a matriz de variância e covariância é estimada a partir da função escore.

6 Advertising

rm(list=ls())
base <- read.csv("C:\\Users\\raphael.fernandez\\Downloads\\Advertising.csv")
base$X = NULL
y=base$sales
x=base$newspaper

7 Ajuste do modelo

set.seed(2019)
model <- function(p,y,x) {
    log_lik <- sum(dgamma(y,p["beta0"]+p["beta1"]*x,sqrt(p["sigma2"]), log = T))
    log_post <- log_lik + dgamma(p["beta0"], 1, 0.1, log = T) + dgamma(p["beta1"], 1, 0.1, log = T)       + dgamma(p["sigma2"],1,0.01, log = T) 
    log_post
}

#Valores iniciais
inits <- c(beta0 = 1, beta1 = 1, sigma2 = 1)

#Ajuste baseado na otimização 
#control = list(fnscale = -1) para maximizar ao invés de minimizar
#hessian = T para retornar a matriz hessiana
fit <- optim(inits, model, control = list(fnscale = -1), hessian = TRUE, y = y,x = x)

#Valores que precisaremos para obter
#Fit$par é o máximo da distribuição de interesse e logo a média da distribuição normal que será utilizada para aproximar
param_mean <- fit$par
#A inversa da negativa da hessiana retorna a matriz de covariância (tem a ver com função score)
param_cov_mat <- solve(-fit$hessian)

round(param_mean, 2)

##  beta0  beta1 sigma2 
##   6.16   0.03   0.25

round(param_cov_mat, 3)

##         beta0  beta1 sigma2
## beta0   0.423 -0.001  0.029
## beta1  -0.001  0.000  0.000
## sigma2  0.029  0.000  0.003

#Gerando amostras da distribuição a posteriori
samples <- rmvnorm(10000, param_mean, param_cov_mat)
beta0    = samples[,1]
beta1    = samples[,2]
sigma2   = samples[,3]

samples = as.mcmc(samples)
summary(samples)

## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean       SD  Naive SE Time-series SE
## beta0  6.15174 0.651095 6.511e-03      6.511e-03
## beta1  0.02707 0.008972 8.972e-05      8.972e-05
## sigma2 0.24747 0.050631 5.063e-04      5.063e-04
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%     75%   97.5%
## beta0  4.87027 5.71505 6.14398 6.59526 7.44012
## beta1  0.00912 0.02099 0.02713 0.03311 0.04449
## sigma2 0.14979 0.21352 0.24667 0.28112 0.34813

plot(beta0,type="l")
abline(-2,0,col="red")

plot(beta1,type="l")
abline(2,0,col="red")

plot(sigma2,type="l")
abline(5,0,col="red")

plot(density(beta0))

plot(density(beta1))

plot(density(sigma2))

8 Modelo final

\(Y_i = 6,15 + 0,03\beta_1\)

Interpreta-se que um acréscimo na veiculação por jornal, aumeneta na média das vendas em \(e^{6,15 + 0,03}\)

9 Árvore de decisão

attach(base)
avr=tree(sales~.,base)
summary(avr)

## 
## Regression tree:
## tree(formula = sales ~ ., data = base)
## Variables actually used in tree construction:
## [1] "TV"    "radio"
## Number of terminal nodes:  8 
## Residual mean deviance:  2.174 = 417.4 / 192 
## Distribution of residuals:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -5.142000 -0.889100  0.007602  0.000000  1.036000  3.600000

plot(avr)
text(avr,pretty=0)

cv.avr=cv.tree(avr)
plot(cv.avr$size,cv.avr$dev,type='b')

prune.avr=prune.tree(avr,best=5)
plot(prune.avr)
text(prune.avr,pretty=0)