Conceitos importantes de Inferência Estatística

Definição 1: Forma Geral da Estatística t

Lembre-se que uma estatística \(t\) tem a forma

\[ t = \frac{\text{valor estimado} - \text{valor hipotetico}}{\text{erro padrao do estimador}} \]

Definição 2: Testando hipóteses sobre \(\beta_1\)

Para testar a hipótese \(H_0: \beta_1 = \beta_{1,0}\) , precisamos realizar os seguintes passos:

1. Calcule o erro padrão de \(\hat{\beta}_1\), \(\text{EP}(\hat{\beta}_1)\)

\[ \text{EP}(\hat{\beta}_1) = \sqrt{ \widehat{Var(\hat{\beta}_1})}, \\ \text{EP}(\hat{\beta}_1) = \frac{1}{n} \times \frac{\frac{1}{n-2} \sum_{i=1}^n (X_i - \overline{X})^2 \hat {u_i}^2 }{ \left[ \frac{1}{n} \sum_{i=1}^n (X_i - \overline{X})^2 \right]^2} \]

Em notação matricial, a matriz de variância e covariância é dada por

\[ \text{Var} (\hat{\beta}) = \sigma^2(X'X)^{-1} \]

onde a diagonal principal da matriz possui os estimadores de variância

2. Calcule a estatística \(t\)

\[ t = \frac{\hat{\beta}_1 - \beta_{1,0}}{ \text{EP}(\hat{\beta}_1) }. \]

3. Dada uma alternativa bilateral (\(H_1:\beta_1 \neq \beta_{1,0}\)) rejeitamos no nível de \(5\%\) se \(|t^{calc}| > 1,96\) ou, equivalentemente, se o \(p-valor\) for inferior a \(0,05\). Lembre-se da definição do \(p-valor\):

\[ \begin{align} p \text{-valor} =& \, \text{Pr}_{H_0} \left[ \left| \frac{ \hat{\beta}_1 - \beta_{1,0} }{ \text{EP}(\hat{\beta}_1) } \right| > \left| \frac{ \hat{\beta}_1^{calc} - \beta_{1,0} }{ \text{EP}(\hat{\beta}_1) } \right| \right] \\ =& \, \text{Pr}_{H_0} (|t| > |t^{calc}|) \\ \approx& \, 2 \cdot \Phi(-|t^{calc} |) \end{align} \]

A última transformação se deve à aproximação normal para amostras grandes.

Teste de hipóteses para Regressão Simples

Considere o modelo que estimamos na Monitoria 1:

\[ \widehat{TestScore} \ = \underset{(9.47)}{698.9} - \underset{(0.49)}{2.28} \times prop\_aluno\_prof \ , \ R^2=0.051 \ , \ SER=18.6. \]

Os resultados podem ser obtidos através dos comandos:

# load the AER package
library(AER)

# load the data set
data(CASchools)   

CASchools = CASchools %>% 
  mutate(prop_aluno_prof = students/teachers,
         score = (read + math)/2)

# estimate the model
linear_model <- lm(score ~ prop_aluno_prof, data = CASchools)   
summary(linear_model)

## 
## Call:
## lm(formula = score ~ prop_aluno_prof, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.727 -14.251   0.483  12.822  48.540 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     698.9329     9.4675  73.825  < 2e-16 ***
## prop_aluno_prof  -2.2798     0.4798  -4.751 2.78e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.58 on 418 degrees of freedom
## Multiple R-squared:  0.05124,    Adjusted R-squared:  0.04897 
## F-statistic: 22.58 on 1 and 418 DF,  p-value: 2.783e-06

Para testar uma hipótese relativa ao parâmetro de inclinação (o coeficiente em \(prop\_aluno\_prof\) ), precisamos de \(\text{EP}(\hat{\beta}_1)\), o erro padrao do respectivo estimador pontual. Como é comum na literatura, os erros padrao são apresentados entre parênteses abaixo das estimativas pontuais.

A definição 2 mostra que é bastante complicado calcular o erro padrao e, portanto, a estatística \(t\) manualmente. A pergunta que você deveria estar se perguntando agora é: podemos obter esses valores com o mínimo esforço usando R? Sim, nós podemos. Vamos primeiro usar summary() para obter um resumo dos coeficientes estimados em linear_model.

# print the summary of the coefficients to the console
summary(linear_model)$coefficients

##                   Estimate Std. Error   t value      Pr(>|t|)
## (Intercept)     698.932949  9.4674911 73.824516 6.569846e-242
## prop_aluno_prof  -2.279808  0.4798255 -4.751327  2.783308e-06

A segunda coluna dos coeficientes relatórios resumidos \(SE(\hat\beta_0)\) e \(SE(\hat\beta_1)\). Além disso, na terceira colunat value, encontramos estatística \(t\) e \(t^{calc}\) adequadas para testes das hipóteses separadas \(H_0: \beta_0=0\) e \(H_0: \beta_1=0\). Além disso, a saída nos fornece os \(p-valor\) correspondentes a ambos os testes em relação às alternativas bilaterais \(H_1:\beta_0\neq0\) respectivamente \(H_1:\beta_1\neq0\) na quarta coluna da tabela.

Vamos dar uma olhada mais de perto no teste de

\[H_0: \beta_1=0 \ \ \ vs. \ \ \ H_1: \beta_1 \neq 0\]

Nós temos

\[ t^{calc} = \frac{-2,279808 - 0}{0,4798255} \approx - 4,75.\]

O que isso nos diz sobre a importância do coeficiente estimado? Rejeitamos a hipótese nula ao nível de significância de \(5\%\) , uma vez que \(|t^{calc}| > 1,96\). Ou seja, a estatística de teste observada cai na região de rejeição como \(p\text{-valor} = 2.78\cdot 10^{-6} < 0.05\) . Concluímos que o coeficiente é significativamente diferente de zero. Em outras palavras, rejeitamos a hipótese de que o tamanho da turma não tem influência nas notas dos alunos nos testes no nível \(5\%\).

Observe que embora a diferença seja insignificante no presente caso, como veremos mais tarde, summary() não realiza a aproximação normal, mas calcula valores \(p\) usando a distribuição \(t\) . Geralmente, os graus de liberdade da distribuição \(t\) assumida são determinados da seguinte maneira:

\[\text{DF} = n - k - 1\]

onde \(n\) é o número de observações utilizadas para estimar o modelo e \(k\) é o número de regressores, excluindo o intercepto. No nosso caso, temos \(n=420\) observações e o único regressor é prop_aluno_prof então \(k=1\) . A maneira mais simples de determinar os graus de liberdade do modelo é

# determine residual degrees of freedom
linear_model$df.residual

## [1] 418

Portanto, para a distribuição amostral assumida de \(\hat\beta_1\) temos

\[ \hat\beta_1 \sim t_{418} \]de modo que o \(p\text{-valor}\) para um teste de significância bilateral possa ser obtido executando o seguinte código:

2 * pt(-4.751327, df = 418)

## [1] 2.78331e-06

Através do seguinte código podemos calcular

O resultado está muito próximo do valor fornecido por summary(). No entanto, como \(n\) é suficientemente grande, também se poderia usar a densidade normal padrao para calcular o \(p\text{-valor}\):

2 * pnorm(-4.751327)

## [1] 2.02086e-06

A diferença é realmente insignificante. Essas descobertas nos dizem que, se \(H_0: \beta_1 = 0\) for verdadeiro e repetirmos todo o processo de coleta de observações e estimativa do modelo, observando um \(\hat\beta_1 \geq |-2.28|\) é muito improvável!

Usando R podemos visualizar como tal afirmação é feita ao usar a aproximação normal. Não deixe que o seguinte trecho de código o detenha: o código é um pouco mais longo que os exemplos usuais e parece desagradável, mas há muita repetição, pois tons de cores e anotações são adicionados em ambas as extremidades da distribuição normal. Recomendamos executar o código passo a passo para ver como o gráfico é aumentado com as anotações.

# Plot the standard normal on the support [-6,6]
t <- seq(-6, 6, 0.01)

plot(x = t, 
     y = dnorm(t, 0, 1), 
     type = "l", 
     col = "steelblue", 
     lwd = 2, 
     yaxs = "i", 
     axes = F, 
     ylab = "", 
     main = expression("Calculating the p-value of a Two-sided Test when"
             ~ t^act ~ "=-4.75"), 
     cex.lab = 0.7,
     cex.main = 1)

tact <- -4.75

axis(1, at = c(0, -1.96, 1.96, -tact, tact), cex.axis = 0.7)

# Shade the critical regions using polygon():

# critical region in left tail
polygon(x = c(-6, seq(-6, -1.96, 0.01), -1.96),
        y = c(0, dnorm(seq(-6, -1.96, 0.01)), 0), 
        col = 'orange')

# critical region in right tail

polygon(x = c(1.96, seq(1.96, 6, 0.01), 6),
        y = c(0, dnorm(seq(1.96, 6, 0.01)), 0), 
        col = 'orange')

# Add arrows and texts indicating critical regions and the p-value
arrows(-3.5, 0.2, -2.5, 0.02, length = 0.1)
arrows(3.5, 0.2, 2.5, 0.02, length = 0.1)

arrows(-5, 0.16, -4.75, 0, length = 0.1)
arrows(5, 0.16, 4.75, 0, length = 0.1)

text(-3.5, 0.22, 
     labels = expression("0.025"~"="~over(alpha, 2)),
     cex = 0.7)
text(3.5, 0.22, 
     labels = expression("0.025"~"="~over(alpha, 2)),
     cex = 0.7)

text(-5, 0.18, 
     labels = expression(paste("-|",t[act],"|")), 
     cex = 0.7)
text(5, 0.18, 
     labels = expression(paste("|",t[act],"|")), 
     cex = 0.7)

# Add ticks indicating critical values at the 0.05-level, t^act and -t^act 
rug(c(-1.96, 1.96), ticksize  = 0.145, lwd = 2, col = "darkred")
rug(c(-tact, tact), ticksize  = -0.0451, lwd = 2, col = "darkgreen")

Imagine que você pudesse extrair todas as amostras aleatórias possíveis de um determinadotamanho. O
intervalo que contém o valor verdadeiro $\beta_i$ em $95\%$ de todas as amostras é dado pela expressão

$$\text{CI}_{0,95}^{\beta_i} = \left[ \hat{\beta}_i - 1,96 \times SE(\hat{\beta}_i) \, , \, \hat{\beta} _i + 1,96 \times SE(\hat{\beta}_i) \right].$$

# set seed for reproducibility
set.seed(4)

# generate and plot the sample data
Y <- rnorm(n = 100, 
           mean = 5, 
           sd = 5)

plot(Y, 
     pch = 19, 
     col = "steelblue")

cbind(CIlower = mean(Y) - 1.96 * 5 / 10, CIupper = mean(Y) + 1.96 * 5 / 10)

##       CIlower  CIupper
## [1,] 4.502625 6.462625

# set seed
set.seed(1)

# initialize vectors of lower and upper interval boundaries
lower <- numeric(10000)
upper <- numeric(10000)

# loop sampling / estimation / CI
for(i in 1:10000) {
  
  Y <- rnorm(100, mean = 5, sd = 5)
  lower[i] <- mean(Y) - 1.96 * 5 / 10
  upper[i] <- mean(Y) + 1.96 * 5 / 10
  
}

# join vectors of interval bounds in a matrix
CIs <- cbind(lower, upper)

mean(CIs[, 1] <= 5 & 5 <= CIs[, 2])

## [1] 0.9487

#> [1] 0.9487

# identify intervals not covering mu
# (4 intervals out of 100)
ID <- which(!(CIs[1:100, 1] <= 5 & 5 <= CIs[1:100, 2]))

# initialize the plot
plot(0, 
     xlim = c(3, 7), 
     ylim = c(1, 100), 
     ylab = "Sample", 
     xlab = expression(mu), 
     main = "Confidence Intervals")

# set up color vector
colors <- rep(gray(0.6), 100)
colors[ID] <- "red"

# draw reference line at mu=5
abline(v = 5, lty = 2)

# add horizontal bars representing the CIs
for(j in 1:100) {
  
  lines(c(CIs[j, 1], CIs[j, 2]), 
        c(j, j), 
        col = colors[j], 
        lwd = 2)
  
}

# compute 95% confidence interval for coefficients in 'linear_model'
confint(linear_model)

##                     2.5 %     97.5 %
## (Intercept)     680.32312 717.542775
## prop_aluno_prof  -3.22298  -1.336636

# compute 95% confidence interval for coefficients in 'linear_model' by hand
lm_sum <- summary(linear_model)

c("lower" = lm_sum$coef[2,1] - qt(0.975, df = lm_sum$df[2]) * lm_sum$coef[2, 2],
  "upper" = lm_sum$coef[2,1] + qt(0.975, df = lm_sum$df[2]) * lm_sum$coef[2, 2])

##     lower     upper 
## -3.222980 -1.336636

Teste de hipóteses para Regressão Múltipla

Discutiremos primeiro como calcular erros padrao, testar hipóteses e construir intervalos de confiança para um único coeficiente de regressão \(\beta_j\) em um modelo de regressão múltipla.

Definição 3: Teste de hipótese \(\beta_j = \beta_{j,0}\) contra a alternativa \(\beta_j \neq \beta{j,0}\)

1. Calcule o erro padrao de \(\hat{\beta_j}\)

2. Calcule a estatística \(t , t^{calc} = \frac{\hat{\beta}_j - \beta_{j,0}} {SE(\hat{\beta_j})}\)

3. Calcule o \(p\text{-valor}\), \(p\text{-valor} = 2 \Phi(-|t^{calc}|)\)

onde \(t^{calc}\) é o valor da estatística \(t\) realmente calculada. Rejeite a hipótese no nível de significância \(5\%\) se o $p\text{-valor}$ for menor que $0,05$ ou, equivalentemente, se \(|t^{calc}| > 1,96\). O erro padrao e (normalmente) a estatística $t$ e o $p$ correspondente para teste $_j = 0$ são calculados automaticamente por funções R adequadas, por exemplo, por `summary`.

Testar uma única hipótese sobre a significância de um coeficiente no modelo de regressão múltipla procede de forma semelhante ao processo no modelo de regressão simples.

Você pode ver isso facilmente inspecionando o resumo dos coeficientes do modelo de regressão

\[TestScore = \beta_0 + \beta_1 \times prop\_aluno\_prof +\beta_2 \times ingles + u\]

model <- lm(score ~ prop_aluno_prof + english, data = CASchools)
coeftest(model, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)     686.032245   8.728225  78.5993  < 2e-16 ***
## prop_aluno_prof  -1.101296   0.432847  -2.5443  0.01131 *  
## english          -0.649777   0.031032 -20.9391  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Você pode verificar se essas quantidades são calculadas como no modelo de regressão simples calculando manualmente as estatísticas \(t\) ou os \(p\text{-valor}\) usando a saída fornecida acima e usando R como calculadora.

Por exemplo, usando a definição do \(p\text{-valor}\) para um teste bilateral conforme dado na Definição 4, podemos confirmar o \(p\text{-valor}\) para um teste da hipótese de que o coeficiente \(\beta_1\), o coeficiente de size, seja aproximadamente zero.

# compute two-sided p-value
2 * (1 - pt(abs(coeftest(model, vcov. = vcovHC, type = "HC1")[2, 3]),
            df = model$df.residual))

## [1] 0.01130921

O cálculo dos intervalos de confiança para coeficientes individuais no modelo de regressão múltipla prossegue como no modelo de regressão simples usando a função confint().

model <- lm(score ~ prop_aluno_prof + english, data = CASchools)
confint(model)

##                       2.5 %      97.5 %
## (Intercept)     671.4640580 700.6004311
## prop_aluno_prof  -1.8487969  -0.3537944
## english          -0.7271113  -0.5724424

Para obter intervalos de confiança em outro nível, digamos \(90\%\) , basta definir o nível do argumento em nossa chamada de confint() de acordo.

confint(model, level = 0.9)

##                         5 %        95 %
## (Intercept)     673.8145793 698.2499098
## prop_aluno_prof  -1.7281904  -0.4744009
## english          -0.7146336  -0.5849200

# scale expenditure to thousands of dollars
CASchools$expenditure <- CASchools$expenditure/1000

# estimate the model
model <- lm(score ~ prop_aluno_prof + english + expenditure, data = CASchools)
coeftest(model)

## 
## t test of coefficients:
## 
##                   Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)     649.577947  15.205717  42.7193 < 2.2e-16 ***
## prop_aluno_prof  -0.286399   0.480523  -0.5960  0.551489    
## english          -0.656023   0.039106 -16.7756 < 2.2e-16 ***
## expenditure       3.867901   1.412122   2.7391  0.006426 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# compute the sample correlation between 'prop_aluno_prof' and 'expenditure'
cor(CASchools$prop_aluno_prof, CASchools$expenditure)

## [1] -0.6199822

linearHypothesis(model, "expenditure=4")

## Linear hypothesis test
## 
## Hypothesis:
## expenditure = 4
## 
## Model 1: restricted model
## Model 2: score ~ prop_aluno_prof + english + expenditure
## 
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1    417 85701                           
## 2    416 85700  1    1.8028 0.0088 0.9255

Teste de hipóteses conjuntas usando a estatística F

O modelo estimado é

\[ \widehat{TestScore} = \underset{(15,21)}{649,58} -\underset{(0,48)}{0,29} \times prop_aluno_prof - \underset{(0,04)}{0,66} \times ingles + \underset{ (1,41)}{3,87} \times despesas. \]

Agora, podemos rejeitar a hipótese de que o coeficiente de tamanho e o coeficiente de despesa são zero? Para responder a isto, temos que recorrer a testes de hipóteses conjuntas. Uma hipótese conjunta impõe restrições aos coeficientes de regressão múltiplos. Isso é diferente de conduzir testes \(t\) individuais, onde uma restrição é imposta a um único coeficiente.

A Estatística-\(F\) é dada por

\[ F = \frac{(SSR_{\text{restrito}} - SSR_{\text{irrestrito}})/q}{SSR_{\text{irrestrito}} / (nk-1)} \]

sendo \(SSR_{restrito}\) a soma dos quadrados dos resíduos da regressão restrita, ou seja, a regressão onde impomos a restrição. \(SSR_{irrestrito}\) é a soma dos resíduos quadrados do modelo completo, q é o número de restrições sob o nulo e \(k\) é o número de regressores na regressão irrestrita.

É fácil realizar testes \(F\) no R. Podemos usar a função linearHypothesis() contida no pacote car.

# estimate the multiple regression model
model <- lm(score ~ prop_aluno_prof + english + expenditure, data = CASchools)

# execute the function on the model object and provide both linear reprop_aluno_profictions 
# to be tested as prop_aluno_profings
linearHypothesis(model, c("prop_aluno_prof=0", "expenditure=0"))

## Linear hypothesis test
## 
## Hypothesis:
## prop_aluno_prof = 0
## expenditure = 0
## 
## Model 1: restricted model
## Model 2: score ~ prop_aluno_prof + english + expenditure
## 
##   Res.Df   RSS Df Sum of Sq      F   Pr(>F)    
## 1    418 89000                                 
## 2    416 85700  2    3300.3 8.0101 0.000386 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A saída revela que a estatística F para este teste de hipótese conjunta é cerca de \(8.01\) e o valor p correspondente é \(0.0004\). Assim, podemos rejeitar a hipótese nula de que ambos os coeficientes são zero em qualquer nível de significância comumente utilizado na prática.

A saída padrao de um resumo do modelo também relata uma estatística \(F\) e o \(p\)-valor correspondente. A hipótese nula pertencente a este teste \(F\) é que todos os coeficientes populacionais no modelo, exceto o intercepto, são zero, então as hipóteses são

\[H_0: \beta_1=0, \ \beta_2 =0, \ \beta_3 =0 \quad \text{vs.} \quad H_1: \beta_j \neq 0 \ \text{para pelo menos um} \ j=1,2,3.\]

Isso também é chamado de estatística de regressão geral \(F\) e a hipótese nula é obviamente diferente do teste se apenas \(\beta_1\) e \(\beta_3\) forem zero.

Verificamos agora se a estatística \(F\) pertencente ao \(p\) listado no resumo do modelo coincide com o resultado relatado por linearHypothesis().

# execute the function on the model object and provide the reprop_aluno_profictions 
# to be tested as a character vector
linearHypothesis(model, c("prop_aluno_prof=0", "english=0", "expenditure=0"))

## Linear hypothesis test
## 
## Hypothesis:
## prop_aluno_prof = 0
## english = 0
## expenditure = 0
## 
## Model 1: restricted model
## Model 2: score ~ prop_aluno_prof + english + expenditure
## 
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    419 152110                                  
## 2    416  85700  3     66410 107.45 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model)$fstatistic

##    value    numdf    dendf 
## 107.4547   3.0000 416.0000

O valor de entrada é a estatística geral \(F\) e é igual ao resultado de linearHypothesis(). O teste \(F\) rejeita a hipótese nula de que o modelo não tem poder para explicar os resultados dos testes.

# draw the 95% confidence set for coefficients on prop_aluno_prof and expenditure
confidenceEllipse(model, 
                  fill = T,
                  lwd = 0,
                  which.coef = c("prop_aluno_prof", "expenditure"),
                  main = "95% Confidence Set",
                  ylab="Coefficients of Expenditure",
                  xlab="Coefficients of prop_aluno_prof")

Especificação do Modelo para Regressão Múltipla

Definição 4: Viés de Variável Omitida

O viés de variável omitida é o viés no estimador OLS que surge quando os regressores se correlacionam com uma variável omitida. Para que surja viés de variável omitida, duas coisas devem ser verdadeiras:

Pelo menos um dos regressores incluídos deve estar correlacionado com a variável omitida.

A variável omitida deve ser um determinante da variável dependente $S$.

Discutiremos agora um exemplo onde enfrentamos um potencial viés de variável omitida em um modelo de regressão múltipla:

Considere novamente a equação de regressão estimada

\[ \widehat{TestScore} = \underset{(8.7)}{686.0} - \underset{(0.43)}{1.10} \times prop\_aluno\_prof - \underset{(0.031)}{0.650} \times ingles. \]

Estamos interessados em estimar o efeito causal do tamanho da turma na pontuação do teste. Pode haver um viés devido à omissão de “oportunidades de aprendizagem externas” da nossa regressão, uma vez que tal medida poderia ser um determinante dos resultados dos testes dos alunos e também poderia ser correlacionada com ambos os regressores já incluídos no modelo.

“Oportunidades de aprendizagem externas” são um conceito complicado e difícil de quantificar. Em vez disso, um substituto que podemos considerar é o contexto económico dos alunos, que provavelmente está fortemente relacionado com oportunidades de aprendizagem externas: pense em pais ricos que são capazes de fornecer tempo e/ou dinheiro para aulas particulares dos seus filhos. Assim, aumentamos o modelo com a variável merenda , a percentagem de alunos que se qualificam para merenda escolar gratuita ou subsidiada devido aos rendimentos familiares abaixo de um determinado limite, e reestimamos o modelo.

# estimate the model and print the summary to console
model <- lm(score ~ prop_aluno_prof + english + lunch, data = CASchools)
coeftest(model, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##                   Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)     700.149957   5.568453 125.7351 < 2.2e-16 ***
## prop_aluno_prof  -0.998309   0.270080  -3.6963 0.0002480 ***
## english          -0.121573   0.032832  -3.7029 0.0002418 ***
## lunch            -0.547345   0.024107 -22.7046 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Assim, a linha de regressão estimada é

\[ \begin{align*} \widehat{TestScore} = & \underset{(5,56)}{700,15} - \underset{(0,27)}{1,00} \times prop\_aluno\_prof - \underset{(0,03)}{0,12} \times ingles \\ &- \underset{(0,02 )}{0,55} \times almoço \end{align*} \]

Não observamos mudanças substanciais na conclusão sobre o efeito do prop_aluno_prof no \(TestScore\): o coeficiente no \(prop\_aluno\_prof\) muda para apenas \(0.1\) e mantém sua significância.

Embora a diferença nos coeficientes estimados não seja grande neste caso, é útil manter o almoço para tornar mais credível a suposição de independência média condicional.

# set seed for reproducibility
set.seed(1)

# generate observations for parking lot space
CASchools$PLS <- c(22 * CASchools$income 
                   - 15 * CASchools$prop_aluno_prof 
                   + 0.2 * CASchools$expenditure
                   + rnorm(nrow(CASchools), sd = 80) + 3000)

# regress test score on PLS
summary(lm(score ~ PLS, data = CASchools))

## 
## Call:
## lm(formula = score ~ PLS, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.372  -9.742   0.592  10.481  36.867 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.575e+02  1.171e+01   39.07   <2e-16 ***
## PLS         6.453e-02  3.836e-03   16.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.73 on 418 degrees of freedom
## Multiple R-squared:  0.4037, Adjusted R-squared:  0.4022 
## F-statistic:   283 on 1 and 418 DF,  p-value: < 2.2e-16

# estimate the correlation between 'calworks' and 'lunch'
cor(CASchools$calworks, CASchools$lunch)

## [1] 0.7394218

# set up arrangement of plots
m <- rbind(c(1, 2), c(3, 0))
graphics::layout(mat = m)

# scatterplots
plot(score ~ english, 
     data = CASchools, 
     col = "steelblue", 
     pch = 20, 
     xlim = c(0, 100),
     cex.main = 0.7,
     xlab="English",
     ylab="Score",
     main = "Percentage of English language learners")

plot(score ~ lunch, 
     data = CASchools, 
     col = "steelblue", 
     pch = 20,
     cex.main = 0.7,
     xlab="Lunch",
     ylab="Score",
     main = "Percentage qualifying for reduced price lunch")

plot(score ~ calworks, 
     data = CASchools, 
     col = "steelblue", 
     pch = 20, 
     xlim = c(0, 100),
     cex.main = 0.7,
     xlab="CalWorks",
     ylab="Score",
     main = "Percentage qualifying for income assistance")

# estimate correlation between student characteristics and test scores
cor(CASchools$score, CASchools$english)

## [1] -0.6441238

cor(CASchools$score, CASchools$lunch)

## [1] -0.868772

cor(CASchools$score, CASchools$calworks)

## [1] -0.6268533

# load the stargazer library
library(stargazer)

# estimate different model specifications
spec1 <- lm(score ~ prop_aluno_prof, data = CASchools)
spec2 <- lm(score ~ prop_aluno_prof + english, data = CASchools)
spec3 <- lm(score ~ prop_aluno_prof + english + lunch, data = CASchools)
spec4 <- lm(score ~ prop_aluno_prof + english + calworks, data = CASchools)
spec5 <- lm(score ~ prop_aluno_prof + english + lunch + calworks, data = CASchools)

# generate a LaTeX table using stargazer
stargazer(spec1, spec2, spec3, spec4, spec5,
          digits = 3,
          header = F,
          column.labels = c("(I)", "(II)", "(III)", "(IV)", "(V)"),
          type = "html")

Demonstração do Teorema FWL

O Teorema de Frisch-Waugh-Lovell (FWL; após a prova inicial de Frisch e Waugh (1933), e posterior generalização por Lovell (1963)) afirma que:

O coeficiente de regressão de qualquer preditor em um modelo multivariado é equivalente ao coeficiente de regressão estimado a partir de um modelo bivariado no qual o resultado residualizado é regredido no componente residualizado do preditor; onde os resíduos são retirados de modelos que regridem o resultado e o preditor em todos os outros preditores na regressão multivariada (separadamente).

Mais formalmente, suponha que temos um modelo de regressão multivariada com \(k\) preditores:

\[ \begin{equation} \hat{y} = \hat{\beta}_{1}x_{1} + ... \hat{\beta}_{k}x_{k} + \epsilon. \end{equation} \]

FWL afirma que cada \(\hat{\beta}_{j}\) na equação é igual a \(\hat{\beta}^{*}_{j}\), o resíduo \(\epsilon = \epsilon^{*}\) em:

\[\begin{equation} \epsilon^{y} = \hat{\beta}^{*}_{j}\epsilon^{x_j} + \epsilon^{*} \end{equation}\]

onde:

\[\begin{aligned} \epsilon^y & = y - \sum_{k \neq j}\hat{\beta}^{y}_{k}x_{k} \\ \epsilon^{ x_{j}} &= x_j - \sum_{k \neq j}\hat{\beta}^{x_{j}}_{k}x_{k}. \end{aligned}\]

onde \(\hat{\beta}^{y}_k\) e \(\hat{\beta}^{x_j}_k\) são os coeficientes de regressão de dois modelos de regressão separados do resultado (omitindo \(x_j\)) e \(x_j\) respectivamente.

Em outras palavras, o FWL afirma que o coeficiente de cada preditor em uma regressão multivariada explica que a variância de \(y\) não é explicado pelo relacionamento dos outros preditores \(k-1\) com o resultado e nem com seu relacionamento com esse preditor, ou seja, o efeito independente de \(x_j\).

## Partial regressions

# Residual of y regressed on x1
y_res <- lm(score ~ prop_aluno_prof, CASchools)$residuals

# Residual of x2 regressed on x1
x_res <- lm(english ~ prop_aluno_prof, CASchools)$residuals

resids <- data.frame(CASchools$score, y_res, x_res)

## Compare the beta values for x2
# Multivariate regression:
summary(lm(score~prop_aluno_prof+english, CASchools))

## 
## Call:
## lm(formula = score ~ prop_aluno_prof + english, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     686.03224    7.41131  92.566  < 2e-16 ***
## prop_aluno_prof  -1.10130    0.38028  -2.896  0.00398 ** 
## english          -0.64978    0.03934 -16.516  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
## F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

# Partials regression
summary(lm(y_res ~ x_res, resids))

## 
## Call:
## lm(formula = y_res ~ x_res, data = resids)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.133e-14  7.049e-01    0.00        1    
## x_res       -6.498e-01  3.930e-02  -16.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.45 on 418 degrees of freedom
## Multiple R-squared:  0.3955, Adjusted R-squared:  0.394 
## F-statistic: 273.4 on 1 and 418 DF,  p-value: < 2.2e-16

summary(lm(y_res ~ x_res, resids))

## 
## Call:
## lm(formula = y_res ~ x_res, data = resids)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.133e-14  7.049e-01    0.00        1    
## x_res       -6.498e-01  3.930e-02  -16.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.45 on 418 degrees of freedom
## Multiple R-squared:  0.3955, Adjusted R-squared:  0.394 
## F-statistic: 273.4 on 1 and 418 DF,  p-value: < 2.2e-16