Processing math: 2%
+ - 0:00:00
Notes for current slide
Notes for next slide

Introdução aos Modelos Lineares

com códigos em R

Setembro de 2024

1 / 61

Agenda

  • O que é e quando usar
  • Parâmetro vs estimador
  • Teste de Hipóteses e valor-p
  • Interpretação dos parâmetros
  • Regressão Linear Múltipla
  • Preditores Categóricos
  • Transformações Não Lineares dos Preditores
2 / 61

Introdução

4 / 61

Motivação

Somos consultores e fomos contratados para dar conselhos para uma empresa aumentar as suas vendas.

Obtivemos o seguinte banco de dados

  • PERGUNTA: Jornal tem mais retorno do que as demais mídias? Quantas vendas terão se eu investir X em jornais?
5 / 61

Motivação

Somos consultores e fomos contratados para dar conselhos para uma empresa aumentar as suas vendas.

Obtivemos o seguinte banco de dados

  • PERGUNTA: Jornal tem mais retorno do que as demais mídias? Quantas vendas terão se eu investir X em jornais?
6 / 61

Modo - Regressão e Classificação

Existem dois principais tipos de problemas em modelagem estatística:

Regressão

Y é uma variável quantitativa contínua

  • Volume de vendas
  • Peso
  • Temperatura
  • Valor de Ações

Classificação

Y é uma variável qualitativa discreta.

  • Fraude/Não Fraude
  • Pegou em dia/Não pagou
  • Cancelou assinatura/Não cancelou
  • Gato/Cachorro/Cavalo/Outro
7 / 61

Exemplo de reta de regressão

8 / 61

Definições e Nomenclaturas

A tabela por trás (do excel, do sql, etc.)

midia investimento vendas
TV 220.3 24.7
newspaper 25.6 5.3
newspaper 38.7 18.3
radio 42.3 25.4
radio 43.9 22.3
TV 139.5 10.3
radio 11.0 7.2
radio 1.6 6.9
9 / 61

Definições e Nomenclaturas

  • X1, X2, ..., Xp: variáveis explicativas (ou variáveis independentes ou features ou preditores).
  • \boldsymbol{X} = {X_1, X_2, \dots, X_p}: conjunto de todas as features.
  • Y: variável resposta (ou variável dependente ou target).
  • Ŷ: valor esperado (ou predição ou estimado ou fitted).
  • f(X) também é conhecida também como "Modelo" ou "Hipótese".

No exemplo:

  • X_1: midia - indicadador de se a propaganda é para jornal, rádio, ou TV.
  • X_2: investimento - valor do orçamento
  • Y: vendas - qtd vendida
10 / 61

Definições e Nomenclaturas

Observado versus Esperado

  • Y é um valor observado (ou verdade ou truth)
  • Ŷ é um valor esperado (ou predição ou estimado ou fitted).
  • Y - Ŷ é o resíduo (ou erro)

Por definição, \hat{Y} = f(x) que é o valor que a função f retorna.

11 / 61

Definições e Nomenclaturas

A tabela por trás depois das predições

midia investimento vendas regressao_linear
TV 220.3 24.7 17.8
newspaper 25.6 5.3 13.8
newspaper 38.7 18.3 14.4
radio 42.3 25.4 15.0
radio 43.9 22.3 15.1
TV 139.5 10.3 13.6
radio 11.0 7.2 13.4
radio 1.6 6.9 12.9
12 / 61

Por que ajustar uma f?

  • Predição
  • Inferência

Predição

Em muitas situações X está disponível facilmente mas, Y não é fácil de descobrir. (Ou mesmo não é possível descobrí-lo). Queremos que \hat{Y} = \hat{f}(X) seja uma boa estimativa (preveja bem o futuro). Neste caso não estamos interessados em como é a estrutura \hat{f} desde que ela apresente predições boas para Y.

Por exemplo:

  • Meu cliente vai atrasar a fatura no mês que vem?
13 / 61

Por que ajustar uma f?

  • Predição
  • Inferência

Inferência

Em inferência estamos mais interessados em entender a relação entre as variáveis explicativas X e a variável resposta Y.

Por exemplo:

  • A droga é eficaz para o tratamento da doença X?
  • Quanto que é o impacto nas vendas para cada real investido em TV?

Neste material focaremos em inferência.

14 / 61

Por que ajustar uma f?

15 / 61

O que é e quando usar

Regressão Linear Simples

y = \beta_0 + \beta_1x

Exemplo:

dist = \beta_0 + \beta_1speed

Ver ISL página 61 (Simple Linear Regression).

16 / 61

O que é e quando usar

Regressão Linear Simples

y = \beta_0 + \beta_1x

Exemplo:

dist = \beta_0 + \beta_1speed

Ver ISL página 61 (Simple Linear Regression).

17 / 61

O que é e quando usar

Regressão Linear Simples

y = \beta_0 + \beta_1x

Exemplo:

dist = \beta_0 + \beta_1speed

18 / 61

O que é e quando usar

Regressão Linear Simples

y = \beta_0 + \beta_1x

# ajuste de uma regressão linear simples no R
melhor_reta <- lm(dist ~ speed, data = cars)
melhor_reta
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Coefficients:
## (Intercept) speed
## -17.579 3.932
19 / 61

Métricas - "Melhor reta" segundo o quê?

Queremos a reta que erre menos.

Exemplo de medida de erro: Root Mean Squared Error.

RMSE = \sqrt{\frac{1}{N}\sum(y_i - \hat{y_i})^2}

20 / 61

Métricas - "Melhor reta" segundo o quê?

Queremos a reta que erre menos.

Exemplo de medida de erro: Root Mean Squared Error.

RMSE = \sqrt{\frac{1}{N}\sum(y_i - (\hat{\beta_0} + \hat{\beta_1}x))^2}

Ou seja, nosso objetivo é

Encontrar \hat{\beta}_0 e \hat{\beta}_1 que nos retorne o ~menor~ RMSE.

21 / 61

Qual o valor ótimo para \beta_0 e \beta_1?

No nosso exemplo, a nossa HIPÓTESE é de que

dist = \beta_0 + \beta_1speed

Então podemos escrever o Erro Quadrático Médio como

EQM = \frac{1}{N}\sum(y_i - \hat{y_i})^2 = \frac{1}{N}\sum(y_i - \color{red}{(\hat{\beta}_0 + \hat{\beta}_1speed)})^2

Com ajuda do Cálculo é possível mostrar que os valores ótimos para \beta_0 e \beta_1 são

\hat{\beta}_1 = \frac{\sum(x_i - \bar{x})(y_i - \bar{y})}{\sum(x_i - \bar{x})^2}

\hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x}

Já que vieram do EQM, eles são chamados de Estimadores de Mínimos Quadrados.

# lembrete: exercício 2 do script!
22 / 61

Depois de estimar...

\hat{y} = \hat{\beta}_0 + \hat{\beta}_1x

Exemplo:

\hat{dist} = \hat{\beta}_0 + \hat{\beta}_1speed

Colocamos um \hat{} em cima dos termos para representar "estimativas". Ou seja, \hat{y}_i é uma estimativa de y_i.

No nosso exemplo,

  • \hat{\beta}_0 é uma estimativa de \beta_0 e vale -17.579.
  • \hat{\beta}_1 é uma estimativa de \beta_1 e vale 3.932.
  • \hat{dist} é uma estimativa de dist e vale -17.579 + 3.932 x speed.
# Exercício: se speed for 15 m/h, quanto que
# seria a distância dist esperada?
23 / 61

Curiosidade - Método numérico

Podemos encontrar a reta que erre menos por meio de algoritmos numéricos.

Exemplo: Modelo de regressão linear f(x) = \beta_0 + \beta_1 x.

24 / 61

Teste de Hipóteses e valor-p

Exemplo: relação entre População Urbana e Assassinatos.

Modelo proposto:

y = \beta_0 + \beta_1 x

Hipótese do pesquisador:

"Assassinatos não estão relacionados com a proporção de população urbana de uma cidade."

Tradução da hipótese em termos matemáticos:

H_0: \beta_1 = 0 \space\space\space\space\space vs \space\space\space\space H_a: \beta_1 \neq 0 Se a hipótese for verdade, então o \beta_1 deveria ser zero. Porém, os dados disseram que \hat{\beta}_1 = 0.02.

0.02 é diferente de 0.00?

25 / 61

0.02 é diferente de 0.00?

Saída do R

##
## Call:
## lm(formula = Murder ~ UrbanPop, data = USArrests)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.537 -3.736 -0.779 3.332 9.728
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.41594 2.90669 2.207 0.0321 *
## UrbanPop 0.02093 0.04333 0.483 0.6312
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.39 on 48 degrees of freedom
## Multiple R-squared: 0.00484, Adjusted R-squared: -0.01589
## F-statistic: 0.2335 on 1 and 48 DF, p-value: 0.6312
26 / 61

0.02 é diferente de 0.00?

Conceito importante: Os estimadores ( \hat{\beta}_0 e \hat{\beta}_1 no nosso caso) têm distribuições de probabilidade.

Simulação de 1000 retas (ajustadas com dados diferentes).

distrib_params

27 / 61

0.02 é diferente de 0.00?

Conceito importante: Os estimadores ( \hat{\beta}_0 e \hat{\beta}_1 no nosso caso) têm distribuições de probabilidade.

A Teoria Assintótica nos fornece o seguinte resultado:

t = \frac{\hat{\beta_1} - \beta_1}{\hat{\sigma}_{\beta_1}} \overset{\text{a}}{\sim} t(N - 2)

Em que

\hat{\sigma}_{\beta_1} = \sqrt{\frac{EQM}{\sum(x_i - \bar{x})^2}}

Usamos essas distribuições assintóticas para testar as hipóteses.

28 / 61

0.02 é diferente de 0.00?

Conceito importante: Os estimadores ( \hat{\beta}_0 e \hat{\beta}_1 no nosso caso) têm distribuições de probabilidade.

A Teoria Assintótica nos fornece o seguinte resultado:

t = \frac{\hat{\beta_1} - \beta_1}{\hat{\sigma}_{\beta_1}} \overset{\text{a}}{\sim} t(N - 2)

Em que

\hat{\sigma}_{\beta_1} = \sqrt{\frac{EQM}{\sum(x_i - \bar{x})^2}}

Usamos essas distribuições assintóticas para testar as hipóteses.

No nosso exemplo, a hipótese é H_0: \beta_1 = 0, então

t = \frac{0.02 - 0}{0.04} = 0.48

29 / 61

0.02 é diferente de 0.00?

Então agora podemos tomar decisão! Se a estimativa cair muito distante da distribuição t da hipótese 0, decidimos por rejeitá-la. Caso contrário, decidimos por não rejeitá-la como verdade.

testes_t_estrelas

## NO R:
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.41594 2.90669 2.207 0.0321 *
## UrbanPop 0.02093 0.04333 0.483 0.6312
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 / 61

Interpretação dos parâmetros

y = \color{darkgblue}{\beta_0} + \color{darkgreen}{\beta_1}x

Interpretações matemáticas

\color{darkgblue}{\beta_0} é o lugar em que a reta cruza o eixo Y.

\color{darkgreen}{\beta_1} é a derivada de Y em relação ao X. É quanto Y varia quando X varia em 1 unidade.

Interpretações estatísticas

\color{darkgblue}{\beta_0} é a distância percorrida esperada quando o carro está parado (X = 0).

\color{darkgreen}{\beta_1} é o efeito médio na distância por variar 1 ml/h na velocidade do carro.

31 / 61

Teste de Hipóteses e valor-p

Exercício 3 do script:

No R, use a função summary(melhor_reta) (ver slide 9) para decidir se speed está associado com dist. Descubra o valor-p associado.

lembrete: o banco de dados se chama cars.

Exercício 4 do script:

Interprete o parâmetro \beta_1.

32 / 61

Intervalo de confiança para \beta_1

Intervalo de 95%

[\hat{\beta} - 1,96 * \hat{\sigma}_{\beta_1}, \hat{\beta} + 1,96 * \hat{\sigma}_{\beta_1}]

Intervalo de 90%

[\hat{\beta} - 1,64 * \hat{\sigma}_{\beta_1}, \hat{\beta} + 1,64 * \hat{\sigma}_{\beta_1}]

Intervalo de 1 - \alpha%

[\hat{\beta}_1 - q_{\alpha} * \hat{\sigma}_{\beta_1}, \hat{\beta} + q_{\alpha} * \hat{\sigma}_{\beta_1}]

em que q_\alpha é o quantil da Normal(0, 1).

Ver ISL página 66 (Assessing the Accuracy of the Model).

33 / 61

O modelo está bom?

EQM e EPR

ERP significa Erro Padrão dos Resíduos e é definido como

EPR = \frac{\sum(y_i - \hat{y_i})^2}{N - 2} = \frac{SQR}{N - 2} O 2 no denominador decorre do fato de termos 2 parâmetros para estimar no modelo.

  • Se y_i = \hat{y}_i \space\space\space \rightarrow \color{green}{EPR = 0 \downarrow}
  • Se y_i >> \hat{y}_i \rightarrow \color{red}{EPR = alto \uparrow}
  • Se y_i << \hat{y}_i \rightarrow \color{red}{EPR = alto \uparrow}

Problema: Como sabemos se o EPR é grande ou pequeno?

Ver ISL página 68 (Assessing the Accuracy of the Model).

34 / 61

O modelo está bom?

R-quadrado ( R^2 )

R^2 = 1 - \frac{\sum(y_i - \color{salmon}{\hat{y_i}})^2}{\sum(y_i - \color{royalblue}{\bar{y}})^2} = 1 - \frac{\color{salmon}{SQR}}{\color{royalblue}{SQT}}

R^2 \approx 1 \rightarrow \color{salmon}{SQR} << \color{royalblue}{SQT}. R^2 \approx 0 \rightarrow \color{salmon}{reta} \text{ em cima da } \color{royalblue}{reta}.

Problema do R^2 é que ele sempre aumenta conforme novos preditores vão sendo incluídos.

Ver ISL página 68 (Assessing the Accuracy of the Model).

35 / 61

O modelo está bom?

R-quadrado ajustado

R^2 = 1 - \frac{\color{salmon}{SQR}}{\color{royalblue}{SQT}}\frac{\color{royalblue}{N-1}}{\color{salmon}{N-p}}

Em que p é o número de parâmetros do modelo (no caso da regressão linear simples, p = 2).

# lembrete: exercícios 5 e 6 do script!
36 / 61

Regressão Linear Múltipla

Regressão Linear Simples

y = \beta_0 + \beta_1x

Exemplo:

dist = \beta_0 + \beta_1speed

### No R:
lm(dist ~ speed, data=cars)

Ver ISL página 61 (Simple Linear Regression).

37 / 61

Regressão Linear Múltipla

Regressão Linear Múltipla

y = \beta_0 + \beta_1x_1 + \dots + \beta_px_p

Exemplo:

mpg = \beta_0 + \beta_1wt + \beta_2disp

### No R:
lm(mpg ~ wt + disp, data=mtcars)
38 / 61

Regressão Linear Múltipla

Regressão Linear Simples

y = \beta_0 + \beta_1x

Regressão Linear Múltipla

y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots + \beta_px_p

# ajuste de uma regressão linear múltipla no R
modelo_boston <- lm(medv ~ lstat + age, data = Boston)
summary(modelo_boston)
# Estimate Std.Error t value Pr(>|t|)
# (Intercept) 33.22 0.73 45.4 < 2e-16 ***
# lstat -1.03 0.04 -21.4 < 2e-16 ***
# age 0.03 0.01 2.8 0.00491 **
39 / 61

Preditores Categóricos

Preditor com apenas 2 categorias

Tamanho das nadadeiras de pinguins são diferentes entre os sexos?

summary(lm(flipper_length_mm ~ sex, data = penguins))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 197.364 1.057 186.792 < 2e-16 ***
# sexmale 7.142 1.488 4.801 2.39e-06 ***
40 / 61

Preditores Categóricos

Preditor com apenas 2 categorias

Tamanho das nadadeiras de pinguins são diferentes entre os sexos?

y_i = \beta_0 + \beta_1x_i \space\space\space\space\space\space \text{em que}\space\space\space\space\space\space x_i = \Bigg\{\begin{array}{ll}1&\text{se a i-ésimo animal for }\texttt{female}\\ 0&\text{se a i-ésimo animal for } \texttt{male}\end{array}

# lembrete: exercícios 8, 9 e 10 do script!

Ver ISL página 84 (Predictors with Only Two Levels).

41 / 61

Preditores Categóricos

Preditor com 3 ou mais categorias

summary(lm(flipper_length_mm ~ species, data = penguins))
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 189.9536 0.5405 351.454 < 2e-16 ***
# speciesChinstrap 5.8699 0.9699 6.052 3.79e-09 ***
# speciesGentoo 27.2333 0.8067 33.760 < 2e-16 ***

Modelo

y_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i}

Em que

x_{1i} = \Bigg \{ \begin{array}{ll} 1 & \text{se for }\texttt{Chinstrap}\\0&\text{caso contrário}\end{array}

x_{2i} = \Bigg \{ \begin{array}{ll} 1 & \text{se for }\texttt{Gentoo}\\0&\text{caso contrário}\end{array}

42 / 61

Preditores Categóricos

Preditor com 3 ou mais categorias

"One hot enconding" ou "Dummies" ou "Indicadores".

species (Intercept) speciesChinstrap speciesGentoo
Adelie 1 0 0
Chinstrap 1 1 0
Gentoo 1 0 1
Adelie 1 0 0
Gentoo 1 0 1
Adelie 1 0 0
Adelie 1 0 0
Chinstrap 1 1 0
43 / 61

Preditores Categóricos

Preditor com 3 ou mais categorias

Interpretação dos parâmetros:

y_{i} = \left\{ \begin{array}{ll} \beta_0 & \text{se for }\texttt{Adelie}\\ \beta_0 + \beta_1&\text{se for } \texttt{Chinstrap}\\ \beta_0 + \beta_2&\text{se for } \texttt{Gentoo}\end{array}\right.

# interprete cada um dos três parâmetros individualmente.
# lembrete: exercício 11 do script!
44 / 61

Transformações Não Lineares dos Preditores

Exemplo: log

Modelo real: y = 10 + 0.5log(x)

Modelo proposto: \small y = \beta_0 + \beta_1log(x)

Outras transformações comuns: raíz quadrada, polinômios, Box-Cox, ...

# lembrete: exercício 13 do script!
45 / 61

Transformações Não Lineares dos Preditores

Exemplo: log

Modelo real: y = 10 + 0.5log(x)

Modelo proposto: \small y = \beta_0 + \beta_1log(x)

Outras transformações comuns: raíz quadrada, polinômios, Box-Cox, ...

# lembrete: exercício 13 do script!
46 / 61

Transformações Não Lineares dos Preditores

Gráfico de Resíduos

47 / 61

Transformações Não Lineares dos Preditores

Exemplo: Regressão Polinomial

Modelo real: y = 500 + 0.4(x-10)^3

Modelo proposto: y = \beta_0 + \beta_1x + \beta_2x^2 + \beta_3x^3

# lembrete: exercício 14 do script!
48 / 61

Transformações Não Lineares dos Preditores

Exemplo: Regressão Polinomial

Modelo real: y = 500 + 0.4(x-10)^3

Modelo proposto: y = \beta_0 + \beta_1x + \beta_2x^2 + \beta_3x^3

y x x2 x3
456.5 5.3 28.2 149.7
492.5 7.4 55.4 412.2
548.4 11.5 131.3 1503.9
758.7 18.2 329.9 5993.0
444.7 4.0 16.3 65.6
748.3 18.0 322.8 5800.8
820.5 18.9 357.0 6744.3
# lembrete: exercício 14 do script!
49 / 61

Transformações Não Lineares dos Preditores

Gráfico de Resíduos

50 / 61

Interações

Interação entre duas variáveis explicativas: Species e Sepal.Length

51 / 61

Interações

Modelo proposto (Matemático): Seja y = Sepal.Width e x = Sepal.Length,

\small \begin{array}{l} y = \beta_0 + \beta_1x\end{array}

Modelo proposto (em R): Sepal.Width ~ Sepal.Length

# lembrete: exercícios 14 ao 17 do script!
52 / 61

Interações

Modelo proposto (Matemático): Seja y = Sepal.Width e x = Sepal.Length,

\small \begin{array}{l} y = \beta_0 + \beta_1x + \beta_2I_{versicolor} + \beta_3I_{virginica}\end{array}

Modelo proposto (em R): Sepal.Width ~ Sepal.Length + Species

# lembrete: exercícios 14 ao 17 do script!
53 / 61

Interações

Modelo proposto (Matemático): Seja y = Sepal.Width e x = Sepal.Length,

\small \begin{array}{l} y = \beta_0 + \beta_1x + \beta_2I_{versicolor} + \beta_3I_{virginica} + \beta_4\color{red}{xI_{versicolor}} + \beta_5\color{red}{xI_{virginica}}\end{array}

Modelo proposto (em R): Sepal.Width ~ Sepal.Length * Species

# lembrete: exercícios 14 ao 17 do script!
54 / 61

Heterocedasticidade

55 / 61

Heterocedasticidade

56 / 61

Heterocedasticidade

57 / 61

Heterocedasticidade

Problema

  • O estimador \hat{\sigma}_{\beta_1} = \sqrt{\frac{EQM}{\sum(x_i - \bar{x})^2}} deixa de ter as melhores propriedades. Poderíamos ter conclusões estranhas para \beta_1.

Diagnóstico

  • Visualização dos resíduos;
  • Testes formais (Breuch-Pagan, White, etc)

Tratamentos

  • Transformações na variável resposta. log(y), sqrt(y), 1/y, etc;
  • Mínimos Quadrados Ponderados/Generalizados
58 / 61

Multicolinearidade

Modelo 1: sem colineares

term estimate std.error statistic p.value
(Intercept) -173.41 43.83 -3.96 0
Limit 0.17 0.01 34.50 0
Age -2.29 0.67 -3.41 0

Modelo 2: com colineares

term estimate std.error statistic p.value
(Intercept) -377.54 45.25 -8.34 0.00
Limit 0.02 0.06 0.38 0.70
Rating 2.20 0.95 2.31 0.02

Problema: Instabilidade numérica, desvios padrão inflados e interpretação comprometida.

Soluções: eliminar uma das variáveis muito correlacionadas ou Consultar o VIF (Variance Inflation Factor)

59 / 61

Multicolinearidade

VIF (Variance Inflation Factor)

Detecta preditores que são combinações lineares de outros preditores.

Procedimento: Para cada preditor X_j,

1) Ajusta regressão linear com as demais: lm(X_j ~ X_1 + ... + X_p).

2) Calcula-se o R-quadrado dessa regressão e aplica a fórmula abaixo

\small VIF(\hat{\beta}_j) = \frac{1}{1 - R^2_{X_j|X_{-j}}}

3) Remova o preditor se VIF maior que 5 (regra de bolso).

# lembrete: exercícios 18 do script!

Ver ISL página 101.

60 / 61

Limitações da Regressão Linear

  • Variável resposta Não Normal
  • Variável resposta Positiva
  • Variável resposta com mais de duas categorias
  • Relação funcional não linear entre X e Y
  • Muitas variáveis disponíveis
  • Muitas interações para testar entre as preditoras
61 / 61

Agenda

  • O que é e quando usar
  • Parâmetro vs estimador
  • Teste de Hipóteses e valor-p
  • Interpretação dos parâmetros
  • Regressão Linear Múltipla
  • Preditores Categóricos
  • Transformações Não Lineares dos Preditores
2 / 61
Paused

Help

Keyboard shortcuts

, , Pg Up, k Go to previous slide
, , Pg Dn, Space, j Go to next slide
Home Go to first slide
End Go to last slide
Number + Return Go to specific slide
b / m / f Toggle blackout / mirrored / fullscreen mode
c Clone slideshow
p Toggle presenter mode
t Restart the presentation timer
?, h Toggle this help
Esc Back to slideshow