Uncanny X-Men: opinião bayesiana na análise do Dr. Silge

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


Leia Também  Lei de Benford: aplicando-se a dados existentes

O caminho bayesiano

A modelagem bayesiana é a prática de atualizar nossas crenças anteriores usando
dados observados para produzir uma distribuição de probabilidade para os valores de
parâmetros desconhecidos. Assim, diferentemente das estimativas pontuais únicas fornecidas por
Abordagens “freqüentistas”, os resultados de uma análise bayesiana são os
distribuições de parâmetros estimados. É por isso que o Dr. Silge
análise de bootstrapping lembrada pela modelagem de regressão bayesiana.

As bibliotecas

o
pacote ‘rstanarm’ foi
usado para ajustar o modelo, e
‘tidybayes’,
‘bayestestR’ e
‘ver’ foram usados ​​para
investigar as estimativas do modelo (‘bayestestR’ e ‘see’ são ambos
de
suíte ‘easystats’ de
pacotes).

library(rstanarm)
library(tidybayes)
library(bayestestR)
library(see)

Montagem do modelo

o stan_glm() função é o equivalente ‘rstanarm’ de glm(). o
apenas argumentos adicionais a serem incluídos são as distribuições anteriores para o
coeficientes preditores e interceptação. Aqui, eu simplifiquei usando
distribuições normais que não eram muito tendenciosas. Uma análise completa
incluir uma seção em que o impacto de diferentes distribuições anteriores
seria avaliado.

bayes_mansion <- stan_glm(
mansion ~ speech + thought + narrative + depicted,
family = binomial(link = "logit"),
data = locations_joined,
prior = normal(location = 0, scale = 0.5),
prior_intercept = normal(location = 0, scale = 3)
)
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000213 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.13 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.22106 seconds (Warm-up)
#> Chain 1: 0.172703 seconds (Sampling)
#> Chain 1: 0.393763 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 6.3e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.188406 seconds (Warm-up)
#> Chain 2: 0.310647 seconds (Sampling)
#> Chain 2: 0.499053 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 5.3e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.53 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 0.159436 seconds (Warm-up)
#> Chain 3: 0.220296 seconds (Sampling)
#> Chain 3: 0.379732 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 1.5e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 0.169659 seconds (Warm-up)
#> Chain 4: 0.129054 seconds (Sampling)
#> Chain 4: 0.298713 seconds (Total)
#> Chain 4:

Avaliação do modelo

Agora que o modelo está em forma, o próximo passo é inspecionar a parte posterior
distribuições dos coeficientes.

Leia Também  SR2 Capítulo 3 Médio | R-bloggers
cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br
plot(bayes_mansion, prob = 0.50, prob_outer = 0.89)

Cada ponto representa a média da distribuição posterior para o
coeficiente, juntamente com os intervalos de densidade de 50% e 89%. Nós podemos ver
que a interceptação é bastante grande e negativa, indicando que, em
Em média, os X-Men tendem a não visitar a X-Mansion. Comparativamente, o
distribuições para os outros coeficientes são muito pequenas e localizadas
próximo a 0. Isso sugere que eles não representam muito mais
informações sobre se os X-Men visitaram ou não a X-Mansion.

Outro gráfico útil na análise bayesiana é da mais alta densidade
Intervalo (IDH), o menor intervalo de valores de parâmetros que contém um dado
densidade da distribuição. Com o IDH de 89% para posterior
distribuição, podemos dizer que, dada a estrutura do modelo e
dados observados, existe uma chance de 89% de que o valor real do parâmetro esteja
dentro do alcance. Este é um método para entender a confiança de
o valor estimado.

plot(bayestestR::hdi(bayes_mansion, ci = c(0.5, 0.75, 0.89, 0.95)))

A partir do IDH mostrado acima, podemos ver que os coeficientes para o
número de balões de fala (speech) e número de vezes que os caracteres
foram retratados depicted em uma questão foram os preditores mais fortes. o
89% do IDH para speech inclui 0, enquanto o IDH 95% para depicted

inclui 0. Portanto, nenhuma dessas distribuições posteriores é
particularmente emocionante, pois são todos muito pequenos (em conjunto com um
interceptação forte) e tem uma boa chance de realmente ser 0.

Duas outras medidas úteis para a análise bayesiana são as
propabilidade da direção (PD) e os região de prática
equivalência
(CORDA). Sem aprofundar muito, o PD é o
probabilidade de um parâmetro ser positivo ou negativo. Varia de 0,5
para 1, um valor de 1 indica que é definitivamente positivo ou negativo
(ou seja, não zero). A corda é um valor semelhante, mas é responsável pelo efeito
tamanho medindo quanto da distribuição posterior está dentro de um
região que nós (analistas) diríamos ser efetivamente zero. Assim, se
o ROPE é alto, é improvável que o valor do parâmetro tenha
muita importância.

A tabela a seguir fornece um resumo das distribuições posteriores
para este modelo.

bayestestR::describe_posterior(bayes_mansion)
#> Possible multicollinearity between depicted and speech (r = 0.7). This might lead to inappropriate results. See 'Details' in '?rope'.
#> # Description of Posterior Distributions
#>
#> Parameter | Median | 89% CI | pd | 89% ROPE | % in ROPE | Rhat | ESS
#> ------------------------------------------------------------------------------------------------
#> (Intercept) | -1.267 | [-2.095, -0.467] | 0.997 | [-0.181, 0.181] | 0 | 1.000 | 5612.356
#> speech | -0.004 | [-0.010, 0.001] | 0.900 | [-0.181, 0.181] | 100 | 1.001 | 2401.035
#> thought | -0.003 | [-0.011, 0.006] | 0.697 | [-0.181, 0.181] | 100 | 1.001 | 2865.077
#> narrative | 0.004 | [-0.005, 0.013] | 0.771 | [-0.181, 0.181] | 100 | 1.000 | 3264.987
#> depicted | 0.007 | [ 0.001, 0.014] | 0.961 | [-0.181, 0.181] | 100 | 1.000 | 2219.555

Podemos ver que, embora o PD para speech e depicted são próximos
a 1,0, indicando que é provável que não seja zero, o "% em ROPE" é 100%
sugerindo que as diferenças não são importantes.

Verificações preditivas posteriores

A última etapa desta análise é fazer previsões usando o modelo.
Primeiro, podemos fazer previsões sobre os dados fornecidos para ver quão bem o
modelo se encaixa nos dados. Segundo, podemos inserir novos pontos de dados que são
É interessante ver como eles afetam as previsões do modelo.

O gráfico abaixo mostra a distribuição de previsões posteriores no
dados originais. Os valores plotados são a probabilidade prevista de que o
X-Mansion foi visitado na edição de quadrinhos, separado por
não a X-Mansion foi realmente visitada. As duas distribuições parecem
quase idêntico. Isso não é surpreendente, porque os coeficientes se ajustam a
as variáveis ​​eram tão pequenas que não fornecem muito mais
informações para a previsão. Portanto, o modelo é principalmente
confiando na interceptação para calcular uma estimativa.

# From the 'rethinking' package.
logistic <- function (x) {
p <- 1/(1 + exp(-x))
p %
mutate(mansion_predict = logistic(predict(bayes_mansion))) %>%
ggplot(aes(x = mansion_predict, color = mansion, fill = mansion)) +
geom_density(size = 1.2, alpha = 0.2) +
geom_vline(xintercept = 0.5, size = 1.2, lty = 2, color = grey) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Set2") +
theme(legend.position = c(0.65, 0.73)) +
labs(x = "predicted probability of being in the X-Mansion",
y = "probability density",
title = "The Bayesian logistic model's posterior predictions",
color = "was in thenX-mansion",
fill = "was in thenX-mansion")

O segundo tipo comum de previsão posterior é criar dados que
varia uma ou duas variáveis ​​e mantém o restante das variáveis
constante. Desde a depicted teve o maior efeito previsto do
coeficientes de não interceptação, decidi realizar esse procedimento posterior
verificação preditiva dos valores dessa variável. Portanto, eu criei
a pred_data quadro de dados que possui 100 valores em toda a faixa de
depicted e apenas os valores médios para o restante das variáveis.
Fazer previsões nesse conjunto de dados artificiais mostrará o efeito de
a depicted variável enquanto mantém as outras variáveis ​​constantes.

Usando o add_fitted_draws() função de 'tidybayes', 200
previsões foram feitas para os dados artificiais.

pred_data %
summarise(across(issue:narrative, mean)) %>%
mutate(depicted = list(modelr::seq_range(locations_joined$depicted, n = 100))) %>%
unnest(depicted) %>%
add_fitted_draws(bayes_mansion, n = 200)
pred_data
#> # A tibble: 20,000 x 10
#> # Groups: issue, speech, thought, narrative, depicted, .row [100]
#> issue speech thought narrative depicted .row .chain .iteration .draw .value
#>          
#> 1 189. 143. 44.2 48.7 16 1 NA NA 1 0.125
#> 2 189. 143. 44.2 48.7 16 1 NA NA 6 0.143
#> 3 189. 143. 44.2 48.7 16 1 NA NA 15 0.117
#> 4 189. 143. 44.2 48.7 16 1 NA NA 47 0.0890
#> 5 189. 143. 44.2 48.7 16 1 NA NA 86 0.272
#> 6 189. 143. 44.2 48.7 16 1 NA NA 118 0.101
#> 7 189. 143. 44.2 48.7 16 1 NA NA 131 0.158
#> 8 189. 143. 44.2 48.7 16 1 NA NA 133 0.279
#> 9 189. 143. 44.2 48.7 16 1 NA NA 140 0.0990
#> 10 189. 143. 44.2 48.7 16 1 NA NA 154 0.109
#> # … with 19,990 more rows

O gráfico a seguir mostra as curvas logísticas para os dados artificiais.
Cada curva representa uma previsão individual no intervalo de
depicted valores. Os dados originais também são plotados na parte superior.

# Just to shift the `mansion` values for plotting purposes.
locations_joined_mod %
mutate(mansion_num = as.numeric(mansion) + ifelse(mansion, -0.1, 0.1))
pred_data %>%
ggplot(aes(x = depicted, y = .value)) +
geom_line(aes(group = .draw), alpha = 0.1) +
geom_jitter(aes(y = mansion_num, color = mansion),
data = locations_joined_mod,
height = 0.08, width = 0,
size = 2.2, alpha = 0.5) +
scale_color_brewer(palette = "Dark2") +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
labs(x = "depicted",
y = "probability of being in the X-mansion",
color = "was in thenX-Mansion",
title = "Posterior predictions of the effect of the numbernof depictions of the main characters",
subtitle = "All other predictors were held constant at their average value.")

Podemos ver que há uma tendência geral para o modelo prever
que o episódio visitou a X-Mansion como o valor representado
aumenta, mas é uma curva sinusoidal muito gradual porque a
a distribuição posterior estava localizada tão perto de zero. Também podemos ver como
a curva é deslocada para mais perto de 0 para a maioria dos valores de depicted. Isto é
por causa do forte valor de interceptação.


Empacotando

No geral, acho que foi uma comparação interessante entre um
abordagem freqüentista na construção de uma distribuição de valores de coeficientes
e o método bayesiano de analisar um modelo. Eu não estou na posição
para fornecer uma comparação teórica entre as duas abordagens, embora eu
diria que, pessoalmente, interpretar o posterior Bayesiano
distribuição é mais intuitiva do que interpretar o bootstrapped
distribuição.

var vglnk = {key: '949efb41171ac6ec1bf7f206d57e90b8'}; (função (d, t) {var s = d.createElement

Para Deixe um comentário para o autor, siga o link e comente no blog: Posts | Joshua Cook.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da informação / dados.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br