2.1 Modelo Linear Generalizado

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)\)

2.2 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.

2.3 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.

2.4 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.

2.5 Advertising

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

2.6 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)

mean = param_mean
var  = param_cov_mat

#Gerando amostras da distribuição a posteriori
samples <- rmvnorm(10000, param_mean, param_cov_mat)
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

2.7 Distribuições a posteriori

par(mfrow=c(1,3))

x <- seq(4, 8.2, length=1000)
y <- dnorm(x, mean=mean[1], sd=sqrt(var[1,1]))
plot(x, y, type="l", lwd=2,main=expression(beta[0]))
abline(v=mean[1],col="red",lwd=2,lty=2)

x <- seq(-0.01, 0.065, length=1000)
y <- dnorm(x, mean=mean[2], sd=sqrt(var[2,2]))
plot(x, y, type="l", lwd=2,main=expression(beta[1]))
abline(v=mean[2],col="red",lwd=2,lty=2)

x <- seq(0, 0.5, length=1000)
y <- dnorm(x, mean=mean[3], sd=sqrt(var[3,3]))
plot(x, y, type="l", lwd=2,main=expression(tau))
abline(v=mean[3],col="red",lwd=2,lty=2)

2.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}\)

2.9 Medida de erro para o modelo Gamma

set.seed(2019)

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

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


mean((yhat-avrt)^2)

## [1] 71.28157

Ressalta-se, fundamental, que o estimador de bayes não é formulado para otimizar a soma do quadrado dos resíduos. Junta-se com o fato do modelo estar apenas trabalhando com uma única variável e o resultado obtido aqui será alto.

3.1 Treino

Tamanho amostral de 120 para o treino

attach(base)
set.seed(2019)
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:  8 
## Residual mean deviance:  2.16 = 198.7 / 92 
## Distribution of residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -4.0140 -0.8115 -0.1137  0.0000  0.8940  4.7860

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

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

3.2 Poda

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

Verifica-se que a melhor predição está associada ao grupo onde os valores de televisão está acima de 156,55 e de rádio está acima de 25,45

3.3 Medições

set.seed(2019)

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

mean((yhat-avrt)^2)

## [1] 3.953593