Importante até agora

• Fazer estatística descritiva
• Entender como ler um modelo
• Rodar um modelo
• Avaliar o modelo

Fazer tudo isso com modelos lineares de um preditor é razoavelmente simples.

Linguagem dos modelos

\[ y_i \backsim Normal(\mu_i, \sigma) \\ \mu_i = \alpha + \beta x_i \\ \beta \backsim Normal(0, 10)\\ \alpha \backsim Normal(0, 10)\\ \sigma \backsim Uniform(0, 1) \]

Interpretação

\[y_i = \alpha + \beta x_i\]

• O que é \(\alpha\)?

• O que \(\beta\) significa?

• Fizemos já um modelo linear de um preditor, mas têm muitos outros tipos de modelos para lidar com vários outros tips de dados, ou hipótese, etc.

• Um teste estatístico simples, e comum, é o t-test.
• O ideia do t-test é simples:
  • temos dois grupos, cada um com uma média;
  • computamos a diferença das médias, e padronizar isso relativa ao desvio padrão dos valores dos grupos.
  • Essa diferença é chamado o “valor t”, ou t-value.
  • Queremos saber se este valor é ‘significamente’ diferente de zero, então comparamos nosso valor t a uma distribuição amostral de valores t (acostumavam ter que fazer com tabelas, mas R faz para nós com qt()).

Vamos fazer um t-test das alturas dos homens e das mulheres na base Howell1.

• Vão precisar dos pacotes rethinking e tidyverse.

• Precisamos arrumar os dados para ter dois dataframes, um de homens, o outro de mulheres.

males <- filter(Howell2, male==1)
females <- filter(Howell2, male==0)

A função no R é t.test. (?t.test para ajuda)

A interpretação do teste t, como tudo no ferquentismo, não é sempre fácil. Temos opções Bayesianas, porém.

library(BEST)  # instalam se não têm

bayes_t <- BESTmcmc(males$height, 
                    females$height, 
                    rnd.seed = 1234)

Pessoalmente, eu acho os gráficos mais fáceis interpretar. Mas, como vimos, a versão Bayesiana demora mais.

Vamos fazer um t-test Bayesiana das médias do peso dos homens e mulheres. Depois façam gráficos e sumários (BESTmcmc tem mais opções do que o que eu mostrei). Comparem com t-test tradicional.

Tem uns outros pacotes que usam MCMC no R, e todos usam um tipo de linguagem mais rápido, como C++ para rodar os modelos. Um pacote que nós vamos usar é o rstanarm, que usa a linguagem de modelagem probabilistica Stan no R. (r + stan e tinha um outro pacote chamado ‘arm’: (Analysis using Regression Models))

# install.packages("rstanarm")
library(rstanarm)

arm é do livro Data Analysis Using Regression and Hierarchical/Multilevel Models por Andrew Gelman e Jennifer Hill, que altamente recomendo.

• rstanarm é uma tentativa de usar a linguagem do lm() para rodar modelos Bayesianos. Já tem modelos pre-compilados (normalmente o código do R tem que ser compilado para C++), então é rápido.

• Com lm() vimos que um modelo linear pode ser rodado com:

model <- lm(formula, data)

• rstanarm põe stan_ na frente do lm:

model <- stan_lm(formula, data, prior)

Vamos rodar um modelo de \(height \backsim weight\) usando stan_lm(), ou stan_glm() com family = "gaussian".

• ?stan_lm()

• ?stan_glm()

stan_model <- stan_lm(height ~ weight, 
                      data = Howell2,
                      prior = R2(location = 0.5, 
                                 what = "mean"), 
                      chains = 4, cores = 4, 
                      seed = 1234, iter = 10000)

#There were 11729 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
#http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems

Este modelo rodou perfeitamente bem com lm(). Agora, contudo, tem avisos de ‘divergent transitions’. Estes são problemas com o sampler, que indicam um modelo fraco. Com rstanarm foi mais fácil ver estes problemas. (Notam que mais iterações vai melhora isso as vezes.)

• Podemos fazer sumários com summary();

• Olhar só às coeficientes com coef();

• Análisar a qualidade do MCMC com stan_trace();

• Fazer vários gráficos com ?plot.stanreg > rstanarm-plots;

• E, melhor de tudo, faz app interativo com launch_shinystan()

Provavelmente a melhor coisa sobre estatística Bayesiana é o fato que você pode usar os resultados do modelo para calcular o posterior predictive distribution, qua usamos para fazer previsões sobre o valor de dados novos usando o nosso modelo. Isso é muito fácil com rstanarm, a função é:

• posterior_predict(modelo, newdata)

Ou para fazer gráficos:

• pp_check(modelo)

• A distribuição prevista é bimodal! Por que será?

• Lembra o t-test? Homens e mulheres tem médias diferentes.

• Para incluir sexo no modelo, precisamos uma outra variável.

• Fazer isso vai nós mostrar duas coisas: multiple regression e regressão com variáveis categoricas.

• Antes de rodar este modelo com mais que um preditor, vamos fazer com só male, deixando peso para o lado por enquanto.

ANOVA

Em nosso modelo \(height_i ~ \alpha + \beta * weight_i\), ambos de height e weight eram variáveis continuas. Um outro teste estatística que é comum é ANOVA, ou análise de variância. Tradicionalmente, todos estes testes tinha nomes diferentes, e as vezes notação diferente, mas com nossa linguagem de modelos que estamos usando vira fácil, pois ANOVA é um modelo linear com preditores categoricas.

• Variável categoricas são também chamado “factors” e “variáveis discretas”. Quando tem a forma 0/1 são “dummies” ou “indicadores”.

• Temos uma variável categorica no Howell1 que podemos utilizar?

• Por que essa variável está já pronto para usar?

• Vamos fazer modelo de height ~ male. Como escrever isso, usando só o variável dependente e depois a regressão linear de \(\mu\)? (Não precisar especificar priors)

\[ height_i \backsim Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta male_i\\ \] rstanarm tem o modelo ANOVA dentro do pacote já. No R, a função é normalmente aov() (podem usar também), e no rstanarm, é stan_aov():

stan_aov(formula, data, prior)

Podem usar stan_aov(), mas sempre dá erro comigo, então eu vou usar stan_lm(), dado que ANOVA é um modelo linear.

Howell_aov <- stan_lm(height ~ male, 
                       data = Howell2, 
                       prior = R2(location = 0.5,
                                  what="mean"),
                       cores=2)

Vamos agora incluir mais uma variável no modelo, que é a variável peso que já usamos. Como eu posso escrever este modelo agora, dado que tem dois preditores?

• \[y_i \backsim Normal(\mu_i, \sigma)\\ \mu_i = \alpha + \beta_1 male_i + \beta_2 weight_i\]

• vamos usar stan_lm() para rodar este modelo (ou lm() se preferem).

multi_stan <- stan_lm(height ~ male + weight, 
                       data = Howell2, 
                       prior = R2(location = 0.5,
                                  what="mean"),
                       cores=2)

• Vamos fazer gráficos e sumários deste modelo como antes. Os objetos que resultam de stan_lm() tem outros métodos disponíveis também (vejam ?stanreg-methods). Por exemplo, pode fazer uma matriz de variância-covariância, que mostra a variância dentro das variáiveis e entre elas, com vcov().

• Como interpretamos as coeficientes agora?

• Qual modelo é melhor:
  • stan_model? \(y_i = \alpha + \beta weight_i\)
  • Howell_aov? \(y_i = \alpha + \beta male_i\)
  • multi_stan? \(y_i = \alpha + \beta_1 male_i + \beta_2 weight_i\)

• Como sabemos? Tem método formal para saber?

Comparação

O rstanarm já tem o pacote loo dentro dele, e podemos acessá-lo com a função loo() e depois, compare():

loo1 <- loo(stan_model)
loo2 <- loo(Howell_aov)
loo3 <- loo(multi_stan)

compare(loo1, loo2, loo3)

• Essa nós dá estes resultados:

> compare(loo1, loo2, loo3)
     looic   se_looic elpd_loo se_elpd_loo p_loo   se_p_loo
loo3  1997.3    35.2   -998.7     17.6         4.7     1.1 
loo1  2116.6    29.5  -1058.3     14.7         3.0     0.4 
loo2  2172.3    30.5  -1086.2     15.2         3.2     0.5 

• O que queremos ver aqui? O looic tem que ser o menor, e o elpd_loo tem que ser o maior.

• Qual é o modelo preferido pelo loo()?

Façam a mesma coisa usando waic(). Os resultados são diferentes?

Prática

• Vamos colocar tudo isso junto num exercicio. Primeiro, vamos carregar uma base que tem no rstanarm:

data(kidiq)
kidiq

?kidiq para mais informações nestes dados.