Teste de Hipótese

Um jeito simples de escolher as variáveis para sua regressão

É interessante que você conheça antes de ler:

• O que são variáveis em estatística
• Os principais tipos de distribuição de variáveis
• Como funciona MQO para regressão
• Aplicação não específica de testes de hipótese

Mas você ainda consegue aproveitar a leitura mesmo que venha apenas com a vontade de aprender!

Photo by Markus Winkler on Unsplash

1 Introdução

Este artigo tem a intenção de viabilizar o aprendizado em temas mais específicos que não estão tão disponíveis em português quanto estão em outras línguas. Espero que este material seja útil e didático o suficiente para suas necessidades com modelos de regressão.

Se você é ou está se preparando para se tornar um cientista de dados, um dos primeiros tipos de modelos que você vai estudar são aqueles de regressão linear ou regressão logística. Embora muitas vezes tenham propósitos diferentes, existe uma simples e fundamentar semelhança entre todos eles: estes modelos são formulados por uma “função de previsão”.

No fim das contas, regressão linear, logística, polinomial e todas as outras são representadas por funções! As regressões são essencialmente uma maneira de ajustar uma forma funcional conhecida (polinomial, linear, logística, etc.) à dispersão de um conjunto de variáveis de interesse, de modo encontrar o melhor modelo de regressão é um trabalho de encontrar as melhores variáveis e o melhor tipo de função para ajustar à elas:

Numa função de previsão, nós vamos sempre nos deparar com uma variável explicada y sendo prevista por variáveis explicativas x₁, x₂,…,xₙ. De modo que, por exemplo:

y = β₀ + β₁x₁ + β₂x₂ + … + βₙxₙ

Formalmente, se usa este símbolo para o valor estimado do coeficiente:

Famoso i-ésimo “beta chapéu”

E este outro símbolo para o valor “verdadeiro populacional” dele:

Famoso i-ésimo “beta”

Neste artigo, vamos usar apenas a última forma, mas para se referir ao coeficiente estimado. O termo u representa os erros de ajuste do modelo.

A equação acima mostra um exemplo mais simples de regressão linear, mudar de modelo, significa mudar o formato da função, mas ainda estaremos lidando com equações! Sobre esta que mostramos acima, podemos dizer que cada um destes βᵢ representam o impacto que cada respectivo xᵢ tem em y. Ou seja, quanto mais distante de zero for o valor de βᵢ, mais relevante deve ser o xᵢ, mas isto pode nos gerar um problema:

1.1 Problema de escala

Digamos que estamos tentando prever um y que oscila entre vinte e cinquenta, e para isso, vamos usar um x₁ que oscila entre oitenta e quatro mil e um x₂ que oscila entre zero e cinco, representando matematicamente, isto seria:

{y ∈ ℝ: 20 ≤ y ≤ 50}
{x₁ ∈ ℝ: 80 ≤ x₁ ≤ 4000}
{x₂ ∈ ℝ: 20 ≤ x₂ ≤ 50}

Uma coisa que podemos dizer sobre estas variáveis, é que possuem uma diferença muito grande entre seus valores mínimos e máximos, sendo 30, 3920 e 5 respectivamente para y, x₁ e x₂. Sendo assim, podemos dizer que elas apresentam intervalo de escala diferentes.

Agora digamos que nós vamos ajustar um modelo de regressão linear com estas variáveis. Pelo Método de Minimização de Quadrados (MMQ) que é usado para obter os valores β₁ e β₂ respectivos para x₁ e x₂, é esperado que β₂ > β₁, ao interpretar os coeficientes, alguém poderia pensar que pelo β₂ ser maior, então x₂ deve ter mais poder preditivo para explicar y, mas isto não é necessariamente verdade. Vamos ver o motivo usando código de R.

dados <- datasets::LifeCycleSavings

# Vamos dar uma olhada nas variáveis de nosso conjunto de dados
# e suas respectivas escalas usando `range()`
colunas <- names(dados)
for (col in seq_along(colunas)) {
  labs <- paste0(colunas[[col]], c("_min", "_max"))
  paste0(labs, ": ", range(dados[[col]])) |> print()
}

# Output:
#    "sr_min: 0.6"      "sr_max: 21.1"
# "pop15_min: 21.44" "pop15_max: 47.64"
# "pop75_min: 0.56"  "pop75_max: 4.7" 
#   "dpi_min: 88.94"   "dpi_max: 4001.89"
#  "ddpi_min: 0.22"   "ddpi_max: 16.71"

Aqui temos um caso parecido com o exemplo que vimos antes, as variáveis pop15, dpi e pop75 apresentam escalas de intervalo consideravelmente distintas. Agora vamos renomear estas variáveis para ter os nomes y, x₁ e x₂, depois vamos ver como esta diferença se apresenta, comparando pop15 visualmente com as demais.

dados <- dados[, c("pop15", "dpi", "pop75")]
names(dados) <- c("y", "x1", "x2")

# Visualizando pop15 vs dpi:
library(ggplot2)
lim <- c(0, 4000)
ggplot(dados, aes(y=y, x=x1)) +
  geom_point() + theme_classic() +
  coord_cartesian(x=lim, y=lim)

Note que a variação só é perceptível ao longo do eixo de x₁, onde a escala é maior. Se quisermos usar esta variável para explicar y, considerando que é uma boa variável para esta tarefa, teríamos um valor de β₁ muito pequeno, porque estamos dando grandes saltos em x₁ para obter uma pequena variação em y.

O coeficiente β₁ teria que traduzir a variação de x₁ para y levando-a à escala de y multiplicando por um número bem próximo de zero. Vamos ajustar uma regressão linear sem intercepto, e obter o valor do coeficiente apenas para exemplificar:

# Coeficiente de y = beta_1.x1
lm(y~x1+0, data=dados)$coefficients
#  x1
#  0.01468412

Agora vamos repetir o processo, mas comparando y com x₂:

# Visualizando pop15 e pop75:
lim <- c(0, 50)
ggplot(dados, aes(y=y, x=x2)) +
  geom_point() + theme_classic() +
  coord_cartesian(x=lim, y=lim)

# Coeficiente de y = beta_2.x2
lm(y~x2+0, data=dados)$coefficients
#  x2
#  10.15053

O mesmo efeito ocorre aqui, mas agora percebemos melhor apenas a variação em y, onde neste caso, a escala é maior. O coeficiente teve um valor muito distante de zero, simplesmente porque a escala de x₂ é consideravelmente inferior à de y.

2 Entendendo o teste de hipótese

Agora vamos nos tornar capazes de generalizar o conceito por trás de qualquer teste de hipótese, conceituando com uma “regressão linear sem inclinação”, isto significa dizer que vamos prever o valor de y usando apenas a média dele mesmo:

ggplot(dados, aes(y=y, x=seq_along(y))) +
  geom_point() + geom_hline(yintercept=mean(dados$y)) +
  theme_classic() + xlab("Índice")

Neste caso particular, parece que a média não é uma boa métrica para prever o valor de y, parece que todos os dados estão distantes da média, num bom ajuste nós veríamos todas as observações se aglomerando ao longo desta linha, mas isto não vai ser um problema, estamos modelando apenas para fins didáticos. Só para deixar claro que ainda estamos falando de uma regressão, aqui está a forma de função deste modelo:

y = β₀

Onde β₀ seria simplesmente o valor médio de y. Os ingredientes que precisaremos para fazer um teste de hipótese seriam apenas:

1. O valor estimado de β₀
2. Uma métrica que controle a escala da variável em sua distribuição, e levando em conta a incerteza da amostra

O primeiro ingrediente já foi obtido, já sabemos a média de y, mas agora precisamos medir a escala da variável com alguma métrica:

2.1 Variância

Vamos dar uma olhada em como a variância se comporta quando medimos a variância em uma variável. Primeiro medimos a dispersão de cada observação, que será simplesmente sua distância em relação ao valor médio da variável:

y <- dados$y
y_bar <- mean(dados$y)
plot(y, pch=20, xlab="Índice")
abline(h=y_bar)
for (i in seq_along(y))
  segments(x0=i, x1=i, y0=y_bar, y1=y[[i]])

Neste caso, a linha horizontal marca a média, enquanto que a linhas verticais marcam a distância entre a média e cada observação individualmente, o somatório destes valores seria o valor total das distâncias. Considerando n a quantidade de observações em nossa variável, e i a i-ésima observação considerada (ao longo do índice), isto seria representado matematicamente desta forma:

D representa o somatório das distâncias entre as observações e a média

Note que se uma observação xᵢ for maior que a média, vamos obter um valor positivo, e no caso contrário, um valor negativo. Fizemos progresso, mas deste jeito as distâncias negativas vão anular as distâncias positivas, vamos aplicar uma alteração para resolver o problema:

D ao quadrado representa o somatório anterior, mas o valor de cada distância é elevado ao quadrado

O maior problema foi resolvido, nenhuma distância será negativa, mas agora surgem mais duas observações:

1. Outliers vão afetar muito mais o valor de D², já que estamos elevando as distâncias ao quadrado;
2. Esta medida vai ser sensível à quantidade de observações, quanto maior o n, mais valores teremos que somar.

O primeiro pode ou não ser um problema, dependendo do seu desafio de regressão, mas o segundo vai ser um problema para qualquer caso, se alterarmos o tamanho dos dados de treinamento, nosso modelo deve sofrer com a mudança da métrica de escala D. Para resolver isso vamos aplicar uma simples mudança:

A variância representada pelo S-quadrado é D ao quadrado dividido pelos graus de liberdade n-1

Temos que dividir por n-1 para levar em consideração os graus de liberdade (uma explicação detalhada do motivo). Agora temos um S² que nada mais é que a média das distâncias à média ao quadrado, esta é a métrica que chamamos de variância. Está quase perfeito para o que precisamos, só falta alguns detalhes:

2.2 Desvio padrão

Conseguimos medir a escala de nossa variável com esta métrica, mas ainda tem um detalhe meio paradoxal para resolver, nós estamos tentando medir a escala de nossa variável, mas estamos usando uma métrica que está fora da escala destas distâncias que medimos. É daí que surge a necessidade de usar outra medida chamada desvio padrão:

O desvio padrão representado po S é representado pela raiz quadrada da variância S-quadrado

2.3 Erro padrão

O último detalhe é levar em consideração ao nível de incerteza que nos deparamos a lidar com um conjunto de dados que representa a população, uma amostra maior vai apresentar menos variabilidade nas propriedades da amostra (média, variância, etc.) quando comparados com a população. Quanto menor a amostra, mais incerteza, por isto nessa situação, usamos o erro padrão para manter isto sob controle:

O erro padrão representado por Se é igual ao desvio padrão dividido pela raiz quadrada do tamanho amostral

Agora quanto maior o valor da amostra, menor será nosso erro obtido, esta métrica cumpre tudo que precisamos para usar no teste de hipótese. Se você ainda não entendeu porque precisamos destes dois ingredientes, não se preocupe, porque vamos visualizar a intuição por trás disto mais tarde no tópico 3.

2.4 Como se faz um teste de hipótese?

Com as etapas anteriores, nós conhecemos os “ingredientes” que precisamos usar para fazer um teste de hipótese e suas respectivas propriedades. Os parâmetros βᵢ são multiplicadores de variáveis em nosso modelo ou constantes que estão presentes, independentemente de qualquer variável, se a possibilidade de seu valor ser igual à zero for estatisticamente significante, quer dizer que temos uma boa certeza estatística de que este parâmetro pode ser removido, ou que devemos mantê-la, caso contrário.

Como exemplo, vamos considerar um modelo polinomial de segundo grau:

y = β₀ + β₁x₁ + β₂x₁² + β₃x₂

Se tivermos uma certeza estatística de que β₂ possa ser igual a zero, podemos usar isso como um argumento sólido para usar um modelo linear em vez de polinomial. Podemos também usar o teste em β₃ para ter um bom argumento para retirar ou manter x₂ das variáveis explicativas, ou em β₀ para incluir ou não um intercepto.

Existem diferentes tipos de testes (bilateral ou unilateral), que dependem da hipótese definida; e também existem diferentes estatísticas de teste (t, z, χ² [chi quadrado], etc.) que dependem da distribuição e do tamanho da amostra que vai ser testada.

A maneira que a MQO e outras metodologias de regressão usam para ajustar a curva de regressão aos dados tem a propriedade de tornar a distribuição dos erros sempre similar à uma normal, independentemente da qualidade deste ajuste. Abaixo vemos um bom ajuste da regressão linear e sua distribuição de erros no exemplo 1:

ggplot(women, aes(y=weight, x=height)) + 
  geom_point() + theme_classic() +
  geom_smooth(method="lm", formula=y~x, se=FALSE, color="black") +
  ggtitle("Exemplo 1: Ajuste da regressão linear")

resids <- lm(weight~height, data=women)$residuals
ggplot(data.frame(resids), aes(resids)) + 
  geom_density(linewidth=1.2) + theme_classic() + 
  ggtitle("Exemplo 1: Distribuição dos erros")

E no exemplo 2, um mal ajuste da regressão e sua distribuição de erros:

lmtest <- lm(qsec~drat, data=mtcars)
ggplot(mtcars, aes(y=drat, x=qsec)) + 
  geom_point() + theme_classic() +
  geom_smooth(method="lm", formula=y~x, se=FALSE, color="black") +
  ggtitle("Exemplo 2: Ajuste da regressão linear")

resids <- lm(qsec~drat, data=mtcars)$residuals
ggplot(data.frame(resids), aes(resids)) + 
 geom_density(linewidth=1.2) + theme_classic() + 
 ggtitle("Exemplo 2: Distribuição dos erros")

Por isso, quando fazemos um teste de hipótese, partimos do pressuposto que os erros estão normalmente distribuídos. Por isso, a estatística de teste que usamos neste caso é:

• T de Student, que é mais adequada para distribuições similares à normal, usamos quando a amostra é pequena, ou;
• Z (normal padrão), usamos quando quando temos uma amostra grande.

O nome da estatística muda, mas a maneira de calcular é a mesma (mesma fórmula). Embora a fórmula da estatística de testes seja sempre a mesma, a fórmula do erro padrão e do valor estimado do parâmetro βᵢ sempre vai mudar, dependendo da forma funcional adotada. Vamos exemplificar com a fórmula para obter o erro padrão de uma regressão linear simples (com uma variável de cada lado da equação):

Esta fórmula de erro padrão na regressão linear é mais complexa do que quando calculamos o erro padrão na regressão constante (com uma variável e uma constante)

ŷᵢ é o valor previsto pela regressão; x̄ é o valor médio da variável x.

Não precisa entender a fórmula de cara, mas fica perceptível que à medida que adicionamos complexidade ao nosso modelo de regressão, o valor dos parâmetros e seus respectivos erros padrão ficam cada vez mais complexos de se obter.

O lado bom é que as estatísticas de teste manterão sempre a mesma fórmula, e não recebem nenhum grau de complexidade enquanto sofisticamos nossos modelos de regressão.

3 Testando

Vamos obter uma estatística T ou Z dependendo do tamanho da amostra, alguns autores usam 30, outros usam 100 observações como o valor máximo para usar a estatística t. Mas a maneira de chegar à ela é a mesma:

A estatística de teste é igual à diferença entre o valor obtido e o testado de beta-i, dividido pelo erro padrão do próprio beta-i

Neste caso, μ é o valor de parâmetro para o qual estamos comparando o βᵢ na hipótese, como estamos sempre testando para βᵢ=0, este termo vai sempre desaparecer na equação. Se quiser testar a hipótese de βᵢ=20, então usaremos a fórmula com 20 no lugar de μ.

O resultado é uma estatística de teste, que não responde a nossa pergunta diretamente, mas nos deixa mais próximos, quando chegamos neste ponto, conhecemos o tamanho amostral n, o coeficiente estimado βᵢ e o erro padrão deste coeficiente Se(βᵢ). Isto nos habilita a formular uma distribuição dos erros deste coeficiente, vamos visualizar esta distribuição revisitando os exemplos que vimos antes:

ggplot(women, aes(y=weight, x=height)) +
  geom_smooth(method="lm", formula=y~x) + geom_point() +
  theme_classic() + ggtitle("Exemplo 1: Distr. coeficientes possíveis")

Perceba que os pontos ficam muito próximos da linha de regressão, e portanto a área cinza fica muito estreita, vamos ver no mal ajuste, esta área se alargar de maneira considerável:

plot <- ggplot(mtcars, aes(y=qsec, x=drat)) +
  geom_smooth(method="lm", formula=y~x) + geom_point() +
  theme_classic() + ggtitle("Exemplo 2: Distr. coeficientes possíveis")
plot

Isto acontece exatamente por causa da variabilidade dos dados em torno da regressão. A área cinza que se forma tem esse formato de ampulheta, pois ela mostra todas as inclinações que possuem alguma significância estatística relevante, junto com os possíveis interceptos, vamos adicionar mais duas linhas a este gráfico para exemplificar:

plot + 
  geom_abline(intercept=12.3793, slope=1.507848, linewidth=1) +
  geom_abline(intercept=21.12573, slope=-.89825, linewidth=1) +
  ggtitle("Exemplo 2: Ajustes possíveis nas distr.")

Dá pra ver pelo código que em uma linha eu diminuí o intercepto e aumentei a inclinação, e fiz o oposto na segunda linha. Como se tratam de distribuições de probabilidade, existem limites superiores inferiores tanto para o intercepto quanto para a inclinação da linha.

Neste exemplo eu exploro valores próximos aos limites estatisticamente significantes (convencionalmente 95%), assim fica mais fácil perceber como a escolha de diferentes valores estatisticamente significantes pode preencher o espaço nesta curva. Se você não entendeu o que eu quis dizer, não precisa se preocupar, porque agora nós vamos estender isto junto com a:

3.1 Estatística de teste

Vamos imaginar uma linha contínua infinita, quando estimamos o valor de um coeficiente de regressão, estamos marcando um ponto nesta linha. Mas como estamos trabalhando com uma amostra para chegar neste valor e não conhecemos a população inteira, não temos total certeza sobre a posição exata do ponto que marcamos nesta linha.

Para nos ajudar com este problema, vamos desenhar um gráfico que representa o nosso nível de certeza:

Repare que, embora acreditemos que outros valores são possíveis, temos mais certeza (parte mais alta da curva) no valor médio (ou valor estimado), quanto mais distante da média, mais improvável nos parece. O critério para decidir o que é estatisticamente significativo é arbitrário, mas convencionalmente se usa 95% da área deste gráfico:

Perceba que esta área cinza marcada possui duas “fronteiras”, e também podemos marcar a localização delas em nossa linha infinita. Consideramos que qualquer valor entre estas duas fronteiras (agora marcadas em vermelho) é estatisticamente significante, e portanto é possível e plausível:

Como vimos no começo, a variância pode mudar drasticamente de uma variável para outra, medir essa área e encontrar suas fronteiras para cada distribuição pode ser complicado, por isso, usamos uma medida padronizada que se ajusta à todas as distribuições possíveis, que é o erro padrão.

Na distribuição z, a fronteira que marcamos acima fica sempre à 1,96 × Se(βᵢ) de distância de βᵢ, já na distribuição T, esta distância depende da quantidade de observações, por isso, geralmente se usa uma tabela para encontrar o valor. Considerando que xᵢ tem n>100, e portanto estamos usando uma curva normal padrão, podemos obter o valor destas fronteiras conhecendo o valor de βᵢ e Se(βᵢ):

Nossa estatística de teste mostra a posição do valor que queremos testar em quantidade de desvios padrão, se encontramos Z=2,5, significa que o zero está a 2,5 × Se(βᵢ) de distância do valor do nosso coeficiente βᵢ estimado. Neste caso, como estamos tratando da distribuição Z e usando o nível de significância convencional, o zero está muito distante do valor estimado para o considerarmos plausível.

Já que a região que consideramos os valores plausíveis não contém o zero neste caso hipotético, dizemos que rejeitamos a hipótese nula de que βᵢ=0, e portanto, escolhemos manter este coeficiente e sua respectiva variável xᵢ em nosso modelo.

3.2 O teste de hipótese é sempre confiável?

Antes de responder a pergunta, vamos lembrar o que são resíduos da regressão: Ao preencher a fórmula da regressão com valores já conhecidos das variáveis explicativas, obtemos um valor estimado da variável resposta ligeiramente diferente do valor observado. Os resíduos são as diferenças entre os valores observados e os valores estimados (usando a notação canônica desta vez):

uᵢ = yᵢ - ŷᵢ

Portanto, a primeira situação em que não conseguimos confiar no teste de hipótese se resume em: “Se ŷ não estiver correlacionado com y, então u estará”. Podemos verificar isto nos dois exemplos que exploramos antes:

# Exemplo 1 (Bom ajuste)
y <- women$weight
u <- resid(lm(weight~height, data=women))
cor(y, u)
#> 0.01...

# Exemplo 2 (Mal ajuste)
y <- mtcars$qsec
u <- resid(lm(qsec~drat, data=mtcars))
cor(y, u)
#> 0.99...

Quando encontramos a situação presente no exemplo 2, significa que nosso modelo não é apropriado para prever y, e portanto, nossos coeficientes estimados não têm razão de o ser, neste caso, nossa estatística não é confiável porque nosso coeficiente não é confiável. A solução pode estar em um ou uma combinação destes pontos:

• Trocar variáveis existentes por outras com maior poder preditivo
• Adicionar mais variáveis
• Mudar a forma funcional da função de regressão

A segunda situação em que não podemos confiar no teste de hipótese é quando não podemos confiar na medida do desvio padrão. Isto pode acontecer quando a variabilidade dos erros não é constante ao longo da curva de regressão, chamamos isto de heterocedasticidade:

Neste caso, nossa estatística de teste pode ser válida apenas para uma porção dos dados, e de maneira geral uma métrica ruim. Mas neste caso, isto não quer dizer necessariamente que temos um mal modelo, mas que teremos que:

• usar outras métricas robustas à heterocedasticidade;
• ou talvez fazer alguma manipulação nos dados para minimizar este efeito (como a transformação de Box-Cox).

Um teste estatístico pode ser usado para detectar a presença deste efeito na sua regressão, alguns exemplos são os de Breusch-Pagan ou o de White, ambos também são testes de hipótese, mas testam a hipótese da presença de heterocedasticidade no modelo a partir de seus resíduos u.

3.3 Recapitulando

Neste artigo aprendemos:

• o conceito de significância estatística;
• o conceito do que é um teste de hipótese;
• sua aplicação ao escolher variáveis para uma regressão;
• os casos em que o teste de hipótese convencional não é confiável;
• maneiras de lidar com esta situação

Espero que este material seja capaz de te ajudar com seu objetivo com modelos de regressão, e que não se limite apenas a este conteúdo. Convido a verificar as referências para mais informações, obrigado pela leitura!

Referências

• Teste de Hipóteses (ufsc.br)
• Distribuição amostral da variância (ufsc.br)
• Ordinary least squares — Wikipedia
• O que é distribuição normal (Z) (rapidtables.org)
• Distribuição t de Student — Wikipédia, a enciclopédia livre (wikipedia.org)
• Heteroscedasticity (r-project.org)
• aulasEstDistrTdeStudent.pdf (professornemer.com)
• Understanding the Standard Error of a Regression Slope — Statology