Probabilidades no R

• Vimos no Dia 1 o exemplo bem conhecido da distribuição binomial, que é utilizado para modelar dados que tem a forma de dois resultados, como ‘sucesso’ ou ‘falha’; ‘sim’ ou ‘não’.

No R, é dbinom:

dbinom(6, size=9, prob=0.5)

E a probabilidade cumulativa é calculado com pbinom

• Exercício com sample.

• Simulem a façam gráficos de jogadas de moedas (‘coin toss’), por 1000 jogadas. Se você ganha a jogada, você ganha $1. Se você perde, você perde $1. Quanto você tem depois 1000 jogadas?
  • funções úteis: cumsum, e plot (da base, type='l').

• vamos repetir este jogo 10,000 vezes e colocar os resultados (o último número do 1000, quer dizer, o dinheiro que tem) numa lista chamado results. Para fazer isso, vão precisar de um for loop.

• depois, façam um histograma mostrando a média.

Como fazer este tipo de coisa no R:

• Antes do for loop, criamos um objeto no que vamos colocar os resultados do loop.
• Normalmente, este vai ser uma lista ou um dataframe. Por enquanto, podemos usar uma lista vazia.

• results <- list()

• Um for loop no R é criado usando:
  • for( sequência ) { faz algo }
• Por exemplo, vou adicionar 2 a cada elemento deste vetor:

X <- seq(from = 200, to = 250, by = 1)

X
 [1] 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
[22] 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
[43] 242 243 244 245 246 247 248 249 250

Tentam criar um for loop que usa um vetor de 50 valores (integers), e, em cada iteração, adiciona 1 ao número no vetor na iteração. Por exemplo, se seu vetor começa com 4, depois o for loop, esse 4 vai ser 5.

Então, precisamos de um for loop que vai de 1 até 10000. Podemos usar uma função útil do dplyr, last, para tirar a última jogada da cada 1000 (que é uma iteração).

Para inserir o resultado de uma iteração i na lista, precisamos colocar este resultado no elemento [[i]] na lista.

• R tem funções assim para tirar amostras de várias distribuições ou calcular estatísticas dessas distribuições.

• O mais comumente usado é provavelmente os ‘norms’: dnorm e rnorm, para a distribuição normal. (?Distributions para mais detalhes.)

• A altura de plantas de milho podem ser representados por uma distribuição normal, com média de 145cm e sd de 22cm.

• Acha a proporção de plantas:
  • maior que 100cm (?pnorm)
  • entre 120cm e 150cm
  • 150cm ou menos

pnorm(100, 145, 22, lower.tail = F)

## [1] 0.979595

pnorm(150, 145, 22) - pnorm(120, 145, 22)

## [1] 0.461992

pnorm(150, 145, 22, lower.tail = F)

## [1] 0.4101058

Estimação

Quadratic Approximation

• O pacote rethinking tem a função map()

  • ‘map’ é ‘Maximum A Posteriori’, quer dizer que a função acha o modo da distribuição.
  • assim é semelhante à ‘Máxima verossimilhança’ (Maximum Likelihood Estimation, MLE), que é provavelmente a técnica mais frequentemente usada para produzir uma estimativa.
  • Os resultados são iguais em vários casos, mas MAP nós dá a distribuição posterior, e MLE não, e também podemos incluir informação prévia.

Modelos lineares

• também chamado ‘modelos normais’
• Muitas fenomenas no mundo são ‘normais’
• É a escolha menos surpreendente e mais conservativa
  • só precisamos definir a média e a variância.
  • a distribuição das 10 mil jogadas era normal.

Por exemplo, a distribuição que usamos aqui é uniform, mas se repetimos mil vezes, a distribuição é normal:

library(rethinking)
hist(replicate(1000, sum(runif(16, -1, 1))),
     main=NULL, xlab=NULL, ylab=NULL)

Como escrever modelos lineares:

• Temos uma variável dependente/outcome (resultado), e temos uma medida por cada i.

\[outcome_i \backsim Normal(\mu_i, \sigma)\]

• Isso quer dizer que o outcome é distribuido normalmente com média \(\mu\) e desvio padrão \(\sigma\).

• O “\(\backsim\)” quer dizer que a relação é probabilistica.

• Num modelo linear, \(\mu_i\) é uma equação, que pode ter uma variável independente/predictor (preditor), ou várias.

\[\mu_i = \alpha + \beta \times predictor_i\]

• \(\beta\) é a coeficiente do preditor, que descreve a relação entre o predictor e o outcome, chamado o “slope”;
• \(\alpha\) é chamado o “intercept”.

• \(\alpha\), \(\sigma\) e \(\beta\) também precisam de ser definadas (os valores aqui não importam):

\[\sigma \backsim Uniform(0, 1)\] \[\beta \backsim Normal(0, 10)\] \[\alpha \backsim Normal(0, 10)\]

• O outcome é frequentemente escrito \(y_i\) e o predictor é escrito \(x_i\).
• O intercept é o valor esperado do \(y_i\) quando \(x_i = 0\).
• O slope é o valor do \(y_i\) quando \(x_i\) muda por uma unidade.

• Tudo junto vai ser algo como:

• ou simplesmente \(y_i = \alpha + \beta x_i\) antes de definir as outras partes.

Vamos usar um exemplo do livro do Rethinking, que é um modelo linear. Este modelo vai dizer que altura é previsto por peso.

• ‘Altura’ é nossa variável dependente, ou o outcome;

• ‘Peso’ é o que?

• Variável independente, ou predictor.

library(rethinking)
data(Howell1)
d <- Howell1

• Fazem dplyr::filter para tirar da base tudo mundo mais jovem que 18.
  • (por quê?)
  • Fazem scatterplot das variáveis que vamos usar.
  • Tem relação entre elas? Que tipo?

- O modelo linear (do livro):

height~dnorm(mu,sigma)
mu <- a + b*weight
a ~ dnorm(155, 20)
b ~ dnorm(0, 10)
sigma ~ dunif(0, 50)

Podemos fazer gráficos dos priors para visualizá-las:

curve(dnorm(x, 155, 20), from=100, to=250)

• O que estamos dizendo sobre altura com um prior assim?

Usando map()

model_1 <- rethinking::map(
  alist(
        height ~ dnorm(mu, sigma),
        mu <- a + b*weight,
        a ~ dnorm(178, 20),
        b ~ dnorm(0, 10),
        sigma ~ dunif(0, 50)
        ),
    data=d)

• O que podemos fazer com essas estimativas do modelo?

• Tem duas estrategias principais:
  • gráficos
  • tabelas

• precis() nós dá um sumário das estimativas numa tabela:

precis(model_1)

        Mean StdDev   5.5%  94.5%
a     114.46   1.90 111.42 117.49
b       0.89   0.04   0.83   0.96
sigma   5.07   0.19   4.77   5.38

Que tipo de relação tem entre altura e peso? Como sabemos?

Eu criei aquele scatterplot usando o R da base, pois usar objetos da map() dá um pouco de trabalho para transformar num dataframe. Mas com ggplot2, podemos fazer regressões lineares simples assim dentro da função que cria o gráfico.

Façam um scatterplot, usando ggplot2, que tem ‘height’ no x-axis e ‘weight’ no y-axis, e que tem uma linha de regressão no gráfico.

• Dica: geom_smooth

• Por que “lm” foi usado na chamada da função geom_smooth?

• lm() é a função básica do R para rodar modelos lineares (Linear Model).

• lm() precisa de uma formula (y ~ x) e os dados.

• Vamos rodar height ~ weight com lm(). Depois vamos usar summary().

Essas tabelas são difíceis, vamos quebrar isso em partes.

Residuals:
     Min       1Q   Median       3Q      Max 
-19.7464  -2.8835   0.0222   3.1424  14.7744 

Ok! Vamos praticar um pouco.

library(nycflights13)

flights

Trabalham juntos para criar hipótese de várias relações entre as variáveis dessa base de dados. Quando têm umas ideias, escrevem seus modelos e os rodam usando lm().

Façam sumários e gráficos dos seus modelos. O que vocês acham dos seus hipóteses agora?

Voltando ao Rethinking, carrega os dados Howell1 de novo, e essa vez seleciona só as observações para quem tem menos que 18 anos.

Roda um modelo linear de height ~ weight usando estes dados (pode usar lm() ou map()). Faz sumários do modelo. É razoável? Ou ruim? Por quê?

T-test

• height when age < 20
• height when age > 20

m1 <- map(alist(
        w ~ dbinom(size=9, prob=p), 
        p ~ dunif(0,1)), data=list(w=6))

precis(m1)
  Mean StdDev 5.5% 94.5%
p 0.67   0.16 0.42  0.92

• Isso quer dizer que a distribuição, que suponhamos é normal, é maximizada ao valor de 0.67, que é a média.
• Dado que temos a distribuição inteira, podemos usá-la para investigar perguntas sobre quantidades de probabilidade em cima ou em baixo de um valor particular (questões sobre intervalos), ou sobre estimativas particulares (point estimates, como a média).
  

O que queremos ver com MLE?