Capítulo 2 - O Modelo de Regressão Linear Múltiplo

Autor

Jessica Kubrusly

Considere a seguinte situação: temos uma variável de interesse, \(Y\), cujo valor é difícil ou impossível de ser observado diretamente, mas acreditamos que essa variável esteja relacionada com outras variáveis para as quais temos fácil acesso: \(X_1\), \(X_2\), \(\ldots\), \(X_k\). Queremos prever o valor de \(Y\) ou, mais precisamente, o valor esperado para \(Y\), a partir da observação dos valores de \(X_1\), \(X_2\), \(\ldots\), \(X_k\). Este é o cenário clássico de um problema de regressão.

2.1 A Formulação do Modelo

O Modelo de Regressão Linear Múltiplo é o modelo que assume que existe uma relação linear entre o valor esperado de \(Y\) e os valores das variáveis observadas \(X_1\), \(X_2\), \(\ldots\), \(X_{k}\). Isto é:

\[E(Y|\mathbf{X}=\mathbf{x}) = \beta_0 + \beta_1 x_1 + \ldots + \beta_k x_k ,\] em que \(Y\) representa a variável de interesse, \(\mathbf{X}^T = (X_1,X_2,\ldots,X_k)\) o vetor de variáveis explicativas, e observadas, e \(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_k)\) o vetor de parâmetros do modelo.

Na prática, considere \(n\) observações das variáveis \(Y\), \(X_1\), \(X_2\), \(\ldots\), \(X_k\). O valor esperado da variável resposta correspondente à \(i\)-ésima observação pode ser expresso como como uma transformação linear dos valores das \(i\)-ésimas observações das variáveis \(X_1\), \(\ldots\), \(X_k\).

\[E(Y|\mathbf{X}=\mathbf{x}_i) = E(Y_i) = \beta_0 + \beta_1 x_{i,1} + \ldots + \beta_k x_{i,k}\]

Considerando uma amostra de tamanho \(n\) para as variáveis \(Y\), \(X_1\), \(\ldots\) \(X_k\), a equação acima pode ser escrita como \[ E(\mathbf{Y}) = X \boldsymbol{{\beta}}, \] sendo \[ \mathbf{Y} = \begin{pmatrix} {Y_1} \\ {Y_2} \\ {Y_3} \\ \vdots \\ {Y_n} \end{pmatrix}_{n \times 1} \quad , \quad X = \begin{pmatrix} 1 & x_{1,1} & x_{1,2} & \ldots & x_{1,k}\\ 1 & x_{2,1} & x_{2,2} & \ldots & x_{2,k}\\ 1 & x_{3,1} & x_{3,2} & \ldots & x_{3,k}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n,1} & x_{n,2} & \ldots & x_{n,k}\\ \end{pmatrix}_{n \times (k+1)} \text{e} \quad \boldsymbol{{\beta}} = \begin{pmatrix} {\beta_0} \\ {\beta_1} \\ \vdots \\ {\beta_k} \end{pmatrix}_{k \times 1}, \]

Importante. O modelo de regressão linear não tem como objetivo estimar exatamente o valor observado \(Y_i\), mas sim o valor esperado de \(Y_i\) condicionado às variáveis explicativas. As diferenças entre os valores observados e os valores previstos serão modeladas por meio dos erros (ou resíduos), que serão introduzidos na próxima seção.

A Figura 1 mostra um exemplo simples, no qual \(Y\) será explicado por uma única variável independente \(X_1\). Os pontos estão dispersos em torno de uma reta, e é possível ver que quanto maior o valor de \(X_1\) maior é o valor de \(Y\). Este é um exemplo onde a hipótese de linearidade entre \(E(Y)\) e \(X_1\) é válida.

Figura 1: Gráfico de dispersão para um problema de regressão com uma única variável independente.

O objetivo consiste em estimar os parâmetros do modelo, isto é, obter estimativas para os coeficientes \(\beta_0\), \(\beta_1\), \(\ldots\), \(\beta_k\), permitindo estimar o valor de \(Y\) a partir das observações de \(\mathbf{X}_i\).

A estimação dos parâmetros será feita a partir de observações tanto para a variável resposta \(Y\) como para as variáveis explicativas \(X_j\). A estimação baseia-se na seguinte ideia: a partir de uma base de dados com \(n\) observações de \(Y\) e de \(X_j\), \(j = 1 \ldots, k\), e, a partir dessas observações, desejamos estimar os parâmetros assumindo a hipótese de linearidade.

\[ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{i,1} + \ldots + \hat{\beta}_k x_{i,k} \]

Considerando uma amostra de tamanho \(n\) das variáveis \(\mathbf{X}\) e considerando valores estimados para o vetor de parâmetros \(\boldsymbol{\beta}\) (depois vamos explorar como essa estimativa pode ser feita), a equação acima pode ser escrita como \[ \mathbf{\hat{y}} = X \boldsymbol{\hat\beta}, \] sendo \[ \mathbf{\hat{y}} = \begin{pmatrix} \hat{y_1} \\ \hat{y_2} \\ \hat{y_3} \\ \vdots \\ \hat{y_n} \end{pmatrix}_{n \times 1} \quad , \quad X = \begin{pmatrix} 1 & x_{1,1} & x_{1,2} & \ldots & x_{1,k}\\ 1 & x_{2,1} & x_{2,2} & \ldots & x_{2,k}\\ 1 & x_{3,1} & x_{3,2} & \ldots & x_{3,k}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n,1} & x_{n,2} & \ldots & x_{n,k}\\ \end{pmatrix}_{n \times (k+1)} \text{e} \quad \boldsymbol{\hat{\beta}} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ \vdots \\ \hat{\beta_k} \end{pmatrix}_{(k+1) \times 1}, \]

A equação acima fornece uma estimativa do valor esperado da variável resposta para a observação \(i\), \(\hat{y}_i\), em termos dos valores observados para as variáveis \(X_j\), \(j=1, \ldots, k\). Observe que \(\hat{y}_i\) representa uma estimativa do valor esperado de \(Y_i\) e não necessariamente o valor observado da variável resposta.

Supondo conhecidos ou estimados os valores de \(\beta\), é possível calcular valores para \(\hat{y}_i\) e a Figura 2 apresenta novamente o caso simplificado de uma única variável independente com o acréscimo da reta de regressão, que definem os valores previstos para \(Y\) para cada valor observado de \(X_1\).

Figura 2: Gráfico de dispersão para um problema de regressão com uma única variável independente.

2.2 Medidas de Erro

Atribuídos valores ao vetor de parâmetros \(\boldsymbol{\beta}\), digamos \(\boldsymbol{\hat\beta}\), o erro na estimativa da observação \(i\) é dados pela diferença entre o valor observado e o valor estimado pelo modelo:

\[ \varepsilon_i = y_i - \hat{y_i} = y_i - ({\hat\beta}_0 + {\hat\beta}_1 x_{i,1} + \ldots + {\hat\beta}_k x_{i,k}) = y_i - \mathbf{l}_i \boldsymbol{{\hat\beta}} \] sendo \(\mathbf{l}_i\) a \(i\)-ésima linha da matriz \(X\). Considerando uma amostra de tamanho \(n\), a equação acima pode ser escrita como \[ \boldsymbol{{\varepsilon}} = \mathbf{y} - \mathbf{\hat{y}} = \mathbf{y} - X\boldsymbol{{\hat\beta}} \] sendo \[ \boldsymbol{{\varepsilon}} = \begin{pmatrix} {\varepsilon_1} \\ {\varepsilon_2}\\ {\varepsilon_3} \\ \vdots \\ {\varepsilon_n} \end{pmatrix}_{n\times 1}. \] Veja que o erro então é função dos valores de \(\beta\). Se mudarmos o valores de \(\beta\), mudamos os valores dos erros.

A Figura 3 acrescenta ao gráfico do exemplo de uma regressão linear simples alguns erros calculados a partir da diferença entre os valores observados e estimados pelo modelo.

Figura 3: Gráfico de dispersão para um problema de regressão com uma única variável independente.

2.3 Função de Custo

A função de custo \(J: \mathbb{R}^{k+1} \to \mathbb{R}\) é uma função que associa cada vetor numérico \(\boldsymbol{\beta} \in \mathbb{R}^{k+1}\) a um custo, isto é, uma medida de erro para as estimativas realizadas a partir do valor atribuído a \(\boldsymbol{\beta}\). Uma possibilidade muito usual nos problemas de regressão, isto é, quando a variável resposta é numérica, é definir \(J\) como a soma dos erros ao quadrado.

\[ \begin{array}{cccl} J: & \mathbb{R}^{k+1} & \to & \mathbb{R}\\ & (\beta_0, \beta_1, \ldots, \beta_{k} ) & \mapsto & \sum_{i=1}^n \left({y_i} - (\beta_0 + x_{i,1}\beta_1 + \ldots + x_{i,k} \beta_{k}) \right)^2 \end{array} \]

ou, na forma matricial,

\[ \begin{array}{cccl} J: & \mathbb{R}^{k+1} & \to & \mathbb{R}\\ & \boldsymbol{\beta} & \mapsto & (\mathbf{y} - X \boldsymbol{\beta})^T(\mathbf{y} - X \boldsymbol{\beta}) \end{array} \]

\[\begin{array}{lcl} J(\boldsymbol{\beta}) & =& (\mathbf{y} - X\boldsymbol{\beta})^T(\mathbf{y} - X\boldsymbol{\beta})\\ & =& (\mathbf{y}^T - \boldsymbol{\beta}^TX^T)(\mathbf{y} - X\boldsymbol{\beta})\\ & =& \mathbf{y}^T\mathbf{y} - \mathbf{y}^TX\boldsymbol{\beta} - \boldsymbol{\beta}^TX^T\mathbf{y} + \boldsymbol{\beta}^TX^TX\boldsymbol{\beta}\\ & =& \mathbf{y}^T\mathbf{y} - \boldsymbol{\beta}^TX^T\mathbf{y} - \boldsymbol{\beta}^TX^T\mathbf{y} + \boldsymbol{\beta}^TX^TX\boldsymbol{\beta}\\ & =& \mathbf{y}^T\mathbf{y} - 2 \boldsymbol{\beta}^TX^T\mathbf{y} + \boldsymbol{\beta}^TX^TX\boldsymbol{\beta}\\ \end{array}\]

2.4 Estimador para \(\boldsymbol{\beta}\) por Mínimos Quadrados

O estimador para \(\boldsymbol{\beta}\) por mínimos quadrados é aquele que minimiza a soma dos quadrados dos erros. Ou seja, aquele que minimiza a função de custo \(J\).

\[\boldsymbol{\hat\beta} = \arg\min J(\boldsymbol{\beta})\]

Para encontrar o ponto de mínimo da função \(J\) podemos usar métodos numéricos, como o gradiente descendente. Mas, diferentemente da maioria dos modelos de aprendizado de máquina, a regressão linear admite uma solução analítica para o problema de minimização.

Para buscar o ponto de mínimo da função \(J\) vamos buscar o ponto onde o seu gradiente é nulo e depois verificar que este de fato é um ponto de mínimo. Para isso precisamos derivar \(J\) em relação a todas as suas variáveis, no caso, em relação a \(\boldsymbol{\beta}\).

\[ \nabla J (\boldsymbol{\beta}) = \dfrac{\partial{J}}{\partial \boldsymbol{\beta}} (\boldsymbol{\beta}) = \dfrac{\partial}{\partial \boldsymbol{\beta}} \left( \mathbf{y}^T\mathbf{y} - 2 \boldsymbol{\beta}^TX^T\mathbf{y} + \boldsymbol{\beta}^TX^TX\boldsymbol{\beta} \right) = - 2 X^T\mathbf{y} + 2 X^TX\boldsymbol{\beta}\] O ponto que anula o gradiente de \(J\) é aquele que satisfaz

\[- 2 X^T\mathbf{y} + 2 X^TX\boldsymbol{\beta} = \mathbf{0} \in \mathbb{R}^{k+1}\]

Veja que a equação acima descreve um sistema linear de \(k\) equações e \(k\) variáveis. A sua solução depende da existência da inversa da matriz \(X^TX\) e é dada por:

\[2 X^TX\boldsymbol{\beta} = 2 X^T\mathbf{y} \Rightarrow \hat{\boldsymbol{\beta}} = (X^TX)^{-1}X^T\mathbf{y}\]

A verificação de que esse ponto é um ponto de mínimo, e não de máximo ou de sela, pode ser feita verificando-se que a matriz Hessiana da função \(J\), avaliada em \(\hat{\boldsymbol{\beta}}\), é positiva definida. No entanto, quando a função \(J\) é dada pela soma dos quadrados dos erros, essa conclusão pode ser obtida de forma mais simples, pois \(J\) é uma função convexa. Como toda função convexa não possui máximos locais nem pontos de sela no interior de seu domínio, qualquer ponto crítico de \(J\) é necessariamente um ponto de mínimo. Com isso,

\[ \boldsymbol{\hat\beta} = (X^TX)^{-1}X^T\mathbf{y} \]

E os valores previstos a partir de \(\boldsymbol{\hat\beta}\) são

\[ \mathbf{\hat{Y}} = X \boldsymbol{\hat\beta} = X (X^TX)^{-1}X^T\mathbf{Y} = H \mathbf{Y} \]

2.5 O Modelo Estatístico

Na prática tanto os estatísticos quando os cientistas de dados trabalham com modelos de regressão linear. A diferença entre a abordagem deles e que, o modelo estatístico faz mais suposições. Além da suposição de que o valor médio de \(Y\) tem uma relação linear com as variáveis explicativas \(\mathbf{X}\), feita anteriormente, o modelo estatístico supõe também que:

  • \(Y_i = Y|\mathbf{X}=\mathbf{x}_i\) segue uma distribuição normal;
  • \(Var(Y_i) = \sigma^2\), ou seja, a variância de \(Y_i = Y|\mathbf{X}=\mathbf{x}_i\) não depende dos valores observados \(\mathbf{x}_i\);
  • \(\mathbf{Y} = (Y_1, Y_2, \ldots, Y_n)^T\) é um vetor aleatório independente.

Nesse caso, com essas novas suposições, podemos escrever:

\[ Y_i = \beta_0 + \beta_1 x_{i,1} + \ldots + \beta_k x_{i,k} + \epsilon_i, \quad \epsilon_i \sim N(0,\sigma^2) \]

sendo \(\epsilon_i\) variáveis aleatórias independentes. Na forma matricial,

\[ \mathbf{Y} = X \boldsymbol{\beta} + \boldsymbol\epsilon, \quad \boldsymbol\epsilon \sim N_p(\mathbf{0},\sigma^2 I) \Rightarrow \mathbf{Y} \sim N_n(X\boldsymbol{\beta},\sigma^2 I) \] sendo \(I\) a matriz identidade \({n \times n}\). Entende-se por \(\boldsymbol\epsilon \sim N_p(\mathbf{0},\sigma^2 I)\) que o vetor aleatório \(\boldsymbol\epsilon\) segue a distribuição normal \(n\)-variada com vetor de médias nulo e matriz de correlação dada por \(\sigma^2 I\).

Veja que mais suposições foram acrescentadas ao modelo. Caso sejam suposições razoáveis, apersar de maior restrição ao modelo, ganhamos inferências sobre os parâmetros, como será visto a seguir. Todos os resultados apresentados referentes ao Modelo Estatístico podem ser aprofundados com os livros Kutner et al. (2005) ou Montgomery et al. (2012).

2.5.1 Estimador para \(\boldsymbol{\beta}\) por Máxima Verossimilhança

Uma vez assumida uma distribuição para a variável aleatória \(Y_i\), \(Y_i \sim N(\mu_i = \mathbf{l}_i\boldsymbol{\beta} \ , \ \sigma^2_i = \sigma^2)\) podemos buscar o estimador de Máxima Verossimilhança para o vetor de parâmetros \(\boldsymbol{\hat\beta}\), isto é, os valores de \(\boldsymbol{\hat\beta}\) que maximizam a função de verossimilhança considerando a amostra observada.

Para isso, primeiro, a Função de Verossimilhança: \[ \begin{array}{rclcl} L(\boldsymbol{\beta} | \mathbf{Y} = \mathbf{y}, \mathbf{X} = \mathbf{x}) & = & \prod_{i=1}^n \dfrac{1}{\sqrt{2 \pi \sigma_i^2}}e^{-\dfrac{1}{2\sigma^2_i}(y_i - \mu_i)^2} & = & \prod_{i=1}^n \dfrac{1}{\sqrt{2 \pi \sigma^2}} e^{-\dfrac{1}{2\sigma^2}(y_i - \mathbf{l}_i\boldsymbol{\beta})^2} \\ & = & \left(\dfrac{1}{\sqrt{2 \pi \sigma^2}}\right)^n e^{-\dfrac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mathbf{l}_i\boldsymbol{\beta})^2} & = & \left(\dfrac{1}{\sqrt{2 \pi \sigma^2}}\right)^n e^{-\dfrac{1}{2\sigma^2} (\mathbf{y} - X\boldsymbol{\beta})^T(\mathbf{y} - X\boldsymbol{\beta})} \end{array} \]

Por definição, o estimador de máxima verossimilhança para o vetor de parâmetros \(\boldsymbol{\beta}\) é dado por:

\[ \boldsymbol{\hat\beta} = \arg\max L(\boldsymbol{\beta} | \mathbf{Y}, \mathbf{X}) \] Veja que \(\max L(\boldsymbol{\beta})\) é o mesmo que \(\min (\mathbf{y} - X\boldsymbol{\beta})^T(\mathbf{y} - X\boldsymbol{\beta})\) e este último problema de minimização é justamente o problema de mínimos quadrados. Ou seja, o ponto que maximiza a função de verossimilhança \(L\) é o mesmo que minimiza a função de custo \(J\). \[ \boldsymbol{\hat\beta} = (X^TX)^{-1}X^T\mathbf{Y} \]

E, novamente, os valores previstos a partir de \(\boldsymbol{\hat\beta}\) são

\[ \mathbf{\hat{Y}} = X \boldsymbol{\hat\beta} = X (X^TX)^{-1}X^T\mathbf{Y} = H \mathbf{Y} \]

2.5.2 Propriedades da Matriz H

Vale ressaltar que estamos considerando \(H = X (X^TX)^{-1}X^T\), sendo \(X\) uma matriz \(n \times p\) tal que existe \((X^TX)^{-1}\). Aqui o número de colunas de \(X\) foi representado por \(p\) e, considerando a notação anterior, \(p = k+1\). Resultados simples de álgebra linear nos mostram que \(X^TX\) e \((X^TX)^{-1}\) são matrizes simétricas.

Proposição 1 \(H\) é uma matriz simétrica, isto é, \(H^T = H\).

Comprovação. \(H^T = (X (X^TX)^{-1}X^T)^T = (X^T)^T ((X^TX)^{-1})^T X^T) = X (X^TX)^{-1} X^T = H\)

Proposição 2 \(H\) é uma matriz idempotente, isto é, \(H^2 = H\).

Comprovação. \(H^2 = (X (X^TX)^{-1}X^T)^2 = X (X^TX)^{-1}X^TX (X^TX)^{-1}X^T = X(X^TX)^{-1}X^T = H\), ou seja, \(H\) é idempotente Toda matriz idempotente define uma projeção.

Proposição 3 Se \(H\) é simétrica e idempotente, então \(I-H\) é simétrica e idempotente.

Comprovação. \((I-H)^T = T^T - H^T = I - H\). \((I-H)^2 = (I-H)(I-H) = I^2 -IH - HI + H^2 = I -2H + H = I - H\).

Proposição 4 \(HX = H\) e \(H\) é a projeção no subespaço definido pelos vetores colunas de \(X\) e o posto de \(H\) é \(p\).

Comprovação. Toda matriz idempotente é uma matriz de projeção e o seu posto é a dimensão do subespaço de projeção.

Como \(X^TX\) é inversível, ela tem posto completo, ou seja, seu posto é \(p\). Com isso concluíamos que o posto de \(H\) só pode ser no máximo \(p\).

Veja que \(HX = X (X^TX)^{-1}X^T X = X\). Ou seja, cada vetor coluna \(X\) é um autovetor de \(H\) associado ao autovalor \(1\). Como \(X^TX\) é inversível sabemos que as colunas de \(X\) são linearmente independente e definem um subespaço de dimensão \(p\). Logo, o posto de \(H\) tem que ser maior ou igual a \(p\). Mas já vimos que o posto de \(H\) é no máximo \(p\) \(\Rightarrow\) posto de \(H\) é igual a \(p\).

2.5.3 Distribuição Amostral de \(\boldsymbol{\hat\beta}\)

Uma vez definido o estimador \(\boldsymbol{\hat\beta}\) como transformação do vetor aleatório \(Y\), cuja distribuição é conhecida, é possível enocontrar a distribuição amostral do estimador \(\boldsymbol{\hat\beta}\).

\[ \boldsymbol{\hat\beta} = (X^TX)^{-1}X^T\mathbf{Y} \] Já foi visto que \(\mathbf{Y} \sim N_n(X\boldsymbol{\beta},\sigma^2 I)\) . Como \(\boldsymbol{\hat\beta}\) é uma combinação linear de \(\mathbf{Y}\), podemos concluir que \(\boldsymbol{\hat\beta}\) também tem distribuição normal multivariada.

\[ E(\boldsymbol{\hat\beta}) = E\left( (X^TX)^{-1}X^T\mathbf{Y} \right) = (X^TX)^{-1}X^T E\left( \mathbf{Y} \right) = (X^TX)^{-1}X^T X\boldsymbol{\beta} = \boldsymbol{\beta} \] Ou seja, \(\boldsymbol{\hat\beta}\) é um estimador não tendencioso para o vetor de parâmetros \(\boldsymbol{\beta}\).

\[ \begin{array}{ll} Var(\boldsymbol{\hat\beta}) &= Var\left( (X^TX)^{-1}X^T\mathbf{Y} \right) \\ &= (X^TX)^{-1}X^T Var\left( \mathbf{Y} \right) \left( (X^TX)^{-1}X^T \right)^T \\ &= (X^TX)^{-1}X^T \sigma^2 I X ( (X^TX)^{-1} )^T \\ &= \sigma^2 (X^TX)^{-1}X^T X ( (X^TX)^{-1} )^T \\ &= \sigma^2 ( (X^TX)^{-1} )^T \\ &= \sigma^2 (X^TX)^{-1} \\ \end{array} \] Dessa forma,

\[ \boldsymbol{\hat\beta} \sim N_n \left(\boldsymbol{\beta} \ , \ \sigma^2 (X^TX)^{-1} \right) \]

2.5.4 O estimador para \(\sigma^2\)

No modelo estatístico é assumido que \(Var(\mathbf{Y}) = \sigma^2I\), sendo \(\sigma^2\) um parâmetro desconhecido. Vejamos que o \(\hat\sigma^2\) apresentado a seguir é um estimador não tendencioso para \(\sigma^2\).

\[ \hat\sigma^2 = \dfrac{SSE}{n-k-1} = \dfrac{\sum_{i=1}^n ( Y_i - \hat{Y}_i)^2}{n - k - 1} \]

Veja que \[ \begin{array}{lll} SSE & = & \sum_{i=1}^n ( Y_i - \hat{Y}_i)^2 \\ & = & (\mathbf{Y} - \mathbf{\hat{Y}})^T (\mathbf{Y} - \mathbf{\hat{Y}})\\ & = & (\mathbf{Y} - H\mathbf{Y})^T (\mathbf{Y} - H\mathbf{{Y}})\\ & = & (X\boldsymbol{\beta} + \boldsymbol{\epsilon}- H(X\boldsymbol{\beta} + \boldsymbol{\epsilon}))^T (X\boldsymbol{\beta} + \boldsymbol{\epsilon} - H(X\boldsymbol{\beta} + \boldsymbol{\epsilon}))\\ & = & (X\boldsymbol{\beta} + \boldsymbol{\epsilon}- HX\boldsymbol{\beta} - H\boldsymbol{\epsilon})^T (X\boldsymbol{\beta} + \boldsymbol{\epsilon} - HX\boldsymbol{\beta} - H\boldsymbol{\epsilon})\\ & = & (X\boldsymbol{\beta} + \boldsymbol{\epsilon}- X\boldsymbol{\beta} - H\boldsymbol{\epsilon})^T (X\boldsymbol{\beta} + \boldsymbol{\epsilon} - X\boldsymbol{\beta} - H\boldsymbol{\epsilon})\\ & = & (\boldsymbol{\epsilon} - H\boldsymbol{\epsilon})^T (\boldsymbol{\epsilon} - H\boldsymbol{\epsilon})\\ & = & ((I - H)\boldsymbol{\epsilon})^T ((I - H)\boldsymbol{\epsilon})\\ & = & \boldsymbol{\epsilon}^T(I - H)^T (I - H)\boldsymbol{\epsilon}\\ & = & \boldsymbol{\epsilon}^T (I - H)\boldsymbol{\epsilon}\\ \end{array} \] Sendo \(\boldsymbol{\epsilon} \sim N_n(\mathbf{0},\sigma^2I)\). Podemos escrever, \[ \dfrac{SSE}{\sigma^2} = \dfrac{1}{\sigma^2}\boldsymbol{\epsilon}^T (I - H)\boldsymbol{\epsilon} = \left(\dfrac{\boldsymbol{\epsilon}}{\sigma} \right)^T (I - H)\dfrac{\boldsymbol{\epsilon}}{\sigma} \] Como \(\dfrac{\boldsymbol{\epsilon}}{\sigma} \sim N_n(\mathbf{0},I)\) e \((I - H)\) é matriz simétrica idempotente de posto \(n-k-1\), podemos concluir que \[ \dfrac{SSE}{\sigma^2} \sim \chi^2_{n - k - 1} \] pois se \(\mathbf{Z} \sim N_n(\mathbf{0},I)\), \(A\) matriz simétrica, \(n \times n\) e idempotente de posto \(r\) então \(\mathbf{Z}^TA\mathbf{Z} \sim \chi^2_r\) (Proposição 5 da Seção 2.9).

Com isso, \[ E(SSE) = E\left( \dfrac{\sigma^2 SSE}{\sigma^2} \right) = \sigma^2 E\left(\dfrac{SSE}{\sigma^2}\right) = \sigma^2 (n - k - 1) \]

\[ E(\hat\sigma^2) = E\left(\dfrac{SSE}{n-k-1}\right) = \dfrac{E(SSE)}{n-k-1} = \dfrac{\sigma^2 (n - k - 1)}{n-k-1} = \sigma^2 \] \(\hat\sigma^2\) é um estimador não tendencioso para \(\sigma^2\).

2.5.5 Quantidade Pivotal para \(\beta_j\)

Vamos mostrar que \(Q\) definida a seguir é uma quantidade pivotal para o parâmetro \(\beta_j\):

\[ Q = \dfrac{\hat\beta_j - \beta_j}{\sqrt{\hat\sigma^2 \ C_{j+1,j+1}} } \sim t_{n-k-1} \]

Considere \[ {Z_j} = \dfrac{{\hat\beta_j} - {\beta}}{\sqrt{\sigma^2C_{j+1,j+1}}} \quad \text{ e } \quad W = \dfrac{(n - k - 1)\hat\sigma^2}{\sigma^2}. \] Veja que \(Z_j\) é função de \(\hat\beta_j\), que por sua vez é função de \(\mathbf{\hat{Y}}\):

\[ \mathbf{\hat{Y}} = X \boldsymbol{\hat\beta} \quad \Rightarrow \quad X^T\mathbf{\hat{Y}} = X^TX \boldsymbol{\hat\beta} \quad \Rightarrow \quad \boldsymbol{\hat\beta} = (X^TX)^{-1}X^T\mathbf{\hat{Y}} \] Veja que \(W\) é função de \(\hat\sigma^2\), que por sua vez é função de \(\mathbf{Y} - \mathbf{\hat{Y}}\). Como \(\mathbf{\hat{Y}}\) e \(\mathbf{Y} - \mathbf{\hat{Y}}\) são vetores aleatórios independentes (Proposição 6 da Seção 2.9), podemos concluir que as variáveis aleatórias \(Z_j\) e \(W\) são variáveis aleatórias independentes.

Uma vez conhecido a distribuição do estimador \(\boldsymbol{\hat\beta}\) podemos verificar que para cada \(\beta_j\), \(j = 0, \ldots, k\), temos \[ \hat\beta_j \sim N(\beta_j \ , \ \sigma^2 C_{j+1,j+1} ) \] sendo \(C_{j+1,j+1}\) o valor na posição \((j+1, j+1)\) da matriz \(C = (X^TX)^{-1}\). Logo, \(Z_j \sim N(0,1)\). Além disso, foi visto que \[ \dfrac{SSE}{\sigma^2} \sim \chi^2_{n-k-1} \ \Rightarrow W = \dfrac{(n - k - 1)\hat\sigma^2}{\sigma^2} \sim \chi^2_{n-k-1} \]

Dessa forma, \[ \dfrac{Z}{\sqrt{W/(n - k - 1)}} \sim t_{n - k - 1} \]

E veja que \[ \begin{array}{lll} \dfrac{Z}{\sqrt{W/(n - k - 1)}} &=& \dfrac{\dfrac{\hat\beta_j - \beta_j}{\sqrt{\sigma^2C_{j+1,j+1}}}}{\sqrt{\dfrac{(n - k - 1)\hat\sigma^2}{\sigma^2 (n - k - 1)}}} = \dfrac{\dfrac{\hat\beta_j - \beta_j}{\sqrt{\sigma^2C_{j+1,j+1}}}}{\sqrt{\dfrac{\hat\sigma^2}{\sigma^2}}}\\ &=& {\dfrac{\hat\beta_j - \beta_j}{\sqrt{\sigma^2C_{j+1,j+1}}}} {\sqrt{\dfrac{\sigma^2}{\hat\sigma^2}}} = {\dfrac{\hat\beta_j - \beta_j}{\sqrt{\hat\sigma^2 \ C_{j+1,j+1}}}} = Q \sim t_{n-k-1} \end{array} \]

2.5.6 Intervalo de Confiança

A partir de \(Q\) fica fácil de derivar os intervalos de confiança para os parâmetros \(\beta_j\):

\[ \begin{array}{lll} Q \sim t_{n-k-1} &\Rightarrow & P(-t_{n-k-1,1-\frac{\alpha}{2}} < Q < t_{n-k-1,1-\frac{\alpha}{2}}) = 1 - \alpha \\ &\Rightarrow & P(-t_{n-k-1,1-\frac{\alpha}{2}} < {\dfrac{\hat\beta_j - \beta_j}{\sqrt{\hat\sigma^2 \ C_{j+1,j+1}}}} < t_{n-k-1,1-\frac{\alpha}{2}}) = 1 - \alpha \\ &\Rightarrow & P(-t_{n-k-1,1-\frac{\alpha}{2}} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}} < \hat\beta_j - \beta_j < t_{n-k-1,1-\frac{\alpha}{2}} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}} ) = 1 - \alpha \\ &\Rightarrow & P(- \hat\beta_j - t_{n-k-1,1-\frac{\alpha}{2}}\sqrt{\hat\sigma^2 \ C_{j+1,j+1}} < - \beta_j < - \hat\beta_j + t_{n-k-1,1-\frac{\alpha}{2}}\sqrt{\hat\sigma^2 \ C_{j+1,j+1}}) = 1 - \alpha \\ &\Rightarrow & P( \hat\beta_j - t_{n-k-1,1-\frac{\alpha}{2}} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}} < \beta_j < \hat\beta_j + t_{n-k-1,1-\frac{\alpha}{2}} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}} ) = 1 - \alpha \\ \end{array} \]

\[IC({\beta_j})_{100(1-\alpha)\%} = \left[ \hat{\beta}_j - t_{1 - \frac{\alpha}{2}, n - k - 1} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}} \ \ , \ \ \hat{\beta}_j + t_{1 - \frac{\alpha}{2}, n - k - 1} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}} \right]\] sendo \(C\) é a matriz \((X^TX)^{-1}\) e \(C_{j+1,j+1}\) o valor na sua posição \((j+1,j+1)\).

2.5.7 Teste de Hipótese

Uma vez deduzido um intervalo de confiança bilateral para \(\beta_j\), podemos criar regras de decisão para testar as hipóteses:

\[H_0: \beta_j = 0 \quad \text{ versus } \quad H_1: \beta_j \neq 0.\]

A regra de decisão:

“Rejeita \(H_0\) se \(0 < \hat{\beta}_j - t_{1 - \frac{\alpha}{2}, n - k - 1} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}}\) ou \(0 > \hat{\beta}_j + t_{1 - \frac{\alpha}{2}, n - k - 1} \sqrt{\hat\sigma^2 \ C_{j+1,j+1}}\), isto é, rejeitamos \(H_0\) se o intervalo de confiança para \(\beta_j\) não contém o valor 0”

resulta em um teste com nível de confiança \(\alpha\).

2.6 Avaliação do Modelo

Em muitas situações nos perguntamos: será que o modelo obtido é bom? Para responder essa pergunta faremos o uso de medidas de qualidade para a avaliação do mdoelo.

Soma dos Erros Quadráticos (SSE - Sum of Squared Errors)

A Soma dos Erros Quadráticos, que já foi definida anteriormente, muitas vezes também é chamada de SSE.

\[ SSE = \sum_{i=1}^n (y_i - \hat{y}_i)^2\] Quanto menor o valor do SSE, melhor será o desempenho do modelo de regressão.

Erro Quadrático Médio (MSE - Mean Squared Error)

O Erro Quadrático Médio é definido por

\[ MSE = \dfrac{1}{N} \sum_{i=1}^N (y_i - \hat{y}_i)^2\] E assim como o SSE, quanto menor o valor do MSE, melhor será o desempenho do modelo de regressão.

Raiz do Erro Quadrático Médio (RMSE - Root Mean Squared Error)

Como o nome já diz, é a raíz da medida anterior, o MSE. A sua vantagem é que quando estraímos a raíz a unidade da medida de erro passa a ser a mesma unidade da veriável resposta \(Y\), o que facilita a sua interpretação.

\[ RMSE = \sqrt{\dfrac{1}{N} \sum_{i=1}^N (y_i - \hat{y}_i)^2}\] Quanto menor o valor do RMSE, melhor será o desempenho do modelo de regressão.

Erro Absoluto Médio (MAE - Mean Absotule Error)

É mais uma medida de erro que considera agora a média dos erros em valores absolutos. ão.

\[ MAE = \dfrac{1}{N} \sum_{i=1}^N |y_i - \hat{y}_i|\] Quanto menor o valor do MAE, melhor será o desempenho do modelo de regressão.

Coeficiente de Determinação \(R^2\)

O Coeficiente de Determinação \(R^2\) é uma medida que compara o desempenho do modelo de regressão obtido com o desempenho de um modelo mais simples possível, aqueçe que prevê o valor de \(Y\) pela média dos valores observados para essa variável. Dessa forma ele nos dá um valor comparativo pela própria natureza da medida, vejamos.

Primeiro definimos a Soma Total dos Quadrado (SST - Total Sum of Squares), que é a soma dos quadrados dos erros do modelo trivial, dado pela média amostral da variável resposta.

\[ SST = \sum_{i=1}^n (y_i - \bar{y})^2 \] O Coeficiente de Determinação \(R^2\) é então definido por:

\[ R^2 = 1 - \dfrac{SSE}{SST}\] Veja algumas possíveis interpretações do \(R^2\):

  • Se \(R^2 = 0 \Rightarrow\) as previsões do modelo são tão boas (ou tão ruins) quanto as do modelo trivial.

  • Se \(R^2 > 0 \Rightarrow\) as previsões do modelo são melhores que as do modelo trivial.

  • Se \(R^2 < 0 \Rightarrow\) as previsões do modelo são piores que as do modelo trivial.

  • Quanto maior o valor do \(R^2\), melhor a qualidade das previsõs.

  • O valor de \(R^2\) é sempre menor que \(1\). Ou seja, quanto mais perto de 1 for o \(R^2\) melhor o modelo de regressão.

2.7 Seleção de Variáveis

Quando dois modelos apresentam capacidade preditiva semelhante, deve-se preferir aquele que utiliza o menor número de variáveis ou possui menor complexidade (princípio da parcimônia). Por isso, uma das etapas importante do modelo de regressão liner é a seleção das variáveis.

O aumento desnecessário da complexidade de um modelo, isto é, o aumentono número de parâmetros sem a melhora significativa da capacidade preditiva, pode trazer problemas como o sobreajuste (overfitting) e a multicolinearidade.

Sobreajuste

O sobreajuste ocorre quando um mdoelo se ajusta aos dados de treinamento sem que tenha entendido a relação entre as variáveis, somente pelo ecesso de parâmetros. Com isso as medidas de qualidade na base de treino são boas, mas quando analisamos a sua capacidade preditiva em novos dados, o seu desempenho é bem pior.

Multicolinearidade

A multicolinearidade ocorre quando duas ou mais variáveis explicativas são muito correlacionadas. A forte correlação entre elas já é uma justificativa para ficar com uma e não todas. Além disso, de forma práticas, manter variáveis muito correlacionadas na base faz com que a matriz \(X^TX\) tenha algum autovalor próximo de zero. Isso gera alguns problemas na estimativa, como por exemplo, um auto valor perto de zero faz com que \(\det(X^TX)\) seja muito pequeno e, logo, a sua inversa com v alores instáveis uma vez que a inversa é produto de \(\dfrac{1}{\det(X^TX)}\) valores grandes. Isso traz valores grandes (e instáveis) não só para a estimativa dos \(\beta\) como também aumenta a amplitude dos seus intervalos de confiança.

Métodos de Seleção

A seguir serão apresentados métodos de seleções de variáveis. Mas, antes mesmo de aplicar os métodos, é possível fazer uma seleção prévia entre as variáveis quantitativas de forma a eliminar o problema de multicolinearidade. Para isso, construa a matriz de correlação entre as variáveis explicativas quantitativas e verifique se existe alguma correlação com valor absoluto maior que 0,7. Caso exista, elimine alguma das variáveis com alta correlação com outras e repita o processo até que a matriz de correlação não apresente valores absolutos maiores que 0,7.

Critério de inclusão ou exclusão de variáveis

Os critérios AIC (Akaike Information Criterion) e BIC (Bayesian Information Criterion) são medidas utilizadas para comparar modelos, equilibrando qualidade do ajuste e complexidade. Ambos penalizam modelos que utilizam muitos parâmetros, seguindo o princípio da parcimônia.

\[AIC = -2\log(L) + 2p \quad \text{ e } \quad BIC = -2\log(L) + \log(n)p\]

sendo \(L\) a função de verossimilhança do modelo; \(p = k+1\) o número de parâmetros estimados (incluindo o intercepto) e \(n\) é o número de observações.

Quanto menor os valores de AIC e/ou BIC, melhor é o modelo de acordo com o critério.

Forward Selection

O método Forward Selection começa com um modelo sem variável alguma. Na primeira etapa são comparados os desempenhos (AIC ou BIC) de todos os modelo simples (aqueles com uma variável) e seleciona-se o com melhor desepenho. Na segunda etapa, o modelo existente é acrescido com uma nova variável e a variável selecionada para entrar é aquela que proporciona melhor melhoria. Em cada etapa, avalia-se qual das variáveis ainda não selecionadas proporciona a maior melhoria no modelo. Essa variável é então adicionada ao modelo. O processo termina quando não há melhoria com o acréscimo de uma variável.

Backward Elimination

O método Backward Elimination começa com o modelo completo, isto é, com todas as variáveis disponíveis. Na primeira etapa são comparados os desempenhos (AIC ou BIC) de todos os modelo gerados a partir do modelo completo com a eliminação de uma variável e seleciona-se aquele o com melhor desepenho. Na segunda etapa, são comparados os modelos gerados a partir do modelo final da etapa anterior com a eliminação de uma nova variável. A variável selecionada para sair é aquela que proporciona melhor melhoria. Em cada etapa, identifica-se a variável menos relevante e ela é removida. O processo continua até que a retirada de uma nova variável não melhora o critério escolhido para comparação.

Stepwise Selection

O Stepwise Selection combina as duas estratégias anteriores.

O procedimento normalmente começa como o Forward Selection, adicionando variáveis uma a uma. Entretanto, após cada inclusão, verifica-se se alguma das variáveis já presentes deixou de ser importante. Se isso ocorrer, essa variável é removida.

2.8 Regularização

A regularização, adota em vários problemas de regressão, não só na regressão linear, são técnicas que buscam diminuir a complexidade do modelo de forma a evitar sobreajuste. Isso pode ser feito de várias maneira, e para os modelos de regressão linear, a ideia principal é penalizar coeficientes grandes no momento de estimativa dos parâmetros.

Ridge

Este é um método de regularização que substituí a função de custo \(J\), antes definida pelo soma dos quadrados dos erros, para

\[J(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^k \beta_j^2\]

sendo \(\lambda > 0\) um parâmetro de regularização que controla o grau de penalização aplicado aos coeficientes: quanto maior o valor de \(\lambda\), maior será essa penalização. Como consequência, os coeficientes são reduzidos em direção a zero (e também uns em direção aos outros) (Hastie 2009).

As soluções obtidas pelo método Ridge não são invariantes à escala das variáveis explicativas. Por esse motivo, é prática comum padronizar as variáveis de entrada antes de resolver o problema de otimização. Além disso, observe que o intercepto \(\beta_0\) foi excluído do termo de penalização.

Lasso

O LASSO é um método de regularização, assim como o Ridge, mas apresenta diferenças sutis e, ao mesmo tempo, importantes. A estimativa dos coeficientes pelo método LASSO é definida por

\[J(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^k |\beta_j|\]

O termo de penalização do Ridge, dado pela soma dos \(\beta\) ao quadrado, reduz a magnitude de todos os coeficientes, aproximando-os de zero, mas não faz com que eles sejam exatamente iguais a zero (exceto no caso extremo em que \(\lambda \to \infty\)) (Huang 2014). Embora isso não represente, em geral, um problema para a capacidade preditiva do modelo, pode dificultar sua interpretação, especialmente quando o número de variáveis explicativas é elevado.

Aí entra a grande vantagem: o LASSO realiza simultaneamente regularização e seleção automática de variáveis. A penalização do LASSO, dada pela soma dos valores absolutos dos parâmetros, tem como efeito forçar alguns deles ao valor zero quando o parâmetro de regularização \(\lambda\) é razoavelmente grande.

2.9 Apêndice

Proposição 5 Seja \(\mathbf{Z} \sim N_n(\mathbf{0},I)\) e \(A\) matriz simétrica, idempotente de posto \(r\). Então a variável aleatória \(W = \mathbf{Z}^T A \mathbf{Z} \sim \chi^2_r\).

Comprovação. Como \(A\) é matriz simétrica, pelo Teorema Espectral, existe uma matriz ortogonal \(P\) tal que \(PP^T = I\) e \(A = P \Lambda P^T\), sendo \(\Lambda = diag(\lambda_1, \lambda_2, \ldots, \lambda_n)\), matriz diagonal com os autovalores de \(A\) na diagonal.

Como \(A\) é idempotente, \(A\) é uma projeção e todos os seus autovalores são 0 ou 1. Isso também poderia ser verificado por: \[ A^2 = A \ \Rightarrow \ P \Lambda P^T P \Lambda P^T = P \Lambda P^T \ \Rightarrow P \Lambda^2 P^T = P \Lambda P^T \ \Rightarrow \Lambda^2 = \Lambda \ \Rightarrow \ \lambda_i = 1 \text{ ou } \lambda_i = 0 \]

O posto de uma matriz é o número de autovalores não nulos da matriz. Como \(A\) só tem autovalores 0 ou 1, se \(posto(A) = r\) isso significa que \(A\) tem \(r\) autovalores iguais a 1 e \(n - r\) autovalores iguais a 0.

Podemos, portanto, escrever \[ \Lambda = \begin{pmatrix} I_r & 0\\ 0 & 0 \end{pmatrix}_{n \times n} \] Sendo \(I_r\) a matriz de identidade \(r \times r\). E \[ W = \mathbf{Z}^T A \mathbf{Z} = \mathbf{Z}^T P \Lambda P^T \mathbf{Z} \] Defina \(\mathbf{Y} = P^T \mathbf{Z}\). \(\mathbf{Y}\) é vetor aleatório com distribuição normal multivariada, já \(\mathbf{Z} \sim N_n(\mathbf{0},I)\). Vejamos seu vetor de médias e matriz de variância.

\[ E(Y) = E(P^T \mathbf{Z}) = P^T E(\mathbf{Z}) = 0 \quad \text{ e } \quad Var(Y) = Var(P^T \mathbf{Z}) = P^T Var(\mathbf{Z}) P = P^T I P = P^TP = I \] Ou seja, \(\mathbf{Y} \sim N_n(\mathbf{0},I)\) e \[ W = \mathbf{Z}^T P \Lambda P^T \mathbf{Z} = \mathbf{Y}^T \Lambda \mathbf{Y} \] Então, \(W\) é soma do quadrado de \(r\) normais padrão independentes, logo, \(W \sim \chi_r^2\).

Proposição 6 Sejam

  • \(\mathbf{Y} = X \boldsymbol{\beta} + \boldsymbol{\epsilon}\), sendo \(X\) uma matriz \(n \times p\) tal que existe \((X^TX)^{-1}\) e \(\boldsymbol{\epsilon} \sim N_n(\mathbf{0}, \sigma^2I)\);

  • \(\mathbf{\hat{Y}} = X(X^TX)^{-1}X^T \mathbf{Y} = H \mathbf{Y}\);

Então os vetores aleatórios \(\mathbf{\hat{Y}}\) e \(\mathbf{Y} - \mathbf{\hat{Y}}\) são independentes.

Comprovação. Primeiro veja que \(\mathbf{Y}\) é uma transformação linear de um vetor aleatório com distribuição normal multivariada, logo, é um vetor aleatório com distribuição normal multivariada.

\[ E(\mathbf{Y}) = X \boldsymbol{\beta} \quad \text{ e } \quad Var(\mathbf{Y}) = \sigma^2I \]

Além disso, \(\mathbf{\hat{Y}} = H\mathbf{Y}\) também é vetor aleatório com distrubição normal multivariada.

\[ E(\mathbf{\hat{Y}}) = HE(\mathbf{Y}) = HX \boldsymbol{\beta} = X\boldsymbol{\beta} \quad \text{ e } \quad Var(\mathbf{\hat{Y}}) = HVar(\mathbf{Y})H^T = H \sigma^2IH^T = \sigma^2H^2 = \sigma^2H \] Considere o vetor \(\mathbf{\tilde{Y}}\) definido por: \[ \mathbf{\tilde{Y}} = \begin{pmatrix} \hat{Y_1}\\ \hat{Y_2}\\ \vdots\\ \hat{Y_n}\\ Y_1 - \hat{Y}_1\\ Y_2 - \hat{Y}_2\\ \vdots\\ Y_n - \hat{Y}_n \end{pmatrix}_{2n \times 1} = \begin{pmatrix} \begin{pmatrix} \\ \\ & \mathbf{\hat{Y}} & \\ \\ \\ \end{pmatrix} \\ \begin{pmatrix} \\ \\ \mathbf{Y-\mathbf{\hat{Y}}} \\ \\ \\ \end{pmatrix} \end{pmatrix} = \begin{pmatrix} \begin{pmatrix} \\ \\ & & & H & & & \\ \\ \\ \end{pmatrix}_{n \times n} \\ \begin{pmatrix} \\ \\ & & I - H & &\\ \\ \\ \end{pmatrix}_{n \times n} \\ \end{pmatrix}_{2n \times n} \begin{matrix} \begin{pmatrix} \\ \\ \mathbf{Y} \\ \\ \\ \end{pmatrix}_{n \times 1} \\ \begin{matrix} \\ \\ \\ \\ \\ \end{matrix} \end{matrix} \]

Como \(\mathbf{\tilde{Y}}=A\mathbf{Y}\), isto é, é escrito como transformação linear de \(\mathbf{{Y}} \sim N_n(X\boldsymbol{\beta},\sigma^2I)\), podemos afirmar que \(\mathbf{\tilde{Y}}\) também tem distribuição normal multivariada tal que \(E(\mathbf{\tilde{Y}})=AE(\mathbf{Y})\) e \(Var(\mathbf{\tilde{Y}}) = AVar(\mathbf{\hat{Y}})A^T\). Além disso, podemos escrever

\[ Var(\mathbf{\tilde{Y}}) = \begin{pmatrix} \begin{pmatrix} \\ & \ \ Var(\mathbf{\hat{Y}}) \ \ & \\ \\ \end{pmatrix} & \begin{pmatrix} \\ Cov(\mathbf{\hat{Y}},\mathbf{Y} - \mathbf{\hat{Y}}) \\ \\ \end{pmatrix} \\ \begin{pmatrix} \\ Cov(\mathbf{Y} - \mathbf{\hat{Y}},\mathbf{\hat{Y}}) \\ \\ \end{pmatrix} & \begin{pmatrix} \\ \ \ Var(\mathbf{Y} - \mathbf{\hat{Y}}) \ \ \\ \\ \end{pmatrix} \end{pmatrix} \] sendo \[ \begin{array}{lll} Cov(\mathbf{Y} - \mathbf{\hat{Y}},\mathbf{\hat{Y}}) &= Cov(\mathbf{Y},\mathbf{\hat{Y}}) - Cov(\mathbf{\hat{Y}},\mathbf{\hat{Y}}) \\ &= Cov(\mathbf{Y},H\mathbf{{Y}}) - Cov(\mathbf{\hat{Y}},\mathbf{\hat{Y}}) \\ &= Cov(\mathbf{Y},\mathbf{{Y}})H ^T - Var(\mathbf{\hat{Y}}) \\ &= \sigma^2H - \sigma^2H = 0 \\ \end{array} \] Como trata-se de um vetor aleatório com distribuição normal multivariada, correlação nula implica em independência. Nesse caso, os vetores aleatórios \(\mathbf{\hat{Y}}\) e \(\mathbf{Y} - \mathbf{\hat{Y}}\) são idependentes.

Referências

Hastie, Trevor. 2009. The elements of statistical learning: data mining, inference, and prediction. Springer.
Huang, Jianhua Z. 2014. An Introduction to Statistical Learning: With Applications in R By Gareth James, Trevor Hastie, Robert Tibshirani, Daniela Witten: Publisher: Springer, 2013. ISBN 978-1-4614-7137-0. Springer.
Kutner, Michael H., Christopher J. Nachtsheim, John Neter, e William Li. 2005. Applied Linear Statistical Models. 5th ed. McGraw-Hill/Irwin.
Montgomery, Douglas C., Elizabeth A. Peck, e G. Geoffrey Vining. 2012. Introduction to Linear Regression Analysis. 5th ed. A John Wiley & Sons.