Uso do Shiny para o ensino-aprendizagem de Estatística Bayesiana

Cristian Villegas
Eduardo E. Ribeiro Jr
Roseli A. Leandro

Autores

  • Cristian Villegas: Doutor em Estatística pela USP (2010). Atualmente é professor doutor na ESALQ-USP e tem experiência na área de modelos simétricos e semi-paramétricos.
    clobos@usp.br

  • Eduardo E. Ribeiro Jr: Bacharel em Estatística pela UFPR (2016). Atualmente é mestrando no PPG em Estatística e Experimentação na ESALQ-USP e tem experiência na área de modelos lineares generalizados e extensões.
    jreduardo@usp.br

  • Roseli A. Leandro: Doutora em Estatística e Experimentação pela ESALQ-USP (1997). Atualmente é professora doutora na ESALQ-USP e tem experiência na área de inferência bayesiana.
    raleandr@usp.br

Disponibilização

Motivação

  • Aumento da utilização da inferência bayesiana para análise de dados por meio de modelos mais fidedignos a realidade;

  • O melhor uso dos métodos disponíveis para inferência bayesiana dependem de um bom entendimento básico da área;

  • Conceitos importantes de inferência bayesiana podem ser melhor ilustrados por meio de recursos interativos.

Agenda

  • Parte I: Introdução ao Shiny;

  • Parte II: Conceitos de inferência bayesiana com Shiny;

  • Parte III: Usando Shiny para diagnóstico de cadeias MCMC.

Parte I

Introdução ao Shiny

O que é o shiny?

http://shiny.rstudio.com/

Pra que serve o shiny?

Criar aplicativos web com códigos R!

Exemplos:

Como construir um Shiny app?

  • Deve-se programar:
    • A interface para inputs do usuário;
    • As instruções que acessam os inputs e geram outputs para interface.

Regras para criação de apps com shiny

  • Crie os outputs

Regras para criação de apps com shiny

  • Renderize os outputs

Regras para criação de apps com shiny

  • Acesse os inputs

Para usuários do RStudio

Para usuários do RStudio

Dicas gerais para programação de apps

  1. Programe com um editor adequado;
    • Indentação;
    • Destaque de delimitadores.
  2. Inspecione o erro;
    • As mensagens de erro são muito informativas;
    • Muitas vezes o erro é devido a ausência de vírgula ou parentêses.
  3. Um passo de cada vez;
    • Comece com exemplos simples;
    • Aumente gradativamente a complexidade.
  4. Consulte a documentação!

Parte II

Conceitos de inferência bayesiana com Shiny

Inferência bayesiana: conceitos básicos

Referências originais

  • LII An essay towards solving a problem in the doctrine of chances. By the late Rev. Mr. Bayes, F. R. S. communicated by Mr. Price, in a letter to John Canton, A. M. F. R. S late Rev. Mr. Bayes, F. R. S.
  • Bayes, Thomas; Price, Mr. (1763). "An Essay towards solving a Problem in the Doctrine of Chances". Philosophical Transactions of the Royal Society of London. 53 (0): 370–418. doi:10.1098/rstl.1763.0053.
  • Barnard, G (1958). "Studies in the History of Probability and Statistics: IX. Thomas Bayes's An Essay Towards Solving a Problem in the Doctrine of Chances". Biometrika. 45 (3–4): 296–315. doi:10.1093/biomet/45.3-4.293.
  • Thomas Bayes "An Essay towards solving a Problem in the Doctrine of Chances". (Bayes' essay in the original notation)

Inferência bayesiana: conceitos básicos

Motivações

  • Richard Price e a existência de Deus;
  • Richard Price descobriu o experimento de Thomas Bayes e seu, agora famoso, teorema em artigos depois de sua morte.
  • Ele acreditava que o teorema de Bayes ajudaria a provar a existência de Deus.

Primeira publicação em periódico científico

  • In 1662: The first volumes of what was the world's first scientific journal were very different from today's journal


Wikipedia: An Essay towards solving a Problem in the Doctrine of Chances Philosophical Transactions − the world's first science journal

O resultado central do trabalho de Bayes

Teorema de Bayes (teorema da probabilidade inversa)

Se \(H\) denota uma hipótese e \(D\) denota os dados o teorema estabelece que:

\[ P(H \mid D) = \dfrac{P(D \mid H) \times P(H)}{ P(D)}, \]

  • com \(P(H)\) representando o conhecimento sobre \(H\) antes da observação dos dados; e
  • \(P(H \mid D)\) representando o conhecimento sobre \(H\) atualizado pelos dados.

Inferência bayesiana

  • O teorema de Bayes mostra como probabilidades mudam à luz de evidências: dados.

  • Tendo especificado \(P(D \mid H)\) e \(P(H)\) o teorema fornece uma solução do problema de como aprender através dos dados.

Inferência bayesiana

Teorema de Bayes, no contexto de inferência estatística

\[ \pi(\boldsymbol{\theta} \mid \boldsymbol{y}) = \dfrac{\mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y}) \pi(\boldsymbol{\theta})}{ \int \mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y}) \pi(\boldsymbol{\theta})d\boldsymbol{\theta}} \propto \mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y}) \pi(\boldsymbol{\theta}), \]

em que,

  • \(\pi(\boldsymbol{\theta})\) é chamada de distribuição a priori;
  • \(\mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y})\) é a função de verossimilhança; e
  • \(\pi(\boldsymbol{\theta} \mid \boldsymbol{y})\) é chamada de distribuição a posteriori.

Distribuição a priori

  • A distribuição a priori, \(\pi(\boldsymbol{\theta})\) que representa o que é conhecido a respeito dos parâmetros desconhecidos antes que os dados sejam observados, desempenha um importante papel na inferência bayesiana.

  • \(\pi(\boldsymbol{\theta})\) pode ser utilizada para representar conhecimento prévio ou ignorância relativa.

  • Classificações e interpretação (diferentes autores com diferentes definições):
    • Objetivas, subjetivas;
    • Informativas, não-informativas;
    • Vagas;
    • De referência.

Inferência bayesiana

  • Do ponto de vista bayesiano, toda inferência a respeito do parâmetro está baseada na distribuição a posteriori, ou seja,

  • Dada a verossimilhança \(\mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y})\) e a distribuição a priori \(\pi(\boldsymbol{\theta})\), o início de qualquer inferência paramétrica ou decisão sobre \(\boldsymbol{\theta}\) está na densidade a posteriori:

Exemplo com a distribuição binomial

Sejam \(\boldsymbol{y} = \{y_{1}, y_{2}, \ldots, y_{n}\}\) uma amostra aleatória de uma distribuição binomial com parâmetro \(\theta\),

  • a função massa de probabilidade é dada por

\[ f(y_i \mid \theta) = {m \choose y_i} \theta^{y_i} (1-\theta)^{m-y_i}, \quad y_i \in \{ 0, 1, 2, \ldots, m\}; \]

  • e a função de verossimilhança por

\[ \begin{aligned} \mathcal{L}(\theta \mid \boldsymbol{y}) &= \prod_{i=1}^{n} {m \choose y_i} \theta^{y_i} (1-\theta)^{n-y_i} \\ &= \theta^{\sum_i^n y_i} (1-\theta)^{n-\sum_i^n y_i} \prod_{i=1}^n {m \choose y_i}, \quad \theta \in [0,1]. \end{aligned} \]

Exemplo com a distribuição binomial

  • Considerando-se a distribuição a priori conjugada para \(\theta\),

\[ \pi(\theta) = \dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1}, \alpha > 0, \beta > 0. \]

  • A distribuição a posteriori é obtida analiticamente,

\[ \begin{gathered} \pi(\theta \mid \boldsymbol{y}) \propto \mathcal{L}(\theta \mid \boldsymbol{y}) \pi(\theta) = \theta^{\sum_i^{n} y_i + \alpha - 1} (1 - \theta)^{nm-\sum_i^n y_i + \beta -1} \\ \theta \mid \boldsymbol{y} \sim \text{Beta} \left( \sum_i^{n} y_i + \alpha,\; nm - \sum_{i=1}^n y_i + \beta \right) \end{gathered} \]

Exemplo com a distribuição binomial

  • A distribuição preditiva

\[ \begin{aligned} f(\tilde{y} \mid \boldsymbol{y}) &= \int_{\mathcal{R}} f(\tilde{y} \mid \theta) \pi(\theta \mid \boldsymbol{y}) d\boldsymbol{\theta} =\\ &= \int_0^1 {\tilde{m} \choose \tilde{y}} \theta^{\tilde{y}} (1-\theta)^{ \tilde{m}-\tilde{y} } \theta^{\sum_i^{n} y_i + \alpha -1}(1 - \theta)^{nm-\sum_i^n y_i + \beta -1} d \theta,\\ &\tilde{y} \in \{ 0, 1, 2, \ldots, \tilde{m}\}. \end{aligned} \]

e portanto, \[ \tilde{y} \mid \boldsymbol{y} \sim \text{Beta-Binomial}( \tilde{m},\; \sum_{i=1}^n y_i + \alpha,\; n\tilde{m}-\sum_i^n y_i + \beta). \]

Exemplo com a distribuição Poisson

Sejam \(\boldsymbol{y}=\{y_{1}, y_{2}, \ldots, y_{n}\}\) uma amostra aleatória de uma distribuição de Poisson com parâmetro \(\theta\),

  • a função massa de probabilidade é dada por

\[ f(y_{i} \mid \theta) = \dfrac{e^{-\theta} \theta^{y_i}}{y_i!}, \quad y_{i}= 0, 1, 2, \ldots; \]

  • e a função de verossimilhança por

\[ \mathcal{L}(\theta \mid \boldsymbol{y}) = \prod_{i=1}^{n} \dfrac{e^{-\theta} \theta^{y_i}}{y_i!} = \dfrac{e^{-n \theta} \theta^{\sum_{i=1}^{n} y_i}}{ \prod_{i=1}^{n}y_{i}!},\quad \theta >0. \]

Exemplo com a distribuição Poisson

  • Considerando-se a distribuição a priori conjugada para \(\theta\),

\[ \pi(\theta) = \dfrac{\beta^\alpha}{ \Gamma(\alpha)}\theta^{\alpha-1} e^{-\beta \theta}. \]

  • A distribuição a posteriori é obtida analiticamente,

\[ \begin{aligned} \pi(\theta \mid \boldsymbol{y}) &= \dfrac{ \mathcal{L}(\theta \mid \boldsymbol{y}) \pi(\theta)}{ \int_{0}^{\infty}\pi(\theta) \mathcal{L}(\theta \mid \boldsymbol{y}) \pi(\theta) d \theta },\\ &= \dfrac{(n+\beta)^{\alpha + n \bar{y} } \theta^{\alpha+n \bar{y}-1} e^{-(n+\beta) \theta}}{\Gamma(\alpha + n \bar{y} )},\\ \theta \mid \boldsymbol{y}&\; \sim \text{Gamma}(\alpha+n \bar{y},\; n + \beta) \end{aligned} \]

Exemplo com a distribuição Poisson

Observe que, desprezando-se a constante normalizadora, o kernel da distribuição a posteriori é dado por:

\[ \theta^{\alpha+n \bar{y}-1} ~ e^{-(n+\beta), \theta} \]

que representa o kernel de uma distribuição gama com parâmetro de forma, \(\alpha+n \bar{y}\) e parâmetro de taxa , \(n + \beta\).

Ou seja, desprezando-se a constante normalizadora, chega-se a mesma conclusão.

Exemplo com a distribuição Poisson

  • A distribuição preditiva, para \(\tilde{y} = 0, 1, \ldots\)

\[ \begin{aligned} f(\tilde{y} \mid \boldsymbol{y}) &= \int_{\mathcal{R}} f(\tilde{y} \mid \theta) \pi(\theta \mid \boldsymbol{y}) d\boldsymbol{\theta} =\\ &= \int_{0}^{\infty} \dfrac{e^{-\theta}\theta^{\tilde{y}}}{\tilde{y}!} \dfrac{\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\tilde{\alpha})} \theta^{\tilde{\alpha}-1} e^{-\tilde{\beta} \theta} d\theta\\ &= \dfrac{\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\tilde{\alpha}) \tilde{y}!} \dfrac{\Gamma(\tilde{\alpha}+\tilde{y}) \Gamma(\tilde{\beta}+1)}{\Gamma(\tilde{\alpha}+\tilde{y} + \tilde{\beta}+1)} \end{aligned} \]

considerando-se \(\tilde{\alpha}=\alpha+n \bar{y}\) e \(\tilde{\beta} = n+\beta\).

Algoritmo Metropolis-Hastings

  • Quando a distribuição a posteriori possuir uma expressão complexa, a obtenção de resumos a posteriori pode não ser simples;

  • Uma forma bastante utilizada na inferência bayesiana para superar essa dificuldade é gerar uma cadeia de Markov cuja distribuição de equilíbrio seja a distribuição de interesse, e assim os resumos a posteriori serão calculados por meio da amostra obtida (sampling based)

  • Um algoritmo bastante utilizado é o algoritmo de Metropolis-Hastings.

Algoritmo Metropolis-Hastings (versão 1)

Seja \(\pi(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:

  1. Propõe um \(\theta^{*} \sim g(\theta \mid \theta^{t})\);
  2. Aceita \(\theta^{t+1} = \theta^{*}\) com probabilidade \(\min\{1,r\}\) em que:

\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{\pi(\theta^{*}\mid y) \;/\; g(\theta^{*} \mid \theta^{t})}{ \pi(\theta^{t} \mid y) \;/\; g(\theta^{t} \mid \theta^{*}) } \]

  1. Caso contrário, faz \(\theta^{t+1} = \theta^{t}\)

Algoritmo Metropolis-Hastings (versão 2)

Seja \(\pi(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:

  1. Propõe um \(\theta^{*} \sim g(\theta)\);
  2. Aceita \(\theta^{t+1} = \theta^{*}\) com probabilidade \(\min\{1,r\}\) em que:

\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{\pi(\theta^{*}\mid y) \;/\; g(\theta^{*})}{ \pi(\theta^{t} \mid y) \;/\; g(\theta^{t}) } \]

  1. Caso contrário, faz \(\theta^{t+1} = \theta^{t}\)

Algoritmo Metropolis-Hastings (versão 3)

Seja \(\pi(\theta \mid y) = k q(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:

  1. Propõe um \(\theta^{*} \sim g(\theta \mid \theta^{t})\);
  2. Aceita \(\theta^{t+1} = \theta^{*}\) com probabilidade \(\min\{1,r\}\) em que:

\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{q(\theta^{*}\mid y) \;/\; g(\theta^{*}\mid \theta^{t})}{ q(\theta^{t} \mid y) \;/\; g(\theta^{t}\mid \theta^{*}) } \]

  1. Caso contrário, faz \(\theta^{t+1} = \theta^{t}\)

Algoritmo Metropolis-Hastings (versão 4)

Seja \(\pi(\theta \mid y) = k q(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:

  1. Propõe um \(\theta^{*} \sim g(\theta)\);
  2. Aceita \(\theta^{t+1} = \theta^{*}\) com probabilidade \(\min\{1,r\}\) em que:

\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{q(\theta^{*}\mid y) \;/\; g(\theta^{*})}{ q(\theta^{t} \mid y) \;/\; g(\theta^{t}) } \]

  1. Caso contrário, faz \(\theta^{t+1} = \theta^{t}\)

Parte III

Usando Shiny para diagnóstico de cadeias MCMC

Inferência bayesiana via MCMC (Softwares)

  • Pacotes que são interfaces para programas standalone:
    • rstan (para o STAN);
    • rjags/R2jags (para o JAGS);
    • R2OpenBUGS/rbugs (para o OpenBUGS);
    • R2WinBUGS (para o WinBUGS).
  • Pacotes R
    • mcmc;
    • MCMCpack;
    • LaplacesDemon.

Pacote mcmctools

Ferramentas utéis para manipulação de objetos MCMC

# install.packages("devtools")
install_github("jreduardo/mcmctools")
help(package = "mcmctools", help_type = "html")

Funções:

  • mcmctools::mct_select
  • mcmctools::mct_burnin
  • mcmctools::mct_thin
  • mcmctools::mct_merge
  • mcmctools::mct_shinyapp

Exemplo de uso

#-------------------------------------------
# (toy) Data set
str(mtcars)
help(mtcars, help_type = "html")
splom(~mtcars, data = mtcars)

#-------------------------------------------
# Define linear model in JAGS
str_model <- "
 model {
   for (i in 1:n) {
     mu[i] <- X[i, ] %*% beta[]
     y[i] ~ dnorm(mu[i], tau)
   }
   for (j in 1:p) {
     beta[j] ~ dnorm(0, 0.001)
   }
   tau ~ dgamma(0.001, 0.001)
 }
"

#-------------------------------------------
# Initialize the JAGS model
X <- model.matrix(
    ~factor(cyl) + factor(am) + hp + drat + qsec,
    data = mtcars)
jags <- jags.model(
    file = textConnection(str_model),
    data = list(
        n = nrow(X),
        p = ncol(X),
        X = X,
        y = mtcars$mpg),
    n.chains = 3)

#-------------------------------------------
# Posterior sampling
samples <- coda.samples(
    model = jags,
    variable.names = c("beta", "tau", "mu"),
    n.iter = 5000)
str(samples)

#-------------------------------------------
mct_shinyapp(samples)

Obrigado

http://www.leg.ufpr.br/~eduardojr/courses/seeb
{clobos | jreduardo | rleandr} @ usp.br