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,03x_i\)

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

9 Medida de erro para o modelo Gamma

train = sample(1:nrow(base), nrow(base)/1.25)

yhat = 6.15 + 0.03*base$newspaper[-train] 
avrt = base[-train,"sales"]

mean((yhat-avrt)^2)

## [1] 80.99163

10 Árvore de decisão

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

## 
## Regression tree:
## tree(formula = sales ~ ., data = base, subset = train)
## Variables actually used in tree construction:
## [1] "TV"    "radio"
## Number of terminal nodes:  9 
## Residual mean deviance:  1.737 = 262.3 / 151 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.08900 -0.90190 -0.01872  0.00000  0.88860  3.66500

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=3)
plot(prune.avr)
text(prune.avr,pretty=0)

11 Medições

yhat=predict(avr,newdata=base[-train,])
  
avrt = base[-train,"sales"]

plot(yhat,avrt)
abline(0,1)

mean((yhat-avrt)^2)

## [1] 2.537941

12 Random Forest

set.seed(02042019)
train = sample(1: nrow(base ), nrow(base )/2)
test=base[-train ,"sales"]  

bag=randomForest(sales~.,data=base ,subset=train, importance =TRUE) 
bag

## 
## Call:
##  randomForest(formula = sales ~ ., data = base, importance = TRUE,      subset = train) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##           Mean of squared residuals: 3.697536
##                     % Var explained: 84.55

yhat.bag = predict(bag,newdata =base[-train ,])
plot(yhat.bag , test)
abline(0,1)

mean(( yhat.bag - test)^2)

## [1] 4.426138

13 Boosting

set.seed(02042019)
boost=gbm(sales~.,data=base[train ,], distribution="gaussian",n.trees =5000 , interaction.depth =4)
summary(boost)

##                 var   rel.inf
## TV               TV 51.003209
## radio         radio 43.415780
## newspaper newspaper  5.581011

plot(boost)

yhat.boost=predict(boost,newdata = base[-train ,],n.trees =5000)
mean((yhat.boost-test)^2)

## [1] 2.151508

14 Bagging manual

#Criando uma lista que contenha todos os valores para treino
train = list()
for (i in 1:10000){
train[[i]]       = sample(1:nrow(base), nrow(base)/2,replace = T)
}

#Criando uma lista que contenha todos os ajustes da árvore baseado no treino
avr = list()
for (i in 1:10000){
avr[[i]]=tree(sales~.,base,subset = train[[i]])
}

#Média para cada árvore
yhat = list()
for (i in 1:10000){
yhat[[i]]=predict(avr[[i]],newdata=base[-train[[i]],])  
}

#Observado para cada árvore
avrt = list()
for (i in 1:10000){
avrt[[i]] = base[-train[[i]],"sales"]
}

#Visualização das 9 iniciais
par(mfrow=c(3,3))
for (i in 1:9){
plot(yhat[[i]],avrt[[i]])
abline(0,1)}

#Medição individual
rss = vector()
for (i in 1:10000){
rss[i] = mean((yhat[[i]]-avrt[[i]])^2)
}

#Medição final 
mean(rss)

## [1] 3.964381

A idéia foi criar um bootstrap manual, que fosse capaz de ralizar uma amostragem com reposição. Logo após, ajustar uma árvore de regressão para cada amostra sorteada. Cada árvore irá gerar uma medida e por fim, a média das medidas será considerada a medida final.

Tem-se que esta medida é equivalente a 3,95. Menor que a obtida pela árvore de regressão, no entando maior que a alcançada pelo método boosting.