Capítulo III. Movendo-se para dimensões superiores


Vimos nos capítulos anteriores como os dispositivos gráficos simples podem ajudar a entender a estrutura e a dependência dos dados. As ferramentas gráficas foram baseadas em representações de dados univariadas (bivariadas) ou em transformações lisadas de informações multivariadas perceptíveis pelo olho humano. A maioria das ferramentas é extremamente útil em uma etapa de modelagem, mas, infelizmente, não fornece a imagem completa do conjunto de dados. Uma razão para isso é que as ferramentas gráficas apresentaram captura apenas em certas dimensões dos dados e não se concentram necessariamente nessas dimensões ou sub-partes dos dados, em análise que carregam informações estruturais máximas. Serão apresentadas ferramentas poderosas para reduzir a dimensão de um conjunto de dados posteriormente. Neste capítulo, como ponto de partida, ferramentas simples e básicas são usadas para descrever a dependência. Elas são construídas a partir de fatos elementares da teoria da probabilidade e das estatísticas introdutórias, por exemplo, covariância e correlação entre duas variáveis.

As Seções III.1 e III.2 mostram como lidar com esses conceitos em uma configuração multivariada e como um teste simples de correlação entre duas variáveis pode ser derivado. Como as relações lineares estão envolvidas nessas medidas, a Seção III.4 apresenta o modelo linear simples para duas variáveis e lembra o teste \(t\) básico para a inclinação. Na Seção III.5, um exemplo simples de análise de variância unifacorial introduz as notações para o teste \(F\) bem conhecido.

Devido ao poder da notação matricial, tudo isso pode ser facilmente estendido a uma configuração multivariada mais geral. A Seção III.3 mostra como as operações matriciais podem ser usadas para definir estatísticas resumidas de um conjunto de dados e para obter os momentos empíricos das transformações lineares dos dados. Esses resultados serão muito úteis na maioria dos próximos capítulos.

Finalmente, a notação matricial nos permite introduzir o modelo linear múltiplo flexível, onde relações mais gerais entre variáveis podem ser analisadas. Na Seção III.6, o ajuste do modelo por mínimos quadrados e as estatísticas usuais de teste são apresentadas com sua interpretação geométrica. Usando essas anotações, o modelo ANOVA é apenas um caso específico do modelo linear múltiplo.


III.1 Covariância


A covariância é uma medida de dependência entre variáveis aleatórias. Dadas duas variáveis aleatórias \(X\) e \(Y\), a covariância é definida por: \[ \sigma_{XY}= \mbox{Cov}(X,Y)=\mbox{E}(XY)-\mbox{E}(X)\mbox{E}(Y)\cdot \]

A definição precisa dos valores esperados foi dada no Capítulo IV. Se \(X\) e \(Y\) forem independentes um do outro, a covariância \(\mbox{Cov}(X,Y)\) é necessariamente igual a zero, consulte o Teorema III.1. O inverso não é verdadeiro. A covariância de \(X\) consigo mesma é a variância: \[ \sigma_{XX}=\mbox{Var}(X)=\mbox{Cov}(X,X)\cdot \]

se \(X\) for uma variável \(p\)-dimensional, ou seja, se \(X=(X_1,X_2,\cdots,X_p)^\top\), então a covariância teórica entre seus elementos pode ser escrita de forma matricial, ou seja, a matriz de covariâncias é da forma \[ \Sigma=\begin{pmatrix} \sigma_{X_1X_1} & \cdots & \sigma_{X_1X_p} \\ \vdots & \ddots & \vdots \\ \sigma_{X_pX_1} & \cdots & \sigma_{X_pX_p}\end{pmatrix}\cdot \]

As propriedades das matrizes de covariância foram detalhadas no capítulo IV. As versões empíricas dessas quantidades são: \[ s_{XY}=\frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})(y_i-\overline{y}) \qquad \mbox{e} \qquad s_{XX}=\frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})^2\cdot \]

Para \(n\) pequeno, digamos \(n\leq 20\), devemos substituir o fator \(1/n\) acima por \(1/(n-1)\) para corrigir um pequeno viés. Para uma variável aleatória \(p\)-dimensional, obtém-se a matriz empírica de covariância: \[ S=\begin{pmatrix} s_{X_1X_1} & \cdots & s_{X_1X_p} \\ \vdots & \ddots & \vdots \\ s_{X_pX_1} & \cdots & s_{X_pX_p}\end{pmatrix}\cdot \]

Para um gráfico de dispersão de duas variáveis, as covariâncias medem o quão perto a dispersão está de uma linha. Os detalhes matemáticos seguem, mas já deve ser entendido aqui que, nesse sentido, a covariância mede apenas a dependência linear.


Exemplo III.1

Se \(X\) for o conjunto de dados bancários inteiro, obtém-se a matriz de covariância, conforme indicado abaixo:

library(mclust)
data(banknote) 
cov(banknote[,2:7])
##               Length        Left       Right     Bottom        Top    Diagonal
## Length    0.14179296  0.03144322  0.02309146 -0.1032462 -0.0185407  0.08430553
## Left      0.03144322  0.13033945  0.10842739  0.2158028  0.1050394 -0.20934196
## Right     0.02309146  0.10842739  0.16327412  0.2841319  0.1299967 -0.24047010
## Bottom   -0.10324623  0.21580276  0.28413191  2.0868781  0.1645389 -1.03699623
## Top      -0.01854070  0.10503945  0.12999673  0.1645389  0.6447234 -0.54961482
## Diagonal  0.08430553 -0.20934196 -0.24047010 -1.0369962 -0.5496148  1.32771633

A covariância empírica entre \(X_4\) (Bottom) e \(X_5\) (Top), isto é, \(s_{X_4 X_5}\) é encontrada na linha 4 e na coluna 5. O valor 0.1645389. É óbvio que esse valor é positivo? No Exercício III.1, discutiremos mais essa questão.

Se \(X_f\) denota as notas falsificadas, obtemos: 

cov(banknote[banknote[,1]=="counterfeit",2:7])
##               Length         Left         Right      Bottom           Top
## Length    0.12401111  0.031515152  0.0240010101 -0.10059596  0.0194353535
## Left      0.03151515  0.065050505  0.0467676768 -0.02404040 -0.0119191919
## Right     0.02400101  0.046767677  0.0889404040 -0.01857576  0.0001323232
## Bottom   -0.10059596 -0.024040404 -0.0185757576  1.28131313 -0.4901919192
## Top       0.01943535 -0.011919192  0.0001323232 -0.49019192  0.4044555556
## Diagonal  0.01156566 -0.005050505  0.0341919192  0.23848485 -0.0220707071
##              Diagonal
## Length    0.011565657
## Left     -0.005050505
## Right     0.034191919
## Bottom    0.238484848
## Top      -0.022070707
## Diagonal  0.311212121

Para \(X_g\), as notas genuínas: 

cov(banknote[banknote[,1]=="genuine",2:7])
##               Length        Left       Right        Bottom         Top
## Length   0.150241414  0.05801313  0.05729293  0.0571262626  0.01445253
## Left     0.058013131  0.13257677  0.08589899  0.0566515152  0.04906667
## Right    0.057292929  0.08589899  0.12626263  0.0581818182  0.03064646
## Bottom   0.057126263  0.05665152  0.05818182  0.4132070707 -0.26347475
## Top      0.014452525  0.04906667  0.03064646 -0.2634747475  0.42118788
## Diagonal 0.005481818 -0.04306162 -0.02377778 -0.0001868687 -0.07530909
##               Diagonal
## Length    0.0054818182
## Left     -0.0430616162
## Right    -0.0237777778
## Bottom   -0.0001868687
## Top      -0.0753090909
## Diagonal  0.1998090909


Observe que a covariância entre \(X_4\) (distância da estrutura até a borda inferior) e \(X_5\) (distância do quadro até a borda superior) é negativa no ssub-grupos anteriores. Por que isso aconteceria? No Exercício III.2, discutiremos essa questão com mais detalhes.

À primeira vista, as matrizes \(S_f\) e \(S_g\) parecem diferentes, mas criam quase os mesmos gráficos de dispersão. Da mesma forma, a análise comum dos componentes principais no Capítulo 11 sugere uma análise conjunta da estrutura de covariância como em Flury and Riedwyl (1988).

Mostramos a seguir o gráfico de dispersão de variáveis de \(X_4\) vs. \(X_5\) do conjunto de dados bancários inteiros assim como nos sub-grupos de notas falsa e genuínas. 

par(mfrow=c(1,3), mar=c(4,4,5,1), pch = 19)
plot(banknote[,5], banknote[,6], xlab = expression(X[4]), ylab = expression(X[5]), 
     xlim = c(7,13), ylim = c(7,13), main="Dispersão no conjunto de dados \n bancários completo")
points(banknote$Bottom[banknote$Status=="genuine"], banknote$Top[banknote$Status=="genuine"], col = "red")
grid()
plot(banknote[banknote[,1]=="counterfeit",5], banknote[banknote[,1]=="counterfeit",6],  
     xlab = expression(X[4]), ylab = expression(X[5]), xlim = c(7,13), 
     ylim = c(7,13), main="Dispersão no conjunto de dados \n de notas bancárias falsas")
grid()
plot(banknote[banknote[,1]=="genuine",5], banknote[banknote[,1]=="genuine",6], xlab = expression(X[4]), 
     ylab = expression(X[5]), xlim = c(7,13), ylim = c(7,13),
     main="Dispersão no conjunto de dados \n de notas bancárias genuínas")
grid()

library(lattice)
xyplot(banknote[,5] ~ banknote[,6] | banknote[,1], xlab = expression(X[4]), ylab = expression(X[5]))

Os gráficos de dispersão acima são inclinadas para cima, mostram variáveis com covariância positiva. Gráficos de dispersão com estrutura de inclinação descendente têm covariância negativa. Observe que no gráfico de dispersão de \(X_4\) vs. \(X_5\) de todo o conjunto de dados bancários a nuvem de pontos está no topo de cima. No entanto, as duas sub-nuvens de notas bancárias falsificadas e genuínas estão inclinadas para baixo.


Exemplo III.2.

Um gerente de loja têxtil está estudando as vendas de pulôveres “Blue Classic” em dez períodos diferentes. Ele observa o número de pulôveres vendidos \(X_1\), variação no preço \(X_2\) em euros, os custos de propaganda nos jornais locais \(X_3\) em euros e a presença de um assistente de vendas \(X_4\) em horas por período. Ao longo dos períodos, ele observa a seguinte matriz de dados:

pullover = read.table("http://leg.ufpr.br/~lucambio/MSM/pullover.txt", header = TRUE, sep = " ")
head(pullover)
##   Sales Price Advert Asst.Hours
## 1   230   125    200        109
## 2   181    99     55        107
## 3   165    97    105         98
## 4   150   115     85         71
## 5    97   120      0         82
## 6   192   100    150        103

Ele está convencido de que o preço deve ter uma grande influência no número de pulôveres vendidos. Então ele faz um gráfico de dispersão de \(X_2\) vs. \(X_1\), mostrado abaixo.

par(mfrow=c(1,1), mar=c(4,4,1,1), pch=19)
plot(pullover$Price, pullover$Sales, ylab = "Vendas", xlab = "Preço")
grid()

Uma impressão aproximada é que a nuvem é um pouco descendente. Um cálculo da covariância empírica gera um valor negativo como esperado.

cov(pullover)
##                 Sales      Price    Advert Asst.Hours
## Sales      1152.45556  -88.91111 1589.6667   301.6000
## Price       -88.91111  244.26667  102.3333  -101.7556
## Advert     1589.66667  102.33333 2915.5556   233.6667
## Asst.Hours  301.60000 -101.75556  233.6667   197.0667


A função de covariância depende da escala. Assim, se os preços deste exemplo estivessem em ienes japoneses (JPY), obteríamos uma resposta diferente (consulte o Exercício III.16). Uma medida de dependência (linear) independente da escala é a correlação, que introduzimos na próxima seção.


III.2 Correlação


A correlação entre \(X\) e \(Y\) é definida da covariância como segue: \[ \rho_{XY}=\dfrac{\mbox{Cov}(X,Y)}{\sqrt{\mbox{Var}(X)\mbox{Var}(Y)}}\cdot \]

A vantagem da correlação é que ela é independente da escala, ou seja, alterar a escala de medição das variáveis não altera o valor da correlação. Portanto, a correlação é mais útil como uma medida de associação entre duas variáveis aleatórias que a covariância. A versão empírica do \(\rho_{XY}\) é a seguinte: \[ r_{XY}=\dfrac{S_{XY}}{\sqrt{S_{XX} S_{YY}}}\cdot \]

A correlação é, em valor absoluto, sempre menor que 1. É zero se a covariância for zero e viceversa. Para vetores \(p\)-dimensionais \((X_1,\cdots,X_p)^\top\) temos a matriz de correlação teórica \[ P=\begin{pmatrix} \rho_{X_1 X_1} & \cdots & \rho_{X_1 X_p} \\ \vdots & \ddots & \vdots \\ \rho_{X_p X_1} & \cdots & \rho_{X_p X_p}\end{pmatrix}, \]

e sua versão empírica, a matriz de correlação empírica que pode ser calculada a partir das observações, \[ R=\begin{pmatrix} r_{X_1 X_1} & \cdots & r_{X_1 X_p} \\ \vdots & \ddots & \vdots \\ r_{X_p X_1} & \cdots & r_{X_p X_p}\end{pmatrix}\cdot \]


Exemplo III.3

Obtemos a seguinte matriz de correlação para os dados das notas bancárias genuínas:

cor(banknote[banknote[,1]=="genuine",2:7])
##              Length       Left      Right        Bottom         Top
## Length   1.00000000  0.4110529  0.4159765  0.2292752146  0.05745277
## Left     0.41105294  1.0000000  0.6639218  0.2420437898  0.20764186
## Right    0.41597649  0.6639218  1.0000000  0.2547217369  0.13289390
## Bottom   0.22927521  0.2420438  0.2547217  1.0000000000 -0.63156375
## Top      0.05745277  0.2076419  0.1328939 -0.6315637468  1.00000000
## Diagonal 0.03163896 -0.2645751 -0.1497015 -0.0006503468 -0.25959830
##               Diagonal
## Length    0.0316389581
## Left     -0.2645751130
## Right    -0.1497015279
## Bottom   -0.0006503468
## Top      -0.2595983041
## Diagonal  1.0000000000

e para as notas bancárias falsificadas:

cor(banknote[banknote[,1]=="counterfeit",2:7])
##               Length        Left         Right      Bottom           Top
## Length    1.00000000  0.35088416  0.2285333950 -0.25236121  0.0867814182
## Left      0.35088416  1.00000000  0.6148524544 -0.08327004 -0.0734828638
## Right     0.22853339  0.61485245  1.0000000000 -0.05502618  0.0006976718
## Bottom   -0.25236121 -0.08327004 -0.0550261778  1.00000000 -0.6809310004
## Top       0.08678142 -0.07348286  0.0006976718 -0.68093100  1.0000000000
## Diagonal  0.05887240 -0.03549615  0.2055160184  0.37766340 -0.0622089101
##             Diagonal
## Length    0.05887240
## Left     -0.03549615
## Right     0.20551602
## Bottom    0.37766340
## Top      -0.06220891
## Diagonal  1.00000000
Como observado anteriormente para \(\mbox{Cov}(X_4,X_5)\), a correlação entre \(X_4\), a distância do quadro até a borda inferior e \(X_5\), a distância da estrutura até a borda superior é negativa. Isso é natural, uma vez que a covariância e a correlação sempre têm o mesmo sinal, veja também o Exercício III.17.


Por que a correlação é uma estatística interessante para estudar? Está relacionado à independência de variáveis aleatórias, que definiremos mais formalmente mais tarde. No momento, podemos pensar na independência como o fato de uma variável não ter influência sobre outra.


Teorema III.1.

Se \(X\) e \(Y\) são independentes, então \(\rho_{X,Y}=\mbox{Cov}(X,Y)=0\).


Exemplo III.4

Considere uma variável aleatória distribuída normal padrão e uma variável aleatória \(Y=X^2\), que certamente não é independente de \(X\). Aqui temos \[ \mbox{Cov}(X,Y)=\mbox{E}(XY)-\mbox{E}(X)\mbox{E}(Y)=\mbox{E}(X^3)=0\cdot \] Portanto \(\rho_{XY}=0\), também. Este exemplo também mostra que correlações e covariâncias medem apenas dependência linear. A dependência quadrática de \(Y=X^2\) em \(X\) não é refletida por essas medidas de dependência.


Para duas variáveis aleatórias normais, o inverso do Teorema III.1 é verdadeiro: a covariância zero para duas variáveis aleatórias normalmente distribuídas implica independência. Isso será mostrado posteriormente no Corolário V.2.

O Teorema III.1 nos permite verificar a independência entre os componentes de uma variável aleatória normal bivariada. Ou seja, podemos usar a correlação e testar se é zero. A distribuição do \(r_{XY}\) para um vetor \((X,Y))\) arbitrário; infelizmente é complicado. A distribuição do \(r_{XY}\) será mais acessível se \((X,Y)\) é conjuntamente normal. Se transformarmos a correlação pela transformação \(Z\) de Fisher, \[ W=\dfrac{1}{2}\left( \dfrac{1+r_{XY}}{1-r_{XY}}\right), \]

obtemos uma variável que possui uma distribuição mais acessível.

Sob a hipótese de que \(\rho=0\), \(W\) tem uma distribuição normal assintótica. As aproximações da esperança e variância de \(W\) são dadas pelo seguinte: \[ \mbox{E}\approx \dfrac{1}{2}\log\left( \dfrac{1+\rho_{XY}}{1-\rho_{XY}}\right), \]

\[ \mbox{Var}(W)\approx \dfrac{1}{n-3}\cdot \]

A distribuição é dada no seguinte teorema.


Teorema III.2.

\[ Z = \dfrac{W-\mbox{E}(W)}{\sqrt{\mbox{Var}(W)}} \overset{L}{\longrightarrow} N(0,1)\cdot \]


O Teorema III.2 nos permite testar diferentes hipóteses sobre correlação. Podemos corrigir o nível de significância \(\alpha\)˛ a probabilidade de rejeitar uma hipótese verdadeira e rejeitar a hipótese se a diferença entre o valor hipotético e o valor calculado de \(Z\) for maior que o valor crítico correspondente da distribuição normal. O exemplo a seguir ilustra o procedimento.


Exemplo III.5

Vamos estudar a correlação entre milhagem (\(X_2\)) e peso (\(X_8\)) para o conjunto de dados dos carros, Exemplo I.2, em que \(n=74\). temos \(r_{X_2 X_8}=-0.823\). Nossas conclusões do BoxPlot, Exemplo I.2, na qual afirmamos que os carros japoneses geralmente têm melhor quilometragem que os outros, precisam ser revisados. Da figura abaixo e do valor de \(r_{X_2 X_8}\), podemos ver que a milhagem está altamente correlacionada com o peso e que os carros japoneses na amostra são de fato mais leves que os outros. 

par(mfrow=c(1,1), pch=19)
car = read.table("http://leg.ufpr.br/~lucambio/MSM/CarData.txt", header = TRUE, sep = " ")
plot(car$M, car$W, xlim = c(10,45), ylim = c(1500,5000), xlab = expression(paste("Milhagem (",X[2],")")),
     ylab = expression(paste("Peso (",X[8],")")), axes = FALSE)
axis(1)
axis(2, las = 1)
Sede.da.empresa = factor(car$C, levels = 1:3, labels = c("U.S.","Japan", "EU"))
points(car$M[Sede.da.empresa=="U.S."], car$W[Sede.da.empresa=="U.S."], col = "red")
points(car$M[Sede.da.empresa=="Japan"], car$W[Sede.da.empresa=="Japan"], col = "blue")
legend(35,4500, legend = c("U.S.","Japan","EU"), col = c("red","blue","black"), lty = 2)
box()
grid()

Se queremos saber se o \(\rho_{X_2 X_8}\) é significativamente diferente de 0, aplicamos a transformação \(Z\) de Fisher. Isso nos dá:

\[ W=\dfrac{1}{2}\log\left( \dfrac{1+r_{X_2 X_8}}{1-r_{X_2 X_8}}\right)=-1.166 \]

e

\[ z = \dfrac{-1.166-0}{\sqrt{1/71}} = -9.825, \]

ou seja, um valor altamente significativo para rejeitar a hipótese de que \(H_0 \ : \, \rho =0\), os quantis de 2.5% e 97.5% da distribuição normal são -1.96 e 1.96, respectivamente.

Se queremos testar a hipótese de que, digamos, \(H_0 \, : \, \rho=-0.75\), obtemos: \[ z=\dfrac{-1.166-(-0.973)}{\sqrt{1/71}}=-1.627\cdot \] Este é um valor não significativo no nível \(\alpha=0.05\) para \(z\), pois é entre os valores críticos no nível de significância de 5%, ou seja, \(-1.96 <z <1.96\).


Exemplo III.6

Vamos considerar novamente o conjunto de dados dos pullovers, Exemplo III.2. Considere a correlação entre a presença dos assistentes de vendas (\(X_4\)) e o número de pulôveres vendidos (\(X_1\)). Aqui calculamos a correlação como \[ r_{X_1 X_4}=0.633 \] A transformação \(Z\) desse valor é \[ w= \dfrac{1}{2}\log\left( \dfrac{1+r_{X_1 X_4}}{1-r_{X_1 X_4}}\right) = 0.746\cdot \]

O tamanho da amostra é \(n=10\), portanto, para a hipótese \(H_0 \, : \, \rho_{X_1 X_4}=0\), a estatística a considerar é: \[ z=\sqrt{7}(0.746-0)=1.974 \] que é apenas estatisticamente significativo no nível de 5%, ou seja, 1.974 é apenas um pouco maior que 1.96. 


par(mfrow=c(1,1), pch=19)
plot(pullover$Asst.Hours, pullover$Sales, xlab = expression(paste("Horas de assistência de vendas (",X[4],")")),
     ylab = expression(paste("Vendas (",X[1],")")))
grid()



As propriedades de normalização e estabilizadora da variância de \(W\) são assintóticas. Além disso, o uso de \(W\) em pequenas amostras, para \(n\leq 25\) é melhorado pela transformação de Hotelling (Hotelling, 1953): \[ W^* = W-\dfrac{3W +\tanh(W)}{4(n-1)}, \qquad \mbox{com} \qquad \mbox{Var}(W^*)=\dfrac{1}{n-1}\cdot \] A variável transformada \(W^*\) é distribuída assintoticamente como uma distribuição normal.


Exemplo III.7

A partir da observação anterior, obtemos \(w^8=0.6663\) e \(\sqrt{10-1} w^* = 1.9989\) para o Exemplo III.6 anterior. Este valor é significativo no nível de 5%.


Observe que a transformação \(Z\) de Fisher é o inverso da função tangente hiperbólica: \(W=\tanh^{-1}(r_{XY}\); equivalentemente \[ r_{XY} = \tanh(W)=\dfrac{e^{2W}-1}{e^{2W}+1}\cdot \]

Sob as suposições de normalidade de \(X\) e \(Y\), podemos testar sua independência \(H_0 \, : \rho_{XY}=0\), usando a distribuição da estatística \(t-student\) exata \[ T = r_{XY} \sqrt{\dfrac{n-2}{1-r_{XY}^2}} \sim t_{n-2}\cdot \]

Definindo a probabilidade do erro tipo I como \(\alpha\), rejeitamos a hipótese nula \(H_0 \, : \rho_{XY}=0\) se \[ |T|\geq t_{1-\alpha/2;n-2}\cdot \] esta é a estatística utilizada nas seguintes situações. 


cor(car$M,car$W)
## [1] -0.8232254
library(ggstatsplot)
ggscatterstats( data = car, x = M, y = W, bf.message = FALSE )


# Teste coeficiente de correlação de Pearson
test <- cor.test(car$M, car$W)
test
## 
##  Pearson's product-moment correlation
## 
## data:  car$M and car$W
## t = -12.045, df = 69, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8863042 -0.7301567
## sample estimates:
##        cor 
## -0.8232254



III.3. Estatísticas de resumo


Esta seção se concentra na representação de estatísticas de resumo básico: médias, covariâncias e correlações, na notação da matricial, pois geralmente aplicamos transformações lineares aos dados. A notação da matricial nos permite derivar instantaneamente as características correspondentes das variáveis transformadas. A transformação Mahalanobis é um exemplo proeminente de tais transformações lineares.

Suponha que observamos \(n\) realizações de uma variável aleatória \(p\)-dimensional; temos uma matriz de dados \(X_{n\times p}\): \[ X=\begin{pmatrix} x_{11} & \cdots & x_{1p} \\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{np} \end{pmatrix}\cdot \] As linhas \(x_i=(x_{i1},\cdots,x_{ip})\in\mathbb{R}^p\) denotam as observações de uma variável aleatória \(p\)-dimensional \(X\in\mathbb{R}^p\).

As estatísticas que foram brevemente introduzidas em seções anteriores podem ser reescritas na forma da matrizes da seguinte forma. O centro de gravidade das \(n\) observações em \(\mathbb{R}^p\) é dado pelo vetor \(\overline{x}\) das médias \(\overline{x}_j\) das \(p\) variáveis: \[ \overline{x}=\begin{pmatrix} \overline{x}_1 \\ \vdots \\ \overline{x}_p\end{pmatrix}= \frac{1}{n}X^\top \pmb{1}_n\cdot \]

A dispersão das \(n\) observações pode ser caracterizada pela matriz de covariâncias das \(p\) variáveis. As covariâncias empíricas definidas são os elementos da seguinte matriz: \[ S=\frac{1}{n}X^\top X-\overline{x}\overline{x}^\top = \frac{1}{n}\Big( X^\top X-\frac{1}{n}X^\top \pmb{1}_n\pmb{1}_n^\top X \Big)\cdot \] Observe que esta matriz é equivalentemente definida por \[ S=\frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})(x_i-\overline{x})^\top\cdot \] A fórmula de covariância pode ser reescrita como \(S=\frac{1}{n}X^\top HX\) com a matriz \(H\), definda como \[ H=\mbox{I}_n-\frac{1}{n}\pmb{1}_n\pmb{1}_n^\top\cdot \]

Observe que a matriz \(H\) é simétrica e idempotente. De fato, \[ H^2 = \Big( \mbox{I}_n-\frac{1}{n}\pmb{1}_n\pmb{1}_n^\top\Big) \Big(\mbox{I}_n-\frac{1}{n}\pmb{1}_n\pmb{1}_n^\top \Big) = \mbox{I}_n-\frac{1}{n}\pmb{1}_n\pmb{1}_n^\top = H\cdot \] Como uma conseqüência, \(S\) é a semidefinita positiva, ou seja, \(S\geq 0\).

Na verdade, para todo \(a\in\mathbb{R}^p\), \[ a^\top S a= \frac{1}{n}a^\top X^\top HX a =\frac{1}{n}\big( a^\top X^\top H^\top\big) \big( HXa\big) \] isto porque \(H^\top H =H\), então \[ a^\top S a = \frac{1}{n}y^\top y = \frac{1}{n}\sum_{i=1}^p y_i^2\geq 0, \]

para \(y=HXa\).

É bem conhecido a partir do caso unidimensional que \(\frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})^2\) é uma estimativa da variância exibindo um viés da ordem \(1/n\) (Breiman, 1973). No caso multidimensional, \(S_u=n/(n-1)S\) é uma estimativa imparcial ou não-viciada da verdadeira covariância.

O coeficiente de correlação amostral entre a \(i\)-ésima variável e a \(j\)-ésima variável é \(r_{X_i X_j}\). Se \(D=\ diag\big( S_{X_i X_j}\big)\), então a matriz de correlação amostral é \[ R=D^{-1/2}SD^{-1/2}, \] onde \(D^{-1/2}\) é a matriz diagonal de elementos \(\big(S_{X_i X_j}\big)^{-1/2}\) na sua diagonal principal.


Exemplo III.8

As covariâncias empíricas são calculadas para o conjunto de dados de pulôver.

O vetor dos meios das quatro variáveis no conjunto de dados é

colMeans(pullover)
##      Sales      Price     Advert Asst.Hours 
##      172.7      104.6      104.0       93.8

A matriz de covariâncias amostrais é

cov(pullover)
##                 Sales      Price    Advert Asst.Hours
## Sales      1152.45556  -88.91111 1589.6667   301.6000
## Price       -88.91111  244.26667  102.3333  -101.7556
## Advert     1589.66667  102.33333 2915.5556   233.6667
## Asst.Hours  301.60000 -101.75556  233.6667   197.0667

A estimativa imparcial da matriz de covariâncias é igual a

(10/9)*cov(pullover)
##                 Sales      Price    Advert Asst.Hours
## Sales      1280.50617  -98.79012 1766.2963   335.1111
## Price       -98.79012  271.40741  113.7037  -113.0617
## Advert     1766.29630  113.70370 3239.5062   259.6296
## Asst.Hours  335.11111 -113.06173  259.6296   218.9630

e a matriz de correlção amostral é

cor(pullover)
##                 Sales      Price    Advert Asst.Hours
## Sales       1.0000000 -0.1675760 0.8672280  0.6328673
## Price      -0.1675760  1.0000000 0.1212619 -0.4637879
## Advert      0.8672280  0.1212619 1.0000000  0.3082688
## Asst.Hours  0.6328673 -0.4637879 0.3082688  1.0000000


Transformações lineares


Em muitas aplicações práticas, precisamos estudar transformações lineares dos dados originais. Isso motiva a questão de como calcular as estatísticas de resumo após essas transformações lineares.

Seja \(A\) uma matriz \(q\times p\) e considere a matriz de dados transformada \[ Y = XA^\top = (y_1,\cdots,y_n)^\top\cdot \]

A linha \(y_d=(y_{i1}\cdots,y_{iq})\in \mathbb{R}^q\) pode ser vista como a observação de uma variável aleatória \(q\)-dimensional \(Y=AX\). Na verdade, temos \(y_i=x_i A^\top\). Obtemos imediatamente a média e a covariância empírica das variáveis (colunas) que formam a matriz de dados \(Y\): \[ \overline{y}=\frac{1}{n}Y^\top \pmb{1}_n =\frac{1}{n}AX^\top \pmb{1}-n = A \overline{x} \] e \[ S_Y=\frac{1}{n}Y^\top HY =\frac{1}{n}A X^\top HXA^\top = A S_X A^\top\cdot \]

Observe que se a transformação linear não for homogênea, isto é, \(y_i=A x_i+b\), onde \(b_{q\times 1}\) somente a média muda: \(\overline{y}=A\overline{x}+b\).


Exemplo III.9

Suponha que \(X\) seja o conjunto de dados dos pulôver. O gerente quer calcular suas despesas médias para anúncio (\(X_3\)) e assistente de vendas (\(X_4\)).

Suponha que o assistente de vendas cobre um salário por hora de 10 euros. Em seguida, o gerente da loja calcula as despesas \(Y\) como \(Y=X_3+10X_4\). As expressões acima dizem que isso é equivalente a definir a matriz \(A_{.\times 1}\) como \(A=(0,0,1,10)\).

Agora, computacionalmente é muito fácil obter a média e variâncias amostrais das despesas gerais:

c(0,0,1,10)%*%colMeans(pullover)
##      [,1]
## [1,] 1042

e

c(0,0,1,10)%*%cov(pullover)%*%c(0,0,1,10)
##          [,1]
## [1,] 27295.56


Transformação de Mahalanobis


Um caso especial dessa transformação linear é \[ z_i=S^{-1/2}(x_i-\overline{x}) \] quando \(i=1,\cdots,n\). Observe que para a matriz de dados transformada \(Z=(z-1,\cdots,z_n)^\top\) \[ S_Z=\frac{1}{n}Z^\top HZ=\mbox{I}_p\cdot \]

Portanto, a transformação de Mahalanobis elimina a correlação entre as variáveis e padroniza a variância de cada variável.


III.4. Modelo linear para duas variáveis


Vimos várias vezes agora em gráficos de dispersão para baixo e para cima. O que o olho define aqui como uma inclinação? Suponha que possamos construir uma linha correspondente à direção geral da nuvem. O sinal da inclinação dessa linha corresponderia às direções para cima e para baixo. Chame a variável no eixo vertical \(Y\) e a do eixo horizontal \(X\). Uma linha de inclinação é uma linear relação entre \(X\) e \(Y\): \[ y_i=\alpha+\beta x_i+\epsilon_i, \qquad i=1,\cdots,n\cdot \] Aqui, \(\alpha\) é o intercepto e \(\beta\) a inclinação da linha. Os erros ou desvios da linha, são indicados como \(\epsilon_i\) e supõe-se que tenham média zero e variância finita \(\sigma^2\). A tarefa de encontrar \((\alpha,\beta)\) é chamada de ajuste linear.

Vamos derivar estimadores para \(\alpha\) e \(\beta\) mais formalmente, além de descrever com precisão o que é um bom estimador. Por enquanto, pode-se tentar encontrar um bom estimador \((\widehat{\alpha},\widehat{\beta})\) via técnicas gráficas. Uma técnica numérica e estatística muito comum é escolher \(\widehat{\alpha}\) e \(\widehat{\beta}\) que minimizam: \[ (\widehat{\alpha},\widehat{\beta})=\arg \min_{(\alpha,\beta)} \sum_{i=1}^n \big( y_i-\alpha-\beta x_i\big)^2\cdot \] A solução para esta tarefa são os estimadores: \[ \widehat{\beta}=\dfrac{S_{XY}}{S_{XX}} \] e \[ \widehat{\alpha}=\overline{y}-\widehat{\beta}\overline{x}\cdot \]

A variância de \(\widehat{\beta}\) é \[ \mbox{Var}(\widehat{\beta})=\dfrac{\sigma^2}{n S_{XX}}\cdot \]

O erro padrão (\(\mbox{SE}\)) desse estimador é \[ \mbox{SE}(\widehat{\beta})=\sqrt{\mbox{Var}(\widehat{\beta})}=\dfrac{\sigma}{\sqrt{n S_{XX}}}\cdot \]

Podemos usar essa fórmula para testar a hipótese de que \(\beta=0\). Em uma aplicação, a variância \(\sigma^2\) deve ser estimada por um estimador \(\widehat{\sigma}^2\) que será apresentado abaixo. Sob uma suposição de normalidade dos erros, o teste \(t\) para a hipótese \(\beta=0\) funciona da seguinte maneira.

Calculamos a estatística \[ t=\dfrac{\widehat{\beta}}{\mbox{SE}(\widehat{\beta})} \] e rejeita-se a hipótese em um nível de significância de 5% se \(|t|\geq t_{0.975;n-2}\), onde o quantil de 97.5% da distribuição do \(t\)-Student é claramente o valor crítico de 95% para o teste bilateral. Para \(n\geq 30\), isso pode ser substituído por 1.96, o quantil 97.5% da distribuição normal. Um estimador \(\widehat{\sigma}^2\) de \(\sigma^2\) será fornecido a seguir.


Exemplo III.10

Vamos aplicar o modelo de regressão linear às vendas de pullovers. O gerente de vendas acredita que há uma forte dependência do número de vendas em função do preço. Ele calcula a linha de regressão como mostrado. 

par(mfrow=c(1,1), pch=19)
plot(pullover$Price, pullover$Sales, xlab = "Preço", ylab = "Vendas")
abline(lm(Sales ~ Price, data = pullover), col = "red", lwd = 2)
grid()


ajuste.pullover = lm(Sales ~ Price, data = pullover)
summary(ajuste.pullover)
## 
## Call:
## lm(formula = Sales ~ Price, data = pullover)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.095  -8.898   2.036   9.805  64.725 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 210.7736    79.9837   2.635   0.0299 *
## Price        -0.3640     0.7571  -0.481   0.6435  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.5 on 8 degrees of freedom
## Multiple R-squared:  0.02808,    Adjusted R-squared:  -0.09341 
## F-statistic: 0.2311 on 1 and 8 DF,  p-value: 0.6435


Quão bom é esse ajuste? Isso pode ser julgado por meio de medidas de bondade de ajuste. Definindo \[ \widehat{y}_i=\widehat{\alpha}+\widehat{\beta} x_i, \]

como os valores preditos de \(y\) como funçao de \(x\). Com \(\widehat{y}\), o gerente da loja têxtil no exemplo acima pode prever as vendas em função dos preços \(x\). A variação na variável de resposta é: \[ n s_{YY}= \sum_{i=1}^n (y_1-\overline{y})^2\cdot \] A variação explicada pela regressão linear com os valores previstos definidos acima é: \[ \sum_{i=1}^n (\widehat{y}_i-\overline{y})^2\cdot \] A soma dos quadrados dos resíduos, é dada por: \[ RSS=\sum_{i=1}^n (y_i-\widehat{y}_i)^2\cdot \] Um estimador imparcial ou não viciado \(\widehat{\sigma}^2\) de \(\sigma^2\) é dado por \(RSS/(n-2)\). A seguinte relação é válida: \[ \underbrace{\sum_{i=1}^n (y_i-\overline{y})^2}_{\mbox{Variação total}} = \underbrace{\sum_{i=1}^n (\widehat{y}_i-\overline{y})^2}_{\mbox{Variação explicada}} + \underbrace{\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{\mbox{Variação não explicada}}\cdot \]

O coeficiente de determinação é: \[ R^2=\dfrac{\sum_{i=1}^n (\widehat{y}_i-\overline{y})^2}{\sum_{i=1}^n (y_i-\overline{y})^2}=\dfrac{\mbox{Variação explicada}}{\mbox{Variação total}}\cdot \]

O coeficiente de determinação aumenta com a proporção de variação explicada pela relação linear. Nos casos extremos em que \(R^2=1\), toda a variação é explicada pela regressão linear. O outro extremo \(R^2=0\), é onde a covariância empírica é \(S_{XY}=0\). o coeficiente de determinação pode ser reescrito como \[ R^2= 1- \dfrac{\sum_{i=1}^n (y_i-\widehat{y}_i)^2}{\sum_{i=1}^n (y_i-\overline{y})^2}=1-\dfrac{\mbox{Variação não explicada}}{\mbox{Variação total}}\cdot \]

Pode-se ver que na regressão linear, \(R^2=r_{XY}^2\), é o quadrado da correlação entre \(X\) e \(Y\).


Exemplo III.11

Para o exemplo de pulôver acima, estimamos \[ \widehat{\alpha}=210.7736 \qquad \mbox{e} \qquad \widehat{\beta}=-0.3640\cdot \] O coeficiente de determinação é \[ R^2= 0.02808\cdot \] O gerente da loja têxtil conclui que as vendas não são muito influenciadas pelo preço, de maneira linear. 


pullover$predicted <- predict(ajuste.pullover)  
pullover$residuals <- residuals(ajuste.pullover)
library(ggplot2)
ggplot(pullover, aes(x = Price, y = Sales)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +
  geom_segment(aes(xend = Price, yend = predicted), alpha = .2) +
  geom_point(aes(color = abs(residuals), size = abs(residuals))) + 
  scale_color_continuous(low = "black", high = "red") +
  guides(color = FALSE, size = FALSE) + geom_point(aes(y = predicted), shape = 1) +
  theme_bw()



Em geral, a regressão de \(Y\) em \(X\) é diferente da de \(X\) em \(Y\). Demonstraremos isso, mais uma vez, usando os dados do Swiss Bank Notes.


Exemplo III.12

Os mínimos quadrados ajustam as variáveis \(X_4\) (x) e \(X_5\) (y) das notas bancárias genuínas. A figura abaixo mostra a linha ajustada se \(X_5\) for aproximada por uma função linear de \(X_4\). Nesse caso, os parâmetros são 


summary(lm(Top ~ Bottom, data = banknote))
## 
## Call:
## lm(formula = Top ~ Bottom, data = banknote)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.02796 -0.53085 -0.00576  0.47463  1.70857 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.90798    0.37252  26.597   <2e-16 ***
## Bottom       0.07884    0.03910   2.016   0.0451 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7968 on 198 degrees of freedom
## Multiple R-squared:  0.02012,    Adjusted R-squared:  0.01517 
## F-statistic: 4.066 on 1 and 198 DF,  p-value: 0.04511


Se prevemos \(X_4\) por uma função de \(X_5\), chegaríamos a um intercepto e inclinação diferentes 


summary(lm(Bottom ~ Top, data = banknote))
## 
## Call:
## lm(formula = Bottom ~ Top, data = banknote)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4129 -1.1562 -0.1856  0.9912  3.6782 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.6994     1.3518   4.956 1.54e-06 ***
## Top           0.2552     0.1266   2.016   0.0451 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.434 on 198 degrees of freedom
## Multiple R-squared:  0.02012,    Adjusted R-squared:  0.01517 
## F-statistic: 4.066 on 1 and 198 DF,  p-value: 0.04511


A regressão linear de \(Y\) em \(X\) é dada minimizando erros verticais. A regressão linear de \(X\) em \(Y\) faz o mesmo, mas aqui os erros a serem minimizados no sentido dos mínimos quadrados são medidos horizontalmente. Como Visto no Exemplo III.12, as duas linhas de mínimos quadrados são diferentes, embora ambas medem, em certo sentido, a inclinação da nuvem de pontos.

Como mostrado no próximo exemplo, ainda há outra maneira de medir a direção principal de uma nuvem de pontos: está relacionada à decomposição espectral de matrizes de covariância.


Suponhamos temos a seguinte matriz de covariâncias: \[ \Sigma = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\cdot \] A figura mostra um gráfico de dispersão de uma amostra de duas variáveis aleatórias normais com uma matriz de covariância com \(\rho=0.8\). 


set.seed(64859)
library(mvtnorm)
amostra = rmvnorm(n = 150, mean = c(0,0), sigma = matrix(c(1,0.8,0.8,1), ncol = 2))
plot(amostra, pch =19, xlab = "X", ylab = "Y")
grid()


Os valores próprios de \(\Sigma\) são soluções de: \[ \left| \begin{pmatrix} 1-\lambda & \rho \\ \rho & 1-\lambda \end{pmatrix}\right| =0\cdot \]

Assim, \(\lambda_1=1+\rho\) e \(\lambda_2=1-\rho\). Portanto \(\Lambda = \mbox{diag}(1+\rho,1-\rho)\). O autovetor correspondente a \(\lambda_1=1_\rho\) pode ser calculado a partir do sistema de equações lineares: \[ \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = (1+\rho)\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \] ou \[ x_1+\rho x_2 = x_1+\rho x_2 \qquad \mbox{e} \qquad \rho x_1+x_2 = x_2+\rho x_2 \] logo \(x_1=x_2\).

O primeiro autovetor padronizado é \[ \gamma_1= \begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix}\cdot \]

A direção deste vetor próprio é a diagonal na figura acima e captura a principal variação nessa direção. Voltaremos a essa interpretação no Capítulo XI. O segundo vetor próprio, ortogonal para \(\gamma_1\), é \[ \gamma_2= \begin{pmatrix} 1/\sqrt{2} \\ -1/\sqrt{2} \end{pmatrix}\cdot \]

Finalmente, \[ \Gamma = (\gamma_1, \gamma_2) = \begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ 1/\sqrt{2} & -1/\sqrt{2} \end{pmatrix} \] e podemos verificar que \[ \Sigma = \Gamma \lambda \Gamma^\top \cdot \]


O primeiro vetor próprio captura a direção principal de uma nuvem de pontos. A regressão linear de \(Y\) em \(X\) e \(X\) em \(Y\) realizou, em certo sentido, a mesma coisa. Em geral, a direção do vetor próprio e da inclinação dos mínimos quadrados são diferentes. A razão é que o estimador de mínimos quadrados minimiza erros verticais ou horizontais, enquanto o primeiro vetor próprio corresponde a uma minimização ortogonal ao vetor próprio.


III.5. Análise de Variâncias simples


Em uma análise de variância simples, isto é, ANOVA com um fator, supõe-se que os valores médios da variável de resposta \(Y\) sejam induzidos por um fator simples. Suponha que esse fator assuma \(p\) valores e que, para cada nível do fator, temos \(m=n/p\) observações.

O objetivo de uma ANOVA simples é analisar a estrutura de observação \[ y_{kl}=\mu_l+\epsilon_{kl} \] para \(k=1,\cdots,m\) e \(l=1,\cdots,p\).

Cada fator tem um valor médio \(\mu_l\). Cada observação \(y_{kl}\) é assumida como sendo uma soma do valor médio do fator correspondente \(\mu_l\) e um erro aleatório \(\epsilon_{kl}\) com média zero. O valor linear do modelo de regressão cai neste esquema com \(m=1\), \(p=n\) e \(\mu_l=\alpha+\beta x_l\), onde \(x_i\) é o valor do nível \(i\) do fator.


Exemplo III.14

A empresa de pulôver analisa o efeito de três Estratégias de marketing
1. Anúncio no jornal local,
2. Presença de assistente de vendas,
3. Apresentação de luxo nas janelas da loja.
Todas essas estratégias são tentadas em dez lojas diferentes. As observações de venda resultantes são apresentadas: 


loja = rep(seq(1:10),3)
marketing = c(rep(1,10), rep(2,10), rep(3,10))
Vendas = c(9,11,10,12,7,11,12,10,11,13,10,15,11,15,15,13,7,15,13,10,18,14,17,19,14,17,16,14,17,15)
simple.pullover = cbind(loja,marketing,Vendas)


Existem \(p=3\) fatores e \(n=mp=30\) observações nos dados. A empresa quer saber se todas as três estratégias de marketing têm o mesmo efeito médio ou se existem diferenças. Ter o mesmo efeito significa que todos \(\mu_l\) são iguais a um valor \(mu\). A hipótese a ser testada é, portanto, \[ H_0 \, : \, \mu_l = \mu \qquad \mbox{para} \qquad l=1,\cdots,p\cdot \] A hipótese alternativa, de que as estratégias de marketing têm efeitos diferentes, podem ser formulados como \[ H_1 \, : \, \mu_l \neq \mu_k \qquad \mbox{para algum} \qquad l\neq k\cdot \] Isso significaria que uma estratégia de marketing é melhor que as outras. 


one.way <- aov(Vendas ~ factor(marketing))
summary(one.way)
##                   Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(marketing)  2  157.3   78.63   16.89 1.75e-05 ***
## Residuals         27  125.7    4.66                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Concluímos que as estratégias de marketing têm efeitos diferentes.


Testando as premissas da ANOVA


Para verificar se o modelo se encaixa na suposição de homoscedasticidade, vamos observar o gráfico: 

par(mfrow=c(2,2))
plot(one.way)

Os gráficos de diagnóstico mostram a variação inexplicada (resíduos) em todo o intervalo dos dados observados. Cada gráfico fornece uma informação específica sobre o ajuste do modelo, mas basta saber que a linha vermelha que representa a média dos resíduos deve ser horizontal e centrada em zero ou em um, no gráfico de localização de escala, o que significa que não há grandes valores discrepantes que possam causar viés no modelo.

O gráfico Q-Q normal traça uma regressão entre os resíduos teóricos de um modelo perfeitamente homocedástico e os resíduos reais do seu modelo, portanto, quanto mais próximo de uma inclinação de 1, melhor. Este gráfico Q-Q é muito próximo, com apenas um pouco de desvio.

A partir desses gráficos diagnósticos podemos dizer que o modelo se encaixa na suposição de homocedasticidade. Se o seu modelo não se encaixar na suposição de homocedasticidade, você pode tentar o teste de Levene.

Embora tenhamos encontrado que há diferenças entre as médias, o resultado da ANOVA só é robusto se as premissas do testes forem satisfeitas.

Uma vez que a ANOVA tem como premissa a homogeneidade dos dados, podemos testar se nossos dados são adequados nesse quesito com o Teste de Levene para homocedasticia. Para isso precisaremos do pacote car: 


library(car)
leveneTest(one.way)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.5064 0.2398
##       27


A hipótese nula do Teste de Levene é de que não há diferença entre as variâncias dos grupos. O \(p\)-valor maior do que 0.05 nos dá uma confiança estatística para afirmar que as variâncias são de fato iguais e portanto nossos dados são homogêneos.

Já premissa da ANOVA referente a normalidade dos resíduos pode ser testada através do teste de Shapiro-Wilk:


shapiro.test(resid(one.way))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(one.way)
## W = 0.94637, p-value = 0.135


A hipótese nula do Teste de Shapiro-Wilk é de que não há diferença entre a nossa distribuição dos dados e a distribuição normal. O \(p\)-valor maior do que 0.05 nos dá uma confiança estatística para afirmar que as distribuição dos nossos resíduos não difere da distribuição normal.

Dessa forma nossos dados satisfazem todas as premissas da ANOVA e portanto, o resultado da nossa ANOVA são válidos.

A ANOVA nos diz se há diferenças entre as médias dos grupos, mas não quais são as diferenças. Para descobrir quais grupos são estatisticamente diferentes um do outro, você pode realizar o teste de Tukey’s Honestly Significant Difference (Tukey’s HSD) para comparações de pares:


tukey.one.way<-TukeyHSD(one.way)
tukey.one.way
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Vendas ~ factor(marketing))
## 
## $`factor(marketing)`
##     diff        lwr      upr     p adj
## 2-1  1.8 -0.5924918 4.192492 0.1681816
## 3-1  5.5  3.1075082 7.892492 0.0000137
## 3-2  3.7  1.3075082 6.092492 0.0019179


Nos resultados dos testes vemos que existem diferenças estatisticamente significativas, ou seja, \(p\)-valor <0.05 entre os tipos de marketing 3 e 1 e entre os tipos de marketing 3 e 2, mas a diferença entre os tipos de marketing 2 e 1 não é estatisticamente significativa. 


plot(tukey.one.way, las = 1)
grid()


As diferenças significativas em grupo estão em que o intervalo de confiança de 95% não inclui zero. Essa é outra maneira de dizer que o \(p\)-valor para essas diferenças pareadas é <0.05. A partir deste gráfico, podemos ver que os tipos de marketing 3 e 1 e os tipos de marketing 3 e 2 são estatísticaente significativos.



III.6. Modelo linear multivariado


O modelo linear simples e o modelo de Análise de Variância podem ser vistos como um caso particular de um modelo linear mais geral onde as variações de uma variável \(Y\) são explicadas por \(p\) variáveis explicativas \(X\) respectivamente. Seja \(Y_{n\times 1}\) um vetor de observações na variável resposta e \(X_{n\times p}\) uma matriz de dados nas \(p\) variáveis explicativas. Uma aplicação importante da teoria desenvolvida é o ajuste de mínimos quadrados. A ideia é aproximar \(y\) por uma combinação linear \(\widehat{y}\) de colunas de \(X\), ou seja, \(\widehat{y}\in C(X)\). O problema é encontrar \(\widehat{\beta}\in\mathbb{R}^p\) tal que \(\widehat{y}=X\widehat{\beta}\) seja o melhor ajuste de \(y\) no sentido de mínimos quadrados. O modelo linear pode ser escrito como \[ Y=X\beta +\epsilon, \]

onde \(\epsilon\) são os erros. A melhor solução de mínimos quadrados é dada por \(\widehat{\beta}\): \[ \widehat{\beta}=\arg \min_\beta (y-X\beta)^\top (y-X\beta)=\arg \min_\beta \epsilon^\top\epsilon\cdot \]

Suponha que \(X^\top X\) seja de posto completo e, portanto, inversível. Minimizando a expressão acima com respecto a \(\beta\), obtemos \[ \widehat{\beta}=(X^\top X)^{-1} X^\top y\cdot \]

Os valores ajustados ou preditos são \(\widehat{y}=X(X^\top X)^{-1}X^\top y=P y\) sendo a projeção de \(y\) em \(C(X)\).

Os resíduos ordinário ou resíduos de mínimos quadrados são \[ e=y-\widehat{y}=y-X\widehat{\beta}=Q y = (\mbox{I}_n-P)y\cdot \] O vetor \(e\) é a projeção de \(y\) no complemento ortogonal de \(C(X)\).


Exemplo III.15

Voltemos ao exemplo dos pulôveres clássicos azuis. No Exemplo III.11, consideramos o ajuste de regressão das vendas \(X_1\) no preço \(X_2\) e concluímos que havia apenas uma pequena influência das vendas alterando os preços. Um modelo linear que incorpora as três variáveis permite aproximar as vendas em função linear do preço \(X_2\), propaganda \(X_3\) e presença de vendedores \(X_4\) simultaneamente. 


reg.pullover <- lm(Sales ~ Price+Advert+Asst.Hours, data= pullover)
summary(reg.pullover)
## 
## Call:
## lm(formula = Sales ~ Price + Advert + Asst.Hours, data = pullover)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.369  -9.406   1.599   5.151  19.729 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 65.66956   57.12507   1.150  0.29407   
## Price       -0.21578    0.32194  -0.670  0.52764   
## Advert       0.48519    0.08678   5.591  0.00139 **
## Asst.Hours   0.84373    0.37400   2.256  0.06491 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.7 on 6 degrees of freedom
## Multiple R-squared:  0.9067, Adjusted R-squared:  0.8601 
## F-statistic: 19.44 on 3 and 6 DF,  p-value: 0.001713
anova(reg.pullover)
## Analysis of Variance Table
## 
## Response: Sales
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Price       1  291.3   291.3  1.8061 0.2275654    
## Advert      1 8292.5  8292.5 51.4203 0.0003715 ***
## Asst.Hours  1  820.7   820.7  5.0893 0.0649088 .  
## Residuals   6  967.6   161.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ou

one.way <- aov(Sales ~ Price + Advert + Asst.Hours, data = pullover)
summary(one.way)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Price        1    291     291   1.806 0.227565    
## Advert       1   8292    8292  51.420 0.000371 ***
## Asst.Hours   1    821     821   5.089 0.064909 .  
## Residuals    6    968     161                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


O coeficiente de determinação é calculado como: \[ R^2 = 1-\dfrac{e^\top e}{\sum_{i=1}^n (y_i-\overline{y})^2}=0.9067\cdot \] Concluímos que a variação de \(X_1\) é bem aproximada pela relação linear.


O coeficiente de determinação é influenciado pelo número de regressores. Para um determinado tamanho de amostra \(n\), o valor de \(R^2\) aumentará adicionando mais regressores ao modelo linear. O valor de \(R^2\) pode, portanto, ser alto mesmo que sejam incluídos regressores possivelmente irrelevantes. Um coeficiente de determinação ajustado para \(p\) regressores e um intercepto constante, ou seja, com número de parâmetros \(p+1\) é \[ R^2_{adj}=R^2-\dfrac{p(1-R^2)}{n-(p+1)}\cdot \]


Exemplo III.16

O coeficiente de determinação corrigido para o Exemplo III.15 é \(R^2_{adj}=0.8601\).


Observe que o modelo linear é muito flexível e pode modelar relações não lineares entre a resposta \(Y\) e as variáveis explicativas \(X\). Por exemplo, uma relação quadrática em uma variável \(x\) pode ser incluída. Então \(y_i=\alpha+\beta_1 x_i+\beta_2x_i^2+\epsilon_i\) poderia ser escrito em notação matricial como \(y =X\beta+\epsilon\) onde \[ X=\begin{pmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ \vdots & \vdots & \vdots \\ 1 & x_n & x_n^2\end{pmatrix}\cdot \]


Propriedades de \(\widehat{\beta}\)


Quando \(y_i\) é a \(i\)-ésima observação de uma variável aleatória \(Y\), os erros também são aleatórios. Sob suposições padrão: independência, média zero e variância constante \(\sigma^2\). A inferência pode ser realizada em \(\widehat{\beta}\). Usando as propriedades: \[ \mbox{E}(\widehat{\beta})=\beta \] e \[ \mbox{Var}(\widehat{\beta})=\sigma^2 (X^\top X)^{-1}\cdot \]

O análogo do teste \(t\) para a situação de regressão linear multivariada é \[ t=\dfrac{\widehat{\beta}_j}{\mbox{SE}(\widehat{\beta}_j)}\cdot \]

O erro padrão de cada coeficient \(\widehat{\beta}_j\) é dado pela raiz quadrada dos elementos diagonais da matriz \(\mbox{Var}(\widehat{\beta})\). Em situações padrão, a variância \(\sigma^2\) do erro não é conhecida. Para o modelo linear com intercepto, pode-se estimá-la por \[ \widehat{\sigma}^2=\dfrac{1}{n-(p+1)}(y-\widehat{y})^\top (y-\widehat{y}), \] onde \(p\) é a dimensão de \(\beta\). Ao testar \(\beta_j= 0\) rejeitamos a hipótese no nível de significância \(\alpha\)˛ se \(|t|\geq t_{1-\alpha/2;n-(p+1)}\).


III.7. Boston Housing


As estatísticas apresentadas até agora podem ser calculadas para o conjunto de dados de habitação de Boston. Observando-se as médias e medianas confirma a assimetria dos componentes.


Boston.dataset = read.table("http://leg.ufpr.br/~lucambio/MSM/Bostondataset.txt", header = TRUE, sep = "")
Mediana = median(Boston.dataset$MEDV)
Boston.dataset$fMEDV <- factor(I(Boston.dataset$MEDV>Mediana), 
                               levels = c("TRUE","FALSE"), labels = c(1,0))
Boston.dataset[,16] = log(Boston.dataset[,1])
Boston.dataset[,17] = Boston.dataset[,2]/10
Boston.dataset[,18] = log(Boston.dataset[,3])
Boston.dataset[,19] = Boston.dataset[,4]
Boston.dataset[,20] = log(Boston.dataset[,5])
Boston.dataset[,21] = log(Boston.dataset[,6])
Boston.dataset[,22] = Boston.dataset[,7]^2.5/10000
Boston.dataset[,23] = log(Boston.dataset[,8])
Boston.dataset[,24] = log(Boston.dataset[,9])
Boston.dataset[,25] = log(Boston.dataset[,10])
Boston.dataset[,26] = exp(0.4*Boston.dataset[,11])/1000
Boston.dataset[,27] = Boston.dataset[,12]/100
Boston.dataset[,28] = sqrt(Boston.dataset[,13])
Boston.dataset[,29] = log(Boston.dataset[,14])
summary(Boston.dataset[,1:14])
##       CRIM                ZN             INDUS            CHAS        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       NOX               RM             AGE              DIS        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       RAD              TAX           PTRATIO            B         
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      LSTAT            MEDV      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00


ou


library(pastecs)
round(stat.desc(Boston.dataset[,1:14]), digits = 2) # arredondado para 2 decimais
##                 CRIM      ZN   INDUS   CHAS    NOX      RM      AGE     DIS
## nbr.val       506.00  506.00  506.00 506.00 506.00  506.00   506.00  506.00
## nbr.null        0.00  372.00    0.00 471.00   0.00    0.00     0.00    0.00
## nbr.na          0.00    0.00    0.00   0.00   0.00    0.00     0.00    0.00
## min             0.01    0.00    0.46   0.00   0.38    3.56     2.90    1.13
## max            88.98  100.00   27.74   1.00   0.87    8.78   100.00   12.13
## range          88.97  100.00   27.28   1.00   0.49    5.22    97.10   11.00
## sum          1828.44 5750.00 5635.21  35.00 280.68 3180.02 34698.90 1920.29
## median          0.26    0.00    9.69   0.00   0.54    6.21    77.50    3.21
## mean            3.61   11.36   11.14   0.07   0.55    6.28    68.57    3.80
## SE.mean         0.38    1.04    0.30   0.01   0.01    0.03     1.25    0.09
## CI.mean.0.95    0.75    2.04    0.60   0.02   0.01    0.06     2.46    0.18
## var            73.99  543.94   47.06   0.06   0.01    0.49   792.36    4.43
## std.dev         8.60   23.32    6.86   0.25   0.12    0.70    28.15    2.11
## coef.var        2.38    2.05    0.62   3.67   0.21    0.11     0.41    0.55
##                  RAD       TAX PTRATIO         B   LSTAT     MEDV
## nbr.val       506.00    506.00  506.00    506.00  506.00   506.00
## nbr.null        0.00      0.00    0.00      0.00    0.00     0.00
## nbr.na          0.00      0.00    0.00      0.00    0.00     0.00
## min             1.00    187.00   12.60      0.32    1.73     5.00
## max            24.00    711.00   22.00    396.90   37.97    50.00
## range          23.00    524.00    9.40    396.58   36.24    45.00
## sum          4832.00 206568.00 9338.50 180477.06 6402.45 11401.60
## median          5.00    330.00   19.05    391.44   11.36    21.20
## mean            9.55    408.24   18.46    356.67   12.65    22.53
## SE.mean         0.39      7.49    0.10      4.06    0.32     0.41
## CI.mean.0.95    0.76     14.72    0.19      7.97    0.62     0.80
## var            75.82  28404.76    4.69   8334.75   50.99    84.59
## std.dev         8.71    168.54    2.16     91.29    7.14     9.20
## coef.var        0.91      0.41    0.12      0.26    0.56     0.41

A matriz de covariância estimada é dada a seguir:

round(cor(Boston.dataset[,1:14]), digits = 2) # arredondado para 2 decimais
##          CRIM    ZN INDUS  CHAS   NOX    RM   AGE   DIS   RAD   TAX PTRATIO
## CRIM     1.00 -0.20  0.41 -0.06  0.42 -0.22  0.35 -0.38  0.63  0.58    0.29
## ZN      -0.20  1.00 -0.53 -0.04 -0.52  0.31 -0.57  0.66 -0.31 -0.31   -0.39
## INDUS    0.41 -0.53  1.00  0.06  0.76 -0.39  0.64 -0.71  0.60  0.72    0.38
## CHAS    -0.06 -0.04  0.06  1.00  0.09  0.09  0.09 -0.10 -0.01 -0.04   -0.12
## NOX      0.42 -0.52  0.76  0.09  1.00 -0.30  0.73 -0.77  0.61  0.67    0.19
## RM      -0.22  0.31 -0.39  0.09 -0.30  1.00 -0.24  0.21 -0.21 -0.29   -0.36
## AGE      0.35 -0.57  0.64  0.09  0.73 -0.24  1.00 -0.75  0.46  0.51    0.26
## DIS     -0.38  0.66 -0.71 -0.10 -0.77  0.21 -0.75  1.00 -0.49 -0.53   -0.23
## RAD      0.63 -0.31  0.60 -0.01  0.61 -0.21  0.46 -0.49  1.00  0.91    0.46
## TAX      0.58 -0.31  0.72 -0.04  0.67 -0.29  0.51 -0.53  0.91  1.00    0.46
## PTRATIO  0.29 -0.39  0.38 -0.12  0.19 -0.36  0.26 -0.23  0.46  0.46    1.00
## B       -0.39  0.18 -0.36  0.05 -0.38  0.13 -0.27  0.29 -0.44 -0.44   -0.18
## LSTAT    0.46 -0.41  0.60 -0.05  0.59 -0.61  0.60 -0.50  0.49  0.54    0.37
## MEDV    -0.39  0.36 -0.48  0.18 -0.43  0.70 -0.38  0.25 -0.38 -0.47   -0.51
##             B LSTAT  MEDV
## CRIM    -0.39  0.46 -0.39
## ZN       0.18 -0.41  0.36
## INDUS   -0.36  0.60 -0.48
## CHAS     0.05 -0.05  0.18
## NOX     -0.38  0.59 -0.43
## RM       0.13 -0.61  0.70
## AGE     -0.27  0.60 -0.38
## DIS      0.29 -0.50  0.25
## RAD     -0.44  0.49 -0.38
## TAX     -0.44  0.54 -0.47
## PTRATIO -0.18  0.37 -0.51
## B        1.00 -0.37  0.33
## LSTAT   -0.37  1.00 -0.74
## MEDV     0.33 -0.74  1.00

ou

library(corrplot)
par(cex = 0.6)
corrplot(cor(Boston.dataset[,1:14]), method = "number", type = "upper") # mostra apenas o lado superior

library(GGally)
par(cex = 0.6)
ggpairs(Boston.dataset[,c(1:4,14)])

par(cex = 0.6)
ggpairs(Boston.dataset[,c(5:9,14)])

par(cex = 0.6)
ggpairs(Boston.dataset[,c(10:14)])

Observemos a correlação entre \(X_{14}\) (MEDV), o valor da casa e todas as outras variáveis, dada pela última linha ou coluna acima. As correlações mais altas (em valores absolutos) estão em ordem decrescente \(X_{13}\) (LSTAT), \(X_6\) (RM), \(X_{11}\) (PTRATIO), \(X_3\) (INDUS), \(X_{10}\) (TAX), etc.

O uso da transformação \(Z\) de Fisher em cada uma das correlações entre \(X_{14}\) e as outras variáveis confirmaria que todas as correlações são significativamente diferentes de zero, exceto a correlação entre \(X_{14}\) e \(X_4\), a variável indicadora para o Rio Charles. Sabemos, no entanto, que a correlação e a transformação \(Z\) de Fisher não são apropriados para a variável binária. 

# testes de correlação para todo o conjunto de dados
library(Hmisc)
res <- rcorr(as.matrix(Boston.dataset[,1:14])) # rcorr() aceita apenas matrizes
# Exibindo os p-valores arredondados para 3 decimais
round(res$P, 3)
##          CRIM    ZN INDUS  CHAS  NOX    RM   AGE   DIS   RAD   TAX PTRATIO
## CRIM       NA 0.000 0.000 0.209 0.00 0.000 0.000 0.000 0.000 0.000   0.000
## ZN      0.000    NA 0.000 0.338 0.00 0.000 0.000 0.000 0.000 0.000   0.000
## INDUS   0.000 0.000    NA 0.157 0.00 0.000 0.000 0.000 0.000 0.000   0.000
## CHAS    0.209 0.338 0.157    NA 0.04 0.040 0.052 0.026 0.869 0.424   0.006
## NOX     0.000 0.000 0.000 0.040   NA 0.000 0.000 0.000 0.000 0.000   0.000
## RM      0.000 0.000 0.000 0.040 0.00    NA 0.000 0.000 0.000 0.000   0.000
## AGE     0.000 0.000 0.000 0.052 0.00 0.000    NA 0.000 0.000 0.000   0.000
## DIS     0.000 0.000 0.000 0.026 0.00 0.000 0.000    NA 0.000 0.000   0.000
## RAD     0.000 0.000 0.000 0.869 0.00 0.000 0.000 0.000    NA 0.000   0.000
## TAX     0.000 0.000 0.000 0.424 0.00 0.000 0.000 0.000 0.000    NA   0.000
## PTRATIO 0.000 0.000 0.000 0.006 0.00 0.000 0.000 0.000 0.000 0.000      NA
## B       0.000 0.000 0.000 0.273 0.00 0.004 0.000 0.000 0.000 0.000   0.000
## LSTAT   0.000 0.000 0.000 0.226 0.00 0.000 0.000 0.000 0.000 0.000   0.000
## MEDV    0.000 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.000   0.000
##             B LSTAT MEDV
## CRIM    0.000 0.000    0
## ZN      0.000 0.000    0
## INDUS   0.000 0.000    0
## CHAS    0.273 0.226    0
## NOX     0.000 0.000    0
## RM      0.004 0.000    0
## AGE     0.000 0.000    0
## DIS     0.000 0.000    0
## RAD     0.000 0.000    0
## TAX     0.000 0.000    0
## PTRATIO 0.000 0.000    0
## B          NA 0.000    0
## LSTAT   0.000    NA    0
## MEDV    0.000 0.000   NA


As mesmas estatísticas descritivas podem ser calculadas para as variáveis transformadas. Como pode ser visto, a maioria das variáveis agora é mais simétrica. Observe que as covariâncias e as correlações são sensíveis a essas transformações não lineares. 

library(corrplot)
par(cex = 0.6)
corrplot(cor(Boston.dataset[,16:29]), method = "number", type = "upper") # mostra apenas o lado superior


res1 <- rcorr(as.matrix(Boston.dataset[,16:29])) 
round(res1$P, 3)
##       V16   V17  V18   V19   V20   V21   V22  V23   V24   V25   V26   V27   V28
## V16    NA 0.000 0.00 0.522 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.000
## V17 0.000    NA 0.00 0.338 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.000
## V18 0.000 0.000   NA 0.070 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.000
## V19 0.522 0.338 0.07    NA 0.062 0.059 0.098 0.05 0.773 0.404 0.002 0.273 0.164
## V20 0.000 0.000 0.00 0.062    NA 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.000
## V21 0.000 0.000 0.00 0.059 0.000    NA 0.000 0.00 0.000 0.000 0.000 0.003 0.000
## V22 0.000 0.000 0.00 0.098 0.000 0.000    NA 0.00 0.000 0.000 0.000 0.000 0.000
## V23 0.000 0.000 0.00 0.050 0.000 0.000 0.000   NA 0.000 0.000 0.000 0.000 0.000
## V24 0.000 0.000 0.00 0.773 0.000 0.000 0.000 0.00    NA 0.000 0.000 0.000 0.000
## V25 0.000 0.000 0.00 0.404 0.000 0.000 0.000 0.00 0.000    NA 0.000 0.000 0.000
## V26 0.000 0.000 0.00 0.002 0.000 0.000 0.000 0.00 0.000 0.000    NA 0.000 0.000
## V27 0.000 0.000 0.00 0.273 0.000 0.003 0.000 0.00 0.000 0.000 0.000    NA 0.000
## V28 0.000 0.000 0.00 0.164 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000    NA
## V29 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000 0.000
##     V29
## V16   0
## V17   0
## V18   0
## V19   0
## V20   0
## V21   0
## V22   0
## V23   0
## V24   0
## V25   0
## V26   0
## V27   0
## V28   0
## V29  NA


Observe que algumas das correlações entre \(\widetilde{X}_{14}\) e as outras variáveis aumentaram. Se queremos explicar as variações do preço \(\widetilde{X}_{14}\) pela variação de todas as outras variáveis \(\widetilde{X}_1,\cdots,\widetilde{X}_{13}\), poderíamos estimar o modelo linear: \[ \widetilde{X}_{14}=\beta_0+\sum_{j=1}^{13} \beta_j \widetilde{X}_j + \epsilon\cdot \] 

modelo.transformado <- lm(Boston.dataset[,29] ~ Boston.dataset[,16]+Boston.dataset[,17]+
                            Boston.dataset[,18]+Boston.dataset[,19]+Boston.dataset[,20]+
                            Boston.dataset[,21]+Boston.dataset[,22]+Boston.dataset[,23]+
                            Boston.dataset[,24]+Boston.dataset[,25]+Boston.dataset[,26]+
                            Boston.dataset[,27]+Boston.dataset[,28])
summary(modelo.transformado)
## 
## Call:
## lm(formula = Boston.dataset[, 29] ~ Boston.dataset[, 16] + Boston.dataset[, 
##     17] + Boston.dataset[, 18] + Boston.dataset[, 19] + Boston.dataset[, 
##     20] + Boston.dataset[, 21] + Boston.dataset[, 22] + Boston.dataset[, 
##     23] + Boston.dataset[, 24] + Boston.dataset[, 25] + Boston.dataset[, 
##     26] + Boston.dataset[, 27] + Boston.dataset[, 28])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9918 -0.1002 -0.0034  0.1117  0.7640 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.176874   0.379017  11.020  < 2e-16 ***
## Boston.dataset[, 16] -0.014606   0.011650  -1.254 0.210527    
## Boston.dataset[, 17]  0.001392   0.005639   0.247 0.805121    
## Boston.dataset[, 18] -0.012709   0.022312  -0.570 0.569195    
## Boston.dataset[, 19]  0.109980   0.036634   3.002 0.002817 ** 
## Boston.dataset[, 20] -0.283112   0.105340  -2.688 0.007441 ** 
## Boston.dataset[, 21]  0.421108   0.110175   3.822 0.000149 ***
## Boston.dataset[, 22]  0.006403   0.004863   1.317 0.188536    
## Boston.dataset[, 23] -0.183154   0.036804  -4.977 8.97e-07 ***
## Boston.dataset[, 24]  0.068362   0.022473   3.042 0.002476 ** 
## Boston.dataset[, 25] -0.201832   0.048432  -4.167 3.64e-05 ***
## Boston.dataset[, 26] -0.040017   0.008091  -4.946 1.04e-06 ***
## Boston.dataset[, 27]  0.044472   0.011456   3.882 0.000118 ***
## Boston.dataset[, 28] -0.262615   0.016091 -16.320  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2008 on 492 degrees of freedom
## Multiple R-squared:  0.765,  Adjusted R-squared:  0.7588 
## F-statistic: 123.2 on 13 and 492 DF,  p-value: < 2.2e-16


O valor de \(R^2=0.765\) e \(R^2_{adj}=0.759\) mostra que a maior parte da variação de \(\widetilde{X}_{14}\) é explicada pelo modelo linear acima. Novamente, vemos que as variações do \(\widetilde{X}_{14}\) sejam explicadas principalmente por: \(\widetilde{X}_{13}\), \(\widetilde{X}_{8}\), \(\widetilde{X}_{11}\), \(\widetilde{X}_{10}\), \(\widetilde{X}_{12}\), \(\widetilde{X}_{6}\), \(\widetilde{X}_{9}\), \(\widetilde{X}_{4}\) e \(\widetilde{X}_{5}\), em ordem decrescente do valor absoluto da estatística \(t\). As outras variáveis \(\widetilde{X}_{1}\), \(\widetilde{X}_{2}\), \(\widetilde{X}_{3}\) e \(\widetilde{X}_{7}\) parecem ter pouca influência nas variações do \(\widetilde{X}_{14}\).


III.8. Exercícios


  1. A covariância amostral \(s_{X_4 X_5}\) entre \(X_4\) e \(X_5\) para o conjunto de todos os dados bancários é positiva. Dadas as definições de \(X_4\) e \(X_5\), esperaríamos uma covariância negativa. Usando uma figura adequada Você pode explicar por que o \(s_{X_4 X_5}\) é positiva?

  2. Considere as duas sub-nuvens de notas bancárias falsificadas e genuínas separadamente. Você ainda espera que a \(s_{X_4 X_5}\), agora calculado separadamente para cada nuvem, seja positiva?

  3. Observamos que, para duas variáveis aleatórias normais, a covariância zero implica independência. Por que essa observação não se aplica ao Exemplo III.4?

  4. Calcule a covariância entre as variáveis
    \(X_2\) = milhas por galão,
    \(X_8\) = Peso do conjunto de dados dos carros (car). Que sinal você espera que a covariância tenha?

  5. Calcule a matriz de correlação das variáveis no Exemplo III.2. Comente o sinal das correlações e teste a hipótese \(\rho_{X_1 X_2}=0\).

  6. Suponha que você tenha observado um conjunto de observações \(\{x_i\}_{i=1}^n\) com \(\overline{x}=0\), \(s_{XX}=1\) e \(\frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})^2\). Defina a variável \(y_i=x_i^2\). Você pode dizer imediatamente se \(r_{XY}\neq 0\)?

  7. Quantas vendas o gerente têxtil espera com um preço de pulôver de \(x=105\)?

  8. Como é um gráfico de dispersão de duas variáveis aleatórias para \(R^2=1\) e \(R^2=0\)?

  9. Mostre que o coeficiente de determinação é o quadrado da correlação simples entre \(X\) e \(Y\).

  10. Em que circunstâncias você obteria os mesmos coeficientes das linhas de regressão linear de \(Y\) em \(x\) e de \(X\) em \(Y\)?

  11. Trate o desenho do planejamento do Exemplo III.14 como se houvesse trinta lojas e não dez. Defina \(x_i\) como o índice da loja, ou seja, \(x_i= i\), \(i=1,2,\cdots,30\). A hipótese nula é uma linha de regressão constante, \(\mbox{E}(Y)=\mu\). Como é a curva de regressão alternativa?

  12. Realize o teste no Exercício III.13 para o exemplo da loja com um nível de significância 0.99. Você ainda rejeita a hipótese de estratégias de marketing iguais?

  13. Calcule um intervalo de confiança aproximado para \(\rho_{X_1 X_4}\) no Exemplo III.2. Para isso, comece de um intervalo de confiança para \(\tanh^{-1}(\rho_{X_1 X_4})\) e depois aplique a transformação inversa.

  14. No Exemplo III.2, usando a taxa de câmbio de 1 Eur = 106 JPY, calcule a mesma covariância empírica usando preços em ienes japoneses e não em euros. Existe uma diferença significativa? Por quê?

  15. Por que a correlação tem o mesmo sinal que a covariância?

  16. Mostre que \(\mbox{posto}(H)=\mbox{tr}(H)=n-1\).

  17. Calcule nos dados de pulôveres a regressão de \(X_1\) em \(X_2, X_3\) e de \(X_1\) em \(X_2, X_4\). Qual delas tem o melhor coeficiente de determinação?

  18. Considere o problema da ANOVA (Seção III.5) novamente. Estabeleça a matriz de restrição \(A\) para testar \(H_0 \, : \, \mu_1=\mu_2\). Teste esta hipótese.