A análise estatística Bayesiana se beneficiou da explosão da computação barata e poderosa nas últimas duas décadas. As técnicas Bayesianas podem agora ser aplicadas a problemas de modelagem complexos onde não poderiam ter sido aplicadas anteriormente. Parece provável que a perspectiva Bayesiana continue a desafiar e talvez subplantar, métodos estatísticos frequentistas tradicionais que dominaram muitas disciplinas da ciência por tanto tempo.

O artigo atual revisará uma função que permite ao usuário conduzir regressão linear e modelagem linear generalizada, ou seja, não gaussiana; por exemplo, Poisson, binomial, etc.. Um modelo bastante simples é especificado, então modelado usando técnicas tradicionais e então modelado com uma abordagem Bayesiana.

Um tratamento completo da perspectiva Bayesiana está além do escopo deste artigo e poderia encher vários livros; e tem. Os leitores interessados podem consultar uma série de textos introdutórios com foco na perspectiva Bayesiana, por exemplo, Berry, 1996; Bolstad, 2004; Gelman, Carlin, Stern, & Rubin, 2004; Hoff, 2009 e outros. De maneira muito geral, a abordagem Bayesiana da inferência estatística difere da inferência freqüentista tradicional por assumir que os dados são fixos e os parâmetros do modelo são aleatórios, o que configura problemas na forma de; qual é a probabilidade de uma hipótese, ou parâmetro, dados os dados disponíveis? Esses tipos de problemas podem ser expressos com símbolos como: \(p(H|D)\). A inferência frequentista tradicional assume que os parâmetros do modelo são fixos, embora desconhecidos e os dados são essencialmente aleatórios; por exemplo, se a hipótese nula for verdadeira, qual é a probabilidade desses dados? Esses tipos de problemas podem ser declarados na forma geral; qual é a probabilidade dos dados dados uma hipótese? Em símbolos, isso se traduz em: \(p (D|H)\).

Os métodos Bayesianos se concentram em cinco elementos essenciais. Primeiro, a incorporação de informações prévias ou a priori, por exemplo, opinião de especialistas, uma revisão completa da literatura das mesmas variáveis ou semelhantes e/ou dados anteriores. A informação prévia é geralmente especificada quantitativamente na forma de uma distribuição, por exemplo, normal/Gaussiana, Poisson, binomial, etc. e representa uma distribuição de probabilidade para um coeficiente; ou seja, a distribuição de valores prováveis para um coeficiente que estamos tentando modelar, por exemplo, um peso \(\beta\).

Pode ajudar pensar no anterior como um melhor palpite educado. Em segundo lugar, o anterior é combinado com uma função de verossimilhança. A função de verossimilhança representa os dados, ou seja, qual é a distribuição da estimativa produzida pelos dados. Terceiro, a combinação da função a priori com a de verossimilhança resulta na criação de uma distribuição posterior dos valores dos coeficientes. Quarto, simulações são extraídas da distribuição posterior para criar uma distribuição empírica de valores prováveis para o parâmetro populacional. Quinto, estatísticas básicas são usadas para resumir a distribuição empírica de simulações a partir da posterior. A moda ou mediana ou média desta distribuição empírica representa a estimativa de máxima verossimilhança do valor populacional do coeficiente verdadeiro, ou seja, parâmetro populacional e intervalos confiáveis podem capturar o valor populacional verdadeiro com probabilidade anexada.

Tenha em mente que os antecedentes devem ser derivados de forma racional e honesta. Eles podem ser fracos ou fortes. Esses termos referem-se à força de crença que temos no(s) anterior(es). A priori fracas resultam quando não temos muitas evidências ou informações prévias nas quais basear a(s) a priori(s). Quando a priori é fraca, a distribuição a priori será ampla, refletindo muitos valores possíveis e a probabilidade será mais influente na criação da distribuição a posteriori. Prioridades fortes, por outro lado, resultam quando temos uma grande quantidade de evidências sobre as quais basear a(s) anterior(es). Quando a priori é forte, a distribuição a priori será estreita, refletindo um intervalo menor de valores possíveis e a probabilidade terá menos influência na criação da a posteriori, a priori fortes influenciarão a posteriori mais do que a verossimilhança. Deve ficar claro que uma característica-chave do anterior é a capacidade de quantificar nossa incerteza. O posterior pode ser pensado como um compromisso entre o anterior e a verossimilhança. Se a priori for fraca, então terá menos influência na criação da a posteriori; se a priori for forte, então será mais influente na criação da a posteriori.


Descrição do conjunto de dados


Conjunto de dados, com 445 observações nas 12 variáveis, descrito a seguir:

Evaluating the Econometric Evaluations of Training Programs with Experimental Data
Author(s): Robert J. LaLonde
Reviewed work(s):
Source: The American Economic Review, Vol. 76, No. 4 (Sep., 1986), pp. 604-620
Published by: American Economic Association
Stable URL: http://www.jstor.org/stable/1806062. 

dados = read.csv(file = "http://leg.ufpr.br/~lucambio/CE231/20221S/Lalonde.csv", sep = ",", header = TRUE)
head(dados)
##   X age educ black hisp married nodegr re74 re75     re78 u74 u75 treat
## 1 1  37   11     1    0       1      1    0    0  9930.05   1   1     1
## 2 2  22    9     0    1       0      1    0    0  3595.89   1   1     1
## 3 3  30   12     1    0       0      0    0    0 24909.50   1   1     1
## 4 4  27   11     1    0       0      1    0    0  7506.15   1   1     1
## 5 5  33    8     1    0       0      1    0    0   289.79   1   1     1
## 6 6  22    9     1    0       0      1    0    0  4056.49   1   1     1


O National Supported Work Demonstration (NSW) foi um programa de emprego temporário projetado para ajudar trabalhadores desfavorecidos sem habilidades básicas de trabalho a entrar no mercado de trabalho, dando-lhes experiência de trabalho e aconselhamento em um ambiente protegido.

Ao contrário de outros programas de emprego e treinamento patrocinados pelo governo federal americano, o programa de NSW designou candidatos qualificados para cargos de treinamento aleatoriamente. Aqueles designados para o grupo de tratamento receberam todos os benefícios do programa NSW, enquanto aqueles designados para o grupo de controle foram deixados à própria sorte.

Durante meados da década de 1970, a Manpower Demonstration Research Corporation (MDRC) operou o programa NSW em dez locais nos Estados Unidos. O MDRC admitiu no programa AFDC mulheres, ex-viciados em drogas, ex-delinquentes e desistentes do ensino médio de ambos os sexos.

Para aqueles designados para o grupo de tratamento, o programa garantiu um emprego por 9 a 18 meses, dependendo do grupo-alvo e do local. O grupo de tratamento foi dividido em equipes de três a cinco participantes que trabalharam juntos e se reuniram frequentemente com um conselheiro de NSW para discutir queixas e desempenho. O programa NSW pagou os membros do grupo de tratamento pelo seu trabalho. A tabela salarial oferecia aos estagiários salários mais baixos do que receberiam em um trabalho regular, mas permitiram que seus ganhos aumentassem para um desempenho e frequência satisfatórios. Os estagiários podiam permanecer em seus empregos apoiados até que seus termos no programa expirassem e fossem forçados a encontrar um emprego regular.

Embora essas diretrizes gerais fossem seguidas em cada local, as agências que operaram o experimento em nível local forneceram aos membros do grupo de tratamento diferentes experiências de trabalho. O tipo de trabalho variava mesmo dentro dos locais. Por exemplo, alguns dos estagiários em Hartford trabalhavam em um posto de gasolina, enquanto outros trabalhavam em uma gráfica. fazer compras.

Em particular, participantes masculinos e femininos frequentemente realizavam diferentes tipos de trabalho. Os participantes do sexo feminino geralmente trabalhavam em ocupações de serviços, enquanto os participantes do sexo masculino tendiam a trabalhar em ocupações de construção. Consequentemente, os custos do programa variaram entre os locais e grupos-alvo. O programa custou US$ 9.100 por participante do AFDC e aproximadamente US$ 6.800 para os estagiários de outros grupos-alvo

O MDRC coletou ganhos e dados demográficos dos membros do grupo de tratamento e controle na linha de base (quando o MDRC atribuiu aleatoriamente os participantes) e a cada nove meses depois disso, realizando até quatro entrevistas pós-linha de base. Muitos participantes não conseguiram completar essas entrevistas, e esse desgaste da amostra potencialmente influencia os resultados experimentais.

Felizmente, a maior fonte de atrito não afeta a integridade do desenho experimental. Em grande parte devido a recursos limitados, os administradores de NSW agendaram uma entrevista no 27º mês para apenas 65% dos participantes e uma entrevista no 36º mês para apenas 24% dos participantes não-AFDC. Nenhum dos participantes do AFDC foi agendado para uma entrevista de 36 meses, mas o ressurgimento do AFDC durante o outono de 1979 entrevistou 75% dessas mulheres entre 27 e 44 meses após a linha de base.

Como os membros do grupo de treinamento e controle foram agendados aleatoriamente para todas essas entrevistas, essa fonte de atrito não influenciou a avaliação experimental do programa NSW.

Naturalmente, os administradores do programa não localizaram todos os participantes agendados para essas entrevistas. A proporção de participantes que não conseguiram completar as entrevistas agendadas variou entre o grupo experimental, o tempo e o grupo-alvo. Enquanto as taxas de resposta foram estatisticamente significativamente maiores para o tratamento em oposição aos membros do grupo de controle, as diferenças nas taxas de resposta foram geralmente apenas alguns pontos percentuais. Para a entrevista do 27º mês, 72% dos tratamentos e 68% dos membros do grupo de controle completaram as entrevistas.

As diferenças nas taxas de resposta foram maiores ao longo do tempo e do grupo-alvo. Por exemplo, 79% dos participantes agendados completaram a entrevista do 9º mês, enquanto 70% completaram a entrevista do 27º mês. Os participantes do AFDC responderam a taxas consistentemente mais altas do que os outros grupos-alvo; 89 por cento dos participantes do AFDC completaram a entrevista do 9º mês em oposição a 76 por cento dos outros participantes. Embora essas taxas de resposta indiquem que os resultados experimentais podem ser tendenciosos, especialmente para os participantes não AFDC, as comparações entre as características basais dos participantes que fizeram e não completaram uma entrevista de 27 meses sugerem que qualquer viés existente pode ser pequeno


Variáveis observadas


Queremos estudar o ganho real em 1978 segundo os selecionados para o tratamento ou não e demais variáveis possivelmente explicativas. Um otro estudo seria identificar padrões entre os grupos tratamento e controle. 

library(lattice)
pairs(dados[dados$hisp == 1,c(2,3,10)], main = "Hispanos")

pairs(dados[dados$black == 1,c(2,3,10)], main = "Negros")

pairs(dados[I(dados$black == 0 & dados$hisp == 0) ,c(2,3,10)], main = "Outros")

ajuste1 = glm(treat ~ age + educ + black + hisp + married + nodegr + I(re78 - re74) + I(re78 - re75), 
              family = binomial(link = "logit"), data = dados)
summary(ajuste1)
## 
## Call:
## glm(formula = treat ~ age + educ + black + hisp + married + nodegr + 
##     I(re78 - re74) + I(re78 - re75), family = binomial(link = "logit"), 
##     data = dados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5187  -1.0076  -0.8565   1.2474   1.8677  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)     1.128e+00  1.060e+00   1.064  0.28754   
## age             3.246e-03  1.439e-02   0.226  0.82151   
## educ           -7.774e-02  7.209e-02  -1.078  0.28085   
## black          -1.951e-01  3.681e-01  -0.530  0.59611   
## hisp           -8.498e-01  5.077e-01  -1.674  0.09415 . 
## married         2.930e-01  2.693e-01   1.088  0.27653   
## nodegr         -8.663e-01  3.145e-01  -2.755  0.00588 **
## I(re78 - re74)  3.053e-05  2.637e-05   1.158  0.24687   
## I(re78 - re75) -9.328e-07  2.960e-05  -0.032  0.97486   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 604.20  on 444  degrees of freedom
## Residual deviance: 583.66  on 436  degrees of freedom
## AIC: 601.66
## 
## Number of Fisher Scoring iterations: 4


Para realizar o GLM Bayesiano, carregue o pacote arm que contém a função bayesglm (Gelman, et al., 2010). Você notará que existem várias dependências. 


library(arm)


Conduza o Modelo Linear Generalizado Bayesiano, aqui famíly = binomial e obtenha o resumo da saída. 


ajuste2 = bayesglm(formula = formula(ajuste1), family = family(ajuste1), data = dados)
summary(ajuste2)
## 
## Call:
## bayesglm(formula = formula(ajuste1), family = family(ajuste1), 
##     data = dados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4996  -1.0082  -0.8586   1.2476   1.8482  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)     1.007e+00  1.038e+00   0.970  0.33191   
## age             3.570e-03  1.426e-02   0.250  0.80232   
## educ           -7.178e-02  7.085e-02  -1.013  0.31104   
## black          -1.644e-01  3.550e-01  -0.463  0.64327   
## hisp           -7.859e-01  4.843e-01  -1.623  0.10469   
## married         2.830e-01  2.658e-01   1.065  0.28700   
## nodegr         -8.383e-01  3.081e-01  -2.721  0.00651 **
## I(re78 - re74)  2.861e-05  2.503e-05   1.143  0.25303   
## I(re78 - re75)  9.373e-07  2.815e-05   0.033  0.97344   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 604.20  on 444  degrees of freedom
## Residual deviance: 583.69  on 436  degrees of freedom
## AIC: 601.69
## 
## Number of Fisher Scoring iterations: 5


Podemos ver que a saída corresponde ao modelo linear tradicional (regressão OLS), bem como ao GLM tradicional. À medida que os tamanhos das amostras aumentam, os resultados devem convergir para os mesmos valores.

Um dos benefícios da perspectiva Bayesiana, para qualquer análise é que ela nos permite fazer declarações de intervalo confiáveis. Intervalos confiáveis são semelhantes aos intervalos de confiança, mas na estrutura Bayesiana, acredita-se que o intervalo realmente contém o verdadeiro parâmetro populacional.

Por exemplo: um intervalo de 95% de credibilidade para um parâmetro sendo estimado é interpretado como; há uma probabilidade de 95% de que o parâmetro real esteja incluído nesse intervalo. Isso porque o intervalo é baseado em informações da distribuição a posteriri de, por exemplo, uma das distribuições a posteriori do coeficiente do preditor, por exemplo, a distribuição a posteriori do coeficiente da variável educ.

A função bayesglm representa uma espécie de atalho da abordagem Bayesiana para inferência. Normalmente, a posteriori não é usada diretamente para fazer inferências. Em vez disso, uma distribuição empírica é construída com base em sorteios da a posteriori e essa distribuição empírica é o que informa a(s) inferência(s). Aqui, estamos usando o bayesglm como proxy para fazer a distribuição empírica adicionada. Com o bayesglm, obtemos uma distribuição de simulados que são usados no lugar de uma distribuição empírica real, que será abordada mais adiante.

Recupere as distribuições a posteriori dos coeficientes. A função head simplesmente lista as primeiras 10 linhas do objeto no qual ela é executada, o head padrão são as 6 primeiras. 


simulados = coef(sim(ajuste2))
head(simulados, 10)
##       (Intercept)          age        educ       black         hisp    married
##  [1,]  1.69926129 -0.006568690 -0.10041613 -0.09926589 -0.571008919  0.2044424
##  [2,]  1.26498036 -0.023023195 -0.07424979  0.08082208 -0.562269583  0.2904215
##  [3,] -0.04289681  0.001311425 -0.01570719  0.02915942 -0.606216253  0.2933196
##  [4,]  0.11035008 -0.014335466 -0.02120320  0.38901516  0.019498702  0.0133695
##  [5,]  0.95737878 -0.001715472 -0.04424906 -0.06782010 -0.615407369  0.6175201
##  [6,]  1.04533998 -0.011270881 -0.04047467 -0.20993521 -0.290127644 -0.0323446
##  [7,]  2.53490824  0.019542500 -0.24503434 -0.13581434 -0.267406317  0.2951474
##  [8,]  1.58768569 -0.009274054 -0.06621900  0.05586426 -0.354422235  0.1337110
##  [9,]  0.39402750  0.015641573 -0.01498250 -0.44847571 -1.109678118  0.6246406
## [10,] -0.21594114  0.033297239 -0.01900837 -0.51922992 -0.003246408  0.6150439
##           nodegr I(re78 - re74) I(re78 - re75)
##  [1,] -0.8069695   9.301838e-06   5.418707e-06
##  [2,] -0.7398508   4.646441e-05   1.149366e-05
##  [3,] -0.3991092   4.673020e-06   1.845875e-05
##  [4,] -0.5523988   1.475403e-05   2.348095e-05
##  [5,] -1.0760254   2.475967e-05   1.085017e-06
##  [6,] -0.7704278   4.621945e-05  -1.569760e-05
##  [7,] -1.1513882   4.725855e-05  -2.243892e-06
##  [8,] -1.4063448   9.038588e-06   1.134518e-05
##  [9,] -0.8187665   1.566632e-05   1.592685e-05
## [10,] -0.5417282   3.456604e-05  -2.637296e-05


Extraia apenas a distribuição a posteriori do coeficiente da variável educ. Novamente, a função head simplesmente lista os 10 primeiros itens do objeto. 


posteriori.educ = simulados[, 3]
head(posteriori.educ, 10)
##  [1] -0.10041613 -0.07424979 -0.01570719 -0.02120320 -0.04424906 -0.04047467
##  [7] -0.24503434 -0.06621900 -0.01498250 -0.01900837


Damos uma olhada na distribuição a posteriori do coeficiente da variável educ. 


densityplot(posteriori.educ, main = "Densidade estimada da variável educ", xlab = "Observações", 
            ylab = "Densidade", pch = 19)


Agora podemos recuperar o intervalo de 95% de credibilidade para o coeficiente da variável aberta. 


quantile(posteriori.educ, c(0.025, 0.975))
##        2.5%       97.5% 
## -0.20660878  0.05186972


Lembre-se, esse intervalo confiável é interpretado como; há uma probabilidade de 95% de que o verdadeiro valor populacional do coeficiente educ esteja entre -0.22 e 0.05. Lembre-se de que esses números flutuarão um pouco com base na natureza iterativa da função.

Para fazer inferências verdadeiramente Bayesianas sobre nossos coeficientes, precisamos fazer a etapa extra de criar a(s) distribuição(ões) empírica(s) mencionada(s) acima. Ir além implica realmente criar uma distribuição empírica baseada em métodos iterativos a partir da a posteriori. A função MCMClogit no pacote MCMCpack (Martin, Quinn, & Park, 2010) nos fornece o método de simulação Markov Chain Monte Carlo para criar a distribuição empírica; o que por si só nos permite calcular as estatísticas descritivas usadas para inferência. Ou seja, a moda, mediana ou média da distribuição empírica de simulações de MCMC é a estimativa de verossimilhança máxima, ou seja, o topo de uma função de densidade do parâmetro populacional. A função MCMClogit também nos dá o intervalo confiável que inclui o valor real do parâmetro da população.

Primeiro, carregue a biblioteca MCMCpack. 

library(MCMCpack)
ajuste3 = MCMClogit(formula = formula(ajuste2), burnin = 1000, mcmc = 10000, data = dados)


par(mar=c(3,2,1,1) + 0.1)
for(i in 1:9) {
  plot(ajuste3[,i], main = colnames(ajuste3)[i]); grid()
  }


summary(ajuste3)
## 
## Iterations = 1001:11000
## 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
## (Intercept)     1.159e+00 1.040e+00 1.040e-02      5.695e-02
## age             3.406e-03 1.379e-02 1.379e-04      7.493e-04
## educ           -7.795e-02 6.959e-02 6.959e-04      3.796e-03
## black          -2.248e-01 3.884e-01 3.884e-03      2.244e-02
## hisp           -9.309e-01 5.544e-01 5.544e-03      3.431e-02
## married         3.135e-01 2.709e-01 2.709e-03      1.499e-02
## nodegr         -8.730e-01 3.021e-01 3.021e-03      1.593e-02
## I(re78 - re74)  3.549e-05 2.671e-05 2.671e-07      1.502e-06
## I(re78 - re75) -5.355e-06 2.988e-05 2.988e-07      1.673e-06
## 
## 2. Quantiles for each variable:
## 
##                      2.5%        25%        50%        75%      97.5%
## (Intercept)    -9.022e-01  4.389e-01  1.220e+00  1.862e+00  3.080e+00
## age            -2.385e-02 -4.911e-03  3.104e-03  1.275e-02  3.081e-02
## educ           -2.146e-01 -1.236e-01 -8.320e-02 -3.063e-02  5.823e-02
## black          -1.006e+00 -4.625e-01 -2.325e-01  2.969e-02  5.483e-01
## hisp           -2.111e+00 -1.297e+00 -9.324e-01 -5.515e-01  7.985e-02
## married        -2.007e-01  1.200e-01  3.265e-01  5.034e-01  8.400e-01
## nodegr         -1.474e+00 -1.061e+00 -8.856e-01 -6.639e-01 -3.104e-01
## I(re78 - re74) -1.394e-05  1.744e-05  3.515e-05  5.166e-05  9.218e-05
## I(re78 - re75) -6.581e-05 -2.455e-05 -5.820e-06  1.449e-05  5.203e-05


Em seguida, aplique a função MCMClogit. Observe que a fórmula do modelo é a mesma, mas aqui temos algumas novas opções. O argumento burnin é usado porque as iterações do MCMC são sensíveis aos seus valores iniciais, portanto, as primeiras, ou seja, 1.000 iterações são descartadas. O mcmc simplesmente emite quantas iterações (pós-combustão) serão usadas para construir a distribuição empírica. O padrão thin é 1 e representa um controle na convergência, de modo que, uma vez alcançada a convergência aproximada, pode ser benéfico manter apenas algumas simulações e descartar o restante para conservar os recursos do computador (Gelman, Carlin, Stern, & Rubin, 2004).

Percebemos que o única variável significatia é nodegr, significa que os indivíduos selecionados no grupo de tratamento menos deles tinham diploma de ensino médio. tecnicamente isso seria um vício na amostragem e, portanto, não seria recomendado a utilização dessa variável para entender o impacto do trenamento nessas pessoas quanto ao ganho no ano 1978.