Análise de Variância (ANOVA)

É uma técnica estatística utilizada para comparar as médias de dois ou mais grupos, ou seja, se as médias de diferentes grupos são realmente diferentes ou se as diferenças observadas podem ser atribuídas ao acaso. É amplamente utilizada em diversas áreas, como ciências sociais, biologia, economia, entre outras, para testar hipóteses sobre a variação entre grupos de dados. Existem diferentes tipos de ANOVA, como ANOVA de um fator (um caminho), ANOVA de dois fatores (dois caminhos) e ANOVA de medidas repetidas. Neste material, revisaremos a ANOVA de um fator e a ANOVA aplica em modelos lineares.

1 Notação

Considere as seguintes notações: \[\begin{eqnarray} Y_{ij} &:& \mbox{ variável resposta para a $j$-ésima unidade do $i$-ésimo grupo}\\ \mu &:& \mbox{ média global (média para toda a população)}\\ \mu_i &:& \mbox{ média do $i$-ésimo grupo}\\ n_i &:& \mbox{ total de unidades no $i$-ésimo grupo} \end{eqnarray}\]

2 Homocedasticidade

Para testarmos se as variâncias são iguais, podemos usar o teste de Bartlett ou o teste de Levene.

O teste de Bartlett é sensível em relação a hipótese de normalidade dos dados. Se rejeitarmos a hipótese de normalidade, é melhor utilizarmos o teste proposto por Levene. Porém, se a hipótese de normalidade não for violada, o teste proposto por Bartlett tem um comportamento melhor que o teste proposto por Levene. O teste F para variâncias somente pode ser usado quando há 2 populações, enquanto os outros 2 podem ser usados para \(K\) populações.

Neste material, a normalidade é um pressuposto comum a todos os testes de média. Sendo assim, usaremos o teste de Bartlett, o teste de Levene e o teste F (quando possível) para compararmos os resultados.

3 Teste t para comparação de médias

Suponha que temos 2 variáveis com distribuição normal e queremos investigar se as médias destas distribuições são iguais. Ou seja, \[\begin{eqnarray} Y_{1j} &\sim& N(\mu_1, \sigma_1^2)\\ Y_{2j} &\sim& N(\mu_2, \sigma_2^2). \end{eqnarray}\] e \[\begin{eqnarray} H_0: && \mu_1 = \mu_2 = \mu\\ H_1: && \mu_1 \neq \mu_2. \end{eqnarray}\]

3.1 Variâncias Iguais (\(\sigma_1^2 = \sigma_2^2\)) e amostras independentes

Pressupostos:

  1. Normalidade.
  2. Independência (Os grupos são independentes e os indivíduos dentro dos grupos também são independentes).
  3. Homocedasticidade (variâncias iguais).

Considerando \(H_0\) verdadeiro, temos que \[T = \frac{\bar{Y}_1 - \bar{Y}_2}{\sqrt{s_p^2\left(\frac{1}{n_1}+\frac{1}{n_2}\right)}} \sim t_{(n_1+n_2-2)},\] sendo \(n_1\) o total de observações do primeiro grupo, \(n_2\) o total de observações do segundo grupo, \(\bar{Y}_i = \frac{\sum_{j=1}^{n_i}{Y_{ij}}}{n_i}\), \(s_p^2 = \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}\), \(s_i^2 = \frac{\sum_{j=1}^{n_i}{(Y_{ij} - \bar{Y_i})^2}}{n_i-1}\), \(i=1,2\).

3.2 Variâncias Diferentes (\(\sigma_1^2 \neq \sigma_2^2\)) e amostras independentes

Pressupostos:

  1. Normalidade.
  2. Independência (Os grupos são independentes e os indivíduos dentro dos grupos também são independentes).
  3. Heterocedasticidade (variâncias diferentes).

Considerando \(H_0\) verdadeiro, temos que \[T = \frac{\bar{Y}_1 - \bar{Y}_2}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}} \sim t_{\nu},\] sendo \(\nu = \frac{(A+B)^2}{\frac{A^2}{n_1-1} + \frac{B^2}{n_2-1}}\), \(A = \frac{s_1^2}{n_1}\) e \(B = \frac{s_2^2}{n_2}\).

3.3 Amostras Dependentes/Pareadas

Pressupostos:

  1. Normalidade.
  2. Dependência (Os grupos são dependentes e os indivíduos dentro dos grupos são independentes).
  3. Heterocedasticidade (variâncias diferentes).

Considerando \(H_0\) verdadeiro, temos que \[T = \frac{\bar{Y}_1 - \bar{Y}_2}{\sqrt{\frac{s_D^2}{n}}} \sim t_{(n-1)},\] sendo \(n = n_1 = n_2\) o total de observações pareadas, \(\bar{Y}_i = \frac{\sum_{j=1}^{n}{Y_{ij}}}{n}\), \(s_D^2 = \frac{\sum_{j=1}^{n}{(D_j - \bar{D})^2}}{n-1}\), \(D_j = Y_{1j} - Y_{2j}\) e \(\bar{D} = \frac{\sum_{j=1}^{n}{D_{j}}}{n}\).

4 Análise de variância de 1 fator (ANOVA)

Suponha que a variável resposta possa ser avaliada em \(K\) grupos e que para cada grupo tenha a seguinte distribuição: \[\begin{eqnarray} Y_{ij} &\stackrel{ind}{\sim}& N(\mu_i, \sigma^2), \quad i=1, \ldots, n_j, \quad j=1, \ldots, K, \end{eqnarray}\] e \(Y_{ij}\) independente de \(Y_{sj}\), para todo \(i \neq s\).

4.1 Supondo 2 grupos

Supondo que há 2 grupos, temos que \(Y_{1j} \sim N(\mu_1, \sigma^2)\) e \(Y_{2j} \sim N(\mu_2, \sigma^2)\) queremos testar se o fator tem influência na variável resposta. Ou seja,

\[\begin{eqnarray} H_0: && \mu_1 = \mu_2 = \mu\\ H_1: && \mu_1 \neq \mu_2. \end{eqnarray}\]

Pressupostos realizados: Para cada grupo, supõe-se uma distribuição normal da seguinte forma: \[Y_{ij} = \mu_i + e_{ij}, \quad e_{ij} \stackrel{iid}{\sim} N(0,\sigma^2).\]

Note que estamos supondo:

  1. Normalidade.
  2. Independência (Os grupos são independentes e os indivíduos dentro dos grupos também são independentes).
  3. Homocedasticidade (variâncias iguais).

Realizar o teste descrito acima, significa comparar os 2 modelos: \(Y_{1j} = \mu_1 + e_{1j}\) e \(Y_{2j} = \mu_2 + e_{2j}\).

4.1.1 Supondo \(H_1\) verdadeiro

Supondo \(H_1\) verdadeiro, temos que \(Y_{ij} = \mu_i + e_{ij}\).

Logo, o estimador de \(\mu_i\) pode ser obtido pelo método dos mínimos quadrados, ou seja, minimizando a soma dos erros quadráticos

\[\sum_{i=1}^2{\sum_{j=1}^{n_i}{e_{ij}^2}} = \sum_{i=1}^2{\sum_{j=1}^{n_i}{ (Y_{ij} - \mu_i)^2}}.\]

Portanto, temos que o estimador de \(\mu_i\) é \[\begin{eqnarray} \hat{\mu}_i &=& \bar{Y}_i = \frac{ \sum_{j=1}^{n_i}{Y_{ij}} }{n_i} \end{eqnarray}\]

Seja \(r_{ij}\) o resíduo do \(i\)-ésimo grupo da \(j\)-ésima unidade. Seja SQD a soma dos quadrados dentro dos grupos (também chamado de SQR), ou seja, a soma dos quadrados dos resíduos, quando \(H_1\) é verdadeiro. Logo, \[\begin{eqnarray} SQD = SQR = \sum_{i=1}^{2}{\sum_{j=1}^{n_i}{r_{ij}^2}} = \sum_{i=1}^2{\sum_{j=1}^{n_i}{(Y_{ij} - \bar{Y}_i)^2}} \end{eqnarray}\]

O estimador não viesado de \(\sigma^2\) é \[\begin{eqnarray} \hat{\sigma}^2 &=& S^2 = \frac{(n_1-1)S_1^2 + (n_2-1)S_2^2}{n_1+n_2-2} \quad \mbox{ sendo } S_i^2 = \frac{\sum_{j=1}^{n_i}{(Y_{ij}-\bar{Y}_i)^2}}{n_i-1} \end{eqnarray}\]

Portanto, \[\begin{eqnarray} SQD = SQR &=& (n_1-1) S_1^2 + (n_2-1) S_2^2 = (n-2) S^2 \end{eqnarray}\]

Note que \(\bar{Y}_i \sim N(\mu_i, \sigma^2/n_i)\). E, portanto, \(\frac{(\bar{Y}_i-\mu_i)\sqrt{n_i}}{S} \sim t_{(n-2)}\). Então, um intervalo de confiança para \(\mu_i\) bilateral com \(\gamma\)% de confiança é \[\left [ \bar{Y}_i - t_{1-\alpha/2} \frac{S}{\sqrt{n_i}} \;,\; \bar{Y}_i + t_{1-\alpha/2} \frac{S}{\sqrt{n_i}}\right]\] sendo \(\alpha = 1-\gamma\) o nível de significância e \(t_{1-\alpha/2}\) o quantil \(1-\alpha/2\) da distribuição t-Student com \(n-2\) graus de liberdade. E, portanto, \[\bar{Y}_1 - \bar{Y}_2 \sim N\left(\mu_1 - \mu_2, \frac{\sigma^2}{n_1} + \frac{\sigma^2}{n_2}\right).\] Sendo assim, \[\frac{\left( \bar{Y}_1 - \bar{Y}_2 \right) - \left( \mu_1 - \mu_2 \right)}{S\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}} \sim t_{(n-2)}.\] Então, um intervalo de confiança para \(\mu_1 - \mu_2\) bilateral com \(\gamma\)% de confiança é \[\left [ \bar{Y}_1 - \bar{Y}_2 - t_{1-\alpha/2} S\sqrt{\frac{1}{n_1} + \frac{1}{n_2}} \;,\; \bar{Y}_1 - \bar{Y}_2 + t_{1-\alpha/2} S\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}\right]\] Dessa forma, rejeitamos \(H_0\), com \(\alpha\)% de significância, quando o zero não estiver incluído no intervalo acima.

4.1.2 Supondo \(H_0\) verdadeiro

Supondo \(H_0\) verdadeiro, temos que \(Y_{ij} = \mu + e_{ij}\).

Logo, o estimador de \(\mu\) pode ser obtido pelo método dos mínimos quadrados, ou seja, minimizando a soma dos erros quadráticos

\[\sum_{i=1}^2{\sum_{j=1}^{n_i}{e_{ij}^2}} = \sum_{i=1}^2{\sum_{j=1}^{n_i}{ (Y_{ij} - \mu)^2}}.\]

Portanto, temos que o estimador de \(\mu\) é \[\begin{eqnarray} \hat{\mu} &=& \bar{Y} = \frac{\sum_{i=1}^2{\sum_{j=1}^{n_i}{Y_{ij}}}}{n} \end{eqnarray}\]

Seja \(r_{ij}\) o resíduo do \(i\)-ésimo grupo da \(j\)-ésima unidade. Seja SQT a soma dos quadrados totais, ou seja, a soma dos quadrados dos resíduos, quando \(H_0\) é verdadeiro. Logo, \[\begin{eqnarray} SQT = \sum_{i=1}^{2}{\sum_{j=1}^{n_i}{r_{ij}^2}} = \sum_{i=1}^2{\sum_{j=1}^{n_i}{(Y_{ij} - \bar{Y})^2}} \end{eqnarray}\]

O estimador não viesado de \(\sigma^2\) é \[\begin{eqnarray} \hat{\sigma}^2 &=& S_T^2 = \frac{1}{n-1}\sum_{i=1}^2{\sum_{j=1}^{n_i}{(Y_{ij}-\bar{Y})^2}} \end{eqnarray}\]

Portanto, \[\begin{eqnarray} SQT = (n-1) S_T^2 \end{eqnarray}\]

4.1.3 Comparando as hipóteses

Seja \(SQE\) a soma dos quadrados entre os grupos. Logo, temos que \[\begin{eqnarray} SQT = SQD + SQE \end{eqnarray}\]

Sendo assim, temos que \[\begin{eqnarray} SQE = \sum_{i=1}^2{n_i(\bar{Y}_i - \bar{Y})^2}. \end{eqnarray}\] A soma acima representa a variabilidade entre as médias amostrais.

Seja \(R^2\) o , ou seja, a proporção da variação explicada pelo modelo. Então \[\begin{eqnarray} R^2 = \frac{SQE}{SQT}. \end{eqnarray}\]

Sob \(H_0\) verdadeiro, temos que \[\begin{eqnarray} T &=& \frac{\bar{Y}_1 - \bar{Y}_2}{S\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}} \sim t_{(n_1 + n_2 - 2)}\\ T^2 &\sim& F_{(1, n_1+n_2-2)} \end{eqnarray}\] sendo \(t_{(\nu)}\) a distribuição t-Student com \(\nu\) graus de liberdade e \(F_{(\nu_1,\nu_2)}\) a distribuição F-Snedecor com \(\nu_1\) e \(\nu_2\) graus de liberdade.

Como \(\bar{Y} = \frac{n_1 \bar{Y}_1 + n_2 \bar{Y}_2}{n_1+n_2}\), temos que \[\begin{eqnarray} SQE &=& \frac{n_1n_2}{n_1+n_2}(\bar{Y}_1-\bar{Y}_2)^2 = \frac{(\bar{Y}_1-\bar{Y}_2)^2}{1/n_1 + 1/n_2} \end{eqnarray}\] Portanto \[\begin{eqnarray} T^2 &=& \frac{(\bar{Y}_1 - \bar{Y}_2)^2}{S^2(1/n_1 + 1/n_2)} = \frac{SQE}{S^2} = F \end{eqnarray}\] Rejeitamos \(H_0\) se \(F>c\) sendo \(c\) o valor crítico da distribuição F-Snedecor com \((1, n-2)\) graus de liberdade cuja probabilidade de uma variável aleatória com esta distribuição assumir um valor acima de \(c\) é \(\alpha\).

Essas informações são agrupadas numa tabela conhecida por ANOVA:

Fontes de Graus de Soma dos Média da soma F
variação liberdade quadrados dos quadrados
Entre 1 SQE MQE = \(\frac{SQE}{1}\) F = \(\frac{MQE}{MQD}\)
Dentro n-2 SQD = SQR MQD = \(\frac{SQD}{n-2}\)
Total n-1 SQT MQT = \(\frac{SQT}{n-1}\)

4.2 Supondo mais de 2 grupos

Seja \(K\) o total de grupos ou fatores. Então \[\begin{eqnarray} Y_{ij} = \mu_i + e_{ij}, \quad i=1, \ldots,K, \quad j=1, \ldots, n_i, \quad n=\sum_{i=1}^{K}{n_i} \end{eqnarray}\]

Facilmente estendemos os resultados encontrados quando temos 2 subpopulações para mais subpopulações. Desta forma, temos que realizar o teste descrito acima, significa comparar \(K\) modelos: \(Y_{ij} = \mu_i + e_{ij}\), \(i=1, \ldots, K\). Sendo assim, temos que

\[\begin{eqnarray} SQE &=& \sum_{i=1}^{K}{n_i(\bar{Y}_i - \bar{Y})^2}\\ SQD &=& \sum_{i=1}^{K}{\sum_{j=1}^{n_i}{(Y_{ij}-\bar{Y}_i)^2}} = \sum_{i=1}^K{(n_i-1)S_i^2} = S^2 (n-K)\\ SQT &=& \sum_{i=1}^{K}{\sum_{j=1}^{n_i}{(Y_{ij}-\bar{Y})^2}} \\ F &=& \frac{MQE}{MQD} = \frac{SQE}{K-1}\frac{n-K}{SQD} = \frac{SQE}{(K-1)S^2} \end{eqnarray}\] sendo \(\bar{Y}_i\) a média amostral do \(i\)-ésimo grupo, \(S_i^2\) a variância amostral do \(i\)-ésimo grupo, \(S^2\) um estimador para a variância populacional quando \(H_1\) é verdadeiro e \(\bar{Y}\) a média amostral de toda a amostra.

Além disso, temos que a tabela da ANOVA é

Fontes de Graus de Soma dos Média da soma F
variação liberdade quadrados dos quadrados
Entre K-1 SQE MQE = \(\frac{SQE}{K-1}\) F = \(\frac{MQE}{MQD}\)
Dentro n-K SQD = SQR MQD = \(\frac{SQD}{n-K}\)
Total n-1 SQT MQT = \(\frac{SQT}{n-1}\)

Rejeitamos \(H_0\) se \(F>c\) sendo \(c\) o valor crítico da distribuição F-Snedecor com \((K-1, n-K)\) graus de liberdade cuja probabilidade de uma variável aleatória com esta distribuição assumir um valor acima de \(c\) é \(\alpha\).

5 Exemplo com 1 fator (Estado) e 2 grupos (RJ e SP)

Exemplo: Suponha que o interesse esteja em avaliar a venda de um determinado produto em 2 estados diferentes (RJ e SP). Deseja-se investigar se o fator Estado influencia na venda. Dessa forma, a população é dividida em 2 grupos, sendo o primeiro composto pelas vendas do RJ e o segundo composto pelas vendas de SP. A variável resposta \(Y_{ij}\) corresponde a \(j\)-ésima venda do \(i\)-ésimo Estado, \(\mu_1\) a venda média do RJ e \(\mu_2\) a venda média de SP. 

#Exemplo de Analise de Variancia
#comparando o lucro recebido pela venda de 1 produto em regioes diferentes
y1 = c(206.32,207.94,206.19,204.45,209.65,203.81,206.75,205.68,204.49,210.86)
y2 = c(217.08,221.43,218.04,224.13,211.82,213.9,221.28,229.43,213.54,214.51)

par(mfrow=c(1,3))
plot(y1,y2, bty="n", lwd=2)
plot(y1-y2, bty="n", lwd=2, xlab="indice")
hist(c(y1, y2), freq=F, ylab="densidade", main="Histograma")

#empilhando os dados
Estados = c("RJ","SP")
fator   = as.factor( rep(Estados,c(length(y1), length(y2)))  )
#fator = as.factor( rep(c(1,2),c(length(y1), length(y2)))  )
z     = c(y1,y2)
data.frame(fator, z)
##    fator      z
## 1     RJ 206.32
## 2     RJ 207.94
## 3     RJ 206.19
## 4     RJ 204.45
## 5     RJ 209.65
## 6     RJ 203.81
## 7     RJ 206.75
## 8     RJ 205.68
## 9     RJ 204.49
## 10    RJ 210.86
## 11    SP 217.08
## 12    SP 221.43
## 13    SP 218.04
## 14    SP 224.13
## 15    SP 211.82
## 16    SP 213.90
## 17    SP 221.28
## 18    SP 229.43
## 19    SP 213.54
## 20    SP 214.51
#calculando medidas para cada grupo
K       = 2
ybarra  = NULL
ni      = NULL
s2i     = NULL
for(i in 1:K)
{
  ybarra[i]   = mean(z[fator==Estados[i]])
  ni[i]       = length(z[fator==Estados[i]])
  s2i[i]      = var(z[fator==Estados[i]])
}

5.1 Ao nível de 5%

#definindo nivel de significancia
alpha = 0.05

#testando independencia
teste = chisq.test(y1, y2)
teste
## 
##  Pearson's Chi-squared test
## 
## data:  y1 and y2
## X-squared = 90, df = 81, p-value = 0.2313
if(teste$p.value<=alpha){print("rejeito independencia")}else{print("nao rejeito independencia")}
## [1] "nao rejeito independencia"
#testando normalidade pelo Shapiro
s     = 0
K     = 2 #total de modelos diferentes
for(i in 1:K)
{   
  teste = shapiro.test(z[fator==Estados[i]]) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
 #testando normalidade pelo Kolmogorov-Smirnov
 s=0
 for(i in 1:K)
 {  
   teste = ks.test(z[fator==Estados[i]],"pnorm",mean(z[fator==Estados[i]]),sd(z[fator==Estados[i]])) 
   if(teste$p.value<=alpha){s=s+1}  
 }
 if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
 #testando homocedasticidade
 teste = bartlett.test(z,fator)
 if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "problema com a homocedasticidade"
 #testando homocedasticidade
 library(lawstat)
 teste = levene.test(z, fator, location = c("median"))
 if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "problema com a homocedasticidade"
 teste = var.test(y1,y2,conf.level=1-alpha)
 if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "problema com a homocedasticidade"
 #Caso considerassemos as variancias diferentes
 A      = s2i[1]/ni[1]
 B      = s2i[2]/ni[2]
 EstatT = (ybarra[1] - ybarra[2])/sqrt(A+B)
 
 nu = (A+B)^2/(A^2/(ni[1]-1) + B^2/(ni[2]-1))
 
 teste = t.test(y1, y2, paired = FALSE, var.equal = FALSE, conf.level=1-alpha)
 c(nu, teste$parameter[1])
##                df 
## 11.99181 11.99181
 c(EstatT, teste$statistic)
##                   t 
## -6.280445 -6.280445
 c(2*pt(EstatT, nu), teste$p.value)
## [1] 4.077261e-05 4.077261e-05
 if(teste$p.value<=alpha){print("Rejeito a igualdade de medias")}else{print("Nao rejeito a igualdade de medias")}
## [1] "Rejeito a igualdade de medias"

Como as variâncias são diferentes, não posso realizar o teste da ANOVA.

5.2 Ao nível de significância de 1%

#definindo nivel de significancia
alpha = 0.01

#testando independencia
teste = chisq.test(y1, y2)
teste
## 
##  Pearson's Chi-squared test
## 
## data:  y1 and y2
## X-squared = 90, df = 81, p-value = 0.2313
if(teste$p.value<=alpha){print("rejeito independencia")} else{print("nao rejeito independencia")}
## [1] "nao rejeito independencia"
#testando normalidade pelo Shapiro
s     = 0
K     = 2 #total de modelos diferentes
for(i in 1:K)
{   
  teste = shapiro.test(z[fator==Estados[i]]) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando normalidade pelo Kolmogorov-Smirnov
s=0
for(i in 1:K)
{   
  teste = ks.test(z[fator==Estados[i]],"pnorm",mean(z[fator==Estados[i]]),sd(z[fator==Estados[i]])) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando homocedasticidade
teste = bartlett.test(z,fator)
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
#testando homocedasticidade
library(lawstat)
teste = levene.test(z, fator, location = c("median"))
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
teste = var.test(y1,y2,conf.level=1-alpha)
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
#Teste t para amostras independentes com variancias iguais
sp2 = ((ni[1]-1)*s2i[1] + (ni[2]-1)*s2i[2])/(ni[1]+ni[2]-2)
EstatT = (ybarra[1] - ybarra[2])/sqrt(sp2*(1/ni[1]+1/ni[2]))
EstatT
## [1] -6.280445
testeT = t.test(y1, y2, paired = FALSE, var.equal = TRUE)
c(ni[1]+ni[2]-2, testeT$parameter[1])
##    df 
## 18 18
c(EstatT, testeT$statistic)
##                   t 
## -6.280445 -6.280445
c(2*pt(EstatT, ni[1]+ni[2]-2), testeT$p.value)
## [1] 6.378389e-06 6.378389e-06
if(testeT$p.value<=alpha){print("Rejeito a igualdade de medias")}else{print("Nao rejeito a igualdade de medias")}
## [1] "Rejeito a igualdade de medias"
#ANOVA (considerando variancias iguais)

tabela = anova(lm(z~fator)) #forma 1
tabela
## Analysis of Variance Table
## 
## Response: z
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## fator      1 708.29  708.29  39.444 6.378e-06 ***
## Residuals 18 323.22   17.96                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res.anova=aov(z~fator)
res.anova
## Call:
##    aov(formula = z ~ fator)
## 
## Terms:
##                    fator Residuals
## Sum of Squares  708.2880  323.2225
## Deg. of Freedom        1        18
## 
## Residual standard error: 4.237547
## Estimated effects may be unbalanced
summary(res.anova) #forma 2
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## fator        1  708.3   708.3   39.44 6.38e-06 ***
## Residuals   18  323.2    18.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Calculando os dados da tabela da ANOVA sem usar a 
#funcao do R

n       = sum(ni)
yt      = mean(z)
SQT     = sum((z-yt)^2)
SQE     = sum(ni*(ybarra-yt)^2)
yb2     = rep(ybarra,ni)
SQD     = sum((z-yb2)^2)
c(SQT, SQD + SQE)
## [1] 1031.511 1031.510
MQE     = SQE/(K-1)
MQD     = SQD/(n-K)
MQT     = SQT/(n-1) #ST^2
s2      = ((ni[1]-1)*s2i[1] + (ni[2]-1)*s2i[2])/(ni[1]+ni[2]-K)
c(SQD, (n-K)*s2)
## [1] 323.2225 323.2225
EstatF      = MQE/MQD
c(SQE/s2, EstatF)
## [1] 39.44399 39.44399
rbind(c(K-1, round(SQE,2), round(MQE,2), round(EstatF,3), 1-pf(EstatF, K-1, n-K)), c(n-K, round(SQD,2), round(MQD,2), NA, NA))
##      [,1]   [,2]   [,3]   [,4]         [,5]
## [1,]    1 708.29 708.29 39.444 6.378389e-06
## [2,]   18 323.22  17.96     NA           NA
tabela 
## Analysis of Variance Table
## 
## Response: z
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## fator      1 708.29  708.29  39.444 6.378e-06 ***
## Residuals 18 323.22   17.96                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
if(tabela$`Pr(>F)`[1] <= alpha){print("Rejeito a hipótese de igualdade de médias" )}else{print("Não rejeito a hipótese de igualdade de médias")}
## [1] "Rejeito a hipótese de igualdade de médias"
#comparando a estatistica de teste da ANOVA com a do teste T
c(testeT$statistic^2, tabela$`F value`[1])
##        t          
## 39.44399 39.44399
#comparando o p-valor da ANOVA com o do teste T
c(testeT$p.value, tabela$`Pr(>F)`[1])
## [1] 6.378389e-06 6.378389e-06

6 Exemplo com 1 fator (Estado) e 3 grupos (RJ, SP e MG)

Exemplo: Suponha que o interesse esteja em avaliar a venda de um determinado produto em 3 estados diferentes (RJ, SP e MG). Deseja-se investigar se o fator Estado influencia na venda. Dessa forma, a população é dividida em 3 grupos, sendo o primeiro composto pelas vendas do RJ, o segundo composto pelas vendas de SP e o terceiro composto pelas vendas de MG. A variável resposta \(Y_{ij}\) corresponde a \(j\)-ésima venda do \(i\)-ésimo Estado, \(\mu_1\) a venda média do RJ, \(\mu_2\) a venda média de SP e \(\mu_3\) a venda média de MG. 

y1 = c(209.9, 209.5, 209.6, 208.8, 209.5, 210.0, 211.7, 211.2, 210.3, 209.7)
y2 = c(250.4 , 249.6 , 247.9 , 249.5 , 250.0 , 250.0 , 249.5 , 249.3 , 251.8, 249.8 , 249.2, 249.0 , 250.4 , 250.6 , 250.4)
y3 = c(298.3 , 299.6 , 299.4 , 298.8 , 299.3 , 300.9 , 299.4 , 300.1 , 299.7, 301.4 , 301.4, 299.9 , 299.9 , 299.9, 299.8, 300.2 , 299.9 , 299.8 , 300.2, 301.7)
y = c(y1, y2, y3)
#y
fator = factor(rep(1:3,c(10,15,20)))

hist(y, freq=F, ylab="densidade", main="Histograma")

#empilhando os dados
z     = c(y1,y2, y3)
#cbind(fator, z)

#testando normalidade pelo Shapiro
alpha = 0.05
s     = 0
K     = 3 #total de modelos diferentes
for(i in 1:K)
{   
  teste = shapiro.test(z[fator==i]) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando normalidade pelo Kolmogorov-Smirnov
s=0
for(i in 1:K)
{   
  teste = ks.test(z[fator==i],"pnorm",mean(z[fator==i]),sd(z[fator==i])) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando homocedasticidade
teste = bartlett.test(z,fator)
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
#testando homocedasticidade
library(lawstat)
teste = levene.test(z, fator, location = c("median"))
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
tabela = anova(lm(z~fator))

res.anova=aov(z~fator)
res.anova
## Call:
##    aov(formula = z ~ fator)
## 
## Terms:
##                    fator Residuals
## Sum of Squares  58018.96     31.10
## Deg. of Freedom        2        42
## 
## Residual standard error: 0.8604724
## Estimated effects may be unbalanced
summary(res.anova)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## fator        2  58019   29009   39180 <2e-16 ***
## Residuals   42     31       1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Calculando os dados da tabela da ANOVA sem usar a 
#funcao do R

ybarra  = NULL
ni      = NULL
s2i     = NULL
for(i in 1:K)
{
  ybarra[i]   = mean(z[fator==i])
  ni[i]       = length(z[fator==i])
  s2i[i]      = var(z[fator==i])
}
n       = sum(ni)
yt      = mean(z)
SQT     = sum((z-yt)^2)
SQE     = sum(ni*(ybarra-yt)^2)
yb2     = rep(ybarra,ni)
SQD     = sum((z-yb2)^2)
c(SQT, SQD + SQE)
## [1] 58050.05 58050.05
MQE     = SQE/(K-1)
MQD     = SQD/(n-K)
MQT     = SQT/(n-1) #ST^2
s2      = ((ni[1]-1)*s2i[1] + (ni[2]-1)*s2i[2] + (ni[3]-1)*s2i[3])/(ni[1]+ni[2]+ni[3]-K) 
c(SQD, (n-K)*s2)
## [1] 31.09733 31.09733
EstatF      = MQE/MQD
c(SQE/(s2*(K-1)), EstatF)
## [1] 39180.15 39180.15
rbind(c(K-1, round(SQE,2), round(MQE,2), round(EstatF,3), 1-pf(EstatF, K-1, n-K)), c(n-K, round(SQD,2), round(MQD,2), NA, NA))
##      [,1]     [,2]     [,3]     [,4] [,5]
## [1,]    2 58018.96 29009.48 39180.15    0
## [2,]   42    31.10     0.74       NA   NA
tabela
## Analysis of Variance Table
## 
## Response: z
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## fator      2  58019 29009.5   39180 < 2.2e-16 ***
## Residuals 42     31     0.7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
if(tabela$`Pr(>F)`[1] <= alpha){print("Rejeito a hipótese de igualdade de médias" )}else{print("Não rejeito a hipótese de igualdade de médias")}
## [1] "Rejeito a hipótese de igualdade de médias"

7 Exemplo com 1 fator e 3 grupos

#Exemplo de Analise de Variancia
y1 = round(rnorm(10, 200, sqrt(1)),1)
y2 = round(rnorm(15, 200, sqrt(1)),1)
y3 = round(rnorm(20, 200, sqrt(1)),1)

y = c(y1, y2, y3)
#y
fator = factor(rep(1:3,c(10,15,20)))

hist(y, freq=F, ylab="densidade", main="Histograma")

#empilhando os dados
z     = c(y1,y2, y3)
#cbind(fator, z)

#testando normalidade pelo Shapiro
alpha = 0.05
s     = 0
K     = 3 #total de modelos diferentes
for(i in 1:K)
{   
  teste = shapiro.test(z[fator==i]) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando normalidade pelo Kolmogorov-Smirnov
s=0
for(i in 1:K)
{   
  teste = ks.test(z[fator==i],"pnorm",mean(z[fator==i]),sd(z[fator==i])) 
  if(teste$p.value<=alpha){s=s+1}   
}
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando homocedasticidade
teste = bartlett.test(z,fator)
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
#testando homocedasticidade
library(lawstat)
teste = levene.test(z, fator, location = c("median"))
if(teste$p.value>alpha){print("homocedasticidade ok")}else{print("problema com a homocedasticidade")}
## [1] "homocedasticidade ok"
tabela = anova(lm(z~fator))

res.anova=aov(z~fator)
res.anova
## Call:
##    aov(formula = z ~ fator)
## 
## Terms:
##                    fator Residuals
## Sum of Squares   2.34817  48.42383
## Deg. of Freedom        2        42
## 
## Residual standard error: 1.073754
## Estimated effects may be unbalanced
summary(res.anova)
##             Df Sum Sq Mean Sq F value Pr(>F)
## fator        2   2.35   1.174   1.018   0.37
## Residuals   42  48.42   1.153
#Calculando os dados da tabela da ANOVA sem usar a 
#funcao do R

ybarra  = NULL
ni      = NULL
s2i     = NULL
for(i in 1:K)
{
  ybarra[i]   = mean(z[fator==i])
  ni[i]       = length(z[fator==i])
  s2i[i]      = var(z[fator==i])
}
n       = sum(ni)
yt      = mean(z)
SQT     = sum((z-yt)^2)
SQE     = sum(ni*(ybarra-yt)^2)
yb2     = rep(ybarra,ni)
SQD     = sum((z-yb2)^2)
c(SQT, SQD + SQE)
## [1] 50.772 50.772
MQE     = SQE/(K-1)
MQD     = SQD/(n-K)
MQT     = SQT/(n-1) #ST^2
s2      = ((ni[1]-1)*s2i[1] + (ni[2]-1)*s2i[2] + (ni[3]-1)*s2i[3])/(ni[1]+ni[2]+ni[3]-K) 
c(SQD, (n-K)*s2)
## [1] 48.42383 48.42383
EstatF      = MQE/MQD
c(SQE/(s2*(K-1)), EstatF)
## [1] 1.018331 1.018331
rbind(c(K-1, round(SQE,2), round(MQE,2), round(EstatF,3), 1-pf(EstatF, K-1, n-K)), c(n-K, round(SQD,2), round(MQD,2), NA, NA))
##      [,1]  [,2] [,3]  [,4]      [,5]
## [1,]    2  2.35 1.17 1.018 0.3699413
## [2,]   42 48.42 1.15    NA        NA
tabela
## Analysis of Variance Table
## 
## Response: z
##           Df Sum Sq Mean Sq F value Pr(>F)
## fator      2  2.348  1.1741  1.0183 0.3699
## Residuals 42 48.424  1.1529
if(tabela$`Pr(>F)`[1] <= alpha){print("Rejeito a hipótese de igualdade de médias" )}else{print("Não rejeito a hipótese de igualdade de médias")}
## [1] "Não rejeito a hipótese de igualdade de médias"

8 Modelo de regressão linear

Seja \(Y_i \stackrel{indep}{\sim} N(\boldsymbol{x}_i^T\boldsymbol{\beta}, \sigma^2)\), para \(i=1, \ldots, n\), \(\boldsymbol{x}_i^T = (x_{i,0}, \ldots, x_{i,p-1})\) e \(\boldsymbol{\beta} = (\beta_0, \ldots, \beta_{p-1})^T\). Podemos estimar os parâmetros desconhecidos \((\boldsymbol{\beta}, \sigma^2)\) através dos seguintes estimadores não viesados:

\[\begin{eqnarray} \boldsymbol{\hat{\beta}} &=& \left(\boldsymbol{x}^T\boldsymbol{x}\right)^{-1}\boldsymbol{x}^T\boldsymbol{Y}\\ \hat{\sigma}^2 &=& \frac{(\boldsymbol{Y}-\boldsymbol{x}\boldsymbol{\beta})^T(\boldsymbol{Y}-\boldsymbol{x}\boldsymbol{\beta})}{n-p}. \end{eqnarray}\]

8.1 ANOVA para o modelo de regressão linear simples

Seja \(Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 x_i, \sigma^2), \; i=1, \ldots, n\). Tem-se que \[\begin{eqnarray} S_{xy} &=& \sum_{i=1}^n{(x_i-\bar{x})(Y_i-\bar{Y})}\\ S_{xx} &=& \sum_{i=1}^n{(x_i-\bar{x})^2}\\ \hat{\beta}_0 &=& \bar{Y}-\hat{\beta}_1\bar{x}\\ \hat{\beta}_1 &=& \frac{S_{xy}}{S_{xx}}\\ \hat{\sigma}^2 &=& \frac{\sum_{i=1}^n{(Y_i-\hat{\beta}_0-\hat{\beta}_1 x_i)^2}}{n-2} = \frac{\sum_{i=1}^n{\hat{e}_i^2}}{n-2} = \frac{SQR}{n-2} \end{eqnarray}\] sendo \(\bar{Y}=\frac{\sum_{i=1}^n{Y_i}}{n}\) e \(\bar{x}=\frac{\sum_{i=1}^n{x_i}}{n}\).

Fontes de Graus de Soma dos Média da soma F
variação liberdade quadrados dos quadrados
Regressão 1 SQReg = \(\frac{S^2_{xy}}{S_{xx}}\) MQReg = \(\frac{SQReg}{1}\) \(F_{MRLS} = \frac{MQReg}{MQR}\)
Residual n-2 SQR=\(\sum{\hat{e}_i^2}\) MQR = \(\frac{SQR}{n-2}\)
Total n-1 SQT=\(\sum(Y_i-\bar{Y})^2\) MQT = \(\frac{SQT}{n-1}\)

\(F_{MRLS} = \frac{MQReg}{MQR} \sim F(1, n-2)\).

8.2 Exemplo MRLS

#MRLS - gerando os dados
n       = 1000
xi      = rnorm(n)
x       = cbind(rep(1,n), xi)
beta    = c(10,1)
sigma2  = 1
y       = rnorm(n, x %*% beta, sqrt(sigma2))

#analisando graficamente os dados gerados
par(mfrow=c(1,2))
hist(y, freq=F, ylab="densidade", main="Histograma")
plot(xi, y, lwd=2, bty="n", xlab="x")
lines(xi, x %*% beta, col="red", lwd=2)

#testando normalidade pelo Shapiro
alpha = 0.05
s     = 0
teste = shapiro.test(y) 
if(teste$p.value<=alpha){s=s+1} 
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#testando normalidade pelo Kolmogorov-Smirnov
s=0
teste = ks.test(y,"pnorm",mean(y),sd(y)) 
if(teste$p.value<=alpha){s=s+1} 
if(s==0){print("normalidade ok")} else{print("problema com a normalidade")}
## [1] "normalidade ok"
#supondo beta0, beta1 e sigma2 desconhecidos
#ajustando um MRLS
ajuste = lm(y~xi)
tabela = anova(ajuste)

res.anova=aov(y~xi)
res.anova
## Call:
##    aov(formula = y ~ xi)
## 
## Terms:
##                       xi Residuals
## Sum of Squares  993.7187  980.1347
## Deg. of Freedom        1       998
## 
## Residual standard error: 0.991009
## Estimated effects may be unbalanced
summary(res.anova)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## xi            1  993.7   993.7    1012 <2e-16 ***
## Residuals   998  980.1     1.0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tabela
## Analysis of Variance Table
## 
## Response: y
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## xi          1 993.72  993.72  1011.8 < 2.2e-16 ***
## Residuals 998 980.13    0.98                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Calculando os dados da tabela da ANOVA sem usar a 
#funcao do R

Sxy           = sum((xi-mean(xi))*(y-mean(y)))
Sxx           = sum((xi-mean(xi))^2)
beta1chapeu   = Sxy/Sxx
beta0chapeu   = mean(y)-beta1chapeu*mean(xi)
SQR           = sum((y-beta0chapeu-beta1chapeu*xi)^2)
MQR           = SQR/(n-2)
sigma2chapeu  = MQR
SQReg         = Sxy^2/Sxx
MQReg         = SQReg/1
SQT           = sum((y-mean(y))^2)
MQT           = SQT/(n-1)
F_MRLS        = MQReg/MQR

round(rbind(c(1, SQReg, MQReg, F_MRLS, 1-pf(F_MRLS, 1, n-2)), c(n-2, SQR, MQR, NA, NA)),2)
##      [,1]   [,2]   [,3]    [,4] [,5]
## [1,]    1 993.72 993.72 1011.83    0
## [2,]  998 980.13   0.98      NA   NA
tabela
## Analysis of Variance Table
## 
## Response: y
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## xi          1 993.72  993.72  1011.8 < 2.2e-16 ***
## Residuals 998 980.13    0.98                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
round(c(beta0chapeu, beta[1]),3)
## [1]  9.988 10.000
round(c(beta1chapeu, beta[2]),3)
## [1] 0.99 1.00
round(c(sigma2chapeu, sigma2),3)
## [1] 0.982 1.000

8.3 ANOVA para o modelo de regressão linear múltiplo

Podemos comparar os seguintes modelos:

\[\begin{eqnarray} H_0: \quad & Y_i &\stackrel{indep.}{\sim} N(\boldsymbol{x}_{M0} \boldsymbol{\beta}_{M0}, \sigma^2_{M0})\\ H_1: \quad & Y_i &\stackrel{indep.}{\sim} N(\boldsymbol{x}_{M1} \boldsymbol{\beta}_{M1}, \sigma^2_{M1}), \end{eqnarray}\] sendo \(q\) a dimensão de \(\boldsymbol{\beta}_{M0}\) e \(p\) a dimensão de \(\boldsymbol{\beta}_{M1}\), com \(q < p\).

Dessa forma, temos que a tabela de análise de variância é:

Fontes de Graus de Soma dos Média da F
variação liberdade quadrados soma dos
quadrados
Modelo com \(\boldsymbol{\beta}_0\) q \(\boldsymbol{\beta_{M0}}^T \boldsymbol{x_{M0}}^T \boldsymbol{y}\)
Melhoria devido ao modelo M1 p-q \(SQReg = \boldsymbol{\beta_{M1}}^T\boldsymbol{x_{M1}}^T\boldsymbol{y} - \boldsymbol{\beta_{M0}}^T \boldsymbol{x_{M0}}^T \boldsymbol{y}\) \(MQReg = \frac{SQReg}{p-q}\) \(F_{MRL} = \frac{MQReg}{MQR}\)
Residual n-p \(SQR = \boldsymbol{Y}^T\boldsymbol{Y} - \boldsymbol{\beta_{M1}}^T\boldsymbol{x_{M1}}^T\boldsymbol{y}\) \(MQR = \frac{SQR}{n-p}\)

\(F_{MRL} = \frac{D_0-D_1}{p-q} \frac{n-p}{D_1} = \frac{MQReg}{MQR} \sim F(p-q, n-p)\)

8.4 Exemplo MRLS (continuação)

#modelo com somente intercepto
q               = 1 #quantidade de parametros
x_M0            = rep(1,n) 
betachapeu_M0   = solve(t(x_M0) %*% x_M0) %*% t(x_M0) %*% y
sigma2chapeu_M0 = t(y-x_M0%*%betachapeu_M0) %*% (y-x_M0%*%betachapeu_M0) / (n-q)
SQ0             = betachapeu_M0*t(x_M0) %*% y
SQ0
##          [,1]
## [1,] 99709.96
#modelo com intercepto + coef. linear
p       = 2 #quantidade de parametros
betachapeu_M1   = solve(t(x) %*% x) %*% t(x) %*% y
rbind(c(beta0chapeu, beta1chapeu), t(betachapeu_M1))
##                      xi
## [1,] 9.988246 0.9898746
## [2,] 9.988246 0.9898746
sigma2chapeu_M1 = t(y-x%*%betachapeu_M1) %*% (y-x%*%betachapeu_M1) / (n-p)
c(sigma2chapeu_M1, sigma2chapeu)
## [1] 0.9820989 0.9820989
SQReg2     = t(betachapeu_M1) %*% t(x) %*% y - betachapeu_M0*t(x_M0)%*% y
c(SQReg2, SQReg)
## [1] 993.7187 993.7187
MQReg2     = SQReg2 / (p-q)
c(MQReg2, MQReg)
## [1] 993.7187 993.7187
SQR2  = t(y)%*%y - t(betachapeu_M1) %*% t(x) %*% y
c(SQR2, SQR)
## [1] 980.1347 980.1347
MQR2  = SQR2/(n-p)
c(MQR2, MQR)
## [1] 0.9820989 0.9820989
SQT2    = t(y)%*%y
SQT2
##          [,1]
## [1,] 101683.8
F_MRL  = MQReg2/MQR2
c(F_MRL, F_MRLS)
## [1] 1011.832 1011.832
#c(q, SQ0, NA, NA, NA)
rbind(c(p-q, SQReg2, MQReg2, F_MRL, 1-pf(F_MRL, p-q, n-p)),
c(n-p, SQR2, MQR2, NA, NA))
##      [,1]     [,2]        [,3]     [,4] [,5]
## [1,]    1 993.7187 993.7186760 1011.832    0
## [2,]  998 980.1347   0.9820989       NA   NA
tabela
## Analysis of Variance Table
## 
## Response: y
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## xi          1 993.72  993.72  1011.8 < 2.2e-16 ***
## Residuals 998 980.13    0.98                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.5 Exemplo MRL múltiplo

#MRLS - gerando os dados
n       = 1000
x1      = rnorm(n)
x2      = rbeta(n,2,3)
x2      = (x2-mean(x2))/sd(x2)
x       = cbind(rep(1,n), x1, x2)
beta    = c(0.5, 1, -2)
sigma2  = 0.05
y       = rnorm(n, x %*% beta, sqrt(sigma2))

#analisando graficamente os dados gerados
par(mfrow=c(2,2))
hist(y, freq=F, ylab="densidade", main="Histograma")
plot(x1, y, lwd=2, bty="n", xlab="x")
points(x1, x %*% beta, col="red", lwd=2)
plot(x2, y-x1*beta[2], lwd=2, bty="n", xlab="x")
points(x2, x %*% beta, col="red", lwd=2)

#modelo com intercepto + uma covariavel apenas
q               = 2 #quantidade de parametros
x_M0            = cbind(rep(1,n), x1)
betachapeu_M0   = solve(t(x_M0) %*% x_M0) %*% t(x_M0) %*% y
sigma2chapeu_M0 = t(y-x_M0%*%betachapeu_M0) %*% (y-x_M0%*%betachapeu_M0) / (n-q)
SQ0             = t(betachapeu_M0) %*% t(x_M0) %*% y

#modelo com intercepto + 2 covariaveis
p               = 3 #quantidade de parametros
x               = cbind(rep(1,n), x1, x2)
betachapeu_M1   = solve(t(x) %*% x) %*% t(x) %*% y
rbind(t(betachapeu_M1), beta)
##                       x1        x2
##      0.4988824 0.9980422 -1.994341
## beta 0.5000000 1.0000000 -2.000000
sigma2chapeu_M1 = t(y-x%*%betachapeu_M1) %*% (y-x%*%betachapeu_M1) / (n-p)
c(sigma2chapeu_M1, sigma2)
## [1] 0.05290923 0.05000000
SQReg           = t(betachapeu_M1) %*% t(x) %*% y - t(betachapeu_M0) %*% t(x_M0)%*% y
MQReg           = SQReg / (p-q)
SQR             = t(y)%*%y - t(betachapeu_M1) %*% t(x) %*% y
MQR             = SQR/(n-p)

F_MRL          = MQR/MQReg

rbind(c(p-q, SQR, MQR, F_MRL, 1-pf(F_MRL, p-q, n-p)),
c(n-p, SQReg, MQReg, NA, NA))
##      [,1]      [,2]         [,3]         [,4]      [,5]
## [1,]    1   52.7505 5.290923e-02 1.334769e-05 0.9970857
## [2,]  997 3963.9248 3.963925e+03           NA        NA