Análise de Dados com o Software R:
Métodos Estatísticos, Computacionais e Econométricos

Prof. Adriano Azevedo Filho (azevedofilho@usp.br)

Análise de regressão I: introdução à regressão linear simples e múltipla

sumário geral | anterior | próximo

Conteúdo do Módulo

 1 - Introdução à análise de regressão
 
 2 - Diagrama de dispersão e superfície de regressão múltipla (2 variáveis)

 3 - Estimativa de parâmetros de regressão pelo método dos quadrados mínimos
 
 4 - Fórmulas e extração de dados de modelos no R 
 
 5 - Modelos lineares vs. modelos não-lineares
 
 6 - Uso da função predict em gráficos de regressões
 
 7 - Uso de variáveis "dummy" ou binárias em regressões com variáveis qualitativas
 
 8 - Intervalo de confiança e intervalo de previsão para a regressão
 
 9 - Notações usadas em modelos de regressão
 
10 - Análise dos resíduos da regressão e diagnóstico de problemas
 
11 - Estudo de caso: mãe fumante vs. peso do bebe / questões sobre o módulo
  

1 - Introdução à análise de regressão

Este tópico tem o objetivo de introduzir notação e procedimentos utilizados na análise de regressão, desenvolvidos no contexto da avaliação de algumas hipóteses. Alguns conjuntos de dados são analisados de forma exploratória, com o intuito principal de exemplificar as técnicas fundamentais usadas em análise de regressão. Não há intenção de se fazer uma análise rigorosa ou se especificar o melhor modelo para cada situação específica, algo que pode ser um exercício interessante para muitos.

1.1 Dados dos resultados dos 100 m nos jogos olímpicos

Carregue o arquivo dados100m.csv no data frame olimp. Esse arquivo contém os resultados, em segundos, da final dos 100 m nas olimpíadas, desde 1900 até 2012, para homens e mulheres, incluindo informações sobre o primeiro, segundo e terceiro colocados. As mulheres só começaram a participar dessa prova (e das provas de atletismo de um modo geral) em 1928.

####
## Arquivo dados100m.csv - resultados dos 100 m nos jogos olimpicos (1900-2012)
####
## ano   - ano dos jogos
## cid   - cidade onde ocorreram os jogos
## elev  - elevação em relação ao nível do mar em m da cidade dos jogos
## mod   - modalidade da prova (100m, não varia no data frame)
## sexo  - sexo do atleta (M - masculino, F - feminino)
## nome  - nome do atleta
## ord   - ordem de classificação na prova (1, 2, 3)
## pais  - código do pais do atleta (3 letras)
## t     - tempo em segundos para completar a prova
####

### Leitura do arquivo
olimp<-read.csv2("http://ihbs.com.br/html/dados100m.csv")
head(olimp) ## listando as primeiras linhas do data frame

##    ano         cid elev  mod sexo               nome ord pais    t
## 1 1928   Amsterdam    2 100m    F Elizabeth Robinson   1  USA 12.2
## 2 1928   Amsterdam    2 100m    F    Fanny Rosenfeld   2  CAN 12.3
## 3 1928   Amsterdam    2 100m    F        Ethel Smith   3  CAN 12.3
## 4 1932 Los Angeles  113 100m    F      S Walasiewicz   1  POL 11.9
## 5 1932 Los Angeles  113 100m    F       Hilda Strike   2  CAN 11.9
## 6 1932 Los Angeles  113 100m    F  Wilhelmina Bremen   3  USA 12.0

Faça um diagrama de dispersão para visualizar os dados dos homens e mulheres, ao longo dos anos, colorindo os pontos com cores rosa (mulheres) e azul (homens).

colvec<-ifelse(olimp$sexo=="M","blue","deeppink")
plot(olimp$ano,olimp$t,col=colvec, xlab="ano",ylab="t")
title(main="Resultados dos 100 m (1900-2012)")
grid()

Identificação de observações no gráfico

Um exame rápido do diagrama de dispersão mostra um dado surpreendente, o record olimpico das mulheres foi estabelecido já há muito tempo e ainda não foi superado. Como encontrar a observação relativa a esse ponto? (é algo útil para identificação de potenciais problemas)

Com o gráfico aberto execute a função identify do R, com as variáveis que definem os eixos \(x\) e \(y\) do gráfico como argumentos:

## execute o comando abaixo com o gráfico aberto
valor<-identify(olimp$ano,olimp$t) 
## agora posicione o ponteiro da tela no ponto desejado
## e pressione o botão esquerdo do mouse (pode clicar em vários pontos)
## 
## para parar: pressione o botão direito do mouse e escolhar "parar" 
##
## os números das linhas no data frame serão colocadas na variável valor 
## e "plotadas" no gráfico ao lado do ponto

Como exercício, identifique a observação que corresponde ao record feminino ainda não superado (ao menos até a olimpiada de 2012). Use o procedimento acima e veja a posição da observação (esse procedimento indicará que a observação está n linha 40) e mostre o conteúdo dessa linha do data frame:

valor

## [1] 40

olimp[valor,] ## todas as colunas da linha 40

##     ano   cid elev  mod sexo              nome ord pais     t
## 40 1988 Seoul   33 100m    F F Griffith-Joyner   1  USA 10.54

As informações apresentadas indicam que esse foi o resultado de Florence Griffith-Joyner, conhecida como Flo-Jo, obtido nos jogos de Seul em 1988, que corresponde ao record mundial até hoje não superado. É também de Flo-Jo o record mundial dos 200m, obtidos nos mesmos jogos e também não superado.

Esses resultados de Flo-Jo são controversos, por estarem muito abaixo dos que seriam razoáveis para atletas femininas de alto nível da época e até hoje. Os exames anti-doping, até 1988, não eram tão rigorosos (com relação ao uso de anabolizantes e outras substâncias). Esses testes não eram aplicados fora dos períodos de competição, algo que incentivava o uso durante esse período sem deixar vestígios das substâncias ilícitas utilizadas durante os testes anti-doping realizados nas épocas de competição. A partir de 1988, até motivado por esse resultado extraordinário de Flo-Jo, os testes anti-doping passaram a ser mais rigorosos, cobrindo o período fora de competições, e isso parece ter tido um efeito nos resultados femininos a partir de então (veja os dados). Esse possível efeito será examinado em um próximo tópico deste módulo.

1.2 Questões que queremos examinar

Discuta as questões descritas a seguir no contexto das provas masculinas

Pergunta 1 - Há alguma evidência forte de que os tempos na prova de 100 m estão caindo ao longo do tempo? Ainda que isso possa parecer óbvio, há algum argumento forte para defender esse ponto de vista? Quanto em média, por ano, o tempo da prova dos 100 m está caindo?

Pergunta 2 - Há evidência de algum impacto da elevação (altitude) da cidade dos jogos nos tempos da prova? Avalie formalmente. Estime, se houver, o efeito da altitude no tempo da prova.

Como a análise é focada inicialmente somente na prova masculina, vamos criar um data frame específico, contendo somente essas observações, denominado olimpm.

## olimpm: data frame só com observações dos homens, todas as variáveis
olimpm<-olimp[olimp$sexo=="M",]  
plot(olimpm$ano,olimpm$t,xlab="ano",ylab="t",col="blue")
title(main="Prova dos 100 m (Homens, 1900-2012)")

1.3 Modelo estatístico e notação

A resposta a questões e testes de hipótese através métodos estatísticos frequentemente se utiliza de:

• modelos estatísticos paramétricos: modelos que envolvem parâmetros que são estimados através de procedimentos específicos. Uma alternativa a esses modelos, são os chamados modelos estatísticos não-paramétricos, que são assunto de desenvolvimentos mais avançados em estatística, não tratados aqui.

Noções importantes nos modelos:

• variável resposta ou variável dependente: representada por \(Y\) ou \(y\), indica a variável que desejamos “explicar”. No caso desejamos explicar a variável t, que representa o tempo no data frame olimpm.

• variáveis explicativas ou explanatórias: representadas por \(x_1\), \(x_2\), \(\ldots\), \(x_m\), indicam as variáveis que estão sendo consideradas para “explicar” a variável resposta. No caso podemos usar variáveis como ano, elev, sexo e outras, no data frame olimpm, como variáveis explicativas.

Modelo linear básico na análise de regressão (definição e premissas)

É um dos modelos mais usados em estatística e muito associado a uma técnica que se chama análise de regressão, que é um dos pilares técnicos da econometria.

Esse modelo, na sua forma mais elementar, é explicitado por

• \(Y_i=b_0+b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im} + \epsilon_i,\;\;\;\text{com}\;\;\;\epsilon_i\sim \text{N}(0,\sigma)\)

• \(i=1,\ldots,n\) representam as observações e no modelo mais elementar assume-se que \(\text{cov}(\epsilon_i,\epsilon_j)=0\), \(\forall i\neq j\), ou (uma premissa mais forte), \(\epsilon_1\), \(\epsilon_2\),\(\ldots\),\(\epsilon_n\) são independentes. Pelas propriedades da distribuição Normal, \(Y_i\) terá distribuição Normal. As variáveis explicativas \(x_i\)s são usualmente assumidas fixas ou, como premissa equivalente (um pouco mais geral), assume-se diretamente que \(\text{cov}(x_j,\epsilon_i)=0\), \(\forall j,i\). Assume-se, também, que não há erros nas medidas das observações.

Quando \(m=1\) temos a regressão linear simples e quando \(m>1\) temos a regressão linear múltipla.

• regressão linear simples (\(m=1\)): * \(Y_i=b_0+b_1 x_{i1}+ \epsilon_i\)

• regressão linear multipla (\(m>1\)): * \(Y_i=b_0+ b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im} + \epsilon_i\)

• no modelo: \(b_0\), \(b_1\), \(\ldots\), \(b_m\) e \(\sigma\) são os parâmetros que devem ser estimados por procedimentos específicos.

• as estimativas desses parâmetros: são representadas por \(\hat b_0\), \(\hat b_1\), \(\ldots\), \(\hat b_m\) e \(\hat \sigma\)

Obs: a relação que os estimadores dos parâmetros \(\hat b_j\)s tem com os (verdadeiros) parâmetros \(b_j\)s é a mesma que um outro estimador mais conhecido, a média amostral \(\bar X\), tem com \(\mu\) a média teórica ou esperança. Nos modelos estatísticos sempre haverá a possibilidade de estimar intervalos de confiança para os \(b_j\)s, da mesma forma que já obtivemos os intervalos de confiança para \(\mu\), a partir dos dados amostrais.

Essa notação utilizada acima assume que \(\epsilon_i\)s são variáveis aleatórias. Os \(Y_i\)s também são variáveis aleatórias, por serem uma função dos \(\epsilon_i\)s e da parte não aleatória do modelo dependente dos \(b_j\)s como \(x_i\)s, que são assumidos constantes.

Temos que, pelas propriedades da esperança matemática e da variância teórica, pela premissa de termos os valores de \(x_{ij}\), fixos (não aleatórios), que:

• \(E[Y_i]=E[b_0+b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im}+\epsilon_i]\)

• \(E[Y_i]=b_0+b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im}\) dado que \(E[\epsilon_i]=0\) (pela premissa)

• \(V[Y_i]=V[b_0+b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im}+\epsilon_i]\)

• \(V[Y_i]=V[\epsilon_i]=\sigma^2\) dado que \(V[\epsilon_i]=\sigma^2\).

Importante: Em muitos desenvolvimentos \(Y_i\) pode ser representado por letra minúscula, \(y_i\), e isso pode ser confuso para iniciantes, pois não deixa claro o status de variáveis aleatórias. É uma tecnicalidade importante para estudos mais avançados. Nos desenvolvimentos em outros tópicos, usaremos \(y_i\) para simplificar a notação.

Definições mais técnicas desse modelo

Para um bom entendimento desse modelo, é importante interpretá-lo, mais técnicamente por (na notação a seguir o índice \(i\) associado à observação foi suprimido, por ser desnescessário nesse contexto)

• \(Y|x_{1},x_{2},\ldots,x_{n} \sim \text{N}(b_0+b_1 x_{1}+\ldots+b_m x_m, \sigma)\)

Ou seja, \(Y|x_1,x_2,\ldots,x_n\) é a distribuição condicional de \(Y\), a variável resposta, dado os valores de \(x_1\), \(x_2\),\(\ldots\), \(x_m\), as variáveis explicativas, e os parâmetros \(b_0\), \(b_1\), \(b_2\),\(\ldots\), \(b_m\) onde

• \(E[Y|x_1,x_2,\ldots,x_n]=b_0+b_1 x_1+\ldots+b_m x_m\) (curva de regressão)
• \(V[Y|x_1,x_2,\ldots,x_n]=\sigma^2\) (variância condicional)

Pela premissa da distribuição de \(\epsilon\sim\text{N}(0,\sigma)\), associada a propriedades da distribuição Normal, podemos concluir que

• \(Y|x_1,x_2,\ldots,x_n\sim \text{N}(b_0+b_1 x_{1}+\ldots+b_m x_m, \sigma)\)

Note que nos modelos lineares desse tipo, o parâmetro \(b_k\) indica a derivada da curva de regressão com relação ao parâmetro \(b_k\), ou seja

• \(b_k=\displaystyle \frac{\partial\, E[Y|x_{1},x_2,\ldots,x_m]}{\partial\, x_k}\;\;\;\;\;\;\;\;\; k=1,\ldots,m\)

• efeito marginal: O parâmetro \(b_k\), \((k=1,\dots,m)\), indica o efeito marginal da variável \(x_k\) associada a ele, ou seja, aproximadamente, quanto varia \(E[Y|x_{1},x_2,\ldots,x_m)]\) para uma variação positiva de 1 unidade na variável associada (no contexto infinitesimal).

1.4 Modelo estatístico 1: regressão linear simples

Motivado pelo próprio comportamento do tempo apresentado no gráfico, vamos considerar um primeiro modelo estatístico, caracterizado por uma regressão linear simples, para representar o fenômeno. Esse modelo é uma caso particular do modelo geral que especificamos no último tópico:

• Modelo 1: \(Y_i=b_0+b_1 x_{i1} + \epsilon_i\) (regressão linear simples)

• onde \(\epsilon_i\sim \text{N}(0,\sigma)\) com \(\text{cov}(\epsilon_i,\epsilon_j)=0,\;\; \forall i\neq j\) e outras premissas usuais (veja na definição geral)

Nesse modelo a variável resposta \(Y_i\) na notação genérica é representada por \(t_i\) (tempo nos 100 m na observação \(i\)) e a variável explicativa \(x_{i1}\) é representada por ano_i (ano dos jogos na observação \(i\)) no data frame olimpm.

Considere agora a hipótese:

• \(\def\H{{\mathbb H}}\H_0\): “os resultados estão dos tempos são constantes ao longo dos anos”

esse hipótese, no contexto desse modelo, no contexto estatístico, seria definida por

• \(\H_0:\; b_1 = 0\)

indicando que não há efeito dos anos (que caracterizam uma evolução técnica) no tempo em segundos para completar a prova.

Na realidade estamos mais interessados nas hipóteses alternativas \(\H_A:\; b_1 \neq 0\) ou, mais especificamente, \(\H_A:\; b_1 < 0\).

Estratégia usual: verificação da possibilidade de rejeição de \(\H_0\). Se houver evidência para rejeição de \(\H_0\), automaticamente \(\H_A\) ficará fortalecida. Parece estranho mas é assim que funciona a “lógica” usual de testes de parâmetros.

A função principal do R para estimação modelos lineares chama-se lm (iniciais de linear model em inglês). Vamos, a seguir, mostrar como estimar os parâmetros desse modelo usando essa função, obtendo com isso evidências úteis para responder às questões que estamos interessados.

## para estimar os resultados e colocá-los na variável modelo1:
modelo1<-lm(t~ano,data=olimpm) ## t = b0+ b1 ano + E
## para ver os resultados
summary(modelo1)

A figura indica os resultados mais importantes em vermelho. A forma de apresentação é típica de programas estatísticos. Nesses resultados temos:

• Estimativas para \(b_0\) e \(b_1\): são indicadas por \(\hat b_0\) e \(\hat b_1\) que tem valores, respectivamente, de 32,538966 e -0,011335. Observe que no início da linha com as estimativas (nos resultados do R) é apresentado nome da variável ao qual o parâmetro está relacionado, ou “Intercept”, para indicar a estimativa para \(b_0\) o termo constante da regressão.

• Estimativa para \(\sigma\): é indicada por \(\hat \sigma\) com valor 0,144.

• \(R^2\): medida (imperfeita) de qualidade de ajuste da regressão (0 é o pior ajuste, 1 é o melhor). Na regressão linear simples, coincide com o quadrado do coeficiente de correlação entre \(Y\) e \(x\), ou seja, as variáveis \(t\) e \(ano\) no data frame olimpm. Um dos problemas dele é que sempre aumenta se aumentarmos adicionarmos novas variáveis à regressão. O valor obtido 0,881, para \(R^2\), mostra um ajuste razoavelmente bom, algo que é comprovado pela análise visual.

• \(R^2\) ajustado: como o nome diz, é uma medida, baseada no \(R^2\) com um ajuste que considera o número de variáveis explicativas na regressão. É um indicador um pouco melhor da qualidade de ajuste que o próprio \(R^2\), para comparação de modelos alternativos.

• valores-p associados às hipóteses \(\H_0\): no final da linha onde aparece cada estimativa \(b_i\), aparece o valor-p associado à hipótese \(\H_0: b_i=0\). No caso, como os valores são próximos de zero, a evidência sugere a rejeição de \(\H_0: b_0=0\) e \(\H_0: b_1=0\). Ao final, aparece um outro valor-p, associado a estatística F, que testa a hipótese \(\H_0: b_1=0,\ldots,b_m=0\) (observe que \(b_0\) não está na lista). Como \(m=1\), esse teste coincide, em resultado, com o valor-p mostrado na linha da estimativa de \(b_1\) (embora sejam testes distintos).

Apresentação gráfica da curva de regressão e intervalos de confiança para parâmetros

A curva de regressão (teórica) nesse caso, na notação genérica será dada por

• \(E[Y|x_{1}] = b_0+b_1 x_{1}\) cuja estimativa será dada por \(\hat E(Y|x_{1})=\hat b_0 + \hat b_1 x_{1}\), numa situação em que \(Y\) corresponde a variável \(t\) e \(x_1\) à variável \(\text{ano}\)

• \(\hat y\) é comumente usado como símbolo para indicar \(\hat E(Y|x_{1},x_2,\ldots,x_m)\) para simplificar a notação. Contudo, é muito importante não esquecer o significado verdadeiro de \(\hat y\) como a média condicional da variável resposta \(Y\), dado o valor da variável explicativa (ou variáveis explicativas, no caso geral).

Assim, com essas considerações e notações, poderiamos definir:

• \(\hat E[Y|x_{1}]=\hat y=\hat b_0 + \hat b_1\,x_1\) ou já substituindo as estimativas \(\hat y=32{,}538966-0{,}011335\,x_1\).

Essa curva de regressão pode ser apresentada no diagrama de dispersão que contém os pontos, com os seguintes comandos:

modelo1<-lm(t~ano,data=olimpm) ## t = b0+ b1 ano
### Próximo comando assume pré-existência do diagrama de dispersão
abline(modelo1)

• intervalos de confiança para os parâmetros: podemos obter o intervalo de confiança para os “verdadeiros” valores dos \(b_0\) e \(b_1\), no contexto desse modelo, através de

## intervalo de confiança a 95% para os parâmetros
confint(modelo1,level=0.95) 

##                   2.5 %      97.5 %
## (Intercept) 30.67841874 34.39951415
## ano         -0.01228458 -0.01038501

Os valores obtidos indicam que há 95% de probabilidade dos “verdadeiros” parâmetros estarem contidos nesse intervalo. No caso de \(b_1\) por exemplo, há 95% de probabilidade do verdadeiro valor estar dentro no intervalo \([-0{,}01228,\;-0{,}01039]\) o que dá uma boa idéia da incerteza associada à estimativa. A interpretação é muito similar à que fizemos anteriormente para o intervalo de confiança da média teórica ou esperança.

Com esses resultados, temos todos os elementos para responder à primeira pergunta colocada no início do tópico 1.2:

• Pergunta 1: Há alguma evidência forte de que os tempos na prova de 100 m estão caindo ao longo do tempo? Ainda que isso possa parecer óbvio, há algum argumento forte para defender esse ponto de vista? Quanto em média, por ano, o tempo da prova dos 100 m está caindo nos últimos anos?

• Resposta: Tanto a evidência relativa ao valor-p da hipótese \(\H_0:\;b_1=0\), próximo de zero, sugerindo rejeição de \(\H_0\), quanto o intervalo de confiança a 95% sugerem os tempos estão caindo de forma estatisticamente significativa. O valor estimado para \(b_1\) está entre \([-0{,}01228,\;-0{,}01039]\). Esse resultado estimado sugere que o tempo nos 100 m tendeu a cair, em média por ano, no período, algo dentro do intervalo \([0{,}01039,\; 0{,}01228]\) segundos com 95% de probabilidade.

1.5 Modelo estatístico 2: regressão linear múltipla

Para responder à pergunta 2 colocada no início do tópico 6.2, envolvendo um possível efeito da elevação ou altitude nos tempos da prova final de 100 m da olimpíada, precisamos utilizar um modelo um pouco mais complexo, envolvendo uma regressão linear múltipla com 2 variáveis explicativas:

• Modelo 2: \(y_i=b_0+b_1 x_{i1} + b_2 x_{i2} +\epsilon_i\) (regressão linear múltipla)

• onde \(\epsilon_i\sim \text{N}(0,\sigma)\) com \(\text{cov}(\epsilon_i,\epsilon_j)=0,\;\; \forall i\neq j\) e outras premissas usuais (veja na definição geral)

Nesse modelo, substituímos a notação \(Y_i\) por \(y_i\), para simplificar, mas com o entendimento de que \(y_i\) é uma variável aleatória, utilizada como a variável resposta, que no conjunto de dados é representada por \(t_i\) (tempo nos 100 m na observação \(i\)). A variável explicativa \(x_{i1}\) é representada por ano_i (ano dos jogos na observação \(i\)). A novidade é a variável \(x_{i2}\) que representa elev_i (elevação em m da cidade da observação \(i\)), do data frame olimpm.

A estimativa dos parâmetros \(b_0\), \(b_1\) e \(b_2\) pode ser feita com muita facilidade através da função lm através de:

## para estimar os resultados e colocá-los na variável modelo1:
modelo2<-lm(t~ano+elev,data=olimpm) ## t = b0+ b1 ano + b2 elev + E
## para ver os resultados
summary(modelo2)

## 
## Call:
## lm(formula = t ~ ano + elev, data = olimpm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34118 -0.07482 -0.00952  0.09358  0.28724 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.228e+01  8.940e-01  36.106  < 2e-16 ***
## ano         -1.119e-02  4.567e-04 -24.510  < 2e-16 ***
## elev        -1.077e-04  3.628e-05  -2.968  0.00403 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1369 on 75 degrees of freedom
## Multiple R-squared:  0.8939, Adjusted R-squared:  0.8911 
## F-statistic: 315.9 on 2 and 75 DF,  p-value: < 2.2e-16

O modelo 2 parece melhor que o modelo 1. O modelo 2 não só apresentou um \(R^2\) ajustado maior (de 0,88 no modelo 1 subiu para 0,891 no modelo 2), como o valor-p associado ao parâmetro \(b_2\) relacionado ao efeito da variável elev é 0,004, o que sugere a rejeição da hipótese \(\H_0: b_2=0\). Isso indica um efeito estatisticamente significativo da variável elev, na explicação do tempo nos 100 m. O valor estimado de \(b_2\) é \(\hat b_2=-1,08e-4\) ou \(\hat b_2=-0,000108\).

O intervalo de confiança para os parâmetros pode ser obtido por

## intervalo de confiança a 95% para os parâmetros
confint(modelo2)

##                    2.5 %        97.5 %
## (Intercept) 30.496169419  3.405789e+01
## ano         -0.012102085 -1.028268e-02
## elev        -0.000179932 -3.538815e-05

para \(b_2\) o intervalo de confiança estimado, a 95%, foi \([-0{,}0000354,\;-0{,}0001799]\). Esse seria o efeito estimado para cada 1 metro adicional de elevação. Para 1000 m esse efeito seria 1000 vezes maior, situando-se entre \([-0{,}0354,\;-0{,}1799]\) segundos com 95% de probabilidade. Ou seja, parece ser vantajoso correr os 100 m em altitudes mais elevadas, em razão do ar ser mais rarefeito, oferecendo menos resistência aos atletas.

Esse modelo 2, que mostrou-se superior ao modelo 1 pelos critérios utilizados, também indica uma nova estimativa para o intervalo de confiança para o \(b_1\), entre \([-0{,}0121,\;-0{,}0108]\) com 95% de probabiliade, muito similar ao obtido no modelo 1.

Assim sendo, podemos oferecer uma possível responder à pergunta 2 colocada no tópico 6.2:

• Pergunta 2: Há evidência de algum impacto da elevação (altitude) da cidade dos jogos nos tempos da prova? Avalie formalmente. Estime, se houver, o efeito da altitude no tempo da prova.

• Resposta: Sim, o resultado da análise sugere que o efeito estimado da elevação no tempo dos 100 m é estatisticamente significativo (valor-p=\(0{,}004\)) e negativo na medida que aumenta a elevação. O efeito esperado no tempo dos 100 m, para uma variação de 1000 m na elevação foi estimado dentro do intervalo de confiança \([-0{,}0354,\;-0{,}1799]\) segundos, com 95% de probabilidade. Num trabalho realizado por Mureika 2008 o efeito reportado é de -0,04 segundos para cada 1000 m, que está dentro do intervalo de confiança estimado, algo que traz suporte adicional para a noção de que o aumento da elevação pode trazer um benefício à velocidade dos atletas nos 100 m.

Referências disponíveis para o uso do R em análise de regressão

As duas referências a seguir são relativamente elementares e disponíveis em pdf, e podendo ajudar no entendimento do uso do R na análise de regressão e econometria introdutória.

• FARAWAY, J. Practical Regression and Anova using R. [nota: versão preliminar do livro Linear Models; os dados e scripts podem se baixados do site do livro]

• FARNSWORTH, G. Econometrics in R (nota: nível básico, focado em técnicas econométricas]

1.6 Problemas sugeridos (já resolvidos na última tarefa)

1. Responda as 2 perguntas colocadas no tópico 1.2 para os resultados das mulheres nos 100 m descritos no data frame olimp. Para facilitar siga o mesmo procedimento utilizado para o caso dos homens, criando um data frame olimpm somente com as informações das mulheres. As questões na tarefa dessa semana, com relação a este tópico, serão relacionadas às estimativas necessárias para resposta a essas 2 perguntas.

2 - Diagrama de dispersão e superfície de regressão múltipla (2 variáveis)

Neste tópico curto apresentaremos recursos do R para a construção de diagramas de dispersão e superfícies de regressão no contexto tridimensional, envolvendo 3 variáveis de interesse. Há várias alternativas de fazer isso no R. A forma utilizada, usando a função scatter3d é uma das mais fáceis e poderosas, requerendo a instalação dos packages rgl e car (para opções mais avançadas do comando scatter3d outros packages podem ser exigidos). Outras possibilidade são descritas aqui no contexto bidimensional e tridimensional.

Os gráficos produzidos pelo comando scatter3d são interativos e podem ser rotacionados com o mouse, facilitando a visualização. As figuras apresentadas junto dos comandos ilustram algumas cenas que poderão ser visualizadas no modo interativo.

Modelo 2 do data frame olimpm (resultados dos homens nos 100 m)

No tópico anterior, consideramos, no modelo 2, que a variável resposta \(t\) (tempo nos 100 m) era explicada pelas variáveis ano (ano dos jogos) e elev (elevação ou altitude da cidade dos jogos) no modelo estatístico:

• \(t_i = b_0+b_1 \text{ano}_i + b_2\,\text{elev}_i + \epsilon_i\)
• com \(\epsilon_i\sim\text{N}(0,\sigma)\), \(\text{cov}(\epsilon_i,\epsilon_j)=0\), para \(i\neq j\), e outras premissas usuais (veja na definição geral)

A equação que define o modelo caracterizará um plano e não mais uma simples curva, algo que ocorreria no caso de uma regressão linear simples. O diagrama de dispersão e o essa superfície de resposta caracterizada pela regressão poderá ser obtida por:

### Instale os packages rgl e car antes executar
require(rgl) 
require(car) 
scatter3d(t ~ ano + elev, data=olimpm)

As figuras evidenciam que há uma situação em que os jogos foram realizados numa cidade com elevação muito acima do padrão geral. Um exame dos dados revelará que foram os jogos olímpicos realizados na cidade do México. De um modo geral, os jogos tendem a ser realizados em cidades de baixa elevação (ou altitude), em muitos casos próxima ao nível do mar.

Data frame Prestige (prestigio ocupacional em função de renda e educação)

O data frame Prestige, examinado a seguir, é relativo a um trabalho antigo e influente, que tenta examinar a relação entre o “prestígio ocupacional”, a renda e a educação em profissões distintas, agrupadas em 3 grupos (profissional, técnico administrativo, não técnico). Uma descrição do problema estudado é apresentado aqui.

O data frame Prestige é automaticamente carregado quando o package car é carregado. Com esse data frame são exploradas algumas outras opções do comando scatter3d (veja mais opções na documentação da função usando ?scatter3d)

require(rgl) 
require(car)
## Uma superfície de regressão para cada "tipo" de ocupação
scatter3d(prestige ~ income + education | type, data=Prestige)
## Diagrama de dispersão em 3 dimensões: raio é proporcional a proporção de mulheres 
scatter3d(prestige ~ income + education | type, radius=women/2, surface=FALSE, data=Prestige)
## mostrando as primeiras linhas do data frame Prestige
head(Prestige)  ## 
### 

##                     education income women prestige census type
## gov.administrators      13.11  12351 11.16     68.8   1113 prof
## general.managers        12.26  25879  4.02     69.1   1130 prof
## accountants             12.77   9271 15.70     63.4   1171 prof
## purchasing.officers     11.42   8865  9.11     56.8   1175 prof
## chemists                14.62   8403 11.68     73.5   2111 prof
## physicists              15.64  11030  5.13     77.6   2113 prof

3 - Estimação de parâmetros de regressão pelo método dos quadrados mínimos

Há vários métodos importantes para obtenção de estimadores em estatística. De longe, os mais utilizados são:

• Método dos quadrados mínimos: mais utilizado no contexto de modelos lineares usuais

• Método da máxima verossimilhança: necessário para modelos mais avançados

Outras possibilidades (também importantes) incluem: método dos momentos, método dos momentos generalizado (GMM) e método Bayesiano

Neste módulo cobriremos o método dos quadrados mínimos.

3.1 Método dos quadrados mínimos (ordinários)

Veja apresentação detalhada sobre o método dos quadrados mínimos nesses slides.

3.2 Obtendo as estimativas do comando lm pela implementação das fórmulas

Vamos agora verificar as fórmulas apresentadas nos slides, para estimativa dos parâmetros do modelo de regressão associado aos dados dos homens nos jogos olímpicos estimado anteriormente.

• modelo1: \(y_i=b_0+b_1 x_{i} + \epsilon_i\) (regressão linear simples)

• onde \(\epsilon_i\sim \text{N}(0,\sigma)\) com \(\text{cov}(\epsilon_i,\epsilon_j)=0,\;\; \forall i\neq j\) e outras premissas usuais (veja na definição geral).

Nesse modelo a variável resposta \(y_i\) na notação genérica é representada por \(t_i\) (tempo nos 100 m na observação \(i\)) e a variável explicativa \(x_{i}\) é representada por ano_i (ano dos jogos na observação \(i\)) no data frame olimpm.

Vimos nos slides que os estimadores de quadrados mínimos \(\hat b_0\) e \(\hat b_1\) de \(b_0\) e \(b_1\) são definidos por

• \(\hat b_1=\displaystyle \frac{\sum_{i=1}^n x_i y_i - n\, \bar x\,\bar y}{\sum_{i=1}^n x_i^2 - n\, \bar x^2}\)

• \(\hat b_0=\displaystyle \bar y - \hat b_1 \bar x\)

No caso, os \(y_i\)s representam as observações da variável t (tempo nos 100 m) e \(x_i\)s as observações da variável ano_i (ano dos jogos na observação \(i\)) no data frame olimpm.

Os valores obtidos pela operacionalização das fórmulas leva a:

x<-olimpm$ano
y<-olimpm$t
n<-length(y)
b1est<-(sum(x*y)-n*mean(x)*mean(y))/(sum(x^2)-n*mean(x)^2)
b0est<-mean(y)-b1est*mean(x)
cat("b0 estimado = ",b0est," b1 estimado = ",b1est,"\n")

## b0 estimado =  32.53897  b1 estimado =  -0.0113348

Podemos também obter os mesmo resultados utilizando a fórmula geral matricial, dada por * \(latex \def\X{{\mathbf X}} \def\y{{\mathbf y}} \def\vecbeta{{\boldsymbol{\beta}}} \hat \vecbeta=(\X^T\X)^{-1}\X^T\y\)

nessa notação, o vetor \(\hat \vecbeta\) conterá o vetor de parâmetros estimados \(\hat b_0\) e \(\hat b_1\). Veja os slides para definição da matriz \(\X\) (observe que essa matriz tem sua primeira coluna com o número 1 em todas as linhas, no caso mais usual):

A seguir a implementação computacional, usando os resultados dos últimos comandos.

## matrix X tem uma coluna com 1 e colunas com as variáveis explicativas
## cbind junta colunas em uma matriz ou data frame
X<-cbind(rep(1,length(x)),x) 
X[1:5,] ## mostrando 5 linhas da matriz X

##           x
## [1,] 1 1900
## [2,] 1 1900
## [3,] 1 1900
## [4,] 1 1904
## [5,] 1 1904

betaest<-solve(t(X)%*%X)%*%t(X)%*%y
betaest

##         [,1]
##   32.5389664
## x -0.0113348

Foram obtidos os mesmos resultados apresentados anteriormente para as estimativas de \(\hat b_0\) e \(\hat b_1\), através do comando lm.

4 - Fórmulas e extração de informações de modelos no R

Este tópico é um material de consulta que serve para explicar como os modelos devem ser definidos nas fórmulas usadas em funções como lm e outras similares (ex. glm que veremos no próximo módulo). Além disso, mostraremos como extrair informações dos modelos estimados, algo que pode ser conveniente em muitas situações.

4.1 Definição de fórmulas

As possibilidades são muitas. Apresentamos abaixo alguns casos comuns na análise de regressão básica. Para especificações mais complexas, deve-se consultar a documentação do R e literatura específica.

modelo	fórmula
\(y=b_0 + \epsilon\)	lm(y~1)
\(y=b_0 + b_1 x_1 + \epsilon\)	lm(y~1+x) lm(y~x)
\(y=b_1 x + \epsilon\)	lm(y~0+x)
\(y=b_0+b_1 x +b_2 x^2 + \epsilon\)	lm(y~x+I(x^2))
\(y=b_0+b_1 x +b_2 x^2 + \epsilon\)	lm(y~poly(x,2,raw=TRUE))
\(y=b_0+b_1 x +b_2 z + b_3 x\, z + \epsilon\)	lm(y~x+z+I(x*z))
\(y=b_0+b_1 x + b_2 z + b_3 x\, z + \epsilon\)	lm(y~x*z)

Algumas opções importantes do comando lm:

• data: permite a definição do data frame no qual as variáveis estão definidas, evitando definições de variável incluindo o nome do data frame.

• ex: lm(t~ano,data=olimp) em lugar de lm(olimp\(y~olimp\)x)

• subset: possibilita usar um subconjunto de observações do data frame, definido por um vetor lógico.

• ex: lm(t~ano,data=olimp,subset=ano>1950) para observações a partir de 1950

4.2 Extração de informações de modelos

Considere o modelo modelo2 que estimamos anteriormente, relativo a resultados dos homens nos 100 m, envolvendo as variáveis ano e elev. Para exemplificar a opção subset, definida no parágrafo final do último tópico, vamos somente considerar os dados relativos aos primeiros colocados. Essa informação está na variável ord, no data frame (o número indica a classificação)

## para estimar os resultados e colocá-los na variável modelo1:
modelo2<-lm(t~ano+elev,data=olimpm,subset=(ord==1)) ## t = b0+ b1 ano + b2 elev + E

Uma vez estimado, podemos extrair muitas informações uteis usando funções como

• summary(modelo): mostra uma síntese dos resultados principais

• coef(modelo): mostra as estimativas dos parâmetros

• confint(modelo): mostra os intervalos de confiança dos parâmetros

• residuals(modelo): mostra os resíduos do modelo (detalhados em outro tópico)

• fitted(modelo) ou predict(modelo): mostra os valores estimados pela regressão para cada observação no conjunto de dados (a função predict é mais geral e serve para outros propósitos)

• names(modelo): mostra todas as informações disponíveis, que podem ser extraídas usando modelo$identificador onde o identificador é o nome usado para identificar a informação específica

As informações do summary podem ser extraídas de uma forma similar. Para ver todas as informações disponíveis para extração do summary, que podem ser extraídas usando summary(modelo)$identificador onde o identificador é o nome usado para identificar a informação específica

• names(summary(modelo)): mostra todas as informações que podem ser extraídas

No contexto do modelo modelo2 o uso dessas opções levará a

## síntese das informações do modelo
summary(modelo2)

## 
## Call:
## lm(formula = t ~ ano + elev, data = olimpm, subset = (ord == 
##     1))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.265239 -0.055924 -0.003714  0.069167  0.215569 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.190e+01  1.418e+00  22.504  < 2e-16 ***
## ano         -1.104e-02  7.241e-04 -15.243 1.63e-13 ***
## elev        -9.295e-05  5.753e-05  -1.616     0.12    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1253 on 23 degrees of freedom
## Multiple R-squared:  0.9135, Adjusted R-squared:  0.9059 
## F-statistic: 121.4 on 2 and 23 DF,  p-value: 5.988e-13

## parâmetros estimados
coef(modelo2)

##   (Intercept)           ano          elev 
## 31.8996357420 -0.0110372121 -0.0000929458

## valores estimados pelo modelo / predict(modelo2) dá o mesmo resultado
fitted(modelo2)

##        60        63        66        69        72        75        78 
## 10.925773 10.871586 10.839334 10.795092 10.706794 10.660880 10.619705 
##        81        84        87        90        93        96        99 
## 10.565239 10.528433 10.397845 10.352674 10.307968 10.265399 10.220971 
##       102       105       108       111       114       117       120 
##  9.970204 10.086015 10.085457 10.034431  9.991304  9.954591  9.912394 
##       123       126       129       132       135 
##  9.846124  9.822423  9.766842  9.731058  9.691464

## resíduos do modelo
residuals(modelo2)

##           60           63           66           69           72 
##  0.074227323  0.128414318 -0.039333897  0.004907897  0.093205594 
##           75           78           81           84           87 
## -0.060879588  0.180294995 -0.265239173 -0.228433043 -0.097845415 
##           90           93           96           99          102 
##  0.047325837  0.192032360 -0.065398870 -0.220971184 -0.020203814 
##          105          108          111          114          117 
##  0.053985306 -0.025457428  0.215569410 -0.001304146 -0.034590962 
##          120          123          126          129          132 
##  0.047606024 -0.006124026  0.047576745  0.083157927 -0.041058347 
##          135 
## -0.061463843

## todas as informações disponíveis
names(modelo2)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

## verificando todas as informações relativas aos parâmetros
modelo2$coefficients

##   (Intercept)           ano          elev 
## 31.8996357420 -0.0110372121 -0.0000929458

Informações relativas ao summary do modelo

## todas as informações que podem ser extraídas do summary
names(summary(modelo2))

##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

## extraíndo o r2 ajustado
summary(modelo2)$adj.r.squared

## [1] 0.9059479

## extraíndo as informações relativas a parâmetros
summary(modelo2)$coefficients

##                  Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept) 31.8996357420 1.417512e+00  22.503956 3.641719e-17
## ano         -0.0110372121 7.240957e-04 -15.242754 1.633515e-13
## elev        -0.0000929458 5.752627e-05  -1.615711 1.197926e-01

## extraíndo somente os informações sobre o valor-p
summary(modelo2)$coefficients[,4]

##  (Intercept)          ano         elev 
## 3.641719e-17 1.633515e-13 1.197926e-01

5 - Modelos lineares vs. modelos não-lineares

Um aspecto que pode parecer restritivo a princípio é o fato dos modelos de regressão mais usuais serem lineares. Os modelos lineares, ao contrário do que indica essa primeira impressão, são muito poderosos e cobrem um amplo espectro de problemas. O contexto usado para o termo linear se refere aos parâmetros do modelo e não às variáveis. Alguns exemplos de modelos lineares e não-lineares são apresentados a seguir a título de ilustração:

• modelos lineares (exemplos):
• \(y=b_0+b_1 x_1+ b_2 x_2+\epsilon\)
• \(y=b_0+b_1 x_1 + b_2 x_1^2 + b_3 \log x_2 + b_4 x_3 ^3+\epsilon\)
• \(\log y=\displaystyle b_0+b_1 \frac{1}{x_1}+\epsilon\)

Note que em todos os casos os modelos são lineares nos parâmetros \(b_j\)s.

• modelos não-lineares (exemplos):
• \(y=\displaystyle \frac{b_0+b_1 x_1}{b_3+ b_4 x_2}\;\; +\epsilon\)
• \(y=b_0+\log (b_1 x_1 + b_2 x_2) + b_3 x_3^3+\epsilon\)
• \(y=\frac{1}{\displaystyle 1+e^{b_0+b_1 x_1}}\;\;\;+\epsilon\)

Alguns modelos não-lineares podem ser estimados pelo método dos quadrados mínimos, mas não com as fórmulas definidas para o caso usual linear, visto no tópico anterior. Deverá haver uma derivação específica que pode ser complicada de operacionalizar. O método da maxima verossimilhança é usado em muitos desses casos.

5.1 Modelos não-lineares que podem ser “linearizados”

Em alguns casos importantes, modelos que são não-lineares à primeira vista podem ser estimados a partir de transformações que criar uma versão linear apropriada para estimativa dos parâmetros.

Um exemplo importante é apresentado neste tópico, operacionalizado no contexto dos dados dos resultados das mulheres nos 100 m nos jogos olimpicos (data frame olimp)

Considere o modelo não-linear dado por:

• \(y^*_i=\displaystyle e^{b_0+b_1 x_i} E_i\)

onde \(E_i\) é uma variável aleatória com distribuição Log-normal (uma distribuição com uma relação estreita com a distribuição Normal). Aplicando logaritmo natural dos dois lados da expressão chegamos a:

• \(\ln y^*_i=b_0+b_1 x_i + \ln E_i\)

Ao troca a notação: \(y_i=\ln y^*_i\) e \(\epsilon_i=\ln E_i\), nosso modelo original resultou num modelo linear nos parâmetros \(b_0\) e \(b_1\):

• \(y_i=b_0+b_1 x_i + \epsilon_i\).

Os parâmetros desse modelo podem ser estimado pelos métodos usuais relacionados à regressão linear. A transformação se beneficia de um resultado associado à distribuição Log-normal: o logarítmo natural de uma variável aleatória Log-normal é Normal. Basta que usar as pressuposições usuais com relação a \(\epsilon_i\) para estimar o modelo e fazer inferências:

• \(\epsilon_i\sim \text{N}(0,\sigma)\) e \(\text{cov}(\epsilon_i,\epsilon_j)=0\) para \(i\neq j\).

Para plotar o modelo, não devemos esquecer de transformar novamente os valores de \(\hat y_i\) em \(\hat y^*_i\) usando:

• \(\hat y^*_i=\displaystyle e^{\displaystyle \hat y_i}\)

5.2 Estimando uma regressão não-linear que pode ser “linearizada”

Considere os dados dos resultados dos 100 m das mulheres no data frame olimp. Para facilitar vamos definir um data frame específico, denominado olimpf, somente com esses dados e fazer um diagrama de dispersão. Vamos também estimar uma regressão linear simples para esse caso, para efeito de comparação:

• \(y_i=b_0+b_1 x_i + \epsilon_i\)
• com \(\epsilon_i\sim \text{N}(0,1)\), \(\text{cov}(\epsilon_i,\epsilon_j)=0\) para \(i\neq j\) e outras premissas usuais (veja na definição geral),

onde \(y_i\) representa t (tempo nos 100 m) e \(x_i\) representa ano (ano dos jogos).

olimpf<-olimp[olimp$sexo=="F",]  
y<-olimpf$t
x<-olimpf$ano
plot(x,y,xlab="ano",ylab="t",col="deeppink")
title(main="Dados da prova 100 m - Mulheres (1928-2012)")
modelo1L<-lm(y~x)
summary(modelo1L)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50928 -0.15007  0.01599  0.16546  0.50862 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.961811   2.280468   18.84   <2e-16 ***
## x           -0.016053   0.001156  -13.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.221 on 57 degrees of freedom
## Multiple R-squared:  0.7718, Adjusted R-squared:  0.7678 
## F-statistic: 192.8 on 1 and 57 DF,  p-value: < 2.2e-16

## colocando a regressão linear no gráfico
abline(modelo1L,col="red") 

## para estimar os resultados e colocá-los na variável modelo1:
modelo2<-lm(t~ano,data=olimpf) ## t = b0+ b1 ano + b2 elev + E
## para ver os resultados
summary(modelo2)

## 
## Call:
## lm(formula = t ~ ano, data = olimpf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50928 -0.15007  0.01599  0.16546  0.50862 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.961811   2.280468   18.84   <2e-16 ***
## ano         -0.016053   0.001156  -13.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.221 on 57 degrees of freedom
## Multiple R-squared:  0.7718, Adjusted R-squared:  0.7678 
## F-statistic: 192.8 on 1 and 57 DF,  p-value: < 2.2e-16

confint(modelo2)

##                  2.5 %      97.5 %
## (Intercept) 38.3952523 47.52836932
## ano         -0.0183677 -0.01373747

Vamos agora estimar um novo modelo, usando os conceitos e premissas estabelecidas no tópico anterior:

• \(y^*_i=e^{b_0+b_1 x_i} E_i\) (mesmas premissas do tópico anterior).

Nesse caso \(y^*_i\) agora representa a t (tempo nos 100 m) e \(x_i\) representa ano (ano dos jogos)

Para linearizar esse modelo não-linear, aplicamos logarítmo natural nos 2 lados, obtendo (usando a mesma notação anterior)

• \(y_i=b_0+b_1 x_i+\epsilon_i\)

onde, \(y_i=\ln y^*_i\). Agora podemos estimar os parâmetros, tomando o cuidado de definir a variável resposta da forma indicada (aplicando o logaritmo):

### transformação da variável resposta original
logy=log(y)
### estimativa dos parâmetros
modelo2NL<-lm(logy~x)
summary(modelo2NL)

## 
## Call:
## lm(formula = logy ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.046888 -0.013493  0.001421  0.014872  0.043177 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.1948900  0.1978539   26.26   <2e-16 ***
## x           -0.0014048  0.0001003  -14.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01918 on 57 degrees of freedom
## Multiple R-squared:  0.7748, Adjusted R-squared:  0.7709 
## F-statistic: 196.2 on 1 and 57 DF,  p-value: < 2.2e-16

Observe que o \(R^2\) ajustado obtido, 0,771, é um pouco mais alto que o anterior obtido no modelo modelo1L, que foi 0,768, o que sugere um melhor ajuste.

Para visualizar essa curva, precisaremos de um procedimento específico, que serve para apresentação de qualquer curva (não será possível utilizar a função abline, que só serve para o caso da regressão linear simples).

A função fundamental para calcular os valores estimados pela regressão (as médias condicionais), chama-se predict. Essa função também será utilizada, em um próximo tópico, para obtenção de intervalos de confiança para as estimativas da regressão (médias condicionais), assim como intervalos de previsão.

A função predict é utilizada, no caso, para obter os valores estimados pela regressão especificada no modelo definido na função, a partir de um conjunto de valores das variáveis explicativas.

## definindo o dominio de x na faixa de anos de interesse
novox<-seq(1928,2012,1)
## estimando os valores da curva de regressão
yest<-predict(modelo2NL,newdata=data.frame(x=novox))
## transformando yest (em logarítmo de segundos) para segundos (tempo)
y<-exp(yest)
## desenhando a nova curva em azul
lines(novox,y,col="blue")

As duas curvas são muito similares, havendo pouca diferença prática entre os 2 modelos (a curva azul está praticamente sobre a curva vermelha).

Para um modelo melhor, que eventualmente capture o efeito de testes anti-doping mais rigorosos após 1988, precisaremos de um modelo mais sofisticado que considera variáveis binárias (próximo tópico)

5.3 Obtenção dos efeitos marginais no modelo não-linear “linearizado”

Um cuidado importante no caso que acabamos de desenvolver, é na interpretação e cálculo do efeito marginal, a partir dos parâmetros obtidos. Observe no caso anterior que

• \(\displaystyle \hat y^*_i =\,e^{\hat b_0+\hat b_1 x_i}\)

e

• \(\displaystyle \frac{d\, \hat y^*_i}{d\, x_i} = \hat b_1\,e^{\hat b_0+\hat b_1 x_i}\)

ou seja, o parâmetro \(\hat b_1\) não será mais, isoladamente, a estimativa do efeito marginal nesse caso. Uma estimativa pode ser obtida pela expressão acima (justificada por conceitos mais avançados), que depende de \(\hat b_0\), \(\hat b_1\) e do próprio valor de \(x_i\).

No caso específico, para \(x_i=2012\), temos

• \(\displaystyle \frac{d\, \hat y^*_i}{d\, x_i}=\displaystyle \hat b_1\,e^{\hat b_0+\hat b_1 2012}.\)

Para facilitar os cálculos podemos utilizar diretamente os valores dos parâmetros já estimados que ficam armazenados dentro do “modelo” estimado, na opção coefficients.

O trecho de código abaixo mostra o procedimento completo para encontrar o efeito marginal associado à variável explicativa \(x\) (no caso ano).

## obtendo os parâmetros estimados no vetor best
best<-coef(modelo2NL)
best

##  (Intercept)            x 
##  5.194890041 -0.001404841

## achando o efeito marginal em x=2012
## lembre que b0 está na posição 1 e b1 na posição 2 do vetor
best[2]*exp(best[1]+best[2]*2012)

##           x 
## -0.01500335

O efeito marginal estimado agora, de -0,015, é muito similar ao estimado anteriormente, -0,014, no modelo1L para mulheres, que considerou uma regressão linear simples.

6 - Uso da funções predict em gráficos de regressões

Um bom entendimento da função predict é central para a elaboração de gráficos, assim como estimativa dos intervalos de confiança e intervalos de previsão relacionados com a variável resposta, condicional aos valores observados das variáveis explicativas (outro tópico). Neste tópico focaremos no primeiro uso (gráficos).

Considere os dados abaixo. Suponha que esses dados se referem a um experimento para caracterizar a curva de resposta (média) \(y\) de um processo de produção, a partir do uso de um insumo \(x\), estimada por uma regressão.

x<-c(0,0,0,50,50,50,100,100,100,150,150,150)
y<-c(70,73,67,80,85,82,87,90,85,90,88,87)

Assuma que o modelo (mod1) relevante é construído com o apoio de uma equação do segundo grau em \(x\) (parábola), para valores de \(x\) entre 0 e 150.

• modelo mod1: \(y_i=b_0+b_1 x_i+b_2 x_i^2 + \epsilon_i\), com \(\epsilon_i\sim \text{N}(0,\sigma)\) e \(\text{cov}(\epsilon_i,\epsilon_j)=0\) para \(i\neq j\) (e outras premissas usuais quanto aos \(x_i\)s)

Problema: estime os parâmetros da regressão e apresente a curva dentro do diagrama de dispersão.

6.1 Alternativa 1: sem uso da função predict (não recomendada)

## estimativa do modelo mod1
## note como introduzir o termo quadrático utilizando I(x^2)
## uma opção equivalente seria: mod1<-lm(y~poly(x,2,raw=TRUE))
mod1<-lm(y~x+I(x^2))
## O sumário dos resultados pode ser acessado por
summary(mod1)

## 
## Call:
## lm(formula = y ~ x + I(x^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1667 -1.3333 -0.1667  1.9167  3.1667 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 70.1666667  1.3219421  53.078 1.50e-12 ***
## x            0.2900000  0.0424584   6.830 7.64e-05 ***
## I(x^2)      -0.0011333  0.0002713  -4.178  0.00238 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.349 on 9 degrees of freedom
## Multiple R-squared:  0.9276, Adjusted R-squared:  0.9115 
## F-statistic: 57.65 on 2 and 9 DF,  p-value: 7.393e-06

Com essas estimativas, a curva de regressão será dada por:

• \(\hat y= 70{,}166667+0{,}29\,x-0{,}001133\, x^2\)

Para facilitar, os parâmetros de um modelo produzido pela função lm ficam armazenados num vetor que pode ser recuperado por

best<-coef(mod1)
best

##  (Intercept)            x       I(x^2) 
## 70.166666667  0.290000000 -0.001133333

Poderíamos desenhar essa curva de regressão usando procedimentos já discutidos para fazer um gráfico de uma função qualquer vista no tópico 4 de módulo anterior, definindo um vetor com pontos do domínio da função, calculando o contradominio através da função desejada e produzindo o gráfico:

## diagrama de dispersão 
plot(x,y,col="blue")
## criando o domínio
novox<-seq(0,150,0.1)
## encontrando os pontos no contradomínio da função
yest<-70.166667+0.29*novox-0.001133 * novox^2
## ou usando diretamente os valores estimados
yest<-best[1]+best[2]*novox+best[3]*novox^2
## desenhando a curva
lines(novox,yest,col="red")

6.2 Alternativa 2: com uso da função predict (recomendada)

O mesmo resultado obtido acima pode ser conseguido com menos esforço, através da função predict, que tem múltiplos usos. O desenvolvimento abaixo assume que o modelo mod1 já está estimado.

## diagrama de dispersão 
plot(x,y,col="blue")
## criando o domínio
novox<-seq(0,150,0.1)
## estimando os valores de yest
yest<-predict(mod1,newdata=data.frame(x=novox))
## desenhando a curva
lines(novox,yest,col="red")

No exemplo, a função utilizou 2 argumentos: o nome do modelo (mod1) e newdata, que foi definido com um data frame contendo os valores das variáveis explicativas, que serão usadas para a obtenção dos valores de \(y\) estimados pela regressão. Toda a complexidade do modelo é “escondida” dentro do próprio modelo estimado (mod1), algo que facilita a operacionalização dos procedimentos.

Compare os 2 procedimentos para entender melhor como utilizar essa função predict que pode facilitar muitos desenvolvimentos mais complexos. No caso de uma regressão linear poderíamos utilizar diretamente abline(mod1) após o modelo ser estimado, para colocar a reta no diagrama de dispersão. É interessante, contudo, conhecer a forma mais geral de se fazer isso, que se aplica a qualquer função, linear ou não-linear.

7. Uso de variáveis “dummy” ou binárias em regressões com variáveis qualitativas

A inclusão de variáveis quantitativas numa regressão, como fizemos no tópico anterior, não é algo complicado. A dificuldade maior está na inclusão de variáveis qualitativas ou categóricas (fatores), como variáveis explicativas ou como variáveis resposta. Examinaremos inicialmente a primeira situação (variáveis explicativas).

7.1 Variáveis “dummy” ou binárias

A inclusão de variáveis qualitativas (também chamadas de fatores) em uma regressão pode se fazer com apoio das chamadas

• variáveis “dummy” ou binárias: são variáveis que assume valores 0 e 1 de acordo com a categoria de uma variável de interesse. Em alguns casos a própria variável categórica já está definida desta forma (como uma variável Bernoulli que já vimos anteriormente)

Para tornar a situação mais concreta, suponha que, na situação dos dados relativos aos resultados das mulheres nos 100 m nos jogos olimpicos (já organizado no data frame olimpf, definido nos outros topicos), desejamos incluir a possibilidade de haver uma mudança de regime no modelo, visando capturar o possível efeito do maior rigor nos exames anti-doping após 1988. Desejamos, inclusive, verificar se esse efeito foi estatisticamente significativo.

Para as análise vamos criar uma nova variável (categórica) denominada contdop para representar o período até 1988 (categoria A) e posterior a 1988 (categoria D):

## parâmetros estimados ficam da opção coefficients do modelo
contdop<-as.factor(ifelse(olimpf$ano<=1988,"A","D"))
contdop

##  [1] A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
## [36] A A A A A A A D D D D D D D D D D D D D D D D D
## Levels: A D

7.2 Análise com regressões separadas (não é boa prática!)

Uma análise menos elaborada, poderia considerar a estimativa de 2 regressões separadas, uma para o período 1928-1988 e outra para o período 1989-2012 e tentar construir o argumento sobre o resultado da análise. Essa primeira alternativa será examinada a seguir. As categorias da variável contdop estariam representadas em cada uma das regressões:

y<-olimpf$t[contdop=="A"]
x<-olimpf$ano[contdop=="A"]
modeloA<-lm(y~x)
novox<-1928:1988
yest<-predict(modeloA,newdata=data.frame(x=novox))
lines(novox,yest,col="blue")
summary(modeloA)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47472 -0.15417  0.02125  0.15394  0.48138 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 53.291804   3.467552   15.37  < 2e-16 ***
## x           -0.021341   0.001769  -12.06 6.59e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2133 on 40 degrees of freedom
## Multiple R-squared:  0.7845, Adjusted R-squared:  0.7791 
## F-statistic: 145.6 on 1 and 40 DF,  p-value: 6.592e-15

### estimando o modelo para o período depois de 1988
y<-olimpf$t[contdop=="D"]
x<-olimpf$ano[contdop=="D"]
modeloD<-lm(y~x)
summary(modeloD)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.132635 -0.109921  0.003536  0.060450  0.256993 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 17.651610   8.317354   2.122   0.0509 .
## x           -0.003364   0.004154  -0.810   0.4307  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1201 on 15 degrees of freedom
## Multiple R-squared:  0.04189,    Adjusted R-squared:  -0.02198 
## F-statistic: 0.6558 on 1 and 15 DF,  p-value: 0.4307

novox<-1989:2012
yest<-predict(modeloD,newdata=data.frame(x=novox))
lines(novox,yest,col="blue")

Note que, no caso da segunda regressão estimada (período 1989-2012), os valores-p associados aos 2 parâmetros estimados foram superiores a 5%, sugerindo a não rejeição das hipóteses \(\H_0: b_0=0\) e \(\H_0: b_1=0\). Isso ocorre, em grande medida, pela razão de que as observações foram divididas, ficando cada regressão com um número menor de observações que as existentes no conjunto de dados como um todo. Um baixo número de observações pode dificultar a rejeição da hipótese da nulidade, se os efeitos não forem muito pronunciados.

7.3 Análise utilizando variável “dummy” (recomendada)

O procedimento ideal, desenvolvido a seguir, não separa os dados em 2 regressões, mas considera um único modelo que comporta as 2 regressões de uma só vez, com o apoio de 1 variável “dummy” que denominaremos z

• \(z\): assume valor 0 para o período A (1928-1988) e valor 1 para o período D (1989-2012)

Considere agora o seguinte modelo:

• \(y_i=b_0+b_1 x_i + b_2 z + b_3 z_i x_i + \epsilon_i\)
• com \(\epsilon_i\sim \text{N}(0,\sigma)\), \(\text{cov}(\epsilon_i,\epsilon_j)\), \(i\neq j\) e outras premissas usuais (veja na definição geral)

No modelo, \(y_i\) é a variável resposta que representa t (tempo em seg), \(x_i\) representa ano e \(z_i\) é a variável “dummy” que assume valor 0 para o período A (1928-1988) e valor 1 para o período D (1989-2012).

Observe como o modelo se especializa na medida que \(z=0\) ou \(z=1\)

• \(z_i=0\) (período A 1928-1988): \(y_i=b_0+b_1 x_i + \epsilon_i\)

• \(z_i=1\) (período D 1989-2012): \(y_i=(b_0+b_2) + (b_1+b_3) x_i + \epsilon_i\)

Observamos então que os parâmetros no modelo geral representam:

• \(b_0\): intercepto da curva de regressão no período A

• \(b_1\): inclinação (efeito marginal) da curva de regressão no período A

• \(b_0+b_2\): intercepto da curva de regressão no período D

• \(b_1+b_3\): inclinação (efeito marginal) da curva de regressão no período D

Logo, nesse modelo, os parâmetros \(b_2\) e \(b_3\) representam os “efeitos”, respectivamente, no intercepto e inclinação da regressão observados ao se passar do período A (1928-1988) para o período D (1989-2012).

A vantagem dessa formulação, com relação à formulação anterior, é a possibilidade de testarmos, diretamente, hipóteses envolvendo \(b_2\) e \(b_3\) utilizando os valores-p associados a esses novos parâmetros, considerando todas as observações disponíveis.

• Só efeito no intercepto: para estimar somente o efeito no intercepto, sem considerar o efeito na inclinação da curva, o modelo seria definido somente por:

• \(y_i=b_0+b_1 x_i + b_2 z + \epsilon_i\)

Operacionalização do procedimento

Inicialmente, vamos criar a variável \(z\) no conjunto de dados, a partir da variável contdop criada anteriormente.

z<-ifelse(contdop=="D",1,0)
z

##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [36] 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Para estimar o modelo, note que precisamos de uma nova variável, associada ao parâmetro \(b_3\), que considera o produto de \(z\) e \(x\). Essa nova variavel é indicada na fórmula do comando lm por \(\text{I}(z * ano)\), dado que ano está associada à variável \(x\). Essa é a notação utilizada para criação de novas variáveis a partir das variáveis existentes.

modelocomp<-lm(t~ano+z+I(ano*z),data=olimpf)
summary(modelocomp)

## 
## Call:
## lm(formula = t ~ ano + z + I(ano * z), data = olimpf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47472 -0.12492  0.01967  0.12821  0.48138 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  53.291804   3.127992  17.037   <2e-16 ***
## ano          -0.021341   0.001596 -13.375   <2e-16 ***
## z           -35.640193  13.687233  -2.604   0.0118 *  
## I(ano * z)    0.017977   0.006844   2.627   0.0111 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1924 on 55 degrees of freedom
## Multiple R-squared:  0.8332, Adjusted R-squared:  0.8241 
## F-statistic: 91.56 on 3 and 55 DF,  p-value: < 2.2e-16

Os resultados mostram que os estimadores dos parâmetros \(b_2\) e \(b_3\), são, respectivamente, \(\hat b_2=-35.64019\) e \(\hat b_3=0.01798\). Os valores-p associados às hipóteses da nulidade \(\H_0: b_2=0\) e \(\H_0: b_3=0\) são próximos de 1%, o que suporta a rejeição dessas hipóteses.

Pelo modelo, as regressões para os períodos A e D seriam dadas por

• Periodo A (1928-1988): \(\hat y= 53.29180-0.02134\, x\)
• Periodo D (1989-2012): \(\hat y= (53.29180-35.64019)+(-0.02134+0.01798)\, x= 17.65161-0.00336\, x\)

A próxima figura mostra essas duas retas dentro do diagrama de dispersão, usando o modelo modelocomp estimado nos parágrafos anteriores.

## diagrama de dispersão
plot(olimpf$ano,olimpf$t,xlab="ano",ylab="t",col="deeppink")
title(main="Resultados das mulheres nos 100 m")
## regressão para o período 1928-1988
novoano=1928:1988
novoz=0 ## valor de z
novoy=predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz))
lines(novoano,novoy,type="l")
## regressão para o período 1989-2012
novoano=1989:2012
novoz=1 ## valor de z
novoy=predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz))
lines(novoano,novoy,type="l")

O resultado obtido oferece uma evidência sugerindo que um impacto de exames anti-doping mais rigorosos introduzidos após 1988 (ex. esteróides anabolizantes, etc.) nos tempos das atletas femininas. Isso indica que, possivelmente, essas substâncias eram utilizadas com maior frequência no período anterior a 1989. Essa evidência também suporta a conclusão de que a própria tendência de redução dos tempos das mulheres ao longo dos anos parece ter sido reduzida substancialmente após 1988.

7.4 Considerando variáveis com mais de 2 categorias

O processo de criação e uso de variáveis “dummy” que mostramos, considerou um caso comum envolvendo uma situação de uma variável qualitativa com somente 2 categorias (variável contdop com categorias: período A e período D).

E se tivessemos 3 categorias para os períodos considerados? por exemplo: periodo A1 (1900-1967), período A2 (1968-1988) e período D (1999-2012). Como os exames anti-doping foram iniciados em 1968, nos jogos olimpicos do México, parece razoável considerar 3 períodos, para que cobririam as situações: sem exame anti-doping (1900-1967), com exames menos rigorosos (1968-1988), com exames muito rigorosos (1989-2012)).

Situação 1: somente efeito no intercepto da regressão

Nesse caso precisaríamos de 2 variáveis “dummy”, \(z_1\) e \(z_2\) para capturar as 3 categorias. Suponha que o periodo A1 será período “base” para comparações e só estamos interessados no efeito no intercepto do modelo.

A codificação de \(z_1\) e \(z_2\) (chamada tecnicamente de contraste de tratamento) podería ser dada por:

• Período A1 (1900-1967): \(z_1=0\) e \(z_2=0\)

• Período A2 (1968-1988): \(z_1=1\) e \(z_2=0\)

• Período D (1989-2012): \(z_1=0\) e \(z_2=1\)

Nosso modelos, nessa situação em que somente estamos interessados no efeito no intercepto, seria definido por:

• \(y_i=b_0+b_1 x_i + b_2 z_1+b_3 z_2 + \epsilon_i\), com \(\epsilon_i\sim \text{N}(0,\sigma)\) e \(\text{cov}(\epsilon_i,\epsilon_j)\), \(i\neq j\) e outras premissas usuais (veja na definição geral)

Com as definições teríamos para cada situação possível de \(z_1\) e \(z_2\) a especialização do modelo para

• Período A1 (\(z_1=0\) e \(z_2=0\)): \(y_i=b_0+b_1 x_i+ \epsilon_i\)

• Período A2 (\(z_1=1\) e \(z_2=0\)): \(y_i=(b_0+b_2)+b_1 x_i + \epsilon_i\)

• Período D (\(z_1=0\) e \(z_2=1\)): \(y_i=(b_0+b_3)+b_1 x_i + \epsilon_i\)

No caso, \(b_2\) e \(b_3\) são os efeitos (no intercepto do período base A1) de passar do período base (A1), para os períodos A2 e D, respectivamente.

Para operacionalizar, inicialmente criaremos as variáveis z1 e z2 de forma consistente com essa formulação:

z1=ifelse(olimpf$ano>=1968&olimpf$ano<=1988,1,0)
z2=ifelse(olimpf$ano>1988,1,0)

O novo modelo estimado será denominado modelocomp2:

## 
## Call:
## lm(formula = t ~ ano + z1 + z2, data = olimpf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56031 -0.09695  0.01781  0.10858  0.45305 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.043605   5.126275   8.202 4.10e-11 ***
## ano         -0.015553   0.002633  -5.907 2.27e-07 ***
## z1          -0.232033   0.101796  -2.279   0.0265 *  
## z2           0.010611   0.157746   0.067   0.9466    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1951 on 55 degrees of freedom
## Multiple R-squared:  0.8284, Adjusted R-squared:  0.8191 
## F-statistic: 88.53 on 3 and 55 DF,  p-value: < 2.2e-16

Situação 2: efeito no intercepto e na inclinação

Para incluir o efeito no na inclinação , precisaríamos utilizar o modelo

• \(y_i=b_0+b_1 x_i + b_2 z_1+b_3 z_2 + b_4 x\, z_1 + b_5 x\, z_2 + \epsilon_i\), com \(\epsilon_i\sim \text{N}(0,\sigma)\) e \(\text{cov}(\epsilon_i,\epsilon_j)\), \(i\neq j\) e outras premissas usuais (veja na definição geral)

Esse novo modelo estimado será denominado modelocomp3 e estimado por

## 
## Call:
## lm(formula = t ~ ano + z1 + z2 + I(z1 * ano) + I(z2 * ano), data = olimpf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52096 -0.10628  0.00524  0.10641  0.45608 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  47.937107   6.018390   7.965 1.27e-10 ***
## ano          -0.018580   0.003091  -6.011 1.74e-07 ***
## z1           -9.338059  14.354244  -0.651   0.5182    
## z2          -30.285497  14.529027  -2.084   0.0420 *  
## I(z1 * ano)   0.004651   0.007277   0.639   0.5255    
## I(z2 * ano)   0.015215   0.007292   2.086   0.0418 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1909 on 53 degrees of freedom
## Multiple R-squared:  0.8417, Adjusted R-squared:  0.8267 
## F-statistic: 56.35 on 5 and 53 DF,  p-value: < 2.2e-16

A representação gráfica do modelo modelocomp3, que considera os 3 períodos, seria dada por

É interessante observar que os efeitos no intercepto e inclinação dados por \(b_2\) e \(b_4\) não foram estatísticamente significativos (valor-p >> 5%), enquanto que \(b_3\) e \(b_5\) foram significativos (valor-p = 4,2%). Apesar de apresentar um R2-ajustado um pouco mais alto (0,827) que na situação em que consideramos apenas 2 períodos (A e D), no modelo modelocomp (R2-ajustado de 0,824): a ausência de significância para \(b_2\) e \(b_3\) pode sugerir uma superioridade do modelo modelocomp sobre o modelo modelocomp3.

Como incluir mais de 3 categorias

A consideração de uma variáveis qualitativa ou categórica (fator) com apoio de variáveis “dummy” requer que tenhamos sempre \(k-1\) variáveis dummy definidas para uma variável com \(k\) categorias (1 a menos). A inclusão de novas categorias pode considerar a extensão do mesmo procedimento utilizado para 2 e 3 categorias, mostrado nos últimos tópicos, definindo-se uma categoria base e definindo de forma apropriada cada variável “dummy”" associada as outras categorias.

7.5 Algumas facilidades para inclusão de variáveis qualitativas no R

Para inclusão de variáveis qualitativas nos modelos, como variáveis explicativas, o R não requer, de forma obrigatória, a criação de variáveis “dummy”, como fizemos nos desenvolvimentos dos últimos tópico. Basta incluir o nome da variável qualitativa ou categórica no modelo, que o R criará (internamente), tantas variáveis dummy quanto necessário para estimar o modelo (num esquema que corresponde ao contraste de tratamentos no contexto mais formal, o qual pode ser alterado).

No exemplo utilizado até o momento, poderíamos diretamente ter considerado o modelo com a variável qualitativa ou categórica (fator) contdop, em lugar de ter criado a variável dummy \(z\), o resultado será o mesmo (veja no tópico 4 o significado da fórmula)

modelocomp4<-lm(t~contdop*ano,data=olimpf)
summary(modelocomp4)

## 
## Call:
## lm(formula = t ~ contdop * ano, data = olimpf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47472 -0.12492  0.01967  0.12821  0.48138 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   53.291804   3.127992  17.037   <2e-16 ***
## contdopD     -35.640193  13.687233  -2.604   0.0118 *  
## ano           -0.021341   0.001596 -13.375   <2e-16 ***
## contdopD:ano   0.017977   0.006844   2.627   0.0111 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1924 on 55 degrees of freedom
## Multiple R-squared:  0.8332, Adjusted R-squared:  0.8241 
## F-statistic: 91.56 on 3 and 55 DF,  p-value: < 2.2e-16

A categoria da variável codificada como 0 (ou base de comparação) será a primeira categoria da variável qualitativa (essa ordem pode ser alterada em caso de interesse). No caso específico será a categoria A, associada ao período (1928-1988). Os efeitos são apresentados da mesma forma que no desenvolvimento anterior. Compare os resultados com os obtidos anteriormente.

Se deseja somente o efeito da variável qualitativa no intercepto, basta usar:

modelocomp5<-lm(t~ano+contdop,data=olimpf)
summary(modelocomp5)

## 
## Call:
## lm(formula = t ~ ano + contdop, data = olimpf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46225 -0.12881  0.01213  0.15640  0.49338 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 51.376350   3.198017  16.065   <2e-16 ***
## ano         -0.020364   0.001631 -12.483   <2e-16 ***
## contdopD     0.311333   0.089656   3.473    0.001 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2023 on 56 degrees of freedom
## Multiple R-squared:  0.8122, Adjusted R-squared:  0.8055 
## F-statistic: 121.1 on 2 and 56 DF,  p-value: < 2.2e-16

Em caso de dúvida, crie as variáveis “dummy” da forma que fizemos anteriormente.

7.6 Aprimoramentos possíveis no modelo

O modelo desenvolvido no último tópico pode ser aprimorado em algumas direções. A variável elev (elevação) pode ser avaliada para inclusão no modelo, pois já vimos que ela pode ser importante. Um outro aspecto possível de melhoria, seria relaxar a premissa sobre o comportamento da variância de \(\epsilon\) que foi considerada fixa (\(\sigma\)) em todo o período. Pode ser que a própria variância seja afetada pelos 2 regimes com relação ao rigor na aplicação dos exames anti-doping (períodos A e D). A própria premissa \(\text{cov}(\epsilon_i,\epsilon_j)=0\) precisará ser examinada, eventualmente sugerindo outras direções de aprimoramento, caso existam indícios de correlação nos erros. Certamente existirão outras possibilidades melhores.

Alguns desses aprimoramentos podem ser facilmente implementados (ex. inclusão da variável elev), outros exigem conhecimento mais avançado. O que tende a ocorrer, a partir de um certo ponto é que a inclusão desse aprimoramentos no modelo pode não levar a ganhos expressivos (ou praticamente nenhum ganho), mas os custos em termos de conhecimento necessário e trabalho podem ser expressivos.

Em um próximo tópico desse módulo discutiremos diagnósticos que tenta avaliar a qualidade das premissas utilizadas, algo que pode dar mais ou menos credibilidade às estimativas realizadas, eventualmente sugerindo novas direções de ajuste dos modelos.

8 - Intervalos de confiança e intervalos de previsão para a regressão

8.1 Intervalo de confiança para a regressão

Da mesma forma que existem intervalos de confiança para a média teórica (esperanca matemática) e para os próprios parâmetros da regressão, também é possível estimarmos o intervalo de confiança no contexto da própria regressão.

Para tanto é necessário relembrar que a regressão estima \(E(Y|x_1,x_2,\ldots,x_m)\) uma média teórica condicional à observação das variáveis explicativas. Sem qualquer condicionamento, a estimativa seria a de \(E(Y)\), a média teórica (esperança) incondicional, que poderia utilizar estimadores como a média e/ou intervalo de confiança da média teórica convencional, para uma inferência mais técnica.

A teoria necessária para obtenção do intervalo de confiança para as estimativas da regressão se fundamenta em princípios muito similares aos utilizados para obtenção do intervalo de confiança para a média teória, a partir de uma quantidade pivotal apropriada. Não entraremos nos detalhes desse procedimento, o qual é implementado de forma apropriada pelo R em funções como predict.

Obtenção do intervalo de confiança para o valor estimado pela regressão

Considere o modelo modelocomp estimado no último tópico, envolvendo os dados femininos da prova dos 100 m, considerando 2 períodos com relação ao rigor na aplicação dos tests anti-doping (A e D).

• Pergunta: qual seria uma estimativa para o tempo médio das atletas na final da próxima olimpíada de 2016?

• Resposta: ara 2016, essa estimativa seria dada pelo ponto na curva de regressão estimada, em 2016, e pelo intervalo de confiança correspondente, obtido pelo uso da função predict, que é demostrado a seguir.

Inicialmente para o valor estimado pela curva de regressão estimada (considerando o modelo modelocomp estimado anteriormente)

novoano<-2016
novoz<-1 ## 0 indica período com controle menos rigido
predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz))

##        1 
## 10.86918

Isso significa que o ponto estimado na curva foi 10,87 seg, para o ano 2016. Mas que precisão teríamos associada a essa estimativa? Esse valor é somente a estimativa do ponto na curva (teórica) de regressão. A resposta leva a noção de intervalo de confiança. Mais especificamente podemos perguntar:

• Pergunta: Qual é o intervalo de confiança para a média teórica condicional do tempo das atletas nos 100 m, estimado para 2016, considerando que estamos no período com exames anti-doping mais rigorosos?

• Resposta: Usando a função predict com algumas opções especificas (interval=“confidence” e level=0.95).

novoano<-2016
novoz<-1 ## 0 indica período com controle mais rígido (após 1988)
predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz), interval="confidence",level=0.95)

##        fit      lwr      upr
## 1 10.86918 10.66174 11.07661

O resultado mostra que a estimativa do ponto na curva foi 10,87 seg, e que da média condicional “teórica” estaria entre 10,66 e 11,08 seg com 95% de probabilidade.

Podemos colocar no gráfico a faixa que define o intervalo de confiança para a regressão como um todo, repetindo esse procedimento para muitos pontos do domínio, conforme exemplificado no trecho de código, apresentado a seguir, que parte do último gráfico elaborado no tópico anterior.

## definição das valores do domínio
novoano=1928:2012
novoz=ifelse(novoano<=1988,0,1)
## definição dos pontos da curva que define o limite inferior do intervalo
liminf=predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz), interval="confidence",level=0.95)[,2]
## definição dos pontos da curva que define o limite superior do intervalo
limsup=predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz), interval="confidence",level=0.95)[,3]
## apresentação das curvas no gráfico
lines(novoano,liminf,col="blue",lty=2)
lines(novoano,limsup,col="blue",lty=2)

Podemos inferir que a “verdadeira” curva de regressão estará contida nesse intervalo com 95% de probabilidade, que na figura é delimitado pela linhas tracejadas azuis (para outra probabilidade, basta alterar a opção level na função predict).

8.2 Intervalo de previsão para a regressão

O intervalo de previsão estima o intervalo dentro do qual observaríamos valores da variável resposta, condicional a valores definidos da variável explicativa, com uma probabilidade definida (ex. 95%).

• Pergunta: qual seria a previsão para o intervalo de valores de tempos na final dos 100 m feminino que podem ocorrer na final olímpica de 2016?

• Resposta: A mais defensável “previsão” estatística é dada exatamente pelo intervalo de previsão para os tempos nos 100 m, para o ano de 2016 e considerando a continuidade de exames anti-doping mais rigorosos. A qualidade da previsão é dependente da qualidade do modelo utilizado na estimativa (que frequentemente pode ser aprimorada).

Para obter o intervalo de previsão para 2016, usamos a mesma função predict com opções específicas (interval=“prediction” e level=0.95):

novoano<-2016
novoz<-1 ## 1 indica período posterior a 1988
predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz), interval="prediction",level=0.95)

##        fit      lwr      upr
## 1 10.86918 10.43133 11.30702

O resultado indica que o valor na curva para 2016 será 10,87 seg, e os valores observados estarão no intervalo 10,43 a 11,31 segundos com 95 % de probabilidade.

Da mesma forma que fizemos para o intervalo de confiança para regressão, podemos estabelecer o intervalo de previsão (linha tracejada vermelha) para toda a regressão:

liminf=predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz), interval="prediction",level=0.95)[,2]
limsup=predict(modelocomp,newdata=data.frame(ano=novoano,z=novoz), interval="prediction",level=0.95)[,3]
lines(novoano,liminf,col="red",lty=2)
lines(novoano,limsup,col="red",lty=2)

9 - Notações utilizadas em modelos de regressão

A literatura que trata de modelos de regressão pode considerar algumas notações particulares, que são equivalentes ao modelo

• regressão linear múltipla: \(y_i=b_0+b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im} + \epsilon_i\) com as premissas:

• \(\epsilon_i\sim \text{N}(0,\sigma)\), para \(\;\;i=1,\ldots,n\) com \(\text{cov}(\epsilon_i,\epsilon_j)=0\), para \(i\neq j\), ou a premissa mais forte, de que os erros \(\epsilon_i\)s são independentes e identicamente distribuídos (i.i.d.)

• \(\text{cov}(\epsilon_i,x_{ij})=0\), \(\forall i,j\), algo que ocorrerá se os \(x_i\) forem considerados “fixos”.

Note que, como \(E[\epsilon_i]=0\), as premissas associadas a \(\epsilon\) são comumente expressas por

• \(E(\epsilon_i\,\epsilon_j)=0\) e \(E(\epsilon_i\,x_{ij})=0\).

[lembrete: para 2 v.a.s \(X\) e \(Y\), \(\text{cov}(X,Y)=E(XY)-E(X)E(Y)\)]

Também é comum o uso da letra grega beta, \(\beta_j\), em lugar de \(b_j\) para identificar parâmetros do modelo, assim como o uso de \(e_i\), em lugar de \(\epsilon_i\).

9.1 Notações vetoriais e matriciais

Para tornar os desenvolvimentos mais compactos, é usual o uso da notação matricial. Na forma mais compacta, o modelo dos parágrafos anteriores, contendo \(n\) observações e \(m\) variáveis explicativas, pode ser definido por

• \(\def\X{{\mathbf X}} \def\x{{\mathbf x}} \def\b{{\mathbf b}} \def\y{{\mathbf y}} \def\I{{\mathbf I}} \def\bfepsilon{{\boldsymbol{\epsilon}}} \y=\X^T\b+\bfepsilon\)

onde \(\y=\left(\begin{array}{c} y_1\\ y_2\\ \vdots\\ y_m\\ \end{array} \right)\; \X=\left(\begin{array}{cccc} 1&x_{11}&\ldots&x_{1m}\\ 1&x_{21}&\ldots&x_{2m}\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_{n1}&\ldots&x_{nm}\\ \end{array}\right) \b=\left(\begin{array}{c} b_0\\ b_1\\ \vdots\\ b_m\\ \end{array} \right)\;\; \text{e}\;\; \bfepsilon=\left(\begin{array}{c} \epsilon_1\\ \epsilon_{2}\\ \vdots\\ \epsilon_n\\ \end{array} \right).\)

A premissa com relação ao erro pode ser expressa por

• \(\bfepsilon\sim \text{N}_n(\boldsymbol{0},\sigma^2\I)\)

que indica que os \(n\) erros tem distribuição Normal-Multivariada (\(\text{N}_n\)), com esperança definida por um vetor com \(n\) zeros, representado por \(\boldsymbol{0}\), e matriz de covariância \(\sigma^2\mathbf{I}\). Nessa representação, \(\mathbf{I}\) representa a matriz identidade (uma matriz com 1 na diagonal e zero nos outros termos) com dimensão \(n\times n\). Essa premissa é técnicamente diferente da anterior, implicando que os erros são independêntes (pela distribuição Normal-Multivariada que definimos). Essa diferença é muito pequena nesse nível introdutório.

Podemos concluir dessa representação, pelas propriedades da Normal-Multivariada que

• \(\y|\X\sim \text{N}_m(\X^T\b,\,\sigma^2\I)\)

Quando queremos particularizar cada observação, num modelo usando a notação vetorial, podemos usar

• $ y_i=^T_i+$ onde \(latex \b=\left(\begin{array}{c} b_0\\ b_1\\ \vdots\\ b_m\\ \end{array} \right)\; \text{e}\; \x_i=\left(\begin{array}{c} 1\\ x_{i1}\\ \vdots\\ x_{im}\\ \end{array} \right)\) para \(\;\;i=1,\ldots,n\).

A notação \(\b^T\) significa o vetor \(\b\) transposto, ou \(\b^T=[b_0,b_1,\ldots,\b_m]\).

Logo, nessa nova notação, como \(\epsilon_i\sim \text{N}(0,\sigma)\), temos

• \(y_i|\x_i\sim \text{N}(\b^T\x_i,\sigma)\) com \(E(y_i|\x_i)=\b^t\x_i\) e \(V[y_i|\x_i]=\sigma^2\), com \(y_i|\x_i\), \(i=1,\ldots,n\) independentes.

Uma consequência importante desse modelo, sob a pressuposição de erros com distribuição Normal-Multivariada, é que \(y_i|\x_i\), \(i=1,\ldots,n\) são condicionalmente independentes (independência frente a observação dos \(\x_i\)). O conceito de independência condicional é relativamente mais avançado e não será explorado nesse momento.

10 - Análise dos resíduos da regressão e diagnóstico de problemas

A qualidade dos resultados da regressão é totalmente dependente da validade das premissas consideradas no modelo selecionado, tendo-se em vista os dados utilizados para estimativa do modelo. Para exame de possíveis problemas, há um grande número de técnicas que são exploradas em detalhe na literatura especializada.

Algumas referências mais detalhadas envolvendo o R sobre o assunto, como Fox, 2009 e Faraway, 2002, capítulos 6 e 7, podem ser úteis para um aprofundamento maior nas muitas técnicas disponíveis. Uma apresentação bem sumarizada e didática, focada nos recursos do R, é apresentada aqui.

Neste tópico apresentaremos somente alguns diagnósticos mais essenciais, visando verificar a existência de problemas mais sérios, com relação às premissas consideradas.

• Alguns problemas importantes examinados num diagnóstico:
• Outliers: valores extremos potencialmente problemáticos.
• Influência exagerada de observações: algumas observações tem influência exagerada sobre as estimativas da regressão.
• Ausência de normalidade: análise dos resíduos sugere rejeição da premissa de normalidade dos erros.
• Correlação serial: análise dos resíduos sugere rejeição da premissa de ausência de correlação dos erros (\(\text{cov}(\epsilon_i,\epsilon_j)=0\) para \(i\neq j\)).
• Existência de multicolinearidade: variáveis explicativas tem alta correlação, o que leva a problemas diversos envolvendo aumento da variância nas estimativas.
• Heteroscedasticia (variância não homogênea): a variância não parece ser constante (como estabelecido nas premissas).
• Correlação não nula entre variáveis explicativas e o erro: pode levar a problemas de tendenciosidade para o estimador de quadrados mínimos ordinários.
• Modelo inadequado: o modelo definido não é adequado para modelar a situação de interesse.

A correção desses problemas, quando possível, pode envolver conhecimentos mais avançados, que estão fora do escopo deste material. As consequência principais de muitos desses problemas, caso não sejam corrigidos de forma apropriada, são pelo menos 3, associados às estimativas de quadrados mínimos (ordinários): (a) estimativas podem apresentar uma precisão acima do que ocorreria se o problema fosse corrigido (ex. correlação dos erros, heteroscedasticia), facilitando a rejeição de \(\H_0\); (b) estimativas não serão não-tendenciosas, caso o problema não seja corrigido (ex. correlação dos erros com a variável explicativa); e (c) inferências podem não se justificar em razão das distribuições usadas nos testes e intervalos de confiança não serem válidas (ex. problema de não-normalidade).

10.1 Caracterização técnica dos resíduos da regressão

Para deixar a discussão mais concreta, considere o modelo comumente utilizado na regressão múltipla:

• \(y_i=b_0+b_1 x_{i1}+b_2 x_{i2}+\ldots+b_m x_{im} + \epsilon_i\)
• com \(\epsilon_i\sim \text{N}(0,\sigma)\), \(\text{cov}(\epsilon_i,\epsilon_j)=0\), \(i\neq j\). Assumimos também que os valores das variáveis explicativas são constantes ou seja, somente os \(y_i\)s e os erros \(\epsilon_i\)s são variáveis aleatórias, levando à premissa \(\text{cov}(x_j,e_i)=0\), \(\forall i,j\).

As estimativas dos pontos na curva de regressão para cada observação \(i\) serão dados por

• \(\hat y_i = \hat b_0 + \hat b_1 x_{i1}+\hat b_2 x_{i2}+\ldots+ \hat b_m x_{im}\).

Dentro desse contexto, definimos:

• resíduo \(i\) da regressão: \(\hat \epsilon_i=y_i - \hat y_i\), para \(i=1,\ldots,n\)

O resíduo \(i\) pode ser entendido como uma estimativa do valor de \(\epsilon_i\), para a observação \(i\), ou seja, o erro na observação \(i\) considerado no modelo. O gráfico a seguir ilustra essa noção: as linhas tracejadas vermelhas definem o resíduo associado a cada observação. A figura mostra em mais detalhe o resíduo de uma observação genérica \(i\).

10.2 Diagnósticos usuais em uma regressão linear

Para as análises que faremos a seguir considere o modelo modeloF que envolve os dados dos resultados femininos nos 100 m, considerando somente a variável explicativa ano. Para esse modelo, consideraremos somente as vencedoras das provas em cada um dos jogos realizados entre 1928 e 2012 (variável ord igual a 1). É utilizado o data frame olimpf, somente com os dados das mulheres, definido anteriormente.

## apresentando o diagrama de dispersão
plot(olimpf$ano[olimpf$ord==1],olimpf$t[olimpf$ord==1],col="blue",xlab="anos",ylab="tempo (seg)")
title(main="Vencedores olímpicos (Mulheres, 100 m)")

## para estimar os resultados do modelo modeloH:
modeloF<-lm(t~ano,data=olimpf, subset=(ord==1)) ## t = b0+ b1 ano + E
## para ver os resultados
summary(modeloF)

## 
## Call:
## lm(formula = t ~ ano, data = olimpf, subset = (ord == 1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43611 -0.10821  0.02756  0.11086  0.34801 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.996768   4.020691   9.948 9.68e-09 ***
## ano         -0.014598   0.002038  -7.163 1.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2273 on 18 degrees of freedom
## Multiple R-squared:  0.7403, Adjusted R-squared:  0.7259 
## F-statistic: 51.31 on 1 and 18 DF,  p-value: 1.139e-06

10.3 Visualização gráfica dos resíduos

O exame gráfico dos resíduos, ordenado pelas observações, ou por alguma variável explicativa de interesse, pode revelar algum padrão suspeito. Em função da premissa da normalidade dos erros e covariância zero entre erros (ou mesmo independência, numa premissa mais rigorosa), não será razoável observar padrões regulares nos resíduos apresentados nos gráficos. Isso certamente irá sugerir que as premissas não são válidas.

Vamos verificar a seguir o que ocorre com os resíduos do modelo modeloF:

## extração dos resíduos do modelo
resid<-residuals(modeloF)
plot(resid,main="Resíduos (modeloF)")
abline(h=0,lty=2)

Um exame visual não revela nenhuma anormalidade mais pronunciada, especialmente em função do número de observações que não é dos maiores. Há uma certa tendência para um leve formato de U na núvem de resíduos, mas nada que pareça anormalmente exagerado, podendo ser obra do acaso. Alguns testes, contudo, podem dar mais segurança a essa análise visual.

As duas figuras abaixo ilustram dois casos suspeitos. A figura da direita (caso suspeito 1) mostra um tipico comportamento de U (o contrário seria um U invertido), que indica uma má especificaçào do modelo (talvez a inclusão de algum componente não-linear na variável explicativa pode solucionar o problema). O segundo (caso suspeito 2) mostra uma situação em que claramente a variância parece estar aumentando com a ordem das observações, algo que conflita com a premissa de variância constante (chamada homoscedasticia), sugerindo heteroscedasticia (variância desigual ou não-homogênea).

10.4 Detecção de outliers e de observações muito influentes

A observação do gráfico com os resíduos pode não revelar claramente as observações que podem ser outliers (valores muito extremos, potencialmente errados ou fruto de alguma anormalidade). Outro problema são as observações que tendem a influenciar substancialmente as estimativas, na medida que essa observação é excluida ou incluida no conjunto de dados. Essas observações são chamadas muito influentes (em alguns casos são também outliers).

Há vários procedimentos para identificação de outliers e valores influentes. Dois procedimentos úteis são o teste de Bonferroni para outliers e a distância de Cook para valores muito influentes. Esses dois procedimentos são operacionalizados a seguir.

Teste de Bonferroni para outliers

Nesse teste, observações com valor-p abaixo de 5% (esse valor pode se alterado) são identificadas como possíveis outliers para que seja possível uma análise mais detalhada da observação de eventuais problemas. Detalhes técnicos são apresentados na documentação da função outlierTest do R (package car)

require(car)
## observações com valor-p menor que 5% são suspeitas
outlierTest(modeloF) 

## 
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferonni p
## 40 -2.187574           0.042961      0.85922

Foi identificada a observação 40 que pode ser observada no data.frame

olimpf[40,] 

##     ano   cid elev  mod sexo              nome ord pais     t
## 40 1988 Seoul   33 100m    F F Griffith-Joyner   1  USA 10.54

Que corresponde ao resultado extraordinário obtido por Flo-Jo, já discutido em outros tópicos.

Distância de Cook para detecção de observações influentes

A distância de Cook é uma medida que tenta caracterizar a importância de cada observação com relação a seu efeito sobre as estimativas, na medida que a observação é retirada ou recolocada no conjunto de dados. Mais detalhes técnicos sobre a distância de Cook são apresentados aqui. Há alguma divergência sobre o valor crítico que traz suspeitas a uma observação: alguns defendem que valores superiores a 1 são suspeitos, outros defendem que valores da distância de Cook superiores a 4/n já seriam suspeitos.

A seguir apresentamos o valor da distância de Cook para cada observação, com as referências com relação a possíveis níveis críticos definidos pelas linhas cinza (4/n) e vermelha (1) para os dados do modeloF.

plot(cooks.distance(modeloF),ylim=c(0,max(1,cooks.distance(modeloF))),
     ylab="Distância de Cook")
title(main="Distância de Cook (modeloF)")
abline(h=1,col="red",lty=2)
abline(h=4/length(olimpf$ano[olimpf$ord==1]),col=gray(0.2),lty=2)

Somente uma das observações, a primeira, revelou-se um pouco influente, pelo menos acima do limite \(4/n\) (mas inferior a 1). A posição da observação contribui para a influência por um efeito “alavanca” sobre a regressão (que é amplificado pelo método de quadrados mínimos). Apesar de não ser indicada como um possível “outlier” seu efeito sobre as estimativas é tido como superior à da própria observação 40, destacada pelo Teste de Bonferroni.

10.5 Detecção de não-normalidade e correlação nos resíduos

A premissa de normalidade é importante dentro do contexto das regressões lineares pois a partir dela são especificados muitos procedimentos de teste e estimativas de intervalos de confiança amplamente usados. Desvios da normalidade podem sugerir uma melhor especificação do modelo ou trasformações nas variáveis visando deixar os residuos mais próximos de “normais” (ex. transformação Potência e Box-Cox)

Outra premissa importante é a ausência de correlação nos erros, que pode ser aferida pela análise do comportamento dos resíduos, especialmente em casos em que as observações são indexadas pelo tempo. Uma grande área de estatística e da econometria que estuda as séries temporais ou séries de tempo trata de modelos que visam corrigem as dificuldades trazidas pela correlação.

Avaliação da aderência da premissa de normalidade

Uma análise gráfica envolvendo quantis teóricos e empíricos pode revelar, visualmente, desvios muito significativos da premissa de normalidade. Uma das várias funções que implementam esse gráfico é a qqPlot do package car, ilustrada a seguir.

## obtenção do qqPlot
require(car)
qqPlot(modeloF, main="QQ Plot")  

A situação desejável, na figura é ver os pontos dentro dos limites das linhas tracejadas que correspondem a um “envelope” que contém os pontos com 95% de probabilidade. No caso nenhum dos pontos está fora dos limites. A probabilidade pode ser alterada para 99%, por exemplo, usando a opção envelope=0.99.

Podemos, também, testar diretamente a normalidade através do teste Anderson-Darling, implementada pela função ad.test do package nortest, já utilizada em outro módulo. Um valor-p igual ou abaixo de 5% nesse teste sugere a rejeição da normalidade:

## obtenção do qqPlot
require(nortest)
ad.test(residuals(modeloF))

## 
##  Anderson-Darling normality test
## 
## data:  residuals(modeloF)
## A = 0.2403, p-value = 0.7423

Como o valor-p é bem superior a 5% o teste não sugere rejeição da hipótese de normalidade dos resíduos.

Diagnóstico de possível correlação nos resíduos

Há várias formas de testar correlação. O teste Durbin-Watson é uma opção tradicional, explicada na maioria dos livros de econometria. A implementação do teste é fácil, com apoio do package car. O teste examina as correlações do resíduo em \(t\) com resíduo em \(t-1\), até \(t-k\), onde \(k\) é definido na opção max.lag.

require(car)
durbinWatsonTest(modeloF,max.lag=5) ## testando até lag 5

##  lag Autocorrelation D-W Statistic p-value
##    1      0.09863548      1.655838   0.284
##    2      0.13262202      1.565803   0.314
##    3     -0.06027882      1.854293   0.948
##    4      0.12840070      1.243086   0.212
##    5     -0.08310435      1.659093   0.974
##  Alternative hypothesis: rho[lag] != 0

Os resultados não sugerem rejeição da hipótese de correlação zero pelo menos até o lag 5, o que suporta a validade das premissas.

Um resultado equivalente, mas apresentado em forma de gráfico, é operacionalizado pela função acf do package stats (já instalado com o R base):

require(stats)
acf(residuals(modeloF))

Um ou mais barras fora das linhas tracejadas (no lag 1 ou superior) sugerem rejeição da hipótese de ausência de correlação ou de covariância zero. Sempre no lag 0 a correlação será sempre igual a 1, não devendo ser alvo de preocupação. Nos outros lags, para o modelo analisado, não se observam valores fora das linhas tracejadas.

10.6 Detecção de heteroscedasticia (variância não constante)

Uma das premissas do modelo de regressão convencional é

• \(\epsilon_i\sim\text{N}(0,\sigma)\), ou seja, os erros tem distribuição Normal com \(E(\epsilon_i)=0\) e \(V(\epsilon_i)=\sigma^2\), para \(i=1,\ldots,n\).

Como \(\sigma^2\) é sempre o mesmo para qualquer \(i\), dizemos que a variância é constante. Essa premissa é atendida em um número grande de situações práticas, mas em outras pode não ser válida exigindo procedimentos específicos para tratar o problema. Em modelos mais avançados, há toda uma modelagem focada especificamente no comportamento da variância.

Para detecção de heterocedasticia, o nome dado para o caso em que a variância não é constante, há testes específicos. Um deles é o teste de Breusch-Pagan, que é implementado na função ncvtest do package car. Veja detalhes sobre esse teste aqui

## implementação do teste de Breusch-Pagan
require(car)
ncvTest(modeloF)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.4762464    Df = 1     p = 0.4901275

O valor-p do teste está associado a hipótese \(\H_0\): variância é constante. Como o valor-p encontrado é bem superior a 5%, não há evidência forte sugerindo a rejeição de \(\H_0\).

9.7 Detecção de multicolinearidade e inflação de variância

Quando há correlação alta dentro das variáveis explicativas isso pode ser trazer problemas sérios à estimação. Nessa situação a matrix \(\X^T\X\) usada na fórmula do estimador de quadrados mínimos pode ficar próxima de singular. Uma das consequências importantes disso é o aumento significativo da variância associada ao estimador (veja mais detalhes sobre isso aqui. Há o que se chama de “inflação de variância” nessa situação. Existem testes que podem auxiliar na detecção da intensidade do problema, medindo diretamente o fator de inflação de variância (fator 1 indica ausência de inflação).

Para exemplificar essa situação, que só será relevante quando tivermos 2 ou mais variáveis explicativas no modelo, considere o modelo modelo2 estimado no inicio do módulo, numa versão mais restrita que só inclui os vencedores da prova dos 100 m:

modelo2<-lm(t~ano+elev,data=olimpm,subset=(ord==1))

Para calcularmos os fatores de inflação de variância nesse caso, que está diretamente ligada à multicolinearidade, podemos utilizar a função vif do package car.

require(car)
vif(modelo2)

##      ano     elev 
## 1.011167 1.011167

Os fatores de inflação de variância calculados são bem próximos de 1, o que indica que o problema é praticamente inexistente nesse caso. Valores acima de 5 já são preocupantes. Acima de 10 são claramente problemáticos.

Para contrastar, considere a situação artificial abaixo, de uma regressão em que há uma correlação próxima de 0,986 entre as variáveis explicativas \(x\) e \(z\) (definida por construção). O fator de inflação de variância chega a mais de 35 nesse caso, capturando o problema de multicolinearidade.

require(car)
n<-100
x<-1:n
z<-x+rnorm(n,0,5)
cor(x,z)

## [1] 0.9858452

y<-2*x+3*z+rnorm(length(x),0,20)
modT<-lm(y~x+z)
vif(modT)

##        x        z 
## 35.57546 35.57546

confint(modT)

##                  2.5 %   97.5 %
## (Intercept) -11.948650 3.047258
## x             1.156548 2.692770
## z             2.331669 3.857179

10.8 Testes gerais de validade das premissas

Pesquisas mais recentes vem propondo metodologias gerais de teste das premissas. Um trabalho importante é

• Pena, E.A.; Slate, E.H. Global Validation of Linear Model Assumptions. Journal of the American Statistical Association (JASA) 101(473):341-354, 2006.

Os procedimentos sugeridos por esses autores estão operacionalizados na função gvlma do package gvlma. O trabalho citado deve ser consultado para detalhes específicos sobre os procedimentos de teste utilizados nesse package.

## extração dos resíduos do modelo
require(gvlma)
summary(gvlma(modeloF)) 

## 
## Call:
## lm(formula = t ~ ano, data = olimpf, subset = (ord == 1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43611 -0.10821  0.02756  0.11086  0.34801 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.996768   4.020691   9.948 9.68e-09 ***
## ano         -0.014598   0.002038  -7.163 1.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2273 on 18 degrees of freedom
## Multiple R-squared:  0.7403, Adjusted R-squared:  0.7259 
## F-statistic: 51.31 on 1 and 18 DF,  p-value: 1.139e-06
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = modeloF) 
## 
##                     Value p-value                   Decision
## Global Stat        6.3190 0.17656    Assumptions acceptable.
## Skewness           0.2152 0.64269    Assumptions acceptable.
## Kurtosis           0.2282 0.63285    Assumptions acceptable.
## Link Function      5.4007 0.02013 Assumptions NOT satisfied!
## Heteroscedasticity 0.4748 0.49079    Assumptions acceptable.

Ainda que o resultado tenha acusado um pequeno problema (não passou no teste “link function”), o modelo passou no “Global stat” e todos os outros testes. Para uma interpretação mais detalhada, veja o trabalho original.

11 - Estudo de caso: mãe fumante vs. peso do bebê / questões sobre o módulo

Considere o conjunto de dados no arquivo bebessemNA.csv correspondente a uma amostra de 1184 gestações na California (EUA) ocorridas há cerca de 30 anos, de mães fumantes e não fumantes, que tiveram um bebê do sexo masculino. Os dados são reais e vieram de um banco de dados de uma empresa gestora de plano de saude.

As variáveis desse conjunto de dados são descritas a seguir:

####
## Arquivo bebessemNA.csv - dados de 1184 gestações de meninos  
##                          de mães fumantes e não fumantes
####
## pesobebe   - peso do bebe no nascimento em kg
## gestacao   - semanas de gestação até o nascimento 
## idade      - idade da mãe no nascimento da criança em anos
## alturamae  - altura da mãe em cm
## pesomae    - peso da mãe em kg
## primipara  - é o primeiro filho (S/N)
## fuma       - mãe é fumante (S/N)
####

O conjunto de dados bebessemNA.csv deve ser “lido” no data frame bebes a partir dos comandos definidos a seguir:

### Leitura do arquivo
bebes<-read.csv2("http://ihbs.com.br/html/bebessemNA.csv")
head(bebes) ## listando as primeiras linhas do data frame

##   pesobebe gestacao idade alturamae pesomae primipara fuma
## 1     3.60    40.57    27       157    45.0         S    N
## 2     3.39    40.29    33       163    60.8         S    N
## 3     3.84    39.86    28       163    51.8         S    S
## 4     3.24    40.29    23       170    56.3         S    S
## 5     4.08    40.86    25       157    41.9         S    N
## 6     4.14    34.86    33       157    80.1         S    N

Questões sobre o data frame bebes

1. Para se familiarizar com os dados, faça uma boa caracterização de todas as variáveis quantitativas do arquivo, uma a uma, usando gráficos e medidas resumo que indiquem centralidade, dispersão, assimetria e curtose. Para as variáveis qualitativas, caracterize-as através de frequências absolutas e relativas.

2. Faça um diagrama de dispersão envolvendo as variáveis pesobebe e gestacao e identifique as 2 observações que são claramente erradas. Qual é o número das linhas dessas 2 observações? Elimine essas 2 observações do data frame bebes. (dica: o comando bebes<-bebes[-c(10,20),] recriaria novamente o data frame bebes a partir do data frame original sem as observações 10 e 20)

3. Utilize a função cor para obter uma matriz de correlação entre as variáveis quantitativas. Isso indicará quais são as variáveis que estão mais correlacionadas. Como a função só aceita variáveis quantitativas, use cor(bebes[,1:5]) para considerar somente as 5 primeiras colunas do data frame, correspondentes às variáveis quantitativas. Quais são as variáveis que parecem ser mais correlacionadas?

4. Estime uma regressão do tipo \(y_i=b_0+b_1 x_i + \epsilon_i\), com \(y_i\) representando a variável pesomae e \(x_i\) representando a alturamae. Represente a curva de regressão dentro do diagrama de dispersão. Com base nos dados dessa regressão, qual seria a média condicional estimada do peso das mães que tem altura 150 cm? e altura 170 cm? Há evidência que suporte a noção de que mães mais altas são mais pesadas? Em média quanto o peso da mãe aumenta para cada 10 cm de aumento na altura? Encontre o intervalo de confiança para esse último resultado, baseado na regressão. Encontre o intervalo de confiança para o peso médio (teórico) e o intervalo de previsão para o peso de mães com 170 cm.

5. Use o comando boxplot(pesobebe~fuma,notch=TRUE,col=“bisque”) para construir um boxplot mostrando a distribuição dos pesos dos bebes para mães fumantes e não fumantes. A opção notch=TRUE mostra o intervalo de confiança para a mediana estimada no boxplot (o intervalo é representado pela abertura que inicia no retângulo externo e finaliza na mediana, como um V invertido de cada lado). O que esse resultado está sugerindo com respeito a um possível efeito do uso do cigarro pela mãe no peso dos bebes?

6. Estime uma regressão tipo \(y_i=b_0+b_1 x_{i1}+b_2 x_{i2} + b_3 x_{i3} + b_4 x_{i4} + \epsilon_i\), com \(y_i\) representando a variável pesobebe, e \(x_{i1}\) representando alturamae, \(x_{i2}\) representando o pesomae, \(x_{i3}\) representando idade, \(x_{i4}\) representando gestacao. Pelo resultado da estimativa dos parâmetros e valores-p, quais variáveis parecem ter efeito estatisticamente significativo sobre o peso do bebe? interprete a estimativa do efeito marginal associado a cada variável, e seu intervalo de confiança.

7. Considerando o mesmo modelo da questão anterior, elimine as variáveis que não apresentaram efeito significativo, e inclua no modelo uma variável “dummy”, representada por \(z\), que assume valor 1 quando a mãe é fumante e 0 quando a mãe é não-fumante, de forma que seu parâmetro indique um efeito no intercepto. Estime os parâmetros da regressão. O efeito da variável \(z\) é estatisticamente significativo? qual é a interpretação do parâmetro associado a \(z\) estimado e seu intervalo de confiança? Use um teste para verificar se a hipótese de normalidade do resíduos se verifica para essa regressão.

8. Com base no modelo estimado na questão 7, qual seria o intervalo de confiança e o intervalo de previsão estimado para o peso de um bebe nascido de uma mãe com 80 kg, 170 cm, com gestação de 42 semanas, fumante? responda também considerando os mesmos dados de altura, peso e gestação, mas para uma mãe não-fumante.

9. Com relação ao modelo utilizado na questão 7, use o diagnóstico de outliers (tópico 9) para identificar outras possíveis observações que deveriam ser eliminadas do data frame. Verifique se os resíduos dessa regressão passam no teste de normalidade. Calcule o fator de inflação de variância associado a cada variável. Há problemas aparentes de multicolinearidade?

Questões sobre outros tópicos:

1. Use o método dos quadrados mínimos para derivar as fórmulas para estimar os parâmetros do modelo: \(y_i=b_0+b_1 x_{i}+b_2 x_{i}^2 + e_i\) (equação do segundo grau). (Dica: repita o procedimento para o caso da regressão linear mostrado nos slides para entender bem o procedimento. Nesse caso existirão 3 condições de primeira ordem).

2. A análise desenvolvida com os dados das mulheres nos 100 m mostrou algum impacto após a adoção de testes mais rigorosos anti-doping depois de 1988? explique usando os resultados do modelo no tutorial que utilizou variáveis dummy.

3. Modelos lineares de regressão podem ser utilizados para situações claramente não-lineares?

4. Quais são as principais premissas da regressão linear múltipla estudada e quais são os diagnósticos principais que servem para verificar se essas premissas são válidas?