Trabalho 1 - Inteligência Computacional

A base de dados

Importando

rm(list=ls())
base <- read.csv("C:/Users/raphael.fernandez/Downloads/Advertising.csv")
base$X = NULL
kable(head(base)) %>%
  kable() %>%
  kable_styling()

Error in as.data.frame.default(x) : 
  cannot coerce class ‘"knitr_kable"’ to a data.frame

Avaliando as variáveis

Como não há informações a priori no quis respeito a natureza das variáveis, faz-se necessária a investigação de algumas características das mesmas.

str(base)

'data.frame':   200 obs. of  4 variables:
 $ TV       : num  230.1 44.5 17.2 151.5 180.8 ...
 $ radio    : num  37.8 39.3 45.9 41.3 10.8 48.9 32.8 19.6 2.1 2.6 ...
 $ newspaper: num  69.2 45.1 69.3 58.5 58.4 75 23.5 11.6 1 21.2 ...
 $ sales    : num  22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...

Verifica-se que todas as variáveis são contínuas e não há evidências empíricas que apontem a presença de observações não-positivas

Visualizacão gráfica

par(mfrow=c(2,2))
plot(density(base$TV),main="Densidade da TV")
plot(density(base$radio),main="Densidade da radio")

plot(density(base$newspaper),main="Densidade do Jornal")
plot(density(base$sales),main="Densidade das Vendas")

A densidade da Televisão e do Rádio apontam para uma mistura de distribuições, mas como é considerado anti ético a observação do conjunto dos dados antes da modelagem a priori, iremos apenas nos ater a variável resposta, que demonstra comportamento assimétrico, com decaimento rápido, comportamento comumente observado em um distribuições da família Gama.

Não-normalidade

jarque.bera.test(base$sales)

    Jarque Bera Test

data:  base$sales
X-squared = 6.9848, df = 2, p-value = 0.03043

hist(base$sales,main="Histograma das Vendas")

qqnorm(base$sales)
qqline(base$sales)

Embora seja possível afrouxar a condição de normalidade, em termos visuais, o teste Jarque Bera elimina qualquer hipótese da mesma, ao nível de 0,05.

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

Verossimilhança

Posto que a distribuição 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ção 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.

Prioris hierárquicas

Em inferência bayesiana, quantidades desconhecidas são atribuídas com distribuição, 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.

Markov Chain Monte Carlo

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 da distribuição a posteriori, desconhecida

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

\(\pi(x)\) é um vetor paramétrico de distribuições a 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 Markov Chain Monte Carlo oferece uma convergência garantida via simulaçao por processos estocásticos.

Simulação da Distribuição a Posteriori

Utilizaremos o pacote R2WinBugs para efeito de c?lculo

y=base$sales
x1=base$newspaper
#Modelo do tipo BUGS
sink("mod1.txt")        
cat("
    MODEL LR1 {
    #Verossimilhan?a 
    for(i in 1:N) {
    y[i] ~ dgamma(mu[i],tau)
    mu[i] <- beta0 + beta1*x1[i]
    }
    #Proris
    beta0 ~ dgamma(1,0.1)
    beta1 ~ dgamma(1,0.1)
    tau ~ dgamma(1,0.1)
    
    #Transforma??o para vari?ncia
    sigma2 <- 1/tau
    } ", fill = TRUE)
sink()
#Obeservação: A dist normal no bugs está definida para a precisão. 
#Tamanho amostral
N = length(y)
#Conjunto de dados observados
data = list("N","y","x1")
#Parâmetros estudados
params = c("beta0", "beta1","tau","sigma2")  
#Valores iniciais
inits <- function () {list(beta0 = rnorm(1),
                           beta1 = rnorm(1),
                           tau = 1)
                           sigma2 = 1}
#Modelo final
result <- bugs(data = data, inits = inits, parameters.to.save = params,
               model.file = "mod1.txt", n.chains = 3, n.iter = 10000,
               n.burnin = 1000, n.thin=9, bugs.directory = "C:/WinBUGS14",
               debug = TRUE, save.history = TRUE, DIC = FALSE)

Foram excluídas da análise as covariáveis onde acreditase, antecipadamente, que estas sejam composições, isto é, soma de outras variáveis. São estas: Radio e TV. Permaneceu no modelo apenas o Jornal como priori, por apresentar um comportamento cujo um modelo mais simples, como o aqui apresentado, é mais adequado.

Todas as atribuições aqui definidas foram retiradas de Gelman et al. 2004

Modelo final

Sumário da distribuição a posteriori

result$summary %>%
  kable() %>%
  kable_styling()

	mean	sd	2.5%	25%	50%	75%	97.5%	Rhat	n.eff
beta0	6.2581520	0.6621779	5.015625	5.7970000	6.243000	6.683000	7.6200500	1.004490	500
beta1	0.0278669	0.0089434	0.010609	0.0219275	0.027775	0.033960	0.0455605	1.000811	3000
tau	0.5070021	0.0507519	0.412195	0.4731750	0.505950	0.539025	0.6136000	1.004720	510
sigma2	1.9924600	0.2029242	1.630000	1.8550000	1.976500	2.113250	2.4260250	1.004718	510

Convergência das cadeias

bugs = result$sims.list
bugs_b0 = bugs$beta0
bugs_b1 = bugs$beta1
bugs_phi = bugs$tau
plot(bugs_b0,type = "l",main=expression(beta[0]),ylab="",
     xlab="",cex.main=2)

plot(bugs_b1,type = "l",main=expression(beta[1]),ylab="",
     xlab="",cex.main=2)

plot(bugs_phi,type = "l",main=expression(phi),ylab="",
     xlab="",cex.main=2)

Verifca-se uma boa convergência, embora usualmente seja feito critérios, aqui não o será. Apenas critérios visuais.

Densidade das distribuições a posteriori

densplot(as.mcmc(bugs_b0),main=expression(beta[0]),xlab="",ylab="",cex.main=2)

densplot(as.mcmc(bugs_b1),main=expression(beta[1]),xlab="",ylab="",cex.main=2)

densplot(as.mcmc(bugs_phi),main=expression(phi),xlab="",ylab="",cex.main=2)

O comportamento “normal” das posterioris simuladas é esperado. Com as prioris impostas, o resultado analático esperado é de uma distribuição Gama para os parâmetros, maginalmente.

Equação do modelo

Na abordagem bayesiana, considera-se alguma função perda para atribuir a estimação de parâmetros. Aqui, foi utilizado a média, isto é, a perda quadrática.

\(Yi = 6,25 + 0,03Xi\)

Interpretação do modelo

Como a função de ligação é a exponencial, um acréscimo de 1 em TV, significa um acréscimo de \(e^{6,25 + 0,03}\) nas vendas