Apenas um lembrete sobre as Distribuições e a Teoria do Limite Central

A teoria do limite central garante que, se a amostra for grande o suficiente (n ≥30) ou se a população é normalmente distribuída, a distribuição amostral das médias tem a média igual à média da população.

Entretanto, precisamos ter um certo nível de confiabildiade de nossas estimativas baseadas em distribuições amostrais (média amostral, coeficientes de correlação, betas de regressões, dentre outros). Neste contexto, surge o conceito de desvio-padrão da distribuição amostral da média que será empregado em diversos momentos na construção do chamado erro padrão da média e seus derivados.

Essa constatação sobre as distribuições amostrais, juntamente com o erro padrão, nos ajuda a investigar o quão confiáveis são as estimativas que encontramos. Veremos esse arcabouço empregado em três técnicas estaísticas:

Intervalo de Confiança
Teste de Hipótese
Análise de Correlação Linear e de Regressão Linear Simples

Em muitas situações da vida real, estamos lidando com amostras, ou seja, o desvio padrão da população (σ) não é conhecido…. neste caso não podemos usar a distribuição normal padrão para caracterizar a distribuição amostral da média.

Para uma variável aleatória que é normalmente distribuída (aproximadamente normalmente distribuída) ou n ≥30, a variável média amostral comporta-se tal qual outra distribuição de probabilidade, a distribuição t-Student. Vimos a distribuição t na Prática 3; relembrando

“A distribuição t é uma família de curvas, cada uma determinada por um parâmetro chamado de graus de liberdade n-k (ou df, degrees of freedom); sendo”n" o tamanho da amostra e “k” o número de coeficientes de análise. Assim, quando usamos a distribuição t para estimar uma média populacional (apenas um coeficiente), os graus de liberdade são iguais ao tamanho da amostra menos um".

ATENÇÃO: a validade de qualquer método de estimativa aumenta quando se utiliza uma estatística amostral que seja não tendenciosa e de variância mínima (aleatória e representativa). Uma estatística é não tendenciosa se não superestima ou subestima o parâmetro populacional.

Intervalo de Confiança

Embora possamos supor uma estimativa pontual, por exemplo calcular uma média amostral e indicá-la como a possível média populacional, existe uma diferença entre a estimativa pontual e o valor real do parâmetro. Essa diferença é chamada de erro de amostragem ou erro amostral.

A boa notícia é que podemos calcular o valor máximo para o erro quando definirmos um nível de confiança e soubermos a distribuição amostral da variável.

O Intervalo de Confiança, portanto, é o procedimento de se calcular uma estimativa pontual, então, usá-la para estabelecer os limites superior e inferior para o parâmetro da população, considerando um nível de confiança, o erro padrão e uma distribuição de probabilidades (no caso a t-Student).

Nível de confiança (c): é a probabilidade de uma estimativa intervalar (não mais um ponto específico, agora um intervalo de pontos) contenha o parâmetro populacional verdadeiro, supondo que o processo de estimação é repetido um grande número de vezes. Geralmente, 90%, 95% ou 99%, o mais comum é usar 95%.

Margem de erro (E): é a maior distância possível entre a estimativa pontual e o valor do parâmetro que ela está estimando. Sua fórmula é composta por:

t_c: t crítico ou t tabelado. Usando a função qt() (vimos na Prática 3)

s: desvio padrão amostral

n: número de elementos na amostra (ou tamanho da amostra)

\[ E=t_c\frac{s}{\sqrt{n}} \]

Vamos fazer um intervalo de confiança para a média da variável consumo de água per capita em m³/ano (cons1). Na Prática 3 plotamos um Qqplot e verificamos que ela possui certa aproximação à distribuição normal. Independente dessa verificação, conseguimos atender ao segundo item obrigatório, ter n ≥30, nosso n = 506.

media <- mean(dados_agua_sna$cons1)
media

[1] 46.62771

s <- sd(dados_agua_sna$cons1)
s

[1] 16.50024

n <- nrow(dados_agua_sna)
n

[1] 506

E <- qt(p=1-0.05/2, df=n-1)*s/sqrt(n)
E

[1] 1.441137

linf <- media - E
linf

[1] 45.18658

lsup <- media + E
lsup

[1] 48.06885

IC95 <- c(linf,lsup)
IC95

[1] 45.18658 48.06885

Interpretação: Com 95% de confiança, podemos dizer que consumo de água per capita em m³/ano médio dos municípios do Estado de São Paulo está entre 45,18 e 48,06 m³/ano. Percebam, usamos uma estimativa pontual (média = 46.63) para encontrar uma estimativa intervalar para a verdadeira média populacional (IC95: 45.19, 48.07).

Teste de Hipótese

Um teste de hipótese é um processo que usa estatísticas amostrais para testar uma afirmação sobre o valor de um parâmetro populacional, ou seja, quando podemos considerar que um resultado é genuíno, portanto, não é fruto do acaso.

Para testar um palpite sobre os dados em análise precisamos definir:

Hipótese nula (H₀): hipótese da igualdade ou a hipótese que prevê que os dados não exibirão nenhum efeito, por exemplo, testar se o parâmetro é igual a 0 (igualdade, sem efeito).

Hipótese alternativa (H_a): é o complemento da hipótese nula. Hipótese da desigualdade ou a hipótese que prevê que os dados exibirão efeito, por exemplo, testar se o parâmetro é diferente de 0 (desigualdade, com efeito).

O teste sempre é feito sobre a hipótese nula, que pode ser rejeitada ou não rejeitada. Caso o teste indique pela rejeição da hipótese nula, então, a hipótese alternativa é aceita como verdadeira sugerindo a possível existência de efeitos.

Como sua decisão de rejeitar ou não rejeitar é baseada em uma amostra, você deve aceitar o fato de que sua decisão pode estar incorreta. O erro mensurável em um teste de hipótese é denominado:

Erro do Tipo I: quando rejeitamos a hipótese nula, mas ela é verdadeira, ou seja, quando não existe efeito e aceitamos que existe.

O valor da probabilidade de cometer um Erro do Tipo I em um teste de hipóteses é conhecido como Nível de Significância e é representado pela letra α. O nível de significância é o complemento do nível de confiança, ou seja, se quero ter 95% de confiança isso significa um nível de erro do tipo I de 5%. Os níveis de significância (α) mais utilizados são de 5%, 1% e 0,1%.

Os testes de hipótese são realizados em função das distribuições de probabilidade, por meio de estatisticas de teste. O valor da estatística de teste será calculado (t calculado) e designado por estatísticas de teste padronizada (como z, t, x², dentre outros) e comparado com um t crítico (ou t tabelado).

Se a estatística teste calculada for maior do que a estatística tabelada na distribuição, então, existem evidências para rejeitar a hipótese nula (t calculado > t crítico).

\[ t_{calculado} = \frac{\bar{x}-\mu}{\frac{s}{\sqrt{n}} } \]

É possível, por fim, saber qual a probabilidade de se obter uma estatística de teste igual ou mais extrema do que aquela observada na amostra (t calculado); a essa probabilidade damos o nome de Valor-p (p-value).

Em outras palavras, de maneira mais simplificada, se estamos testando, por exemplo, a média ser igual a 0 (hipótese nula), então, o valor-p corresponde à probabilidade de que isso seja possível ou simplesmente obra do acaso.
Assim, quanto menor for o valor-p, maior a evidência para rejeitar a hipótese nula, pois, um valor-p muito baixo sugere uma probabilidade muito baixa da ocorrência de H₀. A convenção é usar como referência, para rejeitar H₀, um valor-p menor do que 5%.

ATENÇÃO: a estatística de teste e o valor-p são informações complementares; sem a estatística de teste não conseguimos saber a área de probabilidade para calcular o valor-p. Enquanto o t calculado e o t crítico são pontos na distribuição, o valor-p é a área abaixo da curva da distribuição representando a probabilidade de obter um valor igual ou mais extremo ao t calculado.

Na figura 1, segue o gráfico ilustrativo. O exemplo é unicaudal apenas para facilitar o entendimento.

Figura 1 - Teste de Hipótese

Vamos fazer um exemplo ilustrativo. Suponha que técnicos da agência reguladora do serviço de abastecimento de água afirmam que o consumo de água per capita em m³/ano (cons1) nos municípios paulistas no ano de referência dos dados é de 44 m³/ano. Há evidência suficiente para rejeitar a afirmação da agência considerando o nível de significância α = 0,05?

H₀: cons1 = 44 m³/ano

H_a: cons1 ≠ 44 m³/ano

Usando a distribuição t (bicaudal):

n <- nrow(dados_agua_sna)
n

[1] 506

media <- mean(dados_agua_sna$cons1)
media

[1] 46.62771

mi <- 44
mi

[1] 44

sd <- sd(dados_agua_sna$cons1)
sd

[1] 16.50024

tcalculado <- (media - mi)/(sd/sqrt(n))
tcalculado

[1] 3.582309

tcritico <- qt(p=1-0.05/2, df=n-1)
tcritico

[1] 1.964673

Interpretação: (t calculado > t crítico; 3.582309 > 1.964673) Existe evidência suficiente, ao nivel de significância de 5%, para rejeitar H₀, ou seja, a afirmação de que o consumo de água per capita em m³/ano (cons1) nos municípios paulistas no ano de referência dos dados é de 44 m³/ano.

DICA: Fizemos o passo a passo anteriomente para entender como se calcula, contudo, o R/RStudio possui uma função específica para o teste t, t.test().

Todos os valores que encontramos acima são apresentados na saída do teste. Adicionalmente, temos o valor-p, confirmando a rejeição de H₀, uma vez que a probabilidade de sua ocorrência (valor igual ou maior) é de 0,03735%. Note que temos inclusive a descrição da hipótese alternativa. Outro bônus desta função é o cálculo do intervalo de confiança, exatamente o que encontramos no passo a passo detalhado.

t.test(dados_agua_sna$cons1, mu=44, alternative = "two.sided")

    One Sample t-test

data:  dados_agua_sna$cons1
t = 3.5823, df = 505, p-value = 0.0003735
alternative hypothesis: true mean is not equal to 44
95 percent confidence interval:
 45.18658 48.06885
sample estimates:
mean of x 
 46.62771

Correlação Linear e Regressão Linear Simples

Estrutura dos dados analisados nesta Prática

Corte transversal: um conjunto de observações (em geral amostral) sejam indivíduos, empresas, países, regiões, municípios, dentre outros, analisado em um determinado ponto no tempo. O conjunto de dados usado nesta prática: dados das variáveis consumo de água (y) e PIB per capita (x) para os 506 municípios paulistas no ano de 2010.

Para outras estruturas de dados são necessários outros métodos de estimação, testes e premissas que não serão alvo desta Prática.

Correlação Linear

Uma correlação é uma relação entre duas variáveis. Os dados podem ser representados por pares ordenados (x, y), sendo x a variável independente (ou explanatória) e y a variável dependente (ou resposta). Veremos a relação que se estabelece linearmente.

Antes de calcular o coeficiente de correlação linear vamos plotar o diagrama de dispersão (scatter plot) usando a função plot(). O diagrama de dispersão é um gráfico que exibe, para cada observação, o valor de uma variável contra o valor em outra variável, permitindo a visualização da correlação entre duas variáveis.

Vamos investigar a relação linear entre as variáveis cons1 e PIBpc_rc:

plot(x = dados_agua_sna$PIBpc_rc, 
     y = dados_agua_sna$cons1,
     xlab = "PIB per capita R$ de 2010",
     ylab = "Consumo de água per capita m3/ano")

Figura 2 - Diagrama de dispersão

Interpretação: o diagrama de dispersão parece indicar que existe uma correlação positiva, aparentemente, muito fraca entre as duas variáveis, ou seja, quanto maior a renda, maior o consumo de água.

Um coeficiente de correlação é uma medida padronizada do relacionamento entre duas variáveis. Vamos calcular o coeficiente de correlação de Pearson, que investiga se essa correlação existe de forma linear. A interpretação do coeficiente está no intervalo de -1 e 1:

Figura 3 - Diagrama de dispersão

notação: r para o coeficiente amostral e ρ para o coeficiente populacional.

Usaremos a função cor(), com os seguintes argumentos:

x: PIB per capita (PIBpc_rc)
y: Consumo de Água per capita (cons1)
method = “pearson”: para calcular o coeficiente de correlação de Pearson
use = “complete.obs”: para remover os valores faltantes (NA), no nosso caso esse argumento não interfere, pois já removemos os NAs anteriormente.

cor(x = dados_agua_sna$PIBpc_rc, 
    y = dados_agua_sna$cons1,
    method = "pearson",
    use = "complete.obs")

[1] 0.1220042

Lembre-se de que você esta usando dados amostrais para tomar uma decisão sobre dados populacionais, então, é sempre possível que sua inferência possa estar errada.

Já sabemos a estrutura de um teste de hipóteses, portanto, seguiremos direto para o uso da função cor.test(). A estatística t calculada na função é:

\[ t_{calculado} = \frac{r}{\frac{1-r^2}{n-2} } \]

Usaremos a função cor.test(), com os seguintes argumentos:

x: PIB per capita (PIBpc_rc)
y: Consumo de Água per capita (cons1)
method = “pearson”: para calcular o coeficiente de correlação de Pearson
alternative = “two.sided”: para teste de hipótese bicaudal
conf.level = 0.95: nível de significância como 0,05

cor.test(x = dados_agua_sna$PIBpc_rc, 
         y = dados_agua_sna$cons1,
         method = "pearson",
         alternative = "two.sided",
         conf.level = 0.95)

    Pearson's product-moment correlation

data:  dados_agua_sna$PIBpc_rc and dados_agua_sna$cons1
t = 2.7596, df = 504, p-value = 0.005998
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.03520999 0.20697177
sample estimates:
      cor 
0.1220042

Interpretação: a correlação entre o PIB per capita e o consumo de água per capita nos municípios paulistas no ano de referência dos dados é de 0,12 ou 12%, confirmando a percepção do diagrama de dispersão de uma correlação positiva, porém muito fraca.

o teste t, (t calculado > t crítico; 2.7596 > 1.964673) indica que existe evidência suficiente, ao nivel de significância de 5%, para rejeitar H₀, ou seja, de que pode existir uma correlação entre o PIB per capita e o consumo de água per capita nos municípios paulistas no ano de referência dos dados.

Adicionalmente, temos o valor-p, confirmando a rejeição de H₀, uma vez que a probabilidade de sua ocorrência (valor igual ou maior) é de 0,5998%. Note que temos inclusive a descrição da hipótese alternativa. Outro bônus desta função é o cálculo do intervalo de confiança.

ALERTA: Coeficientes de correlação NÃO dão indicação da causalidade!!!

DICA: pacotes como tidyverse, corrplot e Hmisc permitem explorar matrizes de correlação dentro do conjunto de dados.

Regressão Linear Simples

ALERTA: nesta prática vamos fazer um uso bastante introdutório sobre a regressão linear simples (simples, por ser apenas uma variável X), ou seja, a análise de regressão deve ser alvo de aprofundamento para aqueles que desejam fazer aplicações robustas (quem sabe uma próxima prática!).

Após verificar que a correlação linear entre duas variáveis é significativa, o próximo passo é determinar a equação da reta que melhor modela os dados. Essa reta é chamada de reta de regressão e sua equação pode ser usada para predizer os valores de y para um dado valor de x.

\[ \hat{y} = b + mx \] em que o $\hat{y}$ é o valor previsto de y para um dado valor de x. A inclinação m e o intercepto em y, b, são dados por:

\[ m = \frac{n\sum{xy} - (\sum{x}) (\sum{y})}{n\sum{x^2}-(\sum{x})^2} \]

\[ b = \bar{y} - m\bar{x} = \frac{\sum{y}}{n} - m\frac{\sum{x}}{n} \]

Em que $\bar{y}$ é a média dos valores de y no conjunto de dados, $\bar{x}$ é a média dos valores de x e n é o número de pares de dados. A reta de regressão sempre passa pelo ponto ($\bar{x}$,$\bar{y}$).

Como encontrar a “linha” que melhor se ajusta aos nossos dados? Precisamos considerar a existência do erro.

\[ \hat{y_i} = \beta_0 + \beta_1X_1+\varepsilon_i \]

em que:

$\beta_0$ representa b

$\beta_1$ representa m

$\varepsilon_i$ é o erro do modelo, corresponde à diferença entre o verdadeiro valor de y_i e o y_i estimado, ou seja, $\varepsilon_i = y_i-\hat{y_i}$, também chamado de resíduo da regressão.

Existem alguns métodos para estimar os parámetros Betas ($\beta_0$ e $\beta_1$). Aqui usaremos as funções do R/RStudio que estimam pelo chamado:

Método dos Mínimos Quadrados Ordinários (MQO ou OLS): procedimento matemático/estatístico que minimiza (por meio de cálculo de derivadas, sistema de equações lineares) a soma dos desvios quadrados.

ATENÇÃO: não iremos aprofundar nesta prática a Análise de Variância, Análise dos Resíduos, bem como explicações mais detalhadas sobre as suposições do MQO. Se você pretende aplicar o método com mais robustez é indicado que busque estes conhecimentos.

Suposições importantes na maioria dos métodos de estimação linear (sobretudo no MQO):

Linearidade nos parâmetros (ou linearizáveis)
Amostra Aleatória
Inexistência de Multicolinearidade, ou seja, variabilidade em todas as variáveis explicativas e não existe combinação linear perfeita entre as variáveis explicativas.
Exogeneidade de x_i , ou seja, E(x|u) = 0 (sendo u o resíduo)
Homocesasticidade, ou seja, a variância de u_i dado x_i é a mesma para qualquer x_i, V(u|x) =σ².

Coeficiente de determinação R² e R²_aj

Quando elevamos o coeficiente de correlação de Pearson ao quadrado ele passa a representar a razão entre a variação explicada pelo nosso modelo e a variação total presente no modelo.

Temos, então, o coeficiente de determinação R²; que está entre 0 e 1, sendo que quanto mais próximo de 1, maior o poder explicativo do modelo estimado:

\[ R^2 = \frac{SQM}{SQT} = \frac{SQT-SQR}{SQT} = 1-\frac{SQR}{SQT} \]

pela análise da variância, tem-se que:

SQT: Soma dos Quadrados Totais (desvio total)

SQM: Soma dos Quadrados do Modelo (desvio explicado pelo modelo estimado)

SQR: Soma dos Quadrados dos Resíduos (desvio não explicado pelo modelo estimado)

Como o R² nunca diminui quando adicionamos variáveis ao modelo, podemos fazer uma alteração nele para sabermos se o termo adicional foi útil ou não para aumentar a capacidade de explicação.

Temos, então, o coeficiente de determinação R² ajustado; que está entre 0 e 1, sendo que quanto mais próximo de 1, maior o poder explicativo do modelo estimado:

\[ R^2_{aj} = \frac{SQM/(n-k)}{SQT/(n-1)} \] em que n é o tamanho da amostra e k é o número de coeficientes a serem estimados.

Teste t (bicaudal)

O Teste t (bicaudal) tem por finalidade testar a existência de efeito linear individual de X (via $\beta_1$) e do intercepto ($\beta_0$).

Já sabemos a estrutura de um teste de hipóteses, portanto, veremos na saída da regressão do R/RStudio os valores calculados. A estatística t calculada na função é:

\[ t_{calculado} = \frac{\hat{\beta_1}}{s(\hat{\beta_1}) } \] em que o erro padrão do coeficiente é $s(\hat{\beta_1})$:

\[ s(\hat{\beta_1}) = \sqrt{\frac{SQR/(n-k)} {\sum(x_i-\bar{x})^2}} \]

DICA: outro teste comum em análises de regressão é o Teste F. Ele tem por finalidade testar o efeito conjunto das variáveis explicativas sobre a dependente. Isso significa verificar se, pelo menos, uma das variáveis explicativas do modelo exerce efetivamente alguma influência linear sobre a variável dependente. Sua aplicação ocorre, principalmente, em modelos de regressão múltipla (mais de uma variável explicativa).

A relação linear simples investigada via MQO será:

x: variável independente (indicadora/explicativa/preditora); PIB per capita (PIBpc_rc)
y: variável dependente (resposta/saída); Consumo de Água per capita (cons1)

Usaremos a função lm(), com os seguintes argumentos:

fórmula: que indica a relação linear, seguindo a estrutura “y ~ x”
data: que indica a base de dados
na.action = na.exclude: que indica a exclusão de valores faltantes (NAs)

Para evitar resultados com notação científica usaremos, adicionalmente, a função options().

options(scipen=999)

modelo1 <- lm (formula = cons1 ~ PIBpc_rc, data = dados_agua_sna)
summary(modelo1)

Call:
lm(formula = cons1 ~ PIBpc_rc, data = dados_agua_sna, na.action = na.exclude)

Residuals:
    Min      1Q  Median      3Q     Max 
-46.368  -9.288  -0.895   6.823 143.811 

Coefficients:
               Estimate  Std. Error t value            Pr(>|t|)    
(Intercept) 44.22266474  1.13606748   38.93 <0.0000000000000002 ***
PIBpc_rc     0.00012081  0.00004378    2.76               0.006 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.39 on 504 degrees of freedom
Multiple R-squared:  0.01489,   Adjusted R-squared:  0.01293 
F-statistic: 7.615 on 1 and 504 DF,  p-value: 0.005998

Interpretação:

Reta de regressão estimada por nosso modelo: Y = 44,22 + 0,00012081X.

Caso os testes confirmem que temos coeficientes estatisticamente significativos, então, podemos inferir que, em média, um aumento de R$ 1,00 no PIB per capita está associado a um aumento de 0,000012081 m³/ano no consumo de água (0,000012081 * 1000 = 0,120810 litros/ano).

Testes

Coeficientes de determinação (R² e R² ajustado): o R² é 0,01489 ou 1,489%, podemos interpretar como a capacidade da variável preditora X explicar a variação em Y. O R² ajustado é 0,01293 ou 1,293%, medida alternativa ao R² que penaliza a inclusão de variáveis independentes (X) pouco explicativas. O R² e R² ajustado indicam um modelo com baixíssimo poder explicativo, mostrando ser necessário procurar outras variáveis para contribuir na explicação do consumo de água per capita. Um modelo com bom poder explicativo apresenta coeficientes de determinação acima de 0,6 ou 60%.

Estatística t e p-valor dos coeficientes betas:

Para $\beta_0$: t_calculado = 38,93 (podemos chegar no mesmo valor usando a fórmula do erro padrão do coeficiente = 44,22/1,1360),

dado um t_crítico = 3,309935 (usando a função qt(p=1-0.001/2, df=n-2)para o nível de significância de 0,001, essa informação aparece pelos asteriscos).

t_calculado> t_crítico (38,93 > 3,3099), rejeito H₀, ou seja, o intercepto é estatísticamente significativo, informação reforçada pelo Valor-p abaixo de 5% (0,0000000000000002 ou 0,0%).

Para $\beta_1$: t_calculado = 2,76 (podemos chegar no mesmo valor usando a fórmula do erro padrão do coeficiente =0,00012081/0,00004378),

dado um t_crítico = 2,585619 (usando a função qt(p=1-0.01/2, df=n-2)para o nível de significância de 0,01, essa informação aparece pelos asteriscos).

t_calculado> t_crítico (2,76 > 2,5856), rejeito H₀, ou seja, o coeficiente de inclinação é estatísticamente significativo, informação reforçada pelo Valor-p abaixo de 5% (0.006 ou 0,6%).

Prática 4 - Inferencial

Gislaine Quaglio

2021

Antes de qualquer coisa…

Importar banco de dados e visualizar