Introdução a estatística no `R`

Robert McDonnell

2016-10-03

O que nós vamos fazer neste curso?

revisão do R (tidy data)
introdução às ideias de estatística
como aplicar estas ideias no R
preparação para o curso dos modelos preditivos

O curso é baseado neste livro, que altamente recomendo.

Também acho que R for Data Science, do Hadley Wickham, é útil. Ou simpleR – Using R for Introductory Statistics, que é de graça.

Primeiro:

abra RStudio,
entra $\Rightarrow$ File $\Rightarrow$ New Project.... Vamos criar um novo R Project para este curso.
de New Project..., cliquem no New Directory e Empty Project

Vamos usar vários pacotes, mas por enquanto precisamos desses:

install.packages(c('devtools', 'coda', 
                   'mvtnorm', 'loo',
                   'grid', 'gridExtra',
                   'rstanarm'))
library(devtools)
install_github("rmcelreath/rethinking")
install_github("hadley/tidyverse")

tidyverse tem muitos pacotes que precisamos (dplyr, tidyr, ggplot2, etc.).

O pacote rethinking precisa de Stan, que é um software para fazer sampling. Stan, por sua vez, depende de um compilador de C++. Precisamos disso antes de começar.

https://github.com/stan-dev/rstan/wiki/Installing-RStan-on-Windows

library(tidyverse)
library(rethinking)
library(grid)
library(gridExtra)

O tidyverse é excelente para organizar, arrumar e visualizar dados, mas para fazer estatística, a base R já tem opções ótimas, então nós vamos usar uma mistura dos estilos.

Blueprint de uma análise

uma pergunta: será que os primeiros bebês nascem mais tarde?

e uma outra: como eu começaria responder a uma pergunta assim?

Usando as ferramentas de estatística para responder a essa pergunta:

Coletamos dados confiáveis
Fazemos estatística descritiva: sumarizar, descrever e visualizar os dados, checamos erros;
Construímos modelos;
Testamos hipóteses e fazemos comparação dos modelos;
Produzimos estimativas de valores de interesse.

Desde 1973, o Centers for Disease Control and Prevention tem feito a National Survey of Family Growth (NSFG). Nós vamos usar estes dados para explorar a questão dos bebês e para praticar usar R na análise estatística (descritiva).

Vamos ler os dados da NSFG no R:

preg <- read_csv("2002fempreg.csv")

Quais tipos de variáveis temos?

Usem select() para tirar todas as colunas que não vamos usar no preg.

As que vamos usar são: caseid, prglngth, outcome, pregordr, birthord, birthwgt_lb, birthwgt_oz, agepreg.

(?dplyr::select para quem não lembra.)

library(dplyr)
preg <- select(preg, 
               c(caseid, prglngth, 
               outcome, pregordr, birthord, 
               birthwgt_lb, birthwgt_oz, 
               agepreg))

As variáveis são:

caseid: ID do respondente;
prglngth: duração da gravidez;
outcome: ‘outcome’ da gravidez. 1 indica nascido vivo;
pregordr: número da gravidez (1 para primeiro filho, 2 para o segundo etc.);
birthord: número do nascido vivo (1 para primeiro, etc.);
birthwgt_lb e birthwgt_oz relata o peso dos bebês;
agepreg: idade da mãe no nascimento;

Este survey é cross-sectional, que quer dizer que tem dados de um grupo num momento no tempo. Outros tipos são longitudinal data (e panel data).

É um pouco estranho, mas a variável agepreg é em “centiyears”. Quer dizer, é a idade da mãe multiplicada por 100.
Também birthwgt_lb e birthwgt_oz tem valores que não queremos (e.g. 99 = DON'T KNOW, 98 = REFUSED, 97 = NOT ASCERTAINED). (E tem um erro.)
Vamos usar as ferramentas de tidy data para limpar estes dados, especificamente filter() e mutate().

1: Achem e tirem valores estranhas em birthwgt_lb e birthwgt_oz (tirem da tabela). Lembrem-se que tem 16 ounces por cada libra.
- Dica:
```
sort(unique(preg$birthwgt_lb))
sort(unique(preg$birthwgt_oz))
```
2: Mudem agepreg para ter a idade ‘normal’ das mães;

3: Criem uma nova variável, total_lb, que é birthwgt_lb + birthwgt_oz. Essa variável deve ser em libras (lb).

4: Dado que tratar em pounds e ounces é chato, vamos criar uma nova variável peso_kg que é o peso em quilogramas (1lb = 0.453592 kg). Também vamos usar só outcome = 1.

preg <- preg %>% 
  filter(birthwgt_lb <= 15,
         birthwgt_oz <= 15) %>% 
  mutate(agepreg = agepreg/100) %>% 
  mutate(total_lb = birthwgt_lb + 
           (birthwgt_oz/16)) %>% 
  mutate(peso_kg = total_lb*0.453592) %>% 
  filter(outcome == 1)

Estatística Descritiva

Gráficos ajudam bastante a entender a distribuição dos dados. Um método é simplesmente contar quantas observações há, e as mais comuns são boxplots, barplots e histogramas.

Usamos histogramas com variáveis contínuas (numéricos numa escala) e barplots com variáveis categóricas (pode colocar em categorias ou num ordem). Boxplots são usados por variáveis contínuas visualizadas por níveis de uma variável categórica.

library(ggplot2); library(grid); library(gridExtra)
bar <- ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut))
hist <- ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = carat), binwidth = 0.5)
box <- ggplot(data = diamonds, 
              mapping = aes(x = cut, y = price)) +
  geom_boxplot()
grid.arrange(bar, hist, box, ncol=3)

Vamos usar peso_kg e pregordr para construir um histograma e um barplot. Qual variável é continua e qual é categórica? (na verdade, não é tão óbvio)

> unique(preg$pregordr)
[1]  1  2  3  4  5  6  7  8  9 10 11 12 17

Embora estes valores sejam números, os valores nesta escala representam categorias. pregordr é a ordem dos bebês da mãe: primeiro, segundo, terceiro, etc. (17??!). Variáveis categóricas que têm ordem assim são chamados “ordinais”.

Frequentemente, variáveis ordinais tem uma escala menor, por exemplo numa survey, na qual pessoas respondem sobre sua opinião sobre algo, com 1 sendo ‘ruim’ e 5 ‘excelente’.

Este exemplo não é tão claro, mas normalmente é bastante simples ver qual tipo é uma variável.
Exemplos de variáveis categóricas (também chamado ‘discretas’) são:
- sexo;
- estado civil;
- nacionalidade
Exemplos de variáveis continuas são:
- temperatura;
- renda;
- totais de importação.

Dado que pregordr pode ser considerado discreta, vamos fazer um barplot com ela. Em ggplot2, usamos a sintaxe:

ggplot(preg, aes(pregordr)) + geom_bar()

O que podemos dizer sobre a distribuição dessa variável?
Onde é o valor mais comum (o modo)?

geom_histogram funciona semelhantemente.

Uma coisa importante com histogramas é o número de “bins”, quer dizer, quantas barras quer ter. Destaques importantes dos dados podem ser escondidos usando bins demais ou poucos.

Criem três histogramas com peso_kg, com o número de bins sendo 5, 15 e 50. Coloquem juntos num gráfico só.

aa <- ggplot(preg, aes(peso_kg)) + 
  geom_histogram(bins = 5)
ab <- ggplot(preg, aes(peso_kg)) + 
  geom_histogram(bins = 15)
ac <- ggplot(preg, aes(peso_kg)) + 
  geom_histogram(bins = 50)
grid.arrange(aa, ab, ac, ncol=3)

Como tudo no ggplot2, temos bastante opções:

fi <- ggplot(preg, aes(x=peso_kg)) +
  geom_histogram(binwidth=.5, colour="black", fill="white")

fy <- ggplot(preg, aes(x=peso_kg)) + geom_density()

fo <- ggplot(preg, aes(x=peso_kg)) + 
  geom_histogram(aes(y=..density..),
                 binwidth=.5,
                 colour="black", fill="white") +
  geom_density(alpha=.2, fill="#FF6666")

grid.arrange(fi, fy, fo, ncol=3)  # põe as três juntos

Podemos ter linha para indicar a “média”:

fo + geom_vline(aes(xintercept=mean(rating, na.rm=T)),#Tira NA
  color="red", linetype="dashed", size=1)

Como sumarizar essa distribuição?

A média não é igual ao modo;
A distribuição é quase ‘normal’, ou ‘Gaussiana’;
É “right-skewed”;
Parece que não tem ‘outliers’.

Sumários

Para sumarizar e descrever distribuições, normalmente focamos nas seguintes propriedades (“summary statistics”):

tendência central: Os valores têm ‘cluster’ num ponto específico?

modos: Tem mais que um cluster?

‘spread’: Quanta variabilidade tem nos valores?

‘tails’: Quão rápido os valores caem pra zero quando movemos para longe dos modos?

‘outliers’: Temos valores extremos longe dos modos?

Tendência Central

A média aritmética é a soma de uma coleção de números dividido pela quantidade de números na coleção: \[\bar x = \frac 1 n \sum_{i=1}^n x_i .\]

# No R:
mean(x)

Se você lê livros de estatística, vai ver que o “sample mean”, o “mean” e o “expectation” são usados quase indistintamente. É frequentemente escrito $\mu$.

Pumpkins <- data_frame(Pumpkins= c(1, 3, 3, 591))
ggplot(Pumpkins, aes(x = Pumpkins)) +
         geom_histogram(colour="black", fill="white") + 
         geom_vline(aes(xintercept=mean(Pumpkins, na.rm=T)),
           color="#104E8B", linetype="dashed", size=1)

Variância

A variância é expressada assim:

\[S^2 = \frac 1 n \sum_i (x_i - \bar x)^2,\]

# No R:
var(x)

onde $\bar x$ é a média que acabamos de ver.

Desvio Padrão

Na fórmula da variância, o termo $x_i - \bar x$ é a “deviação da média”, e a variância é isso ao quadrado. A raiz quadrada da variância, S, é o desvio padrão (standard deviation). É frequentemente escrito $\sigma$, e é útil porque é nas mesmas unidades de mensuração do que a média.

# No R:
sd(x)

Vamos olhar a variável agepreg, a idade das mães.

O que é a média da agepreg?
A variância?
O sd?

Podemos também calcular o ‘intervalo interquartil’ (interquartile range), o mínimo, o máximo, a mediana no R usando summary(x).

mean(preg$agepreg)
[1] 38.57059

var(preg$agepreg)
[1] 30.84675

sd(preg$agepreg)
[1] 5.553985

summary(preg$agepreg)
 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
10.83   20.50   24.33   24.94   28.89   44.08

As vezes nossos dados tem muitos ‘outliers’ (valores extremos), e pode ser que e média não seja muito representativo da distribuição; assim, a mediana é melhor, que é o valor no meio da distribuição.

Criem sumários simples (não visuais) da prglngth.
Façam dois histogramas da prglngth no mesmo gráfico.
- Dica: usem facet_grid.
Qual visualização seria boa para ver a diferença entre os dados dos primeiros bebês e dos outros? Façam uma, que mostre essa diferença entre os dois grupos: “primeiros” e “outros”.

# criar outra variável para distinguir entre primeiros 
# bebês e outros:

preg <- preg %>% 
  mutate(grupo = if_else(birthord == 1, 
                         "primeiro", "outro"))

ggplot(preg, aes(x=prglngth, fill=grupo)) +
  geom_histogram(binwidth=.5, colour="black") +
  facet_grid(grupo ~ .)

ggplot(data = preg, 
        mapping = aes(x = grupo, y = prglngth)) +
   geom_boxplot()

Não parecem tão diferentes! Mas para testar a diferença estatisticamente, temos que olhar primeiro à inferência e probabilidade.

relações entre variáveis

scatterplots são populares para ver (de um jeito descritiva) relações entre duas variáveis. Por exemplo (?faithful):

ggplot(data = faithful) + 
  geom_point(mapping = aes(x = eruptions, y = waiting))

Modelos estatísticos

O que é um ‘modelo’?
Por que usamos modelos?

Modelos são andaimes: construímos um andaime, testamos que tudo está sólido, e daí nós construímos o que realmente queremos construir. Depois, tiramos o andaime e a descartamos, mas foi isso que nós permitirem construir algo.
São também como robôs: poderosos, mas só fazem o que nós pedimos; não são ‘verdadeiros’ nem ‘falsos’, são ferramentas para realizar algo.

Por que este estilo de estatística pode ser melhor?
É especializado demais: inovação é difícil;
Muitas vezes, as ferramentas erradas são usadas;
Ferramentas são usadas sem pensar;
Perdemos os avanços na modelagem moderna.
E desenhado para testar hipóteses:
- Mas na maioria dos casos, um hipótese nulo é ‘testado’;
- Não estamos testando o hipótese da pesquisa;
- (Para mais ver o Statistical Rethinking.)

Modelagem moderna

Pode ser um procedimento de teste;
Também podemos:
- medir
- prever
- comparar e debater entre modelos.
  - previsão etc.

Para fazer isso, vamos focar em três áreas:

1: Análise de dados Bayesiana

2: Modelos multinível

3: Comparação de modelos

Estatística frequentista ou clássica cabe como caso especial de estatística Bayesiana, então vamos ver que em muitos casos os resultados são quase iguais. As vezes é mais fácil fazer a versão frequentista, as vezes não. Mas o coração dos dois tipos é o mesmo (o “Likelihood”).

Estatística Bayesiana

\[ Posterior \propto Likelihood \times Prior \]

todos são distribuições;
o likelihood é o número de formas algo pode acontecer segundo nosso pressupostas;
o prior é informação prévia;
o posterior é a distribuição de probabilidade de um elemento ou elementos dos nossos modelos.

modelos Bayesianos ‘aprendem’ cada vez que recebe informação (Bayesian updating). Um exemplo:

source_gist(
  id="c7cf1ab6a9f5e983f07a067a8fe32126", 
            filename="BayesianUpdate.R")

Modelos Multiníveis

Modelos no qual os parâmetros são baseado em outros parâmetros, ou seja, tem modelo(s) dentro do modelo.

por que faz isso?

Principalmente para incluir agrupamento nos modelos (‘clustering’). Se nós sabemos que temos grupos em nossos dados, nossos modelos vão ser melhores se incluímos este fato.

Exemplos:
- Políticos dentro de um partido;
- Países de regiões específicas;
- Crianças em escolas diferentes;
- Medidas feitas em momentos diferentes.

Muitos fenômenos são naturalmente pensados assim. Estes modelos também modelam variância melhor, particularmente a variância entre grupos diferentes.

Comparação de modelos

Baseado em teoria de informação;
Ajuda com overfitting (o modelo é especifico demais à amostra);
Nos permite evitar o uso de uma hipótese nula e ao invés usar outros modelos para testar nossas hipóteses. Isso é importante, pois em muitas situações a hipótese nula não faz sentido ou não podemos formular uma que é realística.

Probabilidade

O likelihood

Tem quatro bolinhas nesta sacola, e elas tem duas cores: azul e branca.

Tem cinco possibilidades:

4 brancas;
3 brancas e uma azul;
2 brancas e duas azuis;
1 brancas e três azuis;
4 azuis.

Essas são nossas conjeturas/hipóteses.

Vamos considerar a conjetura que a sacola tem:

Vamos chamar isso de hipótese 1 (H1). Sob este pressuposto, a primeira retirada tem essas possibilidades:

Nós vamos tirar três bolinhas da sacola. Agora tem essas 64 possibilidades:

Tiramos três bolinhas dessa sacola, colocando a bolinha retirada de volta cada vez:

Estes são nossos dados (Data) agora. O H1 têm 3/64 jeitos de fazer essa retirada.

Umas das outras conjeturas são rejeitadas agora, pois não são logicamente compatíveis com os dados:

As que sobrarem tem as seguintes formas de fazer o Data:

O Prior

João na fábrica de bolinhas te fala que as azuis são mais raras, podemos incluir essa informação prévia:

Esse ‘prior’ é que distinguia estatística Bayesiana, mas hoje em dia a estatística não-Bayesiana também tem formas de incluir essa informação, como penalized likelihoods, por exemplo.

contar assim não é prático, mas podemos usar a matemática de probabilidade para tornar isso mais fácil.

A proporção dos azuis segundo nossa hipótese é p, o valor do parâmetro;

A plausabilidade de H1 depois ver Data é proporcional ao número de formas H1 pode produzir o Data multiplicado pela informação prévia.

“A informação prévia” é chamado o prior, e “o número de formas que H1 pode produzir o Data” é o likelihood.

“A plausabilidade de H1 depois ver Data” é o posterior.

Isso é expressado na linguagem de probabilidade assim:

\[ P(H1|Data) \propto P(Data|H1) \times P(Data) \]

ou seja: