1 Pacotes Utilizados na disciplina

require(matrixcalc)
require(Matrix)
require(ibd)
require(pracma)
require(matlib)
require(MASS)
require(tidyverse)
require(kableExtra)
require(DT)

2 Introdução aos modelos lineares

Vamos aprender sobre equações normais, Estimação e teste de hipóteses sobre parâmetros, modelos de regressão e estatística experimental.

Seja \(\boldsymbol{y}\) uma variável resposta (ou dependente) sendo explicada por um conjunto de variáveis explicativas (independentes), mais um termo de erro: \[ \begin{array}{l} \boldsymbol{y}=f(\boldsymbol{x})+\boldsymbol{\epsilon}=f(x_1,x_2,\dots,x_p)+\boldsymbol{\epsilon}, \mbox{ e considerando a observação do indivíduo }i \mbox{ tem os:}\\ y_i=f(x_{1i},x_{2i},\dots,x_{pi})+\epsilon_i \mbox{ significa que cada indivíduo tem seu erro associado} \end{array} \]

Exemplo 1 Funções lineares e não lineares nos parâmetros:

\[ \begin{array}{|l|l|}\hline \boldsymbol{y}=f(\boldsymbol{x})+\boldsymbol{\epsilon} & \mbox{classificação}\\ \hline y=\beta_0+\beta_1 x_{1}+\beta_2x_2 + \epsilon & \mbox{linear nos parâmetros}\\ y=\beta_0 +\beta_1 x_1+\beta_2 \log(x_2) + \epsilon & \mbox{linear nos parâmetros}\\ y=\beta_0 +\beta_1 x_1^2+\beta_2 \log(x_2) + \epsilon & \mbox{ linear nos parâmetros}\\ y=\beta_0 + x_1^{\beta_1}+\beta_2x_2+\epsilon & \mbox{não linear nos parâmetros}\\ y=\beta_0 + \beta_1^2x_1+\beta_2x_2+\epsilon & \mbox{não linear nos parâmetros}\\\hline \end{array} \]

2.1 O modelo Linear geral

Definições iniciais:

\[ \begin{array}{l} \boldsymbol{y}_{n \times 1} \mbox{ vetor coluna que caracteriza as observações}\\ \boldsymbol{X}_{n \times p} \mbox{ matrix de variáveis independentes, com posto } k \leq min\{n,p\}\\ \boldsymbol{\beta}_{p\times 1} \mbox{vetor de parâmetros desconhecidos a serem estimados}\\ \epsilon_{n \times 1} \mbox{ vetor de erros a serem estimados} \end{array} \] \[ \Downarrow\\ \boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\epsilon \]

2.2 Caracterização de modelos lineares:

\[ \begin{array}{|l|l|} \hline \mbox{forma do modelo} & \mbox{caracterização}\\ \hline I) y_i=\beta_0+\beta_1 x_i+\epsilon_i & \mbox{Regressão linear simples}\\ II) y_i=\beta_0+\beta_1 x_{1i}+ \beta_2 x_{2i}+\ldots+ \beta_k x_{ki}+\epsilon_i & \mbox{Regressão linear múltipla}\\ III) y_i=\mu+\alpha_i +\epsilon_{ij} & \mbox{Experimento inteiramente casualizado}\\ IV) y_{ij}=\mu + \alpha_i + \beta_j + \epsilon_{ij} & \mbox{Experimento em blocos ao acaso}\\ \hline \end{array} \]

2.3 Forma matricial

Para os casos I e II, estende-se a matriz \(\boldsymbol{X}\) com uma coluna de un’s à esquerda, e para os demais casos, cria-se esta matriz para atender a interpretação do experimento, de forma que o modelo geral é \(\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}\)

\[ \begin{array}{l}\hline \mbox{forma do modelo} \\ \hline I) \begin{pmatrix} y_1 \\y_2 \\\dots \\y_n \end{pmatrix} = \begin{pmatrix} 1 & x_1 \\1 & x_2 \\ \dots & \dots \\1 & x_n \end{pmatrix} \begin{pmatrix} \beta_0 \\\beta_1 \end{pmatrix}+\begin{pmatrix} \epsilon_1 \\\epsilon_2 \\\dots \\\epsilon_n \end{pmatrix}\\ II) \begin{pmatrix} y_1 \\y_2 \\\dots \\y_n \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{21} & \dots & x_{k1} \\ 1 & x_{12} & x_{22} & \dots & x_{k2}\\ \vdots & \vdots & \vdots & \vdots\\ 1 & x_{1n} & x_{2n} & \vdots & x_{kn} \end{pmatrix} \begin{pmatrix} \beta_0 \\\beta_1 \\ \vdots \\ \beta_k \end{pmatrix}+\begin{pmatrix} \epsilon_1 \\\epsilon_2 \\\dots \\\epsilon_n \end{pmatrix}\\ III) \begin{pmatrix} y_{11} \\y_{12}\\y_{13}\\y_{21} \\y_{22} \\y_{23} \end{pmatrix}= \begin{pmatrix} 1 & 1 & 0\\ 1 & 1 & 0 \\ 1& 1 & 0 \\ 1 & 0 & 1\\ 1 & 0 & 1 \\ 1& 0&1 \end{pmatrix} \begin{pmatrix} \mu \\ \alpha_1 \\ \alpha_2 \end{pmatrix}+ \begin{pmatrix} \epsilon_{11} \\\epsilon_{12} \\\epsilon_{13} \\ \epsilon_{21} \\ \epsilon_{22} \\ \epsilon_{23} \end{pmatrix}\\ IV) \begin{pmatrix} y_{11} \\y_{12}\\y_{13}\\y_{21} \\y_{22} \\y_{23} \end{pmatrix}= \begin{pmatrix} 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \mu \\ \alpha_1 \\ \alpha_2 \\ \beta_1 \\ \beta_2 \\ \beta_3\end{pmatrix}+ \begin{pmatrix} \epsilon_{11} \\\epsilon_{12} \\\epsilon_{13} \\ \epsilon_{21} \\ \epsilon_{22} \\ \epsilon_{23} \end{pmatrix} \end{array} \]

Exemplo 2 A solução de mínimos quadrados \((\boldsymbol{X}^t\boldsymbol{X})^-\boldsymbol{X}^ty\) para o caso I) se reduz a forma: \[ \begin{array}{l} \hat{\beta_0}=\bar{y}-\hat{\beta_1}\bar{x}\\ \hat{\beta_1}=\frac{\sum \limits_1^n x_iy_i-n\bar{x}\bar{y}}{\sum \limits_1^n x_i^2-n\bar{x}^2} \end{array} \]

Este o caso especial em que \(\boldsymbol{X}^t\boldsymbol{X}\) é uma matriz de posto completo. Por isso é importante que as variáveis explicativas do banco de dados não apresentem colinearidade. \[ (\boldsymbol{X}^t\boldsymbol{X})^-=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\\ \mbox{significa que a inversa generalizada é igual a inversa} \]

3 Modelo Linear de Gauss-Markov

O modelo linear geral com a seguinte caracterização para os erros: \[ \left\{ \begin{array}{l} E(\epsilon_i)=0\\ VAR(\epsilon_i)=\sigma^2, \forall i \end{array} \right. \mbox{ com } Cor(\epsilon_i,\epsilon_j)=0, \forall i\ne j \Rightarrow \left\{ \begin{array}{l} E(\boldsymbol{\epsilon})=\begin{pmatrix} 0 \\0 \\ \dots \\ 0 \end{pmatrix}\\ VAR(\boldsymbol{\epsilon})=\boldsymbol{I} \sigma^2 \end{array} \right. \]

3.1 Propriedades:

\[ \left\{ \begin{array}{l} E(\boldsymbol{y})=E(\boldsymbol{X\theta}+\boldsymbol{\epsilon})=\boldsymbol{X\theta}\\ VAR(\boldsymbol{y})=VAR(\boldsymbol{X\theta}+\boldsymbol{\epsilon})=\boldsymbol{I} \sigma^2\\ \hat{\boldsymbol{\beta}}=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t \boldsymbol{y} \mbox{ é a solução de mínimos quadrados ordinária} \end{array} \right. \]

Casos especiais:

Modelo Linear Ponderado de Gauss Markov: quando \[ \left\{ \begin{array}{l} E(\boldsymbol{\epsilon})=\begin{pmatrix} 0 \\0 \\ \dots \\ 0 \end{pmatrix}\\ VAR(\boldsymbol{\epsilon})=\boldsymbol{D} \sigma^2\mbox{ , com }\boldsymbol{D} \mbox{ matrix diagonal} \\ \Rightarrow \hat{\boldsymbol{\beta}}=(\boldsymbol{X}^t \boldsymbol{D}^ {-1} \boldsymbol{X})^{-1}\boldsymbol{X}^t \boldsymbol{D}^ {-1} \boldsymbol{y} \mbox{ é a solução de mínimos quadrados ponderado} \end{array} \right. \]
Modelo Linear Generalizado de Gauss Markov: quando \[ \left\{ \begin{array}{l} E(\boldsymbol{\epsilon})=\begin{pmatrix} 0 \\0 \\ \dots \\ 0 \end{pmatrix}\\ VAR(\boldsymbol{\epsilon})=\boldsymbol{\Omega} \sigma^2\mbox{ , com }\boldsymbol{\Omega} \mbox{ assumindo forma geral} \\ \Rightarrow \hat{\boldsymbol{\beta}}=(\boldsymbol{X}^t \boldsymbol{\Omega}^ {-1} \boldsymbol{X})^{-1}\boldsymbol{X}^t \boldsymbol{\Omega}^ {-1} \boldsymbol{y} \mbox{ é a solução de mínimos quadrados generalizada} \end{array} \right. \]
Modelo Linear de Gauss Markov Normal: quando \[ \left\{ \begin{array}{l} \boldsymbol{\epsilon} \sim N(\boldsymbol{\emptyset},\boldsymbol{I}\sigma^2) \mbox{(caso ordinário) }\\ \boldsymbol{\epsilon} \sim N(\boldsymbol{\emptyset},\boldsymbol{D}\sigma^2) \mbox{(caso ponderado)} \\ \boldsymbol{\epsilon} \sim N(\boldsymbol{\emptyset},\boldsymbol{\Omega}\sigma^2) \mbox{(caso generalizado) } \end{array} \right. \]

4 Sistema de equações normais

Modelos lineares são utilizados para identificar a relação entre uma variável dependente e uma ou mais variáveis independentes. O objetivo é encontrar o melhor hiperplano (ou linha, em duas dimensões) que se ajusta aos dados. Uma maneira de fazer isso é minimizar o erro quadrático, que é a soma dos quadrados das diferenças entre os valores observados e os valores estimados pelo modelo.

Modelo Linear Geral:

\[ \Downarrow\\ \boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\epsilon \]

Problema de minimizar de funções Consiste em minimizar a função f do erro quadrático: \[ \begin{array}{l} f(\boldsymbol{\beta})= \boldsymbol{\epsilon}^t\boldsymbol{\epsilon}\\ =(\boldsymbol{y}−\boldsymbol{X}\boldsymbol{\beta})^t (\boldsymbol{y}−\boldsymbol{X}\boldsymbol{\beta}) \end{array} \]

Deriva-se f com respeito a \(\boldsymbol{\beta}\) para obter as equações normais: \[ \boldsymbol{X}^t\boldsymbol{X} \boldsymbol{\beta} = \boldsymbol{X}^t \boldsymbol{y} \]

Observações:

Inversão de matriz: Para grandes conjuntos de dados, a inversa de \(\boldsymbol{X}^t\boldsymbol{X}\) é computacionalmente custosa.
Matriz singular: Se as colunas de \(\boldsymbol{X}\) são linearmente dependentes (multicolinearidade), \(\boldsymbol{X}^t\boldsymbol{X}\) não tem inversa.

Exemplo 3: O sistema \(\boldsymbol{X} \boldsymbol{\beta} = \boldsymbol{y}\) abaixo é inconsistente:

\[ \begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 0 \end{pmatrix} \begin{pmatrix} \beta_1\\ \beta_2\\ \beta_3 \end{pmatrix} = \begin{pmatrix} 2\\ 3\\ 5\\ 4 \end{pmatrix} \] \[ \begin{array}{l} \square \mbox{Inversa: não há!}\\ \square \mbox{Inversa generalizada (Penrose):}\\ \boldsymbol{X}^+=\begin{pmatrix} 1/6 & 1/6 & 1/6 & 1/6\\ 1/3 & 1/3 & -1/6 & -1/6\\ -1/6 & -1/6 & 1/3 & 1/3 \end{pmatrix} \end{array} \] \[ \begin{array}{l} \boldsymbol{X} \boldsymbol{X}^+.\boldsymbol{y}\ne\boldsymbol{y} \mbox{ pois:}\\ \begin{pmatrix} 5/2 \\ 5/2 \\ 9/2 \\ 9/2 \end{pmatrix} \ne \begin{pmatrix} 2\\ 3\\ 5\\ 4 \end{pmatrix} \end{array} \] Porém, expandindo o sistema com equações normais do tipo \(\boldsymbol{X}^t\boldsymbol{X} \boldsymbol{\beta} = \boldsymbol{X}^t \boldsymbol{y}\) \[ \begin{pmatrix} 4 & 2 & 2\\ 2 & 2 & 0\\ 2 & 0 & 2 \end{pmatrix} \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2\\ \hat{\beta}_3 \end{pmatrix} = \begin{pmatrix} 14\\ 5\\ 9 \end{pmatrix} \]

\[ \begin{array}{l} \square \mbox{Inversa: não há!}\\ \square \mbox{Inversa generalizada (Penrose):}\\ (\boldsymbol{X}^t.\boldsymbol{X})^+= \begin{pmatrix} 1/9 & 1/18 & 1/18\\ 1/18 & 5/18 & -2/9\\ 1/18 & -2/9 & 5/18 \end{pmatrix}\\ \square \mbox{Consistência:}\\ (\boldsymbol{X}^t.\boldsymbol{X})(\boldsymbol{X}^t.\boldsymbol{X})^{+}\boldsymbol{X}^t\boldsymbol{y}=\boldsymbol{X}^t\boldsymbol{y} \mbox{ pois:}\\ (\boldsymbol{X}^t.\boldsymbol{X})(\boldsymbol{X}^t.\boldsymbol{X})^{+}\boldsymbol{X}^t\boldsymbol{y}=\begin{pmatrix} 14\\ 5\\ 9 \end{pmatrix} \end{array} \]

\[ \begin{array}{l} \square \mbox{Melhor solução aproximada (dentre as várias soluções):} \\ \hat{\boldsymbol{\beta}}=\boldsymbol{X}^+y=\begin{pmatrix} 7/3\\ 1/6\\ 13/6 \end{pmatrix} \end{array} \]

Implementação no R

X=matrix(c(1,1,0,1,1,0,1,0,1,1,0,1),byrow=T,nrow=4)
y=matrix(c(2,3,4,5),byrow=T,nrow=4)
showEqn(X,y,paste0("beta", 1:ncol(X)))

## 1*beta1 + 1*beta2 + 0*beta3  =  2 
## 1*beta1 + 1*beta2 + 0*beta3  =  3 
## 1*beta1 + 0*beta2 + 1*beta3  =  4 
## 1*beta1 + 0*beta2 + 1*beta3  =  5

Inversa generalizada

Xmais=fractions(ginv(X))
Xmais

##      [,1] [,2] [,3] [,4]
## [1,]  1/6  1/6  1/6  1/6
## [2,]  1/3  1/3 -1/6 -1/6
## [3,] -1/6 -1/6  1/3  1/3

Consistência: não há!

fractions(X%*%Xmais%*%y)

##      [,1]
## [1,] 5/2 
## [2,] 5/2 
## [3,] 9/2 
## [4,] 9/2

all.equal(X%*%Xmais%*%y,y)

## [1] "Mean relative difference: 0.1428571"

Sistema de equações normais:

showEqn(t(X)%*%X,t(X)%*%y,paste0("beta_hat", 1:ncol(t(X)%*%X)))

## 4*beta_hat1 + 2*beta_hat2 + 2*beta_hat3  =  14 
## 2*beta_hat1 + 2*beta_hat2 + 0*beta_hat3  =   5 
## 2*beta_hat1 + 0*beta_hat2 + 2*beta_hat3  =   9

Matriz inversa generalizada

try(solve(t(X)%*%X))

## Error in solve.default(t(X) %*% X) : 
##   Lapack routine dgesv: system is exactly singular: U[3,3] = 0

aux=fractions(ginv(t(X)%*%X))
aux

##      [,1] [,2] [,3]
## [1,]  1/9 1/18 1/18
## [2,] 1/18 5/18 -2/9
## [3,] 1/18 -2/9 5/18

Consistência:

aux2=t(X)%*%X
aux3=t(X)%*%y
all.equal(aux2%*%aux%*%aux3,aux3)

## [1] TRUE

Melhor solução aproximada (dentre as várias soluções)

betas=fractions(ginv(X)%*%y)
betas

##      [,1]
## [1,]  7/3
## [2,]  1/6
## [3,] 13/6

Visualiza

pred=X%*%betas
dados=data.frame(x1=X[,1],x2=X[,2],x3=X[,3],y=y,pred=pred)
dados

##   x1 x2 x3 y pred
## 1  1  1  0 2  2.5
## 2  1  1  0 3  2.5
## 3  1  0  1 4  4.5
## 4  1  0  1 5  4.5

Observações:

Este problema é fictício, com poucas observações;
Nos problemas práticos de regressão linear, o banco de dados tem mais observações do que este exemplo, e são feitos estudos para contornar a colinearidade de variáveis independentes, excluindo ou transformando estas variáveis, evitando que a matriz \(\boldsymbol{X}^t.\boldsymbol{X}\) seja singular, e visando um modelo mais parcimonioso (que não inclui muitas variáveis sem necessidade).

5 Aplicações - Análise de Variância

Exemplo 2 (Amarante, 1992), página 82. Os dados abaixo são adaptação dos dados de Vannucchi e colaboradores (1991), referente ao estudo do comportamento dos níveis séricos de aminoácidos em pacientes com pelagra: doença causada pela deficiência de vitamina B3 (niacina), e submetidos a ingestão de vitaminas (níveis baixo, médio ou alto). Este estudo é um estudo típico de análise de variância com um fator (one way), em que a resposta do nível sérico para o indivíduo \(i\) é dada por: \[ y_i = \mu + \alpha_{j} +\epsilon_i, \]

em que \(\mu\) é uma média geral, \(\alpha_1\) e \(\alpha_2\) são os efeitos de ingestão média ou alta de vitaminas, respectivamente e \(\epsilon_i\) é o erro aleatório do indivíduo \(i\).

Input =("
obs;ingestao;valina
1;media;0.0616
2;media;0.1176
3;media;0.1254
4;baixa;0.2173
5;baixa;0.0958
6;baixa;0.1379
7;baixa;0.2072
8;media;0.1064
9;alta;0.2843
10;media;0.1764
11;alta;0.2846
12;alta;0.1751
13;alta;0.2400
14;alta;0.2290
15;baixa;0.1988
")

dados=read.table(textConnection(Input),header=TRUE,sep=";")
ordem_tratamentos=c("baixa","media","alta")
dados$ingestao=factor(dados$ingestao,levels=ordem_tratamentos)
dados = dados %>%
  arrange(ingestao)

datatable(dados,options = list(
            scrollX = TRUE,
            scrollY = "400px",
            pageLength = 5,
            paging = TRUE
          )
)

Visualização gráfica

dados %>% 
  ggplot(aes(x=ingestao,y=valina,fill=ingestao)) +
         geom_boxplot()

Visão matricial:
- para não ter que construir a matriz \(\boldsymbol{X}\), vamos extrair dos resultados diretamente da análise de variância aov no R,
- E vamos exibir o sistema de equações:

anova=aov(valina ~ ingestao , data=dados)
X=as.matrix(model.matrix(anova))
y=matrix(dados$valina,ncol=1)
showEqn(X,y,c("mu","alpha1","alpha2"))

## 1*mu + 0*alpha1 + 0*alpha2  =  0.2173 
## 1*mu + 0*alpha1 + 0*alpha2  =  0.0958 
## 1*mu + 0*alpha1 + 0*alpha2  =  0.1379 
## 1*mu + 0*alpha1 + 0*alpha2  =  0.2072 
## 1*mu + 0*alpha1 + 0*alpha2  =  0.1988 
## 1*mu + 1*alpha1 + 0*alpha2  =  0.0616 
## 1*mu + 1*alpha1 + 0*alpha2  =  0.1176 
## 1*mu + 1*alpha1 + 0*alpha2  =  0.1254 
## 1*mu + 1*alpha1 + 0*alpha2  =  0.1064 
## 1*mu + 1*alpha1 + 0*alpha2  =  0.1764 
## 1*mu + 0*alpha1 + 1*alpha2  =  0.2843 
## 1*mu + 0*alpha1 + 1*alpha2  =  0.2846 
## 1*mu + 0*alpha1 + 1*alpha2  =  0.1751 
## 1*mu + 0*alpha1 + 1*alpha2  =    0.24 
## 1*mu + 0*alpha1 + 1*alpha2  =   0.229

Exibindo o sistema de equações normais:

showEqn(t(X)%*%X,t(X)%*%y,c("mu","alpha1","alpha2"))

## 15*mu + 5*alpha1 + 5*alpha2  =  2.6574 
##  5*mu + 5*alpha1 + 0*alpha2  =  0.5874 
##  5*mu + 0*alpha1 + 5*alpha2  =   1.213

Solução do sistema \(\boldsymbol{\beta}=(\boldsymbol{X}^t.\boldsymbol{X})^{-1}\boldsymbol{X}^t\boldsymbol{y}\)
Não é necessário verificar consistência pois o sistema de equações normais já tem a propriedade de consistência!

sols=try(solve(t(X)%*%X))%*%t(X)%*%y
sols

##                   [,1]
## (Intercept)    0.17140
## ingestaomedia -0.05392
## ingestaoalta   0.07120

Análise de variância diretamente no R com o comando aov - Analysis of Variance

anova$coefficients

##   (Intercept) ingestaomedia  ingestaoalta 
##       0.17140      -0.05392       0.07120

Verificação do efeito de tratamento e comparação de médias para tirar conclusão:
- há diferença entre os níveis médio e alto conforme verificamos na comparação de médias.

summary(anova)

##             Df  Sum Sq  Mean Sq F value  Pr(>F)   
## ingestao     2 0.03939 0.019693   9.087 0.00396 **
## Residuals   12 0.02601 0.002167                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(anova)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = valina ~ ingestao, data = dados)
## 
## $ingestao
##                 diff          lwr        upr     p adj
## media-baixa -0.05392 -0.132470847 0.02463085 0.2012542
## alta-baixa   0.07120 -0.007350847 0.14975085 0.0771658
## alta-media   0.12512  0.046569153 0.20367085 0.0029910

Verificação dos pressupostos da anova:
- Distribuição normal dos resíduos: não rejeitamos a hipótese de normalidade dos resíduos (ao nível de significância \(5\%\))

qqnorm(residuals(anova))
qqline(residuals(anova))

shapiro.test(residuals(anova))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova)
## W = 0.93381, p-value = 0.3108

Verificação da homogeneidade de váriâncias
- não rejeitamos a hipótese de homogeneidade de variâncias (ao nível de significância \(5\%\))

bartlett.test(residuals(anova) ~ ingestao, data = dados)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(anova) by ingestao
## Bartlett's K-squared = 0.21337, df = 2, p-value = 0.8988

Aleatoriedade dos resíduos
- graficamente, verificamos que os resíduos se distribuem de forma aleatória.

# Plot residuals
plot(fitted(anova), residuals(anova),
     ylab="Residuals", xlab="Fitted Values", 
     main="Residuals vs Fitted Values")
abline(h=0, col="red")

6 Aplicações - Análise de Variância - Continuação

Exemplo 3 Leitura e interpretação do material (De Luna & De Olinda, 2008) páginas 112 a 119 notas de curso de modelos Lineares

Exemplo 4 (De Luna & De Olinda, 2015) Os valores observados de uma variável resposta num experimento inteiramente casualizado foram os seguintes:

knitr::include_graphics('figura1_02out.png')

1. Faça uma visualização gráfica e análise de variância usual - pela ferramenta aov do software R, considerando também análise de resíduos ao final
1. Admitindo o modelo \(\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}\), escreva o conjunto de equações relativas a cada observação (sistema com 9 equações)
1. Escreva o sistema de equações do item anterior na forma matricial
1. Escreva o sistema de equações normais
1. Determine a solução das equações normais, e confira com os resultados do item (a)
1. Calcule SQPar = \(\boldsymbol{y}^t\boldsymbol{P}\boldsymbol{y}=\boldsymbol{\hat{\beta}}^t\boldsymbol{X}^t\boldsymbol{y}\), onde \(\boldsymbol{\hat{\beta}}^t\) é o vetor linha com a solução das equações normais
1. Calcule SQRes = \(\boldsymbol{y}^t(\boldsymbol{I}-\boldsymbol{P})\boldsymbol{y}\)=\(\boldsymbol{y}^t\boldsymbol{y}-\boldsymbol{\hat{\beta}}^t\boldsymbol{X}^t\boldsymbol{y}\) e confira com os resultados do item (a)
1. Calcule SQTrat = \(\boldsymbol{y}^t(\boldsymbol{P}-\boldsymbol{P}_1)\boldsymbol{y}\) e confira com os resultados do item (a).

Resolução:

Input =("
trat;resp
1;4
1;3
1;2
2;3
2;4
2;5
3;6
3;7
3;8
")
dados1=read.table(textConnection(Input),header=TRUE,sep=";")%>%
      mutate(trat=as.factor(trat))

1. visualização gráfica: De acordo com o grafico boxplot abaixo, o tratamento 3 apresenta respostas superiores aos demais tratamentos.

dados1 %>% 
  ggplot(aes(x=trat,y=resp,fill=trat)) +
         geom_boxplot()

1. análise de variância usual: De acordo com a tabela de análise de variância, as respostas diferem significativamente de acordo com os tratamentos, ao nível de significância usual de \(5\%\) pois o valor-p da estatística F do teste é igual a 0.00659.

anova1= aov(resp ~ trat, data=dados1)
summary(anova1)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## trat         2     26      13      13 0.00659 **
## Residuals    6      6       1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1. comparação de médias: De acordo com o teste de comparações múltiplas de Tukey, os tratamentos 1 e 2 não diferem, enquanto que o tratamento 3 difere dos demais, ao nível de significância de \(5\%\).

TukeyHSD(anova1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = resp ~ trat, data = dados1)
## 
## $trat
##     diff        lwr      upr     p adj
## 2-1    1 -1.5052356 3.505236 0.4827273
## 3-1    4  1.4947644 6.505236 0.0064937
## 3-2    3  0.4947644 5.505236 0.0242291

1. análise de resíduos: De acordo com a análise de resíduos, rejeita-se a hipótese nula de que os resíduos tem distribuição normal. Pela inspeção visual, verificamos que a normalidade dos resíduos também é prejudicada pelo tamanho reduzido da amostra deste problema.

qqnorm(residuals(anova1))
qqline(residuals(anova1))

shapiro.test(residuals(anova1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(anova1)
## W = 0.82304, p-value = 0.03729

De acordo com o teste de Bartlett, não rejeita-se a hipótese de homogeneidade das variâncias, ao nível de significância de \(5\%\).

bartlett.test(residuals(anova1) ~ trat, data = dados1)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(anova1) by trat
## Bartlett's K-squared = 3.6335e-16, df = 2, p-value = 1

gráficamente, verificamos que os resíduos se distribuem de forma aleatória (nao há presença aparente de autocorrelação nos resíduos). Fica pendente verificar testes de aleatoriedade dos resíduos para confirmar esta afirmação.

plot(fitted(anova1), residuals(anova1),
     ylab="Residuals", xlab="Fitted Values", 
     main="Residuals vs Fitted Values")
abline(h=0, col="red")

1. Como cada observação é dada pela fórmula: \(y_{ij}=\beta_0+\beta_{ij}+\epsilon_{ij}\) (resposta do indivíduo \(i\), no tratamento j), \(i=1,2,\ldots, 9\) e \(j=1,2,3\). \[ \left\{ \begin{array}{l} 4=\beta_0+\beta_1+\epsilon_{11} \\ 3=\beta_0+\beta_1+\epsilon_{21}\\ 2=\beta_0+\beta_1+\epsilon_{31} \\ 3=\beta_0+\beta_2+\epsilon_{42} \\ 4=\beta_0+\beta_2+\epsilon_{52} \\ 5=\beta_0+\beta_2+\epsilon_{62} \\ 6=\beta_0+\beta_3+\epsilon_{73} \\ 7=\beta_0+\beta_3+\epsilon_{83} \\ 8=\beta_0+\beta_3+\epsilon_{93} \end{array} \right. \]

onde \(\beta_0\) é uma média geral, \(\beta_1\), \(\beta_2\) e \(\beta_3\) são os efeitos dos tratamentos 1, 2 e 3, respectivamente, e cada \(\epsilon_{ij}\) é o erro do indivíduo \(i\) no tratamento \(j\) (o indivíduo aparece primeiro no índice).

1. Na forma matricial: \[ \underbrace{ \begin{pmatrix} 4\\ 3\\ 2\\ 3\\ 4\\ 5\\ 6\\ 7\\ 8\\ \end{pmatrix}}_{\boldsymbol{y}}= \underbrace{ \begin{pmatrix} 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1 \end{pmatrix} }_{\boldsymbol{X}} \underbrace{ \begin{pmatrix} \beta_0 \\ \beta_1\\ \beta_2\\ \beta_3 \end{pmatrix} }_{\boldsymbol{\beta}} +\underbrace{ \begin{pmatrix} \epsilon_{11} \\ \epsilon_{21}\\ \epsilon_{31} \\ \epsilon_{42} \\ \epsilon_{52} \\ \epsilon_{62} \\ \epsilon_{73} \\ \epsilon_{83} \\ \epsilon_{93} \end{pmatrix} }_{\boldsymbol{\epsilon}} \]
Alternativamente, é possível escrever no R o conjunto de equações relativas a cada observação: extrai-se a matriz \(\boldsymbol{X}\) no objeto aov.
Considerações importantes: nesta abordagem, o tratamento 1 é considerado referência, e os outros tratamentos serão comparados com o primeiro, sendo assim a matriz \(\boldsymbol{X}\) tem 4 colunas ao invés de 3 colunas, simplificando a visão matricial.

anova1=aov(resp ~ trat , data=dados1)
X=as.matrix(model.matrix(anova1))
X

##   (Intercept) trat2 trat3
## 1           1     0     0
## 2           1     0     0
## 3           1     0     0
## 4           1     1     0
## 5           1     1     0
## 6           1     1     0
## 7           1     0     1
## 8           1     0     1
## 9           1     0     1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$trat
## [1] "contr.treatment"

y=matrix(dados1$resp,ncol=1)
showEqn(X,y,c("beta1","beta2","beta3"))

## 1*beta1 + 0*beta2 + 0*beta3  =  4 
## 1*beta1 + 0*beta2 + 0*beta3  =  3 
## 1*beta1 + 0*beta2 + 0*beta3  =  2 
## 1*beta1 + 1*beta2 + 0*beta3  =  3 
## 1*beta1 + 1*beta2 + 0*beta3  =  4 
## 1*beta1 + 1*beta2 + 0*beta3  =  5 
## 1*beta1 + 0*beta2 + 1*beta3  =  6 
## 1*beta1 + 0*beta2 + 1*beta3  =  7 
## 1*beta1 + 0*beta2 + 1*beta3  =  8

1. Considerando a visão matricial simplificada então o sistema de equações normais \((\boldsymbol{X}^t.\boldsymbol{X}).\boldsymbol{\beta}=\boldsymbol{X}^t\boldsymbol{y}\) pode ser escrito como:

showEqn(t(X)%*%X,t(X)%*%y,c("beta1","beta2","beta3"))

## 9*beta1 + 3*beta2 + 3*beta3  =  42 
## 3*beta1 + 3*beta2 + 0*beta3  =  12 
## 3*beta1 + 0*beta2 + 3*beta3  =  21

Escrevendo no módulo latex: \[ (\boldsymbol{X}^t.\boldsymbol{X}).\boldsymbol{\beta}=\boldsymbol{X}^t\boldsymbol{y} \Rightarrow \left\{ \begin{array}{l} 9\beta_1 + 3\beta_2 + 3\beta_3 = 42 \\ 3\beta_1 + 3\beta_2 = 12 \\ 3\beta_1 + 3\beta_3 = 21 \end{array} \right. \]

onde \(\beta_1\) é o intercepto do modelo, que é igual a média do tratamento 1, \(\beta_2\) é o incremento do tratamento 2 com relação ao tratamento 1 e \(\beta_3\) é o incremento do tratamento 3 com relação ao tratamento 1.

1. Considerando a visão matricial simplificada a solução das equações normais \(\boldsymbol{\hat{\beta}}=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t\boldsymbol{y}\) resulta em \(\hat{\beta_1}=3\) (intercepto), \(\hat{\beta_2}=1\) (incremento do tratamento 2 com relação ao 1) e \(\hat{\beta_3}=4\) (incremento do tratamento 3 com relação ao 1).

beta_hat=solve(t(X)%*%X)%*%t(X)%*%y
beta_hat

##             [,1]
## (Intercept)    3
## trat2          1
## trat3          4

Que confere com os coeficientes da Anova:

anova1$coefficients

## (Intercept)       trat2       trat3 
##           3           1           4

1. SQPar = \(\boldsymbol{Y}^t.\boldsymbol{P}.\boldsymbol{y}\), onde \(\boldsymbol{P}=\boldsymbol{X}(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t \Rightarrow\) SQPar = 222.

P=X%*%solve(t(X)%*%X)%*%t(X)
SQPar=t(y)%*%P%*%y
SQPar

##      [,1]
## [1,]  222

ou equivalentemente: SQPar = \(\boldsymbol{\hat{\beta}}^t.\boldsymbol{X}^t\boldsymbol{y}\), onde \(\boldsymbol{\hat{\beta}}^t\) é vetor linha com a solução das equações normais.

t(beta_hat)%*%t(X)%*%y

##      [,1]
## [1,]  222

1. SQRes=\(\boldsymbol{y}^t(\boldsymbol{I}-\boldsymbol{P})\boldsymbol{y}\), onde \(\boldsymbol{I}\) é uma matriz identidade com a mesma dimensão de \(\boldsymbol{P} \Rightarrow\) SQRes=6

SQRes=t(y)%*%(diag(nrow(P))-P)%*%y
SQRes

##      [,1]
## [1,]    6

ou equivalentemente: SQRes = \(\boldsymbol{y}^t\boldsymbol{y}-\boldsymbol{\hat{\beta}}^t\boldsymbol{X}^t\boldsymbol{y}\)

t(y)%*%y-t(beta_hat)%*%t(X)%*%y

##      [,1]
## [1,]    6

que confere com a Soma de quadrados de resíduos da ANOVA:

sum(anova1$residuals^2)

## [1] 6

1. SQTrat = \(\boldsymbol{y}^t(\boldsymbol{P}-\boldsymbol{P}_1)\boldsymbol{y}\), onde \(\boldsymbol{X}_1\) é formada pela coluna de uns (primeira coluna da matriz \(\boldsymbol{X}\)), e \(\boldsymbol{P}_1=\boldsymbol{X}_1(\boldsymbol{X}_1^t\boldsymbol{X}_1)^{-1}\boldsymbol{X}_1^t \Rightarrow\) SQTrat=26.

X1=X[,1]
P1=X1%*%solve(t(X1)%*%X1)%*%t(X1)
SQTrat=t(y)%*%(P-P1)%*%y
SQTrat

##      [,1]
## [1,]   26

que confere com a Soma de quadrados de tratamentos da ANOVA.

summary(anova1)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## trat         2     26      13      13 0.00659 **
## Residuals    6      6       1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7 Referências

[1] Amarante, A. R.. Um curso em modelos lineares. 1992. 329 f. Dissertação (mestrado) - Universidade Estadual de Campinas, Instituto de Matematica, Estatistica e Ciencia da Computação, Campinas, SP. Disponível em: https://hdl.handle.net/20.500.12733/1576601. Acesso em: 24 set. 2023.

[2] De Luna, J. G., & De Olinda, R. A. (2015). Introdução a Modelos Lineares. Livraria da Física, 1ª edição, 164 p.

[3] De Luna, J. G., & Esteves, D. M. Modelos lineares (notas de curso). Texto elaborado pelos professores, utilizado na disciplina de Modelos lineares do curso de Bacharelado em Estatística da UEPB, Campina Grande - PB, 2008.

[4] Fieller, N. (2016). Basics of Matrix Algebra for Statistics with R. (Chapman & Hall/CRC The R Series), 1ª edição, 244 p..

[5] Reinprecht, J.. (2022). A inversa generalizada de uma matriz e aplicações. Dissertação de mestrado da Universidade Estadual Júlio de Mesquita Filho, Rio Claro, SP.

[6] Vannucchi, H., Moreno, F. S., Amarante, A. R., Oliveira, J. E. D., Marchini, J. S. (1991): Plasma Amino Acid Pattems in Alcoholic Pellagra Patients, Alcohol and Alcoholism, 26,4, 431-436.

Aula de Modelos Lineares - Unidade III