Bastão de Asclépio & Distribuição Normal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
options(warn=-1)
suppressMessages(library(aplpack, warn.conflicts=FALSE))
suppressMessages(library(calculus, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(carData, warn.conflicts=FALSE))
suppressMessages(library(cellWise, warn.conflicts=FALSE))
suppressMessages(library(cluster, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(diptest, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(ellipse, warn.conflicts=FALSE))
suppressMessages(library(emmeans, warn.conflicts=FALSE))
suppressMessages(library(far, warn.conflicts=FALSE))
suppressMessages(library(fmsb, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggfortify, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(HSAUR2, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lawstat, warn.conflicts=FALSE))
suppressMessages(library(lfda, warn.conflicts=FALSE))
suppressMessages(library(lmerTest, warn.conflicts=FALSE))
suppressMessages(library(MANOVA.RM, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(MatchIt, warn.conflicts=FALSE)) # distancia robusta
suppressMessages(library(matlib, warn.conflicts=FALSE))
suppressMessages(library(matrixcalc, warn.conflicts=FALSE))
suppressMessages(library(MatrixModels, warn.conflicts=FALSE))
suppressMessages(library(matrixStats, warn.conflicts=FALSE))
suppressMessages(library(missMethods, warn.conflicts=FALSE))
suppressMessages(library(MomTrunc, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(MuMIn, warn.conflicts=FALSE))
suppressMessages(library(MVA, warn.conflicts=FALSE)) # An Introduction to Applied Multivariate Analysis with R
suppressMessages(library(MVar.pt, warn.conflicts=FALSE))
suppressMessages(library(mvdalab, warn.conflicts=FALSE))
# suppressMessages(library(MVLM, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(MVTests, warn.conflicts=FALSE)) # T^2 robusta
suppressMessages(library(mvtnorm, warn.conflicts=FALSE))
suppressMessages(library(norm, warn.conflicts=FALSE))
suppressMessages(library(phia, warn.conflicts=FALSE))
suppressMessages(library(plotly, warn.conflicts=FALSE))
suppressMessages(library(pracma, warn.conflicts=FALSE))
suppressMessages(library(profileR, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(qcc, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(reticulate, warn.conflicts=FALSE))
suppressMessages(library(rgl, warn.conflicts=FALSE))
suppressMessages(library(rrcov, warn.conflicts=FALSE))
suppressMessages(library(scatterplot3d, warn.conflicts=FALSE))
suppressMessages(library(SHT, warn.conflicts=FALSE))
suppressMessages(library(sjPlot, warn.conflicts=FALSE))
suppressMessages(library(skimr, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
source("summarySEwithin2.R")
fit.lm <- lm(cbind(PC1,PC2)~Breed, data = Dados.org) print(anova.fit <- car::Anova(fit.lm)) print(effectsize::eta_squared(anova.fit, partial = TRUE, ci = 1 - alpha, alternative = "two"))
alphaBonf <- alpha/length(unique(Dados.org$Breed)) fit.lm <- lm(cbind(PC1,PC2)~Breed, data = subset(Dados.org, Breed!="Hereford")) print(anova.fit <- car::Anova(fit.lm)) T2H.fit <- mvdalab::MVComp(subset(Dados.org, Breed=="Angus", select=c(PC1,PC2)), subset(Dados.org, Breed=="Simental", select=c(PC1,PC2)), level = 1-alphaBonf) print(T2H.fit) plot(T2H.fit, Diff2Plot = c(1, 2), include.zero = TRUE) fit.lm <- lm(cbind(PC1,PC2)~Breed, data = subset(Dados.org, Breed!="Angus")) print(anova.fit <- car::Anova(fit.lm)) T2H.fit <- mvdalab::MVComp(subset(Dados.org, Breed=="Hereford", select=c(PC1,PC2)), subset(Dados.org, Breed=="Simental", select=c(PC1,PC2)), level = 1-alphaBonf) print(T2H.fit) plot(T2H.fit, Diff2Plot = c(1, 2), include.zero = TRUE) fit.lm <- lm(cbind(PC1,PC2)~Breed, data = subset(Dados.org, Breed!="Simental")) print(anova.fit <- car::Anova(fit.lm)) T2H.fit <- mvdalab::MVComp(subset(Dados.org, Breed=="Hereford", select=c(PC1,PC2)), subset(Dados.org, Breed=="Angus", select=c(PC1,PC2)), level = 1-alphaBonf) print(T2H.fit) plot(T2H.fit, Diff2Plot = c(1, 2), include.zero = TRUE)
Dados.Mm <- as.matrix(subset(Dados, Sexo=="Masculino", select=c(Estatura, MCT))) Dados.Fm <- as.matrix(subset(Dados, Sexo=="Feminino", select=c(Estatura, MCT))) colMeans(Dados.Mm) colMeans(Dados.Fm) SHT::mean2.1980Johansen(Dados.Mm, Dados.Fm)
mnv <- rrcov::Wilks.test(x=Dados[,c("Estatura","MCT")], grouping=Dados[,"Sexo"], method="mcd",nrep=1e3) print(mnv) result <- MVN::mvn(data = Masc, # R=5e4, # mvnTest = "energy", mvnTest = "hz", univariateTest = "SF", showOutliers=TRUE) print(result$multivariateNormality) print(result$univariateNormality) print(result$multivariateOutliers) print(res <- heplots::boxM(Q6.39[,1:2], Q6.39[,3])) plot(res, main="Q6.39") set.seed(123) dados <- data.frame( fator1 = factor(rep(c("A", "B"), each = 30)), fator2 = factor(rep(c("X", "Y"), times = 30)), medida1 = rnorm(60, mean = 50, sd = 10), medida2 = rnorm(60, mean = 60, sd = 15), medida3 = rnorm(60, mean = 55, sd = 20) ) fit <- manova(cbind(X1,X2)~F1*F2, data=Q6.14b) print(sumario <- summary(fit, test="Pillai")) print(sumario$SS, digits=2) anv <- car::Anova(fit, test="Pillai") eta2 <- effectsize::eta_squared(anv, partial=TRUE, generalized=FALSE, ci=1-alfa/3, alternative="two.sided", verbose=TRUE) eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2) print(eta2, digits=2)
print(a.hat <- solve(Spooled) %*% (centroidif - mu_0), digits = 2) print(a.hat, digits=3)
lda_result <- MASS::lda(Sex ~ Length_ln + Width_ln + Height_ln, data=Q6.18.ln) LDAcoef <- lda_result$scaling print(LDAcoef, digits=3)
print(a.hat/LDAcoef, digits=3) library(reshape2) library(ggplot2) data(iris) Iris.long <- reshape2::melt(iris, id=c('Species'), measured=c('Sepal.Length', 'Sepal.Width', 'Petal.Length', 'Petal.Width')) names(Iris.long) <- c('Group', 'Measure', 'Frequency') IrisBoxplot <- ggplot2::ggplot(Iris.long, ggplot2::aes(Group, Frequency, color = Measure)) + ggplot2::geom_boxplot() + ggplot2::labs(x='Species', y='Size (cm)', color='Measure') plot(IrisBoxplot) rm(Iris.long) setosa_vs_versicolor <- c(1, -1, 0) contrasts(iris$Species) <- setosa_vs_versicolor fit <- lm(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data = iris) summary(fit)
RPubs
“Não é paradoxo dizer que nos nossos momentos de inspiração mais teórica podemos estar o mais próximo possível de nossas aplicações mais práticas.”
WHITEHEAD, AN apud BOYER, CB (1974) História da matemática. São Paulo: Blücher/EDUSP, p. 419.
Análise estatística multivariada tem que produzir valor p para o teste de hipótese nula omnibus.
As ideias desenvolvidas no Capítulo 5 podem ser estendidas para lidar com problemas envolvendo a comparação de vários vetores de média. A teoria é um pouco mais complexa e baseia-se na suposição de distribuições normais multivariadas ou tamanhos de amostra grandes. Da mesma forma, a notação se torna mais pesada. Para contornar esses problemas, frequentemente revisaremos procedimentos univariados para comparar várias médias e depois generalizaremos para os casos multivariados correspondentes por analogia. Os exemplos numéricos que apresentamos ajudarão a consolidar os conceitos.
Como as comparações de médias frequentemente (e deveriam) emanam de experimentos planejados, aproveitamos a oportunidade para discutir alguns dos princípios de uma boa prática experimental. Um delineamento de medidas repetidas (intraparticipantes), útil em estudos comportamentais, é explicitamente considerado, juntamente com modificações necessárias para analisar curvas de crescimento.
Começamos considerando pares de vetores de média. Nas seções posteriores, discutimos várias comparações entre vetores de média organizados de acordo com os níveis de tratamento. As estatísticas de teste correspondentes dependem de uma divisão da variação total em parcelas de variação atribuíveis às fontes de tratamento e erro. Essa divisão é conhecida como análise de variância multivariada (MANOVA).
Uma maneira de classificar os delineamentos é por meio das seguintes características:
\(n\): número de unidades experimentais
\(p\): número de medidas (measures)
\(q\): número de condições experimentais dependentes
\(p\times q\): número de variáveis dependentes
\(g\): número de condições experimentais independentes
UE: unidade experimental
Os arquivos de dados estão, neste texto, sempre no formato wide, mesmo no delineamento de medidas repetidas.
Medições são frequentemente registradas sob diferentes conjuntos de condições experimentais para verificar se as respostas diferem significativamente entre esses conjuntos. Por exemplo, a eficácia de um novo medicamento ou de uma campanha publicitária saturada pode ser determinada ao comparar medições antes do “tratamento” (medicamento ou publicidade) com aquelas após o tratamento. Em outras situações, dois ou mais tratamentos podem ser administrados nas mesmas ou em unidades experimentais semelhantes, e as respostas podem ser comparadas para avaliar os efeitos dos tratamentos.
Uma abordagem racional para comparar dois tratamentos, ou a presença e ausência de um único tratamento, é atribuir ambos os tratamentos à mesma unidade ou unidades idênticas (indivíduos, lojas, lotes de terra, e assim por diante). As respostas em pares podem então ser analisadas calculando suas diferenças, eliminando assim grande parte da variação extrínseca de unidade para unidade.
No caso de uma única resposta (univariada), seja \(X_{i1}\) denotar a resposta ao tratamento 1 (ou a resposta antes do tratamento) e \(X_{i2}\) denotar a resposta ao tratamento 2 (ou a resposta após o tratamento) para a i-ésima tentativa. Ou seja, \((X_{i1},X_{i2})\) são medições registradas na i-ésima unidade ou i-ésimo par de unidades semelhantes. Por delineamento, as \(n\) diferenças
\[ D_i = X_{i1} - X_{i2}, \quad i = 1,2,...,n \tag{6-1} \]
deveriam refletir apenas os efeitos diferenciais dos tratamentos.
Dado que as diferenças \(D_i\) representam observações independentes de uma distribuição \(\mathcal{N}(\delta, \sigma^2_D)\), a variável
\[ t = \dfrac{\bar{D}-\delta}{S_D / \sqrt{n}} \sim t_{n-1} \tag{6-2} \]
em que
\[ \bar{D} = \dfrac{1}{n} \sum_{i=1}^{n} D_i \quad \text{e} \quad S_D^2 = \dfrac{1}{n-1} \sum_{i=1}^{n} (D_i - \bar{D})^2 \tag{6-3} \]
Consequentemente, um teste de nível \(\alpha\) de
\[ H_0: \delta = 0 \quad \text{vs.} \quad H_1: \delta \neq 0 \]
pode ser realizado comparando-se \(|t|\) com \(t_{n-1}(1-\alpha/2)\). Um intervalo de confiança de \(100(1 - \alpha)\%\) para a diferença média \(\delta = \mathbb{E}(X_{i1}-X_{i2})= \mathbb{E}(X_{i1})-\mathbb{E}(X_{i2})\) é fornecido pelo intervalo de confiança:
\[ \text{IC}^{1-\alpha}(\delta)=\left[\bar{D} \pm t_{n-1}(1-\alpha/2) \dfrac{S_D}{\sqrt{n}}\right]\tag{6-4} \]
Uma notação adicional é necessária para a extensão multivariada do procedimento de comparação pareada. É necessário distinguir entre \(p\) respostas, dois tratamentos e \(n\) unidades experimentais. Rotulamos as \(p\) respostas dentro da \(j\)-ésima unidade observacional como:
\[ X_{1j1} = \text{variável 1 sob tratamento 1} \\ X_{1j2} = \text{variável 2 sob tratamento 1} \\ \vdots \\ X_{1jp} = \text{variável } p \text{ sob tratamento 1} \]
\[ \text{---------------} \]
\[ X_{2j1} = \text{variável 1 sob tratamento 2} \\ X_{2j2} = \text{variável 2 sob tratamento 2} \\ \vdots \\ X_{2jp} = \text{variável } p \text{ sob tratamento 2} \]
e as \(p\) variáveis aleatórias de diferença pareada tornam-se:
\[ \begin{align} D_{j1} &= X_{1j1} - X_{2j1} \\ D_{j2} &= X_{1j2} - X_{2j2} \\ \vdots \\ D_{jp} &= X_{1jp} - X_{2jp} \end{align} \tag{6-5} \]
Seja \(\mathbf{D}^{\prime}_j = [D_{j1}\; D_{j2}\; \cdots\; D_{jp}]\) e suponha que, para \(j = 1, 2, \ldots, n\):
\[ \mathbb{E}(\mathbf{D}_j) = \boldsymbol{\delta} = \begin{bmatrix} \delta_1 \\ \delta_2 \\ \vdots\\ \delta_p \end{bmatrix} \]
e
\[ \mathbb{C}(\mathbf{D}_j) = \mathbf{\Sigma}_\mathbf{D} \tag{6-6} \]
Se, além disso, \(\{\mathbf{D}_i\}_{i=1}^{n} \sim \mathcal{N}_p\text{IID}(\boldsymbol{\delta}, \mathbf{\Sigma_\mathbf{D}})\), inferências sobre o vetor das diferenças médias \(\boldsymbol{\delta}\) podem ser baseadas em uma estatística \(T^2\). Especificamente,
\[ T^2 = n(\bar{\mathbf{D}} - \boldsymbol{\delta}_0)^{\prime} \mathbf{S}_\mathbf{D}^{-1} (\bar{\mathbf{D}} - \boldsymbol{\delta}_0) \tag{6-7} \]
em que
\[ \bar{\mathbf{D}} = \dfrac{1}{n} \sum_{i=1}^{n} \mathbf{D}_i \quad \text{e} \quad \mathbf{S}_\mathbf{D} = \frac{1}{n} \sum_{i=1}^{n} (\mathbf{D}_i - \bar{\mathbf{D}})(\mathbf{D}_i - \bar{\mathbf{D}})^{\prime} \tag{6-8} \]
Resultado 6.1. Sejam as diferenças \(\{\mathbf{D}_i\}_{i=1}^{n} \sim \mathcal{N}_p\text{IID}(\boldsymbol{\delta}, \mathbf{\Sigma_\mathbf{D}})\). Então
\[ T^2 = n(\bar{\mathbf{D}} - \boldsymbol{\delta}_0)^{\prime} \mathbf{S}_\mathbf{D}^{-1} (\bar{\mathbf{D}} - \boldsymbol{\delta}_0) \sim \dfrac{(n - 1)p}{n - p}F_{p, n - p}=T^2_{p,n-p} \]
Se \(n - p \gg30\), \(T^2\) segue aproximadamente uma distribuição \(\chi^2_p\), independentemente da forma da população subjacente das diferenças.
Prova. A distribuição exata de \(T^2\) é uma reformulação do resumo em (5-6), com vetores de diferenças para os vetores de observação. A distribuição aproximada de \(T^2\), para \(n - p \gg30\), segue de (4-28).
\[\Diamond\]
A condição \(\boldsymbol{\delta} = \boldsymbol{\delta}_0=\mathbf{0}\) é equivalente a “nenhuma diferença média entre os dois tratamentos”. Para a i-ésima variável, \(\delta_i > 0\) implica que o tratamento 1 é maior, em média, do que o tratamento 2. Em geral, inferências sobre \(\delta\) podem ser feitas usando o Resultado 6.1.
Dadas as diferenças observadas \(\mathbf{d}_j = [d_{j1}, d_{j2}, \ldots, d_{jp}]\), \(j = 1, 2, \ldots, n\), correspondendo às variáveis aleatórias em (6-5), um teste de nível \(\alpha\) de \(H_0: \boldsymbol{\delta} = \mathbf{0}\) versus \(H_1: \boldsymbol{\delta} \ne \mathbf{0}\) para uma população \(\mathcal{N}_p(\boldsymbol{\delta}, \mathbf{\Sigma_\mathbf{D}})\) rejeita \(H_0\) se
\[ T^2 = n\bar{\mathbf{d}}^{\prime} \mathbf{s}_\mathbf{d}^{-1} \bar{\mathbf{d}} > \dfrac{(n - 1)p}{n - p} F_{p, n-p}(1-\alpha) \]
Uma região de confiança de \(100(1 - \alpha)\%\) para \(\boldsymbol{\delta}\) consiste em todos os \(\boldsymbol{\delta}\) tais que:
\[ (\bar{\mathbf{d}}-\boldsymbol{\delta})^{\prime} \mathbf{s}_\mathbf{d}^{-1} (\bar{\mathbf{d}}-\boldsymbol{\delta}) \le \dfrac{(n - 1)p}{n(n - p)} F_{p, n-p}(1-\alpha) \tag{6-9} \]
Além disso, intervalos de confiança simultâneos de \(100(1 - \alpha)\%\) para as diferenças médias individuais \(\delta_i\) são dados por:
\[ \text{IC}^{1-\alpha}(\delta_i)=\left[\bar{d}_i \pm \sqrt{ \dfrac{(n-1)p}{n-p}F_{p,n-p}(1-\alpha)}\sqrt{\dfrac{s_{d_i}^2}{n}}\right]\tag{6-10} \]
Para \(n - p \gg 30\), \(\dfrac{(n - 1)p}{n - p} F_{p,n-p}(1-\alpha) \approx \chi^2_p(1-\alpha)\) e a multinormalidade não precisa ser assumida.
Os intervalos de confiança simultâneos t de Bonferroni de \(100(1 - \alpha)\%\) para as diferenças médias individuais são:
\[ \text{IC}^{1-\alpha}(\delta_i)=\left[\bar{d}_i \pm t_{n-1}\left(1-\dfrac{\alpha}{2p}\right)\sqrt{\dfrac{s_{d_i}^2}{n}}\right]\tag{6-10a} \]
\(p=2, q=2, g=1, n=11\)
UE: amostra retirada de descarregamento num rio ou afluente, homogeneizada e dividida em duas partes iguais; as duas partes são atribuídas aleatoriamente para cada laboratório
Estações de tratamento de águas residuais municipais são obrigadas por lei a monitorar suas descargas em rios e córregos regularmente. Preocupações sobre a confiabilidade dos dados de um desses programas de auto-monitoramento levaram a um estudo no qual amostras de efluentes foram divididas e enviadas para dois laboratórios para teste. Metade de cada amostra foi enviada para o Laboratório Estadual de Higiene de Wisconsin, e a outra metade foi enviada para um laboratório comercial privado rotineiramente usado no programa de monitoramento. Medidas da demanda bioquímica de oxigênio (DBO) e sólidos suspensos (SS) foram obtidas, para \(n = 11\) divisões de amostras, dos dois laboratórios. Os dados estão exibidos na Tabela 6.1.
Tabela 6.1 Dados do Efluente
Amostra \(j\) | Laboratório Comercial (DBO) \(x_{1j1}\) | Laboratório Comercial (SS) \(x_{1j2}\) | Laboratório Estadual de Higiene (DBO) \(x_{2j1}\) | Laboratório Estadual de Higiene (SS) \(x_{2j2}\) |
---|---|---|---|---|
1 | 6 | 27 | 25 | 15 |
2 | 6 | 23 | 28 | 13 |
3 | 18 | 64 | 36 | 22 |
4 | 8 | 44 | 35 | 29 |
5 | 11 | 30 | 15 | 31 |
6 | 34 | 75 | 44 | 64 |
7 | 28 | 26 | 42 | 30 |
8 | 71 | 124 | 54 | 64 |
9 | 43 | 54 | 34 | 56 |
10 | 33 | 30 | 29 | 20 |
11 | 20 | 14 | 39 | 21 |
Os resultados das análises químicas dos dois laboratórios concordam? Se existem diferenças, qual é a sua natureza?
A estatística \(T^2\) para testar \(H_0: \boldsymbol{\delta}^{\prime} = [\delta_1 \;\delta_2] = [0\;0]\) é construída a partir das diferenças de observações pareadas:
\[ \begin{align} d_{j1} &= x_{1j1} - x_{2j1}: -19, -22, -18, -27, -4, -10, -14, 17, 9, 4, -19\\ d_{j2} &= x_{1j2} - x_{2j2}: 12, 10, 42, 15, -1, 11, -4, 60, -2, 10, -7 \end{align} \]
\[ \bar{\mathbf{d}}= \begin{bmatrix} \bar{d}_1\\ \bar{d}_2 \end{bmatrix}= \begin{bmatrix} -9.36\\ 13.27 \end{bmatrix} \]
\[ \mathbf{s}_\mathbf{d}= \begin{bmatrix} 199.26&88.31\\ 88.31&418.61 \end{bmatrix} \]
\[ T^2=\begin{bmatrix}-9.36&13.27\end{bmatrix} \begin{bmatrix}0.0055&-0.0012\\ -0.0012&0.0026\end{bmatrix} \begin{bmatrix}-9.36\\13.27\end{bmatrix}=13.64 \]
Adotando \(\alpha = 0.05\), encontramos que
\[ \dfrac{p(n - 1)}{n - p} F_{p,n-p}(0.95) = \dfrac{2\times10}{9}F_{2,9}(0.95) = 9.46 \]
Como \(T^2 = 13.6 > 9.47\), rejeitamos \(H_0\) e concluímos que existe uma diferença média não nula entre as medições dos dois laboratórios. Ao inspecionar os dados, parece que o laboratório comercial tende a produzir medições de DBO mais baixas e medições de SS mais altas do que o Laboratório Estadual de Higiene. Os intervalos de confiança simultâneos de 95% para as diferenças médias \(\delta_1\) e \(\delta_2\) podem ser calculados usando a equação (6-10). Estes intervalos são:
\[ \text{IC}^{95\%}(\delta_1)=\left[ -9.36 \pm \sqrt{9.47}\sqrt{\dfrac{199.26}{11}} \right]=[-22.45, 3.73] \]
\[ \text{IC}^{95\%}(\delta_2)=\left[ 13.27 \pm \sqrt{9.47}\sqrt{\dfrac{418.61}{11}} \right]=[-5.70, 32.25] \]
Os intervalos de confiança simultâneos de 95% incluem zero, mas a hipótese \(H_0: \boldsymbol{\delta} = \mathbf{0}\) foi rejeitada no nível de 5%. O que devemos concluir?
As evidências apontam para diferenças reais. O ponto \(\boldsymbol{\delta} = \mathbf{0}\) está fora da região de confiança de 95% para \(\boldsymbol{\delta}\) (veja o Exercício 6.1), e esse resultado é consistente com o teste \(T^2\). O nível de confiança simultâneo de 95% se aplica ao conjunto inteiro de intervalos que poderiam ser construídos para todas as combinações lineares possíveis da forma \(a_1 \delta_1 + a_2 \delta_2\). Os intervalos particulares correspondentes às escolhas \((a_1 = 1, a_2 = 0)\) e \((a_1 = 0, a_2 = 1)\) contêm zero. Outras escolhas de \(a_1\) e \(a_2\) produzirão intervalos simultâneos que não contêm zero. (Se a hipótese \(H_0: \boldsymbol{\delta} = \mathbf{0}\) não fosse rejeitada, então todos os intervalos simultâneos incluiriam zero.)
Os intervalos simultâneos de Bonferroni também abrangem zero. (Veja o Exercício 6.2.)
Nossa análise assumiu uma distribuição normal para os \(\mathbf{D}_j\). Na verdade, a situação é ainda mais complexa pela presença de um ou, possivelmente, dois valores discrepantes (outliers). (Veja o Exercício 6.3.) Esses dados podem ser transformados para se tornarem mais próximos de uma distribuição normal, mas com uma amostra tão pequena, é difícil remover os efeitos dos outliers. (Veja o Exercício 6.4.)
Os resultados numéricos deste exemplo ilustram uma circunstância incomum que pode ocorrer ao fazer inferências.
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
x <- read.table("JW6Data/T6-1.dat", quote="\"", comment.char="")
names(x) <- c("BOD_C", "SS_C", "BOD_S", "SS_S")
print.data.frame(x)
BOD_C SS_C BOD_S SS_S
1 6 27 25 15
2 6 23 28 13
3 18 64 36 22
4 8 44 35 29
5 11 30 15 31
6 34 75 44 64
7 28 26 42 30
8 71 124 54 64
9 43 54 34 56
10 33 30 29 20
11 20 14 39 21
alfa <- 0.05
d_BOD <- x$BOD_C-x$BOD_S
d_SS <- x$SS_C-x$SS_S
d <- data.frame(d_BOD,d_SS)
print.data.frame(d)
d_BOD d_SS
1 -19 12
2 -22 10
3 -18 42
4 -27 15
5 -4 -1
6 -10 11
7 -14 -4
8 17 60
9 9 -2
10 4 10
11 -19 -7
$multivariate_normality
Test Statistic p.value Method MVN
1 Henze-Zirkler 0.583 0.079 asymptotic ✓ Normal
$univariate_normality
Test Variable Statistic p.value Normality
1 Shapiro-Wilk d_BOD 0.920 0.322 ✓ Normal
2 Shapiro-Wilk d_SS 0.819 0.017 ✗ Not normal
$descriptives
Variable n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
1 d_BOD 11 -9.364 14.116 -14 -27 17 -19.0 0.0 0.635 2.147
2 d_SS 11 13.273 20.460 10 -7 60 -1.5 13.5 1.294 3.642
$data
d_BOD d_SS
1 -19 12
2 -22 10
3 -18 42
4 -27 15
5 -4 -1
6 -10 11
7 -14 -4
8 17 60
9 9 -2
10 4 10
11 -19 -7
$subset
NULL
$outlierMethod
[1] "none"
attr(,"class")
[1] "mvn"
d_BOD d_SS
-9.363636 13.272727
d_BOD d_SS
d_BOD 199.25455 88.30909
d_SS 88.30909 418.61818
[1] 11
[1] 2
[1] 9.458877
[1] 13.63931
[1] TRUE
F <- T2/((n-1)*p/(n-p))
pv <- 1-pf(F, p, n-p)
cat("\nF(",p,", ",n-p,") = ", F,", p = ",pv, "\n", sep="")
F(2, 9) = 6.13769, p = 0.02082779
a.hat_abs <- abs(solve(s)%*%centroidif)
colnames(a.hat_abs) <- "Importância"
rownames(a.hat_abs) <- c("BOD", "SS")
print(proportions(a.hat_abs), digits=2)
Importância
BOD 0.59
SS 0.41
me1 <- sqrt(T2crit*s[1,1]/n)
cat("IC95(delta1) = [", centroidif[1]-me1,", ", centroidif[1]+me1,"]", sep="")
IC95(delta1) = [-22.45327, 3.726]
me2 <- sqrt(T2crit*s[2,2]/n)
cat("IC95(delta2) = [", centroidif[2]-me2,", ", centroidif[2]+me2,"]", sep="")
IC95(delta2) = [-5.700119, 32.24557]
c <- sqrt(T2crit/n)
car::ellipse(center=centroidif,
shape=s,
radius=c,
fill=TRUE,
fill.alpha=0.1,
grid=FALSE,
col="black",
add=FALSE,
xlab=expression(mu[BOD[C]]-mu[BOD[S]]),
ylab = expression(mu[SS[C]]-mu[SS[S]]),
main="Região elíptica de confiança de 95%")
points(H0[1], H0[2], pch=3, col="black")
text(H0[1], H0[2], pos=1, expression(H[0]))
[,1] [,2]
d_BOD -22.453272 3.72600
d_SS -5.700119 32.24557
One Sample Hotelling T Square Test
Hotelling T Sqaure Statistic = 13.63931
F value = 6.138 , df1 = 2 , df2 = 9 , p-value: 0.0208
Descriptive Statistics
d_BOD d_SS
N 11.000000 11.00000
Means -9.363636 13.27273
Sd 14.115755 20.46016
Detection important variable(s)
Lower Upper Mu0 Important Variables?
d_BOD -22.453272 3.72600 0 FALSE
d_SS -5.700119 32.24557 0 FALSE
Multivariate Paired Hotelling T Square Test
Hotelling T Sqaure Statistic = 13.63931
F value = 6.138 , df1 = 2 , df2 = 9 , p-value: 0.0208
Descriptive Statistics (The First Treatment)
BOD_C SS_C
Means 25.27273 46.45455
Sd 19.68294 31.84451
Descriptive Statistics (The Second Treatment)
BOD_S SS_S
Means 34.63636 33.18182
Sd 10.45249 19.07259
Descriptive Statistics (The Differences)
BOD_S SS_S
Means 9.363636 -13.27273
Sd 14.115755 20.46016
# Usando matriz de contraste
xmean <- colMeans(x)
sx <- cov(x)
I <- diag(1,p)
C <- cbind(I,-I)
T2_C <- as.numeric(n*t(C%*%xmean)%*%solve(C%*%sx%*%t(C))%*%(C%*%xmean)) # 6-15
T2 # 6-7
[1] 13.63931
# # Bulut, H (2021). A robust Hotelling test statistic for one sample case
# # in high dimensional data,
# # Communication in Statistics: Theory and Methods.
# B <- 1e3
# dqsim <- MVTests::simRHT2(n, p, nrep=B)
# MVTests::RHT2(data=d,
# mu0=H0,
# d=dqsim$d,
# q=dqsim$q)
#
# Comparison <- mvdalab::MVComp(x[,1:2], x[,3:4], level = 1-alfa)
# Comparison
# plot(Comparison, Diff2Plot = c(1,2), include.zero = TRUE)
# MANOVA.RM::multRM: non-normal error terms and/or heteroscedastic
Dados.long <- readxl::read_xlsx("JW6Data/T6-1_long.xlsx")
Dados.long$Amostra <- factor(Dados.long$Amostra)
Dados.long$Lab <- factor(Dados.long$Lab)
print.data.frame(Dados.long)
Amostra Lab BOD SS
1 1 Com 6 27
2 2 Com 6 23
3 3 Com 18 64
4 4 Com 8 44
5 5 Com 11 30
6 6 Com 34 75
7 7 Com 28 26
8 8 Com 71 124
9 9 Com 43 54
10 10 Com 33 30
11 11 Com 20 14
12 1 State 25 15
13 2 State 28 13
14 3 State 36 22
15 4 State 35 29
16 5 State 15 31
17 6 State 44 64
18 7 State 42 30
19 8 State 54 64
20 9 State 34 56
21 10 State 29 20
22 11 State 39 21
alfa <- 0.05
fit <- MANOVA.RM::multRM(cbind(BOD,SS) ~ Lab,
data=Dados.long,
alpha=alfa,
within="Lab",
subject="Amostra")
summary(fit)
Call:
cbind(BOD, SS) ~ Lab
A multivariate repeated measures analysis with 1 within-subject factor(s) ( Lab )and 0 between-subject factor(s).
Descriptive:
Lab n BOD SS
1 Com 11 25.273 46.455
2 State 11 34.636 33.182
Wald-Type Statistic (WTS):
Test statistic df p-value
Lab "13.639" "2" "0.001"
modified ANOVA-Type Statistic (MATS):
Test statistic
Lab 3.348
p-values resampling:
paramBS (WTS) paramBS (MATS)
Lab "0.022" "0.032"
\[\Diamond\]
O pesquisador no Exemplo 6.1 realmente dividiu uma amostra primeiro agitando-a e depois despejando-a rapidamente para frente e para trás em duas garrafas para análise química. Isso foi prudente porque uma simples divisão da amostra em duas partes, obtida despejando a metade superior em uma garrafa e o restante em outra, poderia resultar em mais sólidos suspensos na metade inferior devido à sedimentação. Os dois laboratórios então não estariam trabalhando com as mesmas unidades experimentais, ou mesmo unidades semelhantes, e as conclusões não se relacionariam com a competência do laboratório, técnicas de medição e assim por diante.
Sempre que um investigador pode controlar a atribuição de tratamentos a unidades experimentais, um emparelhamento apropriado de unidades e uma atribuição randomizada de tratamentos podem aprimorar a análise estatística. Diferenças, se houver, entre unidades supostamente idênticas devem ser identificadas e as unidades mais semelhantes emparelhadas. Além disso, uma atribuição aleatória do tratamento 1 a uma unidade e do tratamento 2 à outra unidade ajudará a eliminar os efeitos sistemáticos de fontes não controladas de variação. A randomização pode ser implementada jogando uma moeda para determinar se a primeira unidade de um par recebe o tratamento 1 (cara) ou tratamento 2 (coroa). O tratamento restante é então atribuído à outra unidade. Uma randomização independente separada é realizada para cada par. Pode-se conceber o processo da seguinte forma:
Concluímos nossa discussão sobre comparações pareadas observando que \(\bar{\mathbf{d}}\) e \(\mathbf{s}_\mathbf{d}\), consequentemente, \(T_1\), podem ser calculados a partir das quantidades da amostra completa \(\bar{\mathbf{x}}\) e \(\mathbf{s}\). Aqui, \(\bar{\mathbf{x}}\) é o vetor \(2p \times 1\) das médias amostrais para as \(p\) variáveis nos dois tratamentos, dado por:
\[ \bar{\mathbf{x}}^{\prime} = [\bar{x}_{11}\; \bar{x}_{12}\; \cdots\;\bar{x}_{1p}\; \bar{x}_{21}\; \bar{x}_{22}\; \cdots\; \bar{x}_{2p}] \tag{6-11} \]
\(\mathbf{s}\) é a matriz \(2p \times 2p\) de variâncias e covariâncias amostrais organizada como:
\[ \mathbf{s}= \begin{bmatrix} \mathbf{s}_{11} & \mathbf{s}_{12} \\ \mathbf{s}_{21} & \mathbf{s}_{22} \\ \end{bmatrix} \tag{6-12} \]
em que \(\mathbf{s}_{11}, \mathbf{s}_{12}, \mathbf{s}_{21}, \mathbf{s}_{22}\) são matrizes \(p \times p\).
A matriz \(\mathbf{s}_{11}\) contém as variâncias e covariâncias amostrais para as \(p\) variáveis no tratamento 1. Da mesma forma, \(\mathbf{s}_{22}\) contém as variâncias e covariâncias amostrais calculadas para as \(p\) variáveis no tratamento 2. Finalmente, \(\mathbf{s}_{12} = \mathbf{s}_{21}^{\prime}\) são as matrizes de covariâncias amostrais calculadas a partir das observações em pares das variáveis do tratamento 1 e tratamento 2.
Definindo a matriz:
\[ \underset{p\times 2p}{\mathbf{C}}=\left[\underset{p\times p}{\mathbf{I}}\;-\underset{p\times p}{\mathbf{I}}\right] \]
Podemos verificar (veja o Exercício 6.9) que:
\[ \begin{align} \mathbf{D}_i &= \mathbf{C}\mathbf{X}_{i}\\ i &= 1, 2, \ldots, n \\\\ \bar{\mathbf{D}}&=\mathbf{C}\bar{\mathbf{X}}\\ \mathbf{S}_{\mathbf{D}}&=\mathbf{C}\mathbf{S}\mathbf{C}^{\prime} \end{align} \tag{6-14} \]
Assim,
\[ T^2 = n\left(\mathbf{C}\bar{\mathbf{X}}\right)^{\prime}(\mathbf{C}\mathbf{S}\mathbf{C}^{\prime})^{-1} \mathbf{C}\bar{\mathbf{X}} \tag{6-15} \]
E não é necessário primeiro calcular as diferenças \(\mathbf{D}_1, \mathbf{D}_2, \ldots, \mathbf{D}_{n}\). Por outro lado, é aconselhável calcular essas diferenças para verificar a normalidade e a suposição de uma amostra aleatória.
Cada linha \(\mathbf{c}^{\prime}_i\) da matriz \(\mathbf{C}\) na Equação (6-13) é um vetor de contraste, porque seus elementos somam zero. A atenção geralmente é centrada em contrastes ao comparar tratamentos. Cada contraste é perpendicular ao vetor \(\mathbf{1}^{\prime} = [1\; 1\;\ldots\; 1]\) uma vez que \(\mathbf{c}^{\prime}_i \mathbf{1} = 0\). O componente \(\mathbf{1}^{\prime} \mathbf{x}_i\), representando a soma total do tratamento, é ignorado pela estatística de teste \(T^2\) apresentada nesta seção.
# Usando matriz de contraste
xmean <- colMeans(x)
sx <- cov(x)
I <- diag(1,p)
C <- cbind(I,-I)
T2_C <- as.numeric(n*t(C%*%xmean)%*%solve(C%*%sx%*%t(C))%*%(C%*%xmean)) # 6-15
T2_C
[1] 13.63931
[1] 13.63931
Outra generalização da estatística t pareada univariada surge em situações onde \(q\) tratamentos são comparados em relação a uma única variável resposta. Cada sujeito ou unidade experimental recebe cada tratamento uma vez ao longo de períodos sucessivos de tempo. A i-ésima observação é
\[ \mathbf{X}_{i}= \begin{bmatrix} X_{i1}\\X_{i2}\\\vdots\\X_{iq} \end{bmatrix}\\ i = 1, 2, \ldots ,n \]
em que \(X_{ij}\) é a resposta ao \(j\)-ésimo tratamento na \(i\)-ésima unidade. O nome medidas repetidas deriva do fato de que todos os tratamentos são administrados a cada unidade.
Para fins comparativos, consideramos os contrastes dos componentes de
\[ \boldsymbol{\mu}_i = \mathbb{E}(\mathbf{X}_{i}) \]
Esses podem ser:
\[ \begin{bmatrix} \mu_1 - \mu_2\\ \mu_1 - \mu_3\\ \vdots\\ \mu_1 - \mu_q \end{bmatrix}=\begin{bmatrix} 1&-1&0&\cdots&0&0\\ 1&0&-1&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 1&0&0&\cdots&0&-1\\ \end{bmatrix}= [\mathbf{1}\;-\mathbf{I}] \begin{bmatrix} \mu_{1}\\\mu_{2}\\\vdots\\\mu_{q} \end{bmatrix}= \mathbf{C}_1\boldsymbol{\mu} \]
ou
\[ \begin{bmatrix} \mu_2 - \mu_1\\ \mu_3 - \mu_2\\ \vdots\\ \mu_{q-1} - \mu_q \end{bmatrix}= \begin{bmatrix} -1&1&0&\cdots&0&0\\ 0&-1&1&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&-1&1\\ \end{bmatrix} \begin{bmatrix} \mu_{1}\\\mu_{2}\\\vdots\\\mu_{q} \end{bmatrix}=\mathbf{C}_2\boldsymbol{\mu} \]
Ambos \(\mathbf{C}_1\) e \(\mathbf{C}_2\) são chamadas de matrizes de contraste, porque suas \(q-1\) linhas são linearmente independentes e cada uma é um vetor de contraste (soma nula e ortogonais ao vetor unitário).
A natureza do desenho elimina grande parte da influência da variação de unidade para unidade nas comparações de tratamento. Claro, o pesquisador deve randomizar a ordem na qual os tratamentos são apresentados a cada sujeito.
Quando as médias dos tratamentos são iguais, \(\mathbf{C}_1\boldsymbol{\mu} = \mathbf{C}_2\boldsymbol{\mu} = \mathbf{0}\). Em geral, a hipótese de que não há diferenças nos tratamentos (médias de tratamento iguais) torna-se \(\mathbf{C}\boldsymbol{\mu} = \mathbf{0}\) para qualquer escolha da matriz de contraste \(\mathbf{C}\).
Consequentemente, com base nos contrastes \(\mathbf{C}\mathbf{x}_i\) nas observações, temos médias \(\mathbf{C}\bar{\mathbf{x}}\) e matriz de covariância \(\mathbf{C}\mathbf{S}\mathbf{C}^{\prime}\), e testamos \(\mathbf{C}\boldsymbol{\mu} = \mathbf{0}\) usando a estatística \(T^2\):
\[ T^2 = n\left(\mathbf{C}\bar{\mathbf{X}}\right)^{\prime}(\mathbf{C}\mathbf{S}\mathbf{C}^{\prime})^{-1} \mathbf{C}\bar{\mathbf{X}} \]
Considere uma população \(\mathcal{N}_q(\boldsymbol{\mu}, \mathbf{\Sigma})\) e seja \(\mathbf{C}\) uma matriz de contraste. Um teste de nível \(\alpha\) de \(H_0: \mathbf{C}\boldsymbol{\mu} = \mathbf{0}\) vs. \(H_1: \boldsymbol{\mu} \ne \mathbf{0}\) é o seguinte:
Rejeite \(H_0\) se:
\[ \begin{align} T^2 &= n\left(\mathbf{C}\bar{\mathbf{x}}\right)^{\prime}(\mathbf{C}\mathbf{s}\mathbf{C}^{\prime})^{-1} \mathbf{C}\bar{\mathbf{x}} \\ T^2 &> \dfrac{(n-1)(q-1)}{n-q+1}F_{q-1, n-q+1}(1-\alpha)=T^2_{n-1,q-1}(1-\alpha) \end{align} \tag{6-16} \]
Pode-se mostrar que \(T^2\) não depende da escolha particular de \(\mathbf{C}\).
Qualquer par de matrizes de contraste \(\mathbf{C}_1\) e \(\mathbf{C}_2\) deve estar relacionado por \(\mathbf{C}_1 = \mathbf{B}\mathbf{C}_2\), com \(\mathbf{B}\) não singular. Isso ocorre porque cada \(\mathbf{C}\) tem o maior número possível, \(q - 1\), de linhas linearmente independentes, todas perpendiculares ao vetor 1. Então, temos:
\[ \begin{align} \mathbf{C}_1^{\prime}(\mathbf{C}_1\mathbf{S}\mathbf{C}_1^{\prime})^{-1}\mathbf{C}_1&=(\mathbf{B}\mathbf{C}_2)^{\prime}(\mathbf{B}\mathbf{C}_2\mathbf{S}\mathbf{C}_2^{\prime}\mathbf{B}^{\prime})^{-1}(\mathbf{B}\mathbf{C}_2) \\ &= \mathbf{C}_2^{\prime}\mathbf{B}^{\prime}(\mathbf{B}^{\prime})^{-1}(\mathbf{C}_2\mathbf{S}\mathbf{C}_2^{\prime})^{-1}\mathbf{B}^{-1}\mathbf{B}\mathbf{C}_2 \\ \mathbf{C}_1^{\prime}(\mathbf{C}_1\mathbf{S}\mathbf{C}_1^{\prime})^{-1}\mathbf{C}_1&= \mathbf{C}_2^{\prime}(\mathbf{C}_2\mathbf{S}\mathbf{C}_2^{\prime})^{-1}\mathbf{C}_2 \end{align} \]
Assim, \(T^2\) calculado com \(\mathbf{C}_2\) ou \(\mathbf{C}_1 = \mathbf{B}\mathbf{C}_2\) produz o mesmo resultado.
Uma região de confiança para os contrastes \(\boldsymbol{C}\boldsymbol{\mu}\), com \(\boldsymbol{\mu}\) sendo a média de uma população normal, é determinada pelo conjunto de todos \(\boldsymbol{C}\boldsymbol{\mu}\) tais que:
\[ n(\mathbf{C}\bar{\mathbf{x}} - \mathbf{C}\boldsymbol{\mu})^{\prime}(\mathbf{C}\mathbf{s}\mathbf{C}^{\prime})^{-1}(\mathbf{C}\bar{\mathbf{x}} - \mathbf{C}\boldsymbol{\mu}) \leq \dfrac{(n-1)(q - 1)}{n - q+1}F_{q-1, n-q+1}(1-\alpha) \tag{6-17} \]
em que \(\bar{\mathbf{x}}\) e \(\mathbf{s}\) são definidos como na Equação (6-16). Portanto, intervalos de confiança simultâneos de \(100(1 - \alpha)\%\) para contrastes simples \(\mathbf{c}^{\prime}\boldsymbol{\mu}\) (e.g., cada linha de \(\mathbf{C}\)) para quaisquer vetores de contraste de interesse são dados por (veja o Resultado 5A.1):
\[ \text{IC}^{1-\alpha}(\mathbf{c}^{\prime}\boldsymbol{\mu})=\left[\mathbf{c}^{\prime}\bar{\mathbf{x}} \pm \sqrt{\dfrac{(n-1)(q - 1)}{n - q+1}F_{q-1, n-q+1}(1-\alpha)} \sqrt{\dfrac{\mathbf{c}\mathbf{s}\mathbf{c}^{\prime}}{n}} \right] \tag{6-18} \]
\(p=1, q=4, g=1, n=19\)
UE: cão sob pentabarbitol
Anestésicos melhorados são frequentemente desenvolvidos primeiro
estudando seus efeitos em animais não-humanos. Em um estudo, 19 cães
foram inicialmente administrados com o medicamento
pentobarbitol*
. Cada cão foi então administrado com dióxido
de carbono (CO\(_2\)) em dois níveis de
pressão. Em seguida, halotano (H) foi adicionado, e a administração de
CO\(_2\) foi repetida. A resposta, em
milissegundos entre batimentos cardíacos, foi medida para as quatro
combinações de tratamento:
*
: Em contextos reais, o pentobarbitol é um barbitúrico
que pode ser usado como anestésico ou sedativo. Em experimentos com
animais, pode ser usado para sedar ou anestesiar os animais antes de
administrar outros tratamentos ou realizar procedimentos, garantindo
assim que os animais estejam calmos e não sintam dor ou desconforto.
Neste caso, é possível que o pentobarbitol tenha sido usado para sedar
ou anestesiar os cães antes da administração de CO2 e
halotano, de modo a controlar variáveis externas, reduzir o estresse ou
desconforto dos animais e garantir a segurança do procedimento.
A Tabela 6.2 contém as quatro medições para cada um dos 19 cães, onde:
Analisaremos os efeitos anestésicos da pressão CO2 e do halotano a partir deste delineamento de medidas repetidas.
Existem três contrastes de tratamento que podem ser de interesse no experimento. Sejam \(\boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \boldsymbol{\mu}_3\) e \(\boldsymbol{\mu}_4\) as respostas médias para os tratamentos 1, 2, 3 e 4, respectivamente. Então:
\[ (\boldsymbol{\mu}_3 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2) \]
\[ (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_3) - (\boldsymbol{\mu}_2 + \boldsymbol{\mu}_4) \]
\[ (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_2 + \boldsymbol{\mu}_3) \]
ANOVA bifatorial relacionada (ou para medidas repetidas) (rmANOVA)
não-balanceada é possível em arquivo de dados no formato long
se a identificação da unidade experimental é modelada como efeito
aleatório. No entanto, esta ANOVA bifatorial relacionada usando
lmerTest::lmer(heartbeats ~ H*CO2 + (1|UE), data=Dados)
,
car::Anova(modelo,test.statistic="F")
e
summary(modelo, correl=FALSE)
não produz o valor p
omnibus, não sendo assim uma autêntica análise multivariada. A
maior vantagem desta análise é permitir desbalanceamento num
delineamento intraparticipantes. A modelagem multivariada com \(T^2\) usa arquivo de dados no formato
wide ou apenas o vetor de média e matriz de covariância
amostrais (que pode ser obtida mesmo com dados faltantes, em geral, de
maneira inadequada; FIML necessita de valores faltantes aleatoriamente
(MAR: Missing At Random)).
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
x <- read.table("JW6Data/T6-2.dat", quote="\"", comment.char="")
names(x) <- c("CO2high_Hwithout", "CO2low_Hwithout",
"CO2high_Hwith", "CO2low_Hwith")
print.data.frame(x)
CO2high_Hwithout CO2low_Hwithout CO2high_Hwith CO2low_Hwith
1 426 609 556 600
2 253 236 392 395
3 359 433 349 357
4 432 431 522 600
5 405 426 513 513
6 324 438 507 539
7 310 312 410 456
8 326 326 350 504
9 375 447 547 548
10 286 286 403 422
11 349 382 473 497
12 429 410 488 547
13 348 377 447 514
14 412 473 472 446
15 347 326 455 468
16 434 458 637 524
17 364 367 432 469
18 420 395 508 531
19 397 556 645 625
[,1] [,2] [,3] [,4]
[1,] -1 -1 1 1
[2,] 1 -1 1 -1
[3,] 1 -1 -1 1
CO2high_Hwithout CO2low_Hwithout CO2high_Hwith CO2low_Hwith
368.2105 404.6316 479.2632 502.8947
CO2high_Hwithout CO2low_Hwithout CO2high_Hwith CO2low_Hwith
CO2high_Hwithout 2819.287 3568.415 2943.497 2295.357
CO2low_Hwithout 3568.415 7963.135 5303.991 4065.459
CO2high_Hwith 2943.497 5303.991 6851.316 4499.640
CO2low_Hwith 2295.357 4065.459 4499.640 4878.988
[1] 19
[1] 4
[1] 10.93119
[1] 116.0163
[1] TRUE
F <- T2/((n-1)*(q-1)/(n-q+1))
pv <- 1-pf(F, q-1, n-q+1)
cat("\nF(",q-1,", ",n-q+1,") = ", F,", p = ",pv, "\n", sep="")
F(3, 16) = 34.37521, p = 3.317767e-07
a.hat_abs <- abs(solve(C%*%s%*%t(C))%*%(C%*%mean))
colnames(a.hat_abs) <- "Importância"
rownames(a.hat_abs) <- c("H", "CO2p", "H:CO2p")
print(proportions(a.hat_abs), digits=2)
Importância
H 0.562
CO2p 0.375
H:CO2p 0.063
# Interação
c <- C[3,]
me <- sqrt(T2crit*t(c)%*%s%*%c/n)
cat("IC95%(interação) = [",t(c)%*%mean-me,", ",
t(c)%*%mean+me, "]", sep="")
IC95%(interação) = [-78.72858, 53.14964]
# Halotano
c <- C[1,]
me <- sqrt(T2crit*t(c)%*%s%*%c/n)
cat("IC95%(Halotano) = [",t(c)%*%mean-me,", ",
t(c)%*%mean+me, "]", sep="")
IC95%(Halotano) = [135.6503, 282.9813]
# Pressão de CO2
c <- C[2,]
me <- sqrt(T2crit*t(c)%*%s%*%c/n)
cat("IC95%(Pressão de CO2) = [",t(c)%*%mean-me,", ",
t(c)%*%mean+me, "]", sep="")
IC95%(Pressão de CO2) = [-114.7271, -5.37818]
x$UE <- factor(1:nrow(x))
dados_wide <- x
# Transformar os dados de formato wide para longo
Dados_long <- tidyr::gather(dados_wide, key = "variavel",
value = "FCLatencia", -UE)
# Separar as variáveis CO2 e H
Dados_long$CO2 <- ifelse(grepl("high", Dados_long$variavel),
"High", "Low")
Dados_long$H <- ifelse(grepl("without", Dados_long$variavel,
ignore.case = TRUE),
"Without", "With")
# Remover a coluna variável
Dados_long$variavel <- NULL
print(ftable(Dados_long$UE, Dados_long$H, Dados_long$CO2))
High Low
1 With 1 1
Without 1 1
2 With 1 1
Without 1 1
3 With 1 1
Without 1 1
4 With 1 1
Without 1 1
5 With 1 1
Without 1 1
6 With 1 1
Without 1 1
7 With 1 1
Without 1 1
8 With 1 1
Without 1 1
9 With 1 1
Without 1 1
10 With 1 1
Without 1 1
11 With 1 1
Without 1 1
12 With 1 1
Without 1 1
13 With 1 1
Without 1 1
14 With 1 1
Without 1 1
15 With 1 1
Without 1 1
16 With 1 1
Without 1 1
17 With 1 1
Without 1 1
18 With 1 1
Without 1 1
19 With 1 1
Without 1 1
'data.frame': 76 obs. of 4 variables:
$ UE : Factor w/ 19 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ FCLatencia: int 426 253 359 432 405 324 310 326 375 286 ...
$ CO2 : chr "High" "High" "High" "High" ...
$ H : chr "Without" "Without" "Without" "Without" ...
source("summarySEwithin2.R")
alfa <- 0.05
print(psych::describe(FCLatencia ~ H*CO2,
data=Dados_long), digits=2)
Descriptive statistics by group
H: With
CO2: High
vars n mean sd median trimmed mad min max range skew
FCLatencia 1 19 479.26 82.77 473 477.18 72.65 349 645 296 0.32
kurtosis se
FCLatencia -0.62 18.99
------------------------------------------------------------
H: Without
CO2: High
vars n mean sd median trimmed mad min max range skew
FCLatencia 1 19 368.21 53.1 364 371.12 60.79 253 434 181 -0.43
kurtosis se
FCLatencia -0.92 12.18
------------------------------------------------------------
H: With
CO2: Low
vars n mean sd median trimmed mad min max range skew
FCLatencia 1 19 502.89 69.85 513 504.29 65.23 357 625 268 -0.21
kurtosis se
FCLatencia -0.66 16.02
------------------------------------------------------------
H: Without
CO2: Low
vars n mean sd median trimmed mad min max range skew
FCLatencia 1 19 404.63 89.24 410 402.53 63.75 236 609 373 0.32
kurtosis se
FCLatencia -0.12 20.47
alfaBonf <- alfa/(length(unique(Dados_long$H))*
length(unique(Dados_long$CO2)))
ic <- summarySEwithin2(Dados_long,
measurevar="FCLatencia",
withinvars=c("H","CO2"),
idvar="UE",
na.rm=TRUE,
conf.interval=1-alfaBonf)
Automatically converting the following non-factors to factors: H, CO2
Anexando pacote: 'data.table'
O seguinte objeto é mascarado por 'package:reshape':
melt
Os seguintes objetos são mascarados por 'package:dplyr':
between, first, last
O seguinte objeto é mascarado por 'package:DescTools':
%like%
H CO2 N FCLatencia FCLatenciaNormed sd se ci
1 With High 19 479.2632 479.2632 41.53206 9.528107 26.43601
2 With Low 19 502.8947 502.8947 40.83881 9.369066 25.99474
3 Without High 19 368.2105 368.2105 40.78560 9.356859 25.96087
4 Without Low 19 404.6316 404.6316 48.36286 11.095201 30.78396
grf <- ggplot2::ggplot(ic,
ggplot2::aes(x=H,
y=FCLatencia,
colour=CO2)) +
ggplot2::geom_errorbar(position=ggplot2::position_dodge(.9),
width=.1,
ggplot2::aes(ymin=FCLatencia-ci,
ymax=FCLatencia+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white",
position=ggplot2::position_dodge(.9)) +
ggplot2::ylab("FCLatencia") +
ggplot2::ggtitle("H & CO2: FCLatencia\nWithin-subject CI95% Bonferroni") +
ggplot2::theme_bw()
print(grf)
ANOVA
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: FCLatencia
F Df Df.res Pr(>F)
H 112.5668 1 54 8.071e-15 ***
CO2 9.2655 1 54 0.003604 **
H:CO2 0.4203 1 54 0.519558
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: FCLatencia ~ H * CO2 + (1 | UE)
Data: Dados_long
REML criterion at convergence: 797.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.06165 -0.42869 -0.00709 0.41075 2.49423
Random effects:
Groups Name Variance Std.Dev.
UE (Intercept) 3779 61.48
Residual 1849 43.00
Number of obs: 76, groups: UE, 19
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 479.26 17.21 30.60 27.846 < 2e-16 ***
HWithout -111.05 13.95 54.00 -7.961 1.14e-10 ***
CO2Low 23.63 13.95 54.00 1.694 0.096 .
HWithout:CO2Low 12.79 19.73 54.00 0.648 0.520
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect size analysis
eta2g <- as.numeric(MuMIn::r.squaredGLMM(modelo)[2])
cat("\nTamanho de efeito: eta^2 omnibus =", eta2g)
Tamanho de efeito: eta^2 omnibus = 0.7860638
es <- effectsize::interpret_eta_squared(eta2g)
names(es) <- c("Tamanho de efeito omnibus: estimativa pontual")
print(es)
Tamanho de efeito omnibus: estimativa pontual
"large"
(Rules: field2013)
eta2 <- effectsize::eta_squared(anv,
partial=FALSE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=6)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 98.3333% CI | interpret
--------------------------------------------------------------
H | 0.675806 | [0.488697, 0.784409] | large
CO2 | 0.146455 | [0.005033, 0.361202] | large
H:CO2 | 0.007722 | [0.000000, 0.146028] | very small
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=6)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 98.3333% CI | interpret
--------------------------------------------------------------
H | 0.675806 | [0.488697, 0.784409] | large
CO2 | 0.146455 | [0.005033, 0.361202] | large
H:CO2 | 0.007722 | [0.000000, 0.146028] | very small
# Grafico de perfis de medias
o.par <- par()
fit.means <- phia::interactionMeans(modelo)
plot(fit.means,
errorbar=paste0("ci",
round((1-alfaBonf)*100,4)),
abbrev.levels=FALSE)
par(o.par)
plot(effects::effect(c("H"), modelo,
confidence.level=1-alfa/length(unique(Dados_long$H))),
ci.style = "bars")
NOTE: H is not a high-order term in the model
plot(effects::effect(c("CO2"), modelo,
confidence.level=1-alfa/length(unique(Dados_long$CO2))),
ci.style = "bars")
NOTE: CO2 is not a high-order term in the model
plot(effects::effect(c("H", "CO2"), modelo, confidence.level=1-alfaBonf),
multiline = TRUE, ci.style = "bars")
NOTE: HCO2 is not a high-order term in the model
Post hoc tests
H
EMM.A <- emmeans::emmeans(modelo,
specs=pairwise~"H",
adjust="holm",
level=1-alfa,
lmer.df="satterthwaite",
lmerTest.limit=nrow(Dados_long))
NOTE: Results may be misleading due to involvement in interactions
H emmean SE df lower.CL upper.CL
With 491 15.7 22.1 458 524
Without 386 15.7 22.1 354 419
Results are averaged over the levels of: CO2
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
contrast estimate SE df lower.CL upper.CL t.ratio p.value
With - Without 105 9.86 54 84.9 124 10.610 <.0001
Results are averaged over the levels of: CO2
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
print(multcomp::cld(object=EMM.A$emmeans,
level=1-alfa,
adjust="holm",
Letters=letters,
alpha=alfa))
H emmean SE df lower.CL upper.CL .group
Without 386 15.7 22.1 349 424 a
With 491 15.7 22.1 453 529 b
Results are averaged over the levels of: CO2
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 2 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
CO2
EMM.C <- emmeans::emmeans(modelo,
specs=pairwise~"CO2",
adjust="holm",
level=1-alfa,
lmer.df="satterthwaite",
lmerTest.limit=nrow(Dados_long))
NOTE: Results may be misleading due to involvement in interactions
CO2 emmean SE df lower.CL upper.CL
High 424 15.7 22.1 391 456
Low 454 15.7 22.1 421 486
Results are averaged over the levels of: H
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
contrast estimate SE df lower.CL upper.CL t.ratio p.value
High - Low -30 9.86 54 -49.8 -10.2 -3.044 0.0036
Results are averaged over the levels of: H
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
print(multcomp::cld(object=EMM.C$emmeans,
level=1-alfa,
adjust="holm",
Letters=letters,
alpha=alfa))
CO2 emmean SE df lower.CL upper.CL .group
High 424 15.7 22.1 386 462 a
Low 454 15.7 22.1 416 492 b
Results are averaged over the levels of: H
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 2 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
A partir da Equação (6-16), temos \(T^2 = 116 > 10.93\), e rejeitamos \(H_0: \mathbf{C}\boldsymbol{\mu} = \mathbf{0}\) (sem efeitos de tratamento). Para determinar quais dos contrastes são responsáveis pela rejeição de \(H_0\), construímos intervalos de confiança simultâneos de 95% para esses contrastes. A partir da Equação (6-18), o contraste
\[ \mathbf{c}^{\prime}\boldsymbol{\mu} = (\boldsymbol{\mu}_3 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2) = \text{influência do halotano} \]
é estimado pelo intervalo
\[ \begin{align} \text{IC}^{95\%}((\boldsymbol{\mu}_3 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2)) &=\left[(\bar{x}_3 + \bar{x}_4) - (\bar{x}_1 + \bar{x}_2) \pm \sqrt{\dfrac{18\times3}{16} F_{3,16}(0.95)}\sqrt{\dfrac{\mathbf{c}_1^{\prime}\mathbf{s}\mathbf{c}_1}{19}}\right]\\ &= \left[209.31 \pm \sqrt{10.94}\sqrt{\dfrac{9432.32}{16}}\right] \\ &= \left[209.31 \pm 73.70\right]\\ \text{IC}^{95\%}((\boldsymbol{\mu}_3 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2))&=[135.61,283.01] \end{align} \]
em que \(\mathbf{c}_1\) é a primeira linha de \(\mathbf{C}\). Da mesma forma, os contrastes remanescentes são estimados por:
Influência da pressão CO2:
\[ \text{IC}^{95\%}((\boldsymbol{\mu}_1 + \boldsymbol{\mu}_3) - (\boldsymbol{\mu}_2 + \boldsymbol{\mu}_4))=[5.35,114.75] \]
Interação H-CO2: \[ \text{IC}^{95\%}((\boldsymbol{\mu}_1 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_2 + \boldsymbol{\mu}_3))=[-78.76,53.18] \]
A análise dos efeitos principais e de interação é hierárquica. Primeiramente, é analisado o efeito de interação. Se a hipótese nula de ausência de efeito de interação não é rejeitada. Os efeitos principais são diretamente analisáveis pelos respectivos intervalos de confiança.
O contraste de interação H-CO2, \((\boldsymbol{\mu}_1 + \boldsymbol{\mu}_4) - (\boldsymbol{\mu}_2 + \boldsymbol{\mu}_3)\) não é significativamente diferente de zero (veja o terceiro intervalo de confiança).
Portanto, o primeiro intervalo de confiança implica que há um efeito do halotano. A presença de halotano produz tempos mais longos entre batimentos cardíacos.
O segundo intervalo de confiança indica que há um efeito devido à pressão CO2: a menor pressão CO2 produz tempos mais longos entre batimentos cardíacos.
Deve-se ter cautela na interpretação dos resultados, pois os ensaios com halotano devem seguir aqueles sem halotano. O aparente efeito do H pode ser devido a uma tendência temporal. (Idealmente, a ordem temporal de todos os tratamentos deveria ser determinada aleatoriamente.)
\[\Diamond\]
O teste na Equação (6-16) é apropriado quando a matriz de covariância, \(\mathbb{C}(\mathbf{X}) = \mathbf{\Sigma}\), não pode ser assumida com nenhuma estrutura especial. Se for razoável supor que \(\mathbf{\Sigma}\) tenha uma estrutura particular, os testes projetados com essa estrutura em mente têm maior poder do que o da Equação (6-16). (Para \(\mathbf{\Sigma}\) com a estrutura de correlação igual (8-14), consulte uma discussão sobre o delineamento “bloco randomizado” em [17] ou [22].)
Um teste \(T^2\) para testar a igualdade de médias vetoriais de duas populações multivariadas pode ser desenvolvido por analogia com o procedimento univariado. (Consulte [11] para uma discussão sobre o caso univariado.) Esta estatística \(T^2\) é apropriada para comparar respostas de um conjunto de condições experimentais (população 1) com respostas independentes de outro conjunto de condições experimentais (população 2). A comparação pode ser feita sem controlar explicitamente a variabilidade de unidade para unidade, como no caso de comparação pareada.
Se possível, as unidades experimentais devem ser atribuídas aleatoriamente aos conjuntos de condições experimentais. A randomização irá, até certo ponto, mitigar o efeito da variabilidade de unidade para unidade em uma subsequente comparação de tratamentos. Embora alguma precisão seja perdida em relação às comparações pareadas, as inferências no caso de duas populações são, normalmente, aplicáveis a uma coleção mais geral de unidades experimentais, simplesmente porque a homogeneidade da unidade não é exigida.
Considere uma amostra aleatória de tamanho \(n_1\) da população 1 e uma amostra de tamanho \(n_2\) da população 2. As observações sobre \(p\) variáveis podem ser organizadas da seguinte forma:
Amostra (População 1)
\[ \mathbf{x}_{11}, \mathbf{x}_{12}, \ldots, \mathbf{x}_{1n_1} \]
\[ \bar{\mathbf{x}}_1 = \dfrac{1}{n_1} \sum_{i=1}^{n_1} \mathbf{x}_{1i} \quad \mathbf{s}_1 = \dfrac{1}{n_1} \sum_{i=1}^{n_1} (\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)(\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)^{\prime} \]
Amostra (População 2)
\[ \mathbf{x}_{21}, \mathbf{x}_{22}, \ldots, \mathbf{x}_{2n_2} \]
\[ \bar{\mathbf{x}}_2 = \dfrac{1}{n_2} \sum_{i=1}^{n_2} \mathbf{x}_{2i} \quad \mathbf{s}_2 = \dfrac{1}{n_2} \sum_{i=1}^{n_2} (\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)(\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)^{\prime} \]
Nesta notação, o primeiro subscrito — 1 ou 2 — denota a população.
Queremos fazer inferências sobre
\[ \text{vetor médio da população 1} - \text{vetor médio da população 2} = \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 \]
Por exemplo, queremos responder à questão: \(\boldsymbol{\mu}_1 = \boldsymbol{\mu}_2\) (ou, equivalentemente, \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = \mathbf{0}\))? Além disso, se \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 \neq \mathbf{0}\), quais médias componentes são diferentes?
Com algumas suposições, somos capazes de fornecer respostas para essas perguntas.
A amostra \(\mathbf{X}_{11}, \mathbf{X}_{12}, \ldots, \mathbf{X}_{1n_1}\) é uma amostra aleatória de tamanho \(n_1\) de uma população p-variada com vetor médio \(\boldsymbol{\mu}_1\) e matriz de covariância \(\mathbf{\Sigma}_1\).
A amostra \(\mathbf{X}_{21}, \mathbf{X}_{22}, \ldots, \mathbf{X}_{2n_2}\) é uma amostra aleatória de tamanho \(n_2\) de uma população p-variada com vetor médio \(\boldsymbol{\mu}_2\) e matriz de covariância \(\mathbf{\Sigma}_2\).
Além disso, \(\mathbf{X}_{11}, \mathbf{X}_{12}, \ldots, \mathbf{X}_{1n_1}\) são independentes de \(\mathbf{X}_{21}, \mathbf{X}_{22}, \ldots, \mathbf{X}_{2n_2}\). (6-19)
Veremos adiante, para amostra grande, essa estrutura é suficiente para fazer inferências sobre o vetor \(p \times 1\) \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\). No entanto, quando os tamanhos das amostras \(n_1\) e \(n_2\) são pequenos, são necessárias mais suposições.
Ambas as populações são multivariadas normais.
Além disso, \(\mathbf{\Sigma}_1 = \mathbf{\Sigma}_2\) (homocedasticidade multivariada). (6-20)
A segunda suposição, de que \(\mathbf{\Sigma}_1 = \mathbf{\Sigma}_2\), é muito mais forte que sua contraparte univariada. Aqui estamos supondo que vários pares de variâncias e covariâncias são iguais.
Quando \(\mathbf{\Sigma}_1 = \mathbf{\Sigma}_2 = \mathbf{\Sigma}\),
\[ \sum_{i=1}^{n_1} (\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)(\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)^{\prime} \]
é uma estimativa de \((n_1 - 1)\mathbf{\Sigma}\) e
\[ \sum_{i=1}^{n_2} (\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)(\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)^{\prime} \]
é uma estimativa de \((n_2 - 1)\mathbf{\Sigma}\). Consequentemente, podemos agrupar as informações em ambas as amostras para estimar a covariância comum \(\mathbf{\Sigma}\).
Definimos
\[ \mathbf{S}_{\text{comb}} = \dfrac{1}{n_1 + n_2 - 2} \left( \sum_{i=1}^{n_1} (\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)(\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)^{\prime} + \sum_{i=1}^{n_2} (\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)(\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)^{\prime} \right) \]
ou, de forma equivalente,
\[ \mathbf{S}_{\text{comb}} = \dfrac{n_1 - 1}{n_1 + n_2 - 2} \mathbf{S}_1 + \dfrac{n_2 - 1}{n_1 + n_2 - 2} \mathbf{S}_2 \tag{6-21} \]
Dado que \(\sum_{i=1}^{n_1} (\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)(\mathbf{x}_{1i} - \bar{\mathbf{x}}_1)^{\prime}\) tem \(n_1 - 1\) graus de liberdade e \(\sum_{i=1}^{n_2} (\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)(\mathbf{x}_{2i} - \bar{\mathbf{x}}_2)^{\prime}\) tem \(n_2 - 1\) graus de liberdade, o divisor \((n_1 - 1) + (n_2 - 1)\) em (6-21) é obtido combinando os dois graus de liberdade componentes. [Veja (4-24)]. Suporte adicional para o procedimento de agrupamento vem da consideração da verossimilhança normal multivariada. (Veja o Exercício 6.11.)
Para testar a hipótese de que \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = \boldsymbol{\delta}_0\), um vetor especificado, consideramos a distância estatística ao quadrado de \(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\) para \(\boldsymbol{\delta}_0\). Agora,
\[ \begin{align} \mathbb{E}(\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2) &= \mathbb{E}(\bar{\mathbf{X}}_1) - \mathbb{E}(\bar{\mathbf{X}}_2) \\ \mathbb{E}(\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2)&= \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 \end{align} \]
Dado que a suposição de independência em (6-19) implica que \(\bar{\mathbf{X}}_1\) e \(\bar{\mathbf{X}}_2\) são independentes e, portanto, \(\mathbb{C}(\bar{\mathbf{X}}_1, \bar{\mathbf{X}}_2) = \mathbf{0}\) (veja o Resultado 4.5), por (3-9), segue-se que:
\[ \begin{align} \mathbb{C}(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2) &= \mathbb{C}(\bar{\mathbf{X}}_1) + \mathbb{C}(\bar{\mathbf{X}}_2) \\ &= \dfrac{1}{n_1}\mathbf{\Sigma} + \dfrac{1}{n_2}\mathbf{\Sigma} \\ \mathbb{C}(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2) &= \left( \dfrac{1}{n_1} + \dfrac{1}{n_2} \right)\mathbf{\Sigma} \end{align} \tag{6-22} \]
Como \(\mathbf{S}_{\text{comb}}\) estima \(\mathbf{\Sigma}\), vemos que
\[ \hat{\mathbb{C}}(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2)=\left( \dfrac{1}{n_1} + \dfrac{1}{n_2} \right)\mathbf{S}_{\text{comb}} \]
é um estimador de \(\mathbb{C}(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2)\).
\[ \hat{\mathbb{C}}(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2)=\dfrac{\mathbf{S}_{\text{comb}}}{\bar{n}_h/2} \]
em que \(\bar{n}_h=\dfrac{2}{\dfrac{1}{n_1}+\dfrac{1}{n_2}}\) é a média harmônica dos tamanhos de amostra das duas condições independentes.
O teste da razão de verossimilhança de:
\[ H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = \boldsymbol{\delta}_0 \]
é baseado no quadrado da distância estatística, \(T^2\), e é dado por (consulte [1]).
Rejeite \(H_0\) se:
\[ T^2 = (\bar{\mathbf{x}}_1- \bar{\mathbf{x}}_2 - \boldsymbol{\delta}_0)^{\prime} \left(\left( \dfrac{1}{n_1} + \dfrac{1}{n_2} \right)\mathbf{s}_{\text{comb}}\right)^{-1} (\bar{\mathbf{x}}_1- \bar{\mathbf{x}}_2 - \boldsymbol{\delta}_0) > c^2 \tag{6-23} \]
em que a distância crítica \(c^2\) é determinada a partir da distribuição da estatística \(T^2\) de duas condições independentes.
Resultado 6.2. Se \(\{\mathbf{X}_{1i}\}_{i=1}^{n_1} \sim \mathcal{N}_p\text{IID}(\boldsymbol{\mu}_1, \mathbf{\Sigma})\) e ${{2i}}{i=1}^{n_2} _p(_2, ) $, então:
\[ \begin{align} T^2 &= \left(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2 - (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)\right)^{\prime} \left(\left( \dfrac{1}{n_1} + \dfrac{1}{n_2}\right) \mathbf{S}_{\text{comb}}\right)^{-1} \left(\bar{\mathbf{X}}_1- \bar{\mathbf{X}}_2 - (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)\right) \\&\sim \dfrac{(n_1 + n_2 - 2)p}{n_1 + n_2 - p - 1}F_{p, n_1 + n_2 - p - 1}=T^2_{p,n_1 + n_2 - 2} \end{align} \tag{6-24} \]
\[\Diamond\]
Estamos principalmente interessados em regiões de confiança para \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\). A partir de (6-24), concluímos que todos os \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\) dentro da distância estatística ao quadrado \(c^2\) de \(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\) constituem a região de confiança. Esta região é um elipsoide centrado na diferença observada \(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\) e cujos eixos são determinados pelos autovalores e autovetores de \(\mathbf{S}_{\text{comb}}\) (ou \(\mathbf{S}_{\text{comb}}^{-1}\)).
Cinquenta barras de sabão são fabricadas de duas maneiras diferentes. Duas características, \(X_1\) (espumosidade) e \(X_2\) (suavidade), são medidas. As estatísticas resumidas para barras produzidas pelos métodos 1 e 2 são:
\[ \bar{\mathbf{x}}_1 = \begin{bmatrix} 10.2 \\ 3.9 \end{bmatrix} \quad \mathbf{s}_1 = \begin{bmatrix} 2 & 1 \\ 1 & 6 \end{bmatrix} \] \[ \bar{\mathbf{x}}_2 = \begin{bmatrix} 8.3 \\ 4.1 \end{bmatrix} \quad \mathbf{s}_2 = \begin{bmatrix} 2 & 1 \\ 1 & 4 \end{bmatrix} \]
Obtenha uma região de confiança de 95% para \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\).
Primeiramente, observamos que \(\mathbf{s}_1\) e \(\mathbf{s}_2\) são aproximadamente iguais, então é razoável agrupá-los. Assim, a partir de (6-21), temos:
\[ \mathbf{s}_{\text{comb}} = \dfrac{49}{98} (\mathbf{s}_1 + \mathbf{s}_2) \]
Além disso,
\[ \bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2 = \begin{bmatrix} -1.9 \\ 0.2 \end{bmatrix} \]
então a elipse de confiança é centrada em \([-1.9, 0.2]^{\prime}\). Os autovalores e autovetores de \(\mathbf{s}_{\text{comb}}\) são obtidos a partir da equação:
\[ \det(\mathbf{s}_{\text{comb}} - \lambda \mathbf{I}) = 0 \]
o que nos dá \(\lambda^2 - 7\lambda + 9 = 0\). Daí, \(\lambda_1 = 5.303\) e \(\lambda_2 = 1.697\). Os autovetores correspondentes, \(\mathbf{e}_1\) e \(\mathbf{e}_2\), determinados a partir de \(\mathbf{S}_{\text{comb}} \mathbf{e}_i = \lambda_i \mathbf{e}_i\), são:
\[ \mathbf{e}_1 = \begin{bmatrix} 0.290 \\ 0.957 \end{bmatrix} \quad \text{e} \quad \mathbf{e}_2 = \begin{bmatrix} 0.957\\ -0.290\end{bmatrix} \]
Pelo Resultado 6.2,
\[ \left(\dfrac{1}{n_1}+\dfrac{1}{n_2}\right)c^2=\dfrac{2}{50} \dfrac{98\times 2}{97} F_{2, 97}(0.95) = 0.25 \]
já que \(F_{2,97}(0.95) = 3.1\). A elipse de confiança se estende por \(\sqrt{\lambda_1}\sqrt{\left(\dfrac{1}{n_1}+\dfrac{1}{n_2}\right)c^2}=\sqrt{\lambda_1}\sqrt{0.25}\) se estende no autovetor \(\mathbf{e}_i\), ou \(1.15\) unidades na direção do autovetor \(\mathbf{e}_1\) e \(0.65\) unidades na direção \(\mathbf{e}_2\). A elipse de confiança de 95% é mostrada na Figura 6.1. Claramente, \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = \mathbf{0}\) não está na elipse, e concluímos que os dois métodos de fabricação de sabão produzem resultados diferentes. Parece que os dois processos produzem barras de sabão com aproximadamente a mesma suavidade (\(X_2\)), mas aquelas do segundo processo têm mais espuma (\(X_1\)).
Figura 6.1 Elipse de confiança de 95% para diferença de duas médias populacionais.
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
centroidif <- c(-1.9, 0.2)
sp <- matrix(c(2, 1,
1, 5),
nrow=2, byrow=TRUE)
p <- 2
n1 <- n2 <- 50
alfa <- 0.05
T2crit <- ((n1+n2-2)*p/(n1+n2-p-1))*qf(1-alfa, p, n1+n2-p-1)
T2crit
[1] 6.244089
[1] 52.47222
F <- T2/((n1+n2-2)*p/(n1+n2-p-1))
pv <- 1-pf(F, p, n1+n2-p-1)
cat("\nF(",p,", ",n1+n2-p-1,") = ", F,", p = ",pv, "\n", sep="")
F(2, 97) = 25.9684, p = 9.286081e-10
c <- sqrt(T2crit*(1/n1+1/n2))
car::ellipse(center=centroidif,
shape=sp,
radius=c,
fill=TRUE,
fill.alpha=0.1,
grid=FALSE,
col="black",
add=FALSE,
xlab=expression(mu[11] - mu[21]),
ylab=expression(mu[12] - mu[22]),
xlim=c(-3,1),
ylim=c(-1,2),
main="Região elíptica de confiança de 95%\nlather vs. mildness")
abline(v=0,h=0,lty=2)
points(H0[1], H0[2], pch=9, col="black")
text(H0[1], H0[2], pos=1, expression(H[0]))
a.hat_abs <- abs(solve(sp)%*%centroidif)
colnames(a.hat_abs) <- "Importância"
rownames(a.hat_abs) <- c("lather", "mildness")
print(proportions(a.hat_abs), digits=2)
Importância
lather 0.81
mildness 0.19
É possível derivar intervalos de confiança simultâneos para os componentes do vetor \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\). Esses intervalos de confiança são desenvolvidos a partir de uma consideração de todas as possíveis combinações lineares das diferenças nos vetores médios. Assume-se que as populações multivariadas parentais são normais com uma covariância comum \(\mathbf{\Sigma}\).
Resultado 6.3
\[ \text{IC}^{1 - \alpha}(\boldsymbol{\mu}_{11} - \boldsymbol{\mu}_{21})=\left[\mathbf{a}^{\prime} (\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2) \pm \\\sqrt{\dfrac{(n_1 + n_2 - 2)p}{n_1 + n_2 - p - 1} F_{p,n_1+n_2-p-1}(1-\alpha)}\sqrt{\mathbf{a}^{\prime}\left(\dfrac{1}{n_1} + \dfrac{1}{n_2}\right) \mathbf{S}_{\text{comb}}\mathbf{a}}\right] \]
Em particular,
\[ \text{IC}^{1 - \alpha}(\mu_{1i}- \mu_{2i})=\left[\bar{X}_{1i} - \bar{X}_{2i} \pm \\\sqrt{\dfrac{(n_1 + n_2 - 2)p}{n_1 + n_2 - p - 1} F_{p,n_1+n_2-p-1}(1-\alpha)}\sqrt{\left(\dfrac{1}{n_1} + \dfrac{1}{n_2}\right) S_{ii,\text{comb}}}\right]\\ i = 1,2,...,p \]
Observação. Para testar \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = 0\), a combinação linear \(\hat{\mathbf{a}}^{\prime} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)\), com vetor de coeficiente \(\hat{\mathbf{a}} \propto \mathbf{S}_{\text{comb}}^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)\), quantifica a maior diferença populacional. Ou seja, se \(T^2\) rejeitar \(H_0\), então \(\hat{\mathbf{a}}^{\prime} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)\) terá uma média diferente de zero. Frequentemente, tentamos interpretar os componentes dessa combinação linear tanto para a importância do assunto quanto para a importância estatística.
\(p = 2 \;(\text{on-, off-peak}), q = 1, g = 2 \;(\text{com, sem ar-condicionado}):n_1=45, n_2=55\)
“On-peak” refere-se aos períodos do dia em que a demanda por eletricidade é maior, geralmente durante as horas mais quentes do dia quando muitas pessoas estão usando ar condicionado, eletrodomésticos, iluminação, etc. Durante esses períodos, as tarifas de eletricidade podem ser mais altas devido à maior demanda.
“Off-peak” refere-se aos períodos de menor demanda, geralmente à noite ou nas primeiras horas da manhã, quando menos dispositivos e aparelhos estão sendo usados. Durante esses períodos, as tarifas de eletricidade podem ser mais baixas.
No contexto fornecido, o texto está comparando o consumo elétrico “on-peak” e “off-peak” entre os proprietários de casas com ar condicionado e aqueles sem ar condicionado. Isso sugere que aqueles com ar condicionado podem consumir mais eletricidade durante os períodos de pico devido ao uso do aparelho de ar condicionado.
Foram coletadas amostras de tamanhos \(n_1 = 45\) e \(n_2 = 55\) de proprietários de casas em Wisconsin com e sem ar-condicionado, respectivamente. (Dados cortesia do Laboratório Estatístico, Universidade de Wisconsin.) Foram consideradas duas medições do uso elétrico (em quilowatt-hora). A primeira é uma medida do consumo total no pico (\(X_1\)) durante julho e a segunda é uma medida do consumo total fora do pico (\(X_2\)) durante julho. As estatísticas resumidas resultantes são:
Com ar-condicionado:
\[ \bar{\mathbf{x}}_1 = \begin{bmatrix} 204.4 \\ 556.6 \end{bmatrix} \quad \mathbf{s}_1 = \begin{bmatrix} 13825.3 & 23823.4 \\ 23823.4 & 73107.4 \end{bmatrix} \quad n_1 = 45 \]
Sem ar-condicionado:
\[ \bar{\mathbf{x}}_2 = \begin{bmatrix} 130.0 \\ 355.0 \end{bmatrix} \quad \mathbf{s}_2 = \begin{bmatrix} 8632.0 & 19616.7 \\ 19616.7 & 55964.5 \end{bmatrix} \quad n_2 = 55 \]
(O consumo fora do pico é maior que o consumo no pico porque há mais horas fora do pico em um mês.)
Vamos encontrar intervalos de confiança simultâneos de 95% para as diferenças nos componentes do vetor de média.
Embora haja alguma discrepância nas variâncias amostrais, para fins ilustrativos, vamos calcular a matriz de covariância da amostra combinada. Aqui:
\[ \mathbf{s}_{\text{comb}} = \dfrac{(n_1 - 1)\mathbf{s}_1 + (n_2 - 1)\mathbf{s}_2}{n_1 + n_2 - 2}= \begin{bmatrix} 10963.7&21505.5\\ 21505.5&63661.3 \end{bmatrix} \]
e
\[ \begin{align} c^2 &= \dfrac{(n_1 + n_2 - 2)p}{n_1 + n_2 - p - 1} F_{p,n_1+n_2-p-1}(1-\alpha)\\ &=\dfrac{98\times2}{97}F_{2,97}(0.95)\\ &=2.02\times3.1\\ c^2 &=6.26 \end{align} \]
Com \(\boldsymbol{\mu}_1^{\prime} - \boldsymbol{\mu}_2^{\prime} = [\mu_{11} - \mu_{21}, \mu_{12} - \mu_{22}]\), os intervalos de confiança simultâneos de 95% para as diferenças populacionais são:
On-peak:
\[ \text{IC}^{95\%}(\mu_{11}- \mu_{21})=\left[ 204.4 - 130.0 \pm \sqrt{6.26} \sqrt{\left(\dfrac{1}{45} + \frac{1}{55}\right)10963.7}\right]= [21.7,127.1] \]
Off-peak:
\[ \text{IC}^{95\%}(\mu_{12}- \mu_{22})=\left[ 556.6 - 355.0 \pm \sqrt{6.26} \sqrt{\left(\dfrac{1}{45} + \frac{1}{55}\right)63661.3}\right]= [74.7,328.5] \]
Concluímos que há uma diferença no consumo elétrico entre aqueles com ar-condicionado e aqueles sem. Esta diferença é evidente tanto no consumo no pico quanto fora do pico.
A região elíptica de confiança de 95% para \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\) é determinada a partir dos pares de autovalores e autovetores:
\[ \lambda_1 = 71323.5, \mathbf{e}_1 = \begin{bmatrix} 0.336 \\ 0.942 \end{bmatrix} \]
e
\[ \lambda_2 = 3301.5, \mathbf{e}_2 = \begin{bmatrix} 0.942 \\ -0.336 \end{bmatrix} \]
Dado que:
\[ \sqrt{\lambda_1 \left( \dfrac{1}{45} + \dfrac{1}{55} \right) c^2} = \sqrt{71323.5 \left( \frac{1}{45} + \frac{1}{55} \right) 6.26} = 134.3 \]
e
\[ \sqrt{\lambda_2 \left( \dfrac{1}{45} + \dfrac{1}{55} \right) c^2} = \sqrt{3301.5 \left( \dfrac{1}{45} + \dfrac{1}{55} \right) 6.26} = 28.9 \]
Obtemos a região elíptica de confiança de 95% para \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\), ilustrada na Figura 6.2.
Como a região elíptica de confiança para a diferença nas médias não cobre \(\mathbf{0}^{\prime} = [0, 0]\), a estatística de teste \(T^2\) rejeita \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = 0\) no nível de 5%.
Figura 6.2: Elipse de confiança de 95% para diferença dos vetores de médias populacionais.
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
centroidif <- c(204.4 - 130.0,
556.6 - 355.0)
sp <- matrix(c(10963.7, 21505.5,
21505.5, 63661.3),
nrow=2, byrow=TRUE)
p <- 2
n1 <- 45
n2 <- 55
alfa <- 0.05
T2crit <- ((n1+n2-2)*p/(n1+n2-p-1))*qf(1-alfa, p, n1+n2-p-1)
T2crit
[1] 6.244089
[1] 16.0662
[1] TRUE
F <- T2/((n1+n2-2)*p/(n1+n2-p-1))
pv <- 1-pf(F, p, n1+n2-p-1)
cat("\nF(",p,", ",n1+n2-p-1,") = ", F,", p = ",pv, "\n", sep="")
F(2, 97) = 7.951131, p = 0.000634382
c <- sqrt(T2crit*(1/n1+1/n2))
car::ellipse(center=centroidif,
shape=sp,
radius=c,
fill=TRUE,
fill.alpha=0.1,
grid=FALSE,
col="black",
add=FALSE,
xlab=expression(mu[11] - mu[21]),
ylab=expression(mu[12] - mu[22]),
xlim=c(-5,330),
ylim=c(-5,330),
main="Região elíptica de confiança de 95%\n on- vs. off-peak")
abline(v=0,h=0,lty=2)
points(H0[1], H0[2], pch=9, col="black")
text(H0[1], H0[2], pos=3, expression(H[0]))
a.hat_abs <- abs(solve(sp)%*%centroidif)
colnames(a.hat_abs) <- "Importância"
rownames(a.hat_abs) <- c("on-peak", "off-peak")
print(proportions(a.hat_abs), digits=2)
Importância
on-peak 0.4
off-peak 0.6
\[\Diamond\]
Os intervalos de confiança simultâneos t de Bonferroni de \(100(1 - \alpha)\)% para as \(p\) diferenças das médias populacionais são
\[ \text{IC}^{1 - \alpha}(\mu_{1i}- \mu_{2i})=\left[\bar{X}_{1i} - \bar{X}_{2i} \pm t_{n_1+n_2-2}\left(\dfrac{\alpha}{2p}\right)\sqrt{\left(\dfrac{1}{n_1} + \dfrac{1}{n_2}\right) S_{ii,\text{comb}}}\right]\\ i = 1,2,...,p \]
Quando \(\mathbf{\Sigma}_1 \neq \mathbf{\Sigma}_2\), somos incapazes de encontrar uma medida de “distância” como \(T^2\), cuja distribuição não depende dos desconhecidos \(\mathbf{\Sigma}_1\) e \(\mathbf{\Sigma}_2\). O teste de Bartlett é usado para testar a igualdade de \(\mathbf{\Sigma}_1\) e \(\mathbf{\Sigma}_2\) em termos de variâncias generalizadas. Infelizmente, as conclusões podem ser seriamente enganosas quando as populações são não normais. Não normalidade e covariâncias desiguais não podem ser separadas com o teste de Bartlett. (Veja também a Seção 6.6.) Um método de teste da igualdade de duas matrizes de covariância que é menos sensível à suposição de multinormalidade foi proposto por Tiku e Balakrishnan. No entanto, mais experiência prática é necessária com este teste antes de podermos recomendá-lo incondicionalmente.
Sugerimos, sem muito suporte factual, que qualquer discrepância da ordem de \(\mathbf{\sigma}_{1,ii} = 4\mathbf{\sigma}_{2,ii}\) ou vice-versa, provavelmente é séria. Isso é verdade no caso univariado. O tamanho das discrepâncias que são críticas na situação multivariada provavelmente depende, em grande parte, do número de variáveis \(p\).
Transformação não linear pode melhorar as coisas quando as variâncias marginais são bastante diferentes. No entanto, para \(n_1 - p-g\ge30\) e \(n_2 -p-g\ge30\), podemos evitar as complexidades devido a matrizes de covariância desiguais.
Resultado 6.4. Seja o tamanho da amostra tal que \(n_1 - p - g\ge30\) e \(n_2 - p - g\ge30\). Então, a região elíptica de confiança aproximada de \(100(1- \alpha)\%\) para \(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2\) é dado por todos \(\boldsymbol{\mu}\) que satisfazem
\[ (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2 - (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2))^{\prime} \left( \dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2 \right)^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2 - (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)) \leq \chi^2_p(1-\alpha) \]
Além disso, intervalos de confiança simultâneos de \(100(1 - \alpha)\%\) para todas as combinações lineares \(\mathbf{a}^{\prime}(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)\) são fornecidos por
\[ \text{IC}^{1-\alpha}(\mathbf{a}^{\prime}(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2))= \left[\mathbf{a}^{\prime}(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2) \pm \sqrt{\chi^2_p(1-\alpha)} \sqrt{\mathbf{a}^{\prime} \left( \dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2 \right) \mathbf{a}}\right] \]
Prova.
A partir das equações (6-22) e (3-9), temos:
\[ \mathbb{E}(\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2) = \mathbb{E}(\bar{\mathbf{X}}_1) - \mathbb{E}(\bar{\mathbf{X}}_2)= \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 \]
\[ \mathbb{C}(\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2) = \mathbb{C}(\bar{\mathbf{X}}_1) + \mathbb{C}(\bar{\mathbf{X}}_2) = \dfrac{1}{n_1} \mathbf{S}_1 + \dfrac{1}{n_2} \mathbf{S}_2 \]
Pelo teorema do limite central, \(\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2 \underset{a}{\sim} \mathcal{N}_p\left(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2, \dfrac{1}{n_1} \mathbf{\Sigma}_1 + \dfrac{1}{n_2} \mathbf{\Sigma}_2\right)\).
\[\Diamond\]
Observação: Se \(n_1 = n_2 = n\), então \(\dfrac{n - 1}{n + n - 2} = \dfrac{1}{2}\), assim
\[ \begin{align} \dfrac{1}{n_1} \mathbf{S}_1 + \dfrac{1}{n_2} \mathbf{S}_2 &= \dfrac{1}{2} (\mathbf{S}_1 + \mathbf{S}_2)\\ &= \dfrac{(n - 1)\mathbf{S}_1 + (n - 1)\mathbf{S}_2}{n + n - 2}\left(\dfrac{1}{n}+\dfrac{1}{n}\right) \\ \dfrac{1}{n_1} \mathbf{S}_1 + \dfrac{1}{n_2} \mathbf{S}_2 &= \mathbf{S}_\text{comb}\left(\dfrac{1}{n}+\dfrac{1}{n}\right) \end{align} \]
Com tamanhos de amostra iguais, o procedimento para amostra grande é essencialmente o mesmo que o procedimento baseado na matriz de covariância agrupada. (Veja o Resultado 6.2.) Em uma dimensão, é bem conhecido que o efeito de heterocedasticidade é menor quando há balanceamento, i.e., \(n_1 = n_2\), e é maior quanto maior é o desbalanceamento.
Vamos analisar os dados de consumo elétrico discutidos no Exemplo 6.4 usando a abordagem de amostras grandes. Primeiro, calculamos:
\[ \begin{align} \dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2 &= \dfrac{1}{45} \begin{bmatrix} 13825.3 & 23823.4 \\ 23823.4 & 73107.4 \\ \end{bmatrix} + \dfrac{1}{55} \begin{bmatrix} 8632.0 & 19616.7 \\ 19616.7 & 55964.5 \\ \end{bmatrix}\\ \dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2&= \begin{bmatrix} 464.17 & 886.08 \\ 886.08 & 2642.15 \\ \end{bmatrix} \end{align} \]
Os intervalos de confiança simultâneos de 95% para as combinações lineares
\[ \mathbf{a}^{\prime}_1(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2) = [1\;0]\begin{bmatrix} \mu_{11}-\mu_{21}\\ \mu_{12}-\mu_{22} \end{bmatrix} =\mu_{11} - \mu_{21} \]
e
\[ \mathbf{a}^{\prime}_2(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2) = [0\;1]\begin{bmatrix} \mu_{11}-\mu_{21}\\ \mu_{12}-\mu_{22} \end{bmatrix}=\mu_{12} - \mu_{22} \]
são (veja o Resultado 6.4):
\[ \begin{align} \text{IC}^{95\%}(\mu_{11} - \mu_{21}) &=\left[74.4 \pm \sqrt{5.99} \sqrt{464.17}\right]\\ \text{IC}^{95\%}(\mu_{11} - \mu_{21})&=[21.7,127.1]\\\\ \text{IC}^{95\%}(\mu_{12} - \mu_{22}) &=\left[201.6 \pm \sqrt{5.99} \sqrt{2642.15}\right]\\ \text{IC}^{95\%}(\mu_{12} - \mu_{22})&=[75.8,327.4] \end{align} \]
Note que esses intervalos diferem desprezivelmente dos intervalos no Exemplo 6.4, em que o procedimento de agrupamento foi empregado. A estatística \(T^2\) para testar \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = 0\) é:
\[ T^2 = (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)^{\prime} \left( \dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2 \right)^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2) = 16.07 \]
Para \(\alpha = .05\), o valor crítico é \(\chi^2_2(0.95) = 5.99\) e, como \(T^2 = 16.07 > \chi^2_2(0.95) = 5.99\), rejeitamos \(H_0\).
A combinação linear mais crítica que leva à rejeição de \(H_0\) tem o vetor de coeficiente
\[ \hat{\mathbf{a}} \propto \left( \dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2 \right)^{-1} (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2) = \begin{bmatrix} 0.04 \\ 0.06 \\ \end{bmatrix} \]
A diferença no consumo elétrico “on-peak” entre aqueles com ar condicionado e aqueles sem contribui mais do que a diferença correspondente no consumo “off-peak” para a rejeição de \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = \mathbf{0}\).
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
centroidif <- c(204.4 - 130.0,
556.6 - 355.0)
sdif <- matrix(c(464.17, 886.08,
886.08, 2642.15),
nrow=2, byrow=TRUE)
p <- 2
n1 <- 45
n2 <- 55
alfa <- 0.05
X2crit <- qchisq(1-alfa, p)
X2crit
[1] 5.991465
[1] 15.6585
[1] TRUE
pv <- formatC(1-pchisq(T2, p),
format="e", digits=2)
cat("\nX^2(",p,") = ", T2,", p = ",pv, "\n", sep="")
X^2(2) = 15.6585, p = 3.98e-04
c <- sqrt(X2crit)
car::ellipse(center=centroidif,
shape=sdif,
radius=c,
fill=TRUE,
fill.alpha=0.1,
grid=FALSE,
col="black",
add=FALSE,
xlab=expression(mu[11] - mu[21]),
ylab=expression(mu[12] - mu[22]),
xlim=c(-5,130),
ylim=c(-5,330),
main="Região elíptica de confiança de 95%\n on- vs. off-peak")
abline(v=0,h=0,lty=2)
points(H0[1], H0[2], pch=9, col="black")
text(H0[1], H0[2], pos=3, expression(H[0]))
a.hat_abs <- abs(solve(sp)%*%centroidif)
colnames(a.hat_abs) <- "Importância"
rownames(a.hat_abs) <- c("on-peak", "off-peak")
print(proportions(a.hat_abs), digits=4)
Importância
on-peak 0.3965
off-peak 0.6035
\[\Diamond\]
A estatística similar ao \(T^2\) que é menos sensível a observações atípicas para amostras pequenas e de tamanho moderado foi desenvolvida por Tiku e Singh [24]. No entanto, se o tamanho da amostra for moderado a grande, \(n_1-p-g>12\) e \(n_2-p-g>12\), \(T^2\) de Hotelling é notavelmente não afetado por pequenos desvios da normalidade e/ou pela presença de alguns outliers.
É possível testar \(H_0: \boldsymbol{\mu}_1 - \boldsymbol{\mu}_2 = \mathbf{0}\) quando as matrizes de covariância da população são desiguais, mesmo que os dois tamanhos de amostra não sejam grandes, desde que as duas populações sejam multivariadas normais. Essa situação é frequentemente chamada de problema multivariado de Behrens-Fisher. O resultado exige que \(n_1 >p\) e \(n_2 >p\). A abordagem depende de uma aproximação para a distribuição da estatística
\[ T^2 = (\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2-(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2))^{\prime} \left( \dfrac{1}{n_1} \mathbf{S}_1 + \dfrac{1}{n_2} \mathbf{S}_2 \right)^{-1} (\bar{\mathbf{X}}_1 - \bar{\mathbf{X}}_2-(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2))\underset{a}{\sim}T^2_{p,\nu}\tag{6-27} \]
que é idêntica à estatística de amostra grande no Resultado 6.4. No entanto, em vez de usar a aproximação qui-quadrado para obter o valor crítico para testar \(H_0\), a aproximação recomendada para amostras menores (veja [15] e [19]) é dada por
\[ T^2_{p,\nu} = \dfrac{\nu p}{\nu-p+1}F_{p,\nu-p+1}\tag{6-28} \]
em que o número de graus de liberdade \(\nu\) de \(T^2\) é estimado a partir das matrizes de covariância amostrais usando a relação:
\[ \hat{\nu}=\dfrac{p(1+p)}{ \sum_{i=1}^{2}{\dfrac{1}{n_i}\left( \text{tr}\left(\left(\dfrac{1}{n_i}\mathbf{S}_i\left( \dfrac{1}{n_1}\mathbf{S}_1 + \dfrac{1}{n_2}\mathbf{S}_2 \right)^{-1}\right)^2\right) + \left(\text{tr}\left(\dfrac{1}{n_i}\mathbf{S}_i\left( \dfrac{1}{n_1}\mathbf{S}_1 + \dfrac{1}{n_2}\mathbf{S}_2 \right)^{-1}\right)\right)^2\right)}} \tag{6-29} \]
em que \(\min(n_1, n_2) \le \nu \le n_1+ n_2\). Esta aproximação reduz-se à solução usual de Welch para o problema de Behrens-Fisher (heterocedasticidade) no caso univariado (\(p = 1\)).
Para populações normais, a aproximação para a distribuição de \(T^2\) dada por (6-28) e (6-29) geralmente fornece resultados razoáveis.
O teste \(T^2\) de Johansen é robusto à falta de homogeneidade de inclinação (delineamento pré-pós) e homocedasticidade multivariada.
# Outros métodos robustos
smallX <- matrix(rnorm(10*3),ncol=3)
smallY <- matrix(rnorm(10*3),ncol=3)
SHT::mean2.1980Johansen(smallX, smallY)
Two-sample Test for Multivariate Means by Johansen (1980)
data: smallX and smallY
T2 = 0.39345, p-value = 0.9506
alternative hypothesis: true means are different.
Two-sample Test for Multivariate Means by Yao (1965).
data: smallX and smallY
T2 = 0.39345, p-value = 0.9495
alternative hypothesis: true means are different.
Hotelling's T-squared Test for Independent Samples with Equal
Covariance Assumption.
data: smallX and smallY
T2 = 0.39345, p-value = 0.9491
alternative hypothesis: true means are different.
Johansen, S (1980) The Welch-James approximation to the distribution of the residual sum of squares in a weighted linear regression. Biometrika 67(1). doi:10.1093/biomet/67.1.85.
J. Zhang, “An Approximate Hotelling T2-Test for Heteroscedastic One-Way MANOVA,” Open Journal of Statistics, Vol. 2 No. 1, 2012, pp. 1-11. doi: 10.4236/ojs.2012.21001. https://www.scirp.org/journal/paperinformation.aspx?paperid=17138
Gokpınar, E et al. (?) A Computational Approach Test for the Equality of Two Multivariate Normal Mean Vectors under Heterogeneity of Covariance Matrices: https://www.ine.pt/revstat/pdf/AComputationalApproachTest.pdf
Embora os tamanhos das amostras sejam bastante grandes para os dados de consumo elétrico no Exemplo 6.4, usamos esses dados e os cálculos do Exemplo 6.5 para ilustrar os cálculos que levam à distribuição aproximada de \(T^2\) quando as matrizes de covariância da população são diferentes.
Primeiro, calculamos:
\[ \dfrac{1}{n_1} \mathbf{s}_1 = \dfrac{1}{45} \begin{bmatrix} 13825.2 & 23823.4 \\ 23823.4 & 73107.4 \\ \end{bmatrix} = \begin{bmatrix} 307.227 & 529.409 \\ 529.409 & 1624.609 \\ \end{bmatrix} \]
e
\[ \dfrac{1}{n_2} \mathbf{s}_2 = \dfrac{1}{55} \begin{bmatrix} 8632.0 & 19616.7 \\ 19616.7 & 55964.5 \\ \end{bmatrix} = \begin{bmatrix} 156.945 & 356.667 \\ 356.667 & 1017.536 \\ \end{bmatrix} \]
E usando um resultado do Exemplo 6.5,
\[ \left(\dfrac{1}{n_1} \mathbf{s}_1 + \dfrac{1}{n_2} \mathbf{s}_2\right)^{-1} = 10^{-4} \begin{bmatrix} 59.874 & -20.080 \\ -20.080 & 10.519 \\ \end{bmatrix} \]
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
centroidif <- c(204.4 - 130.0,
556.6 - 355.0)
s1 <- matrix(c(13825.2, 23823.4,
23823.4, 73107.4),
nrow=2, byrow=TRUE)
s2 <- matrix(c(8632.0, 19616.7,
19616.7, 55964.5),
nrow=2, byrow=TRUE)
p <- 2
n1 <- 45
n2 <- 55
alfa <- 0.05
compute_nu_hat <- function(s1, s2, n1, n2, p) {
common_part <- (1/n1)*s1 + (1/n2)*s2
sum_term <- function(n, s) {
A <- (1/n)*s %*% solve(common_part)
tr_A <- sum(diag(A))
return(tr_A^2 + (sum(diag(A %*% A))) )
}
sum_val <- (1/n1)*sum_term(n1, s1) + (1/n2)*sum_term(n2, s2)
nu_hat <- (p*(1+p)) / sum_val
return(nu_hat)
}
nu_hat <- compute_nu_hat(s1, s2, n1, n2, p)
T2crit <- (nu_hat*p/(nu_hat-p+1))*qf(1-alfa, p, nu_hat-p+1)
T2crit
[1] 6.313406
H0 <- rep(0, p)
T2 <- as.numeric(t(centroidif-H0)%*%solve((1/n1)*s1 + (1/n2)*s2)
%*%(centroidif-H0))
T2
[1] 15.65853
[1] TRUE
F <- T2/(nu_hat*p/(nu_hat-p+1))
pv <- 1-pf(F, p, nu_hat-p+1)
cat("\nF(",p,", ",nu_hat-p+1,") = ", F,", p = ",pv, "\n", sep="")
F(2, 76.59217) = 7.728364, p = 0.0008763188
c <- sqrt(T2crit)
car::ellipse(center=centroidif,
shape=(1/n1)*s1 + (1/n2)*s2,
radius=c,
fill=TRUE,
fill.alpha=0.1,
grid=FALSE,
col="black",
add=FALSE,
xlab=expression(mu[11] - mu[21]),
ylab=expression(mu[12] - mu[22]),
xlim=c(-5,130),
ylim=c(-5,330),
main="Região elíptica de confiança de 95% (Welch)\n on- vs. off-peak")
abline(v=0,h=0,lty=2)
points(H0[1], H0[2], pch=9, col="black")
text(H0[1], H0[2], pos=3, expression(H[0]))
a.hat_abs <- abs(solve((1/n1)*s1 + (1/n2)*s2)%*%centroidif)
colnames(a.hat_abs) <- "Importância"
rownames(a.hat_abs) <- c("on-peak", "off-peak")
print(proportions(a.hat_abs), digits=2)
Importância
on-peak 0.39
off-peak 0.61
Frequentemente, é necessário comparar mais do que duas populações. Amostras aleatórias coletadas de cada uma das \(g\) populações (condições experimentais independentes) são organizadas como:
População 1: \(\mathbf{X}_{11}, \mathbf{X}_{12}, \ldots, \mathbf{X}_{1n_1}\)
População 2: \(\mathbf{X}_{21}, \mathbf{X}_{22}, \ldots, \mathbf{X}_{2n_2}\)
População g: \(\mathbf{X}_{g1}, \mathbf{X}_{g2}, \ldots, \mathbf{X}_{gn_g}\)
(6-31)
MANOVA unifatorial independente (one-way MANOVA) é usado primeiramente para investigar se os vetores de média populacional são os mesmos e, se não forem, quais componentes do vetores de média diferem significantemente.
A condição 3 pode ser relaxada recorrendo ao teorema do limite central (Resultado 4.13) quando os tamanhos de amostra $n_$ são grandes, i.e.,
\[ n_\ell -pq-g\ge 30 \\ \ell = 1, 2, \ldots, g \]
Na situação univariada, as suposições são de que \(\left\{X_{\ell i}\right\}_{i=1}^{n_\ell } \sim \mathcal{N}\text{IID}\left(\mu_{\ell },\sigma^2\right)\), \(\ell = 1, 2, \ldots, g\), e que as \(g\) amostras aleatórias \(\left\{\left\{X_{\ell i}\right\}_{i=1}^{n_\ell }\right\}_{\ell =1}^{g}\) são independentes. Embora a hipótese nula de igualdade de médias possa ser formulada como \(\mu_{\ell } = \mu_2 = \cdots = \mu_g\), é costumeiro considerar $_$ como a soma de um componente de média global, como \(\mu\), e um componente devido à população específica. Por exemplo, podemos escrever \(\mu_{\ell } = \mu + (\mu_{\ell } - \mu)\) ou \(\mu_{\ell } = \mu + \tau_{\ell }\), em que \(\tau_{\ell } = \mu_{\ell } - \mu\).
As populações geralmente correspondem a diferentes conjuntos de condições experimentais e, portanto, é conveniente investigar os desvios \(\tau_{\ell }\) associados à $$-ésima população (tratamento).
A reparametrização
\[ \mu_{\ell } = \mu + \tau_{\ell } \]
em que:
leva a uma reformulação da hipótese de igualdade das médias. A hipótese nula torna-se
\[ H_0: \tau_1 = \tau_2 = \cdots = \tau_g = 0 \]
A resposta \(X_{\ell j} \sim \mathcal{N}\left(\mu + \tau_{\ell }, \sigma^2\right)\), pode ser expressa da seguinte forma forma
\[ X_{\ell j} = \mu + \tau_{\ell } + e_{\ell j} \]
em que:
em que os \(e_{\ell j} \sim \mathcal{N}\text{IID}\left(0,\sigma^2\right)\), \(j = 1, 2, \ldots, n_{\ell}\) e \(\ell = 1, 2, \ldots, g\) são independentes. Para definir de forma única os parâmetros do modelo e suas estimativas de mínimos quadrados, é costumeiro impor a restrição \(\sum_{\ell =1}^g n_\ell \tau_{\ell } = 0\).
Motivado pela decomposição anterior, a análise de variância é baseada em uma decomposição análoga das observações,
\[ x_{\ell j} = \bar{x} + (\bar{x}_{\ell } - \bar{x}) + (x_{\ell j} - \bar{x}_{\ell }) \]
em que:
Considere as seguintes amostras independentes:
População 1: 9, 6, 9
População 2: 0, 2
População 3: 3, 1, 2
Por exemplo, dado que
\[ \bar{x}_3 = \dfrac{3 + 1 + 2}{3} = 2 \quad \text{e} \quad \bar{x} = \dfrac{9 + 6 + 9 + 0 + 2 + 3 + 1 + 2}{8} = 4, \]
encontramos que
\[ \begin{align} 3 &= x_{31} \\ &= \bar{x} + (\bar{x}_3 - \bar{x}) + (x_{31} - \bar{x}_3) \\ &= 4 + (2 - 4) + (3 - 2) \\ 3&= 4 + (-2) + 1 \end{align} \]
Repetindo esta operação para cada observação, obtemos a seguinte igualdade matricial:
\[ \begin{bmatrix} 9&6&9\\0&2&&\\3&1&2 \end{bmatrix}= \begin{bmatrix} 4&4&4\\4&4&&\\4&4&4 \end{bmatrix}+ \begin{bmatrix} 4&4&4\\-3&-3&&\\-2&-2&-2 \end{bmatrix}+ \begin{bmatrix} 1&-2&1\\-1&1&&\\1&-1&0 \end{bmatrix} \]
\[ \begin{align} x_{\ell j} &= \bar{x} + (\bar{x}_{\ell } - \bar{x}) + (x_{\ell j} - \bar{x}_{\ell })\\ x_{\ell j} &= \hat{\mu}+\hat{\tau}_{\ell }+\hat{e}_{\ell j} \end{align} \]
observação = média + efeito de tratamento + resíduo
A questão da igualdade das médias é respondida avaliando se a contribuição do vetor de tratamento é grande em relação aos resíduos.
Nossas estimativas \(\hat{\tau}_{\ell }=\bar{x}_{\ell } - \bar{x}\) sempre satisfazem \(\sum_{\ell =1}^g n_\ell \hat{\tau}_{\ell } = 0\).
Sob \(H_0\), cada \(\hat{\tau}_{\ell }\) é uma estimativa nula. Se a contribuição do tratamento for grande, \(H_0\) deve ser rejeitada. O tamanho de um vetor é quantificado transformando as linhas do vetor em um único vetor e calculando seu comprimento ao quadrado. Esta quantidade é chamada de soma dos quadrados (SS). Para as observações, construímos o vetor \(\mathbf{y}^{\prime} = [9\; 6\; 9\; 0\; 2\; 3\; 1\; 2]\). Seu comprimento ao quadrado é:
\[ \text{SS}_{\text{obs}} = 9^2 + 6^2 + 9^2 + 0^2 + 2^2 + 3^2 + 1^2 + 2^2 = 216 \]
Similarmente,
\[ \begin{align} \text{SS}_{\text{mean}} &= 4^2 + 4^2 + 4^2 + 4^2 + 4^2 + 4^2 + 4^2 + 4^2 = 8(4^2) = 128\\\\ \text{SS}_{\text{trt}} &= 4^2 + 4^2 + 4^2 + (-3)^2 + (-3)^2 + (-2)^2 + (-2)^2 + (-2)^2 \\ \text{SS}_{\text{trt}} &= 3(4^2) + 2(-3)^2 + 3(-2)^2 = 78 \end{align} \]
e a soma dos quadrados residual é
\[ \text{SS}_{\text{res}} = 1^2 + (-2)^2 + 1^2 + (-1)^2 + 1^2 + 1^2 + (-1)^2 + 0^2 = 10 \]
As somas dos quadrados satisfazem a mesma decomposição, (6-34), que as observações. Portanto,
\[ \text{SS}_{\text{obs}} = \text{SS}_{\text{mean}} + \text{SS}_{\text{trt}} + \text{SS}_{\text{res}} \]
ou
\[ 216 = 128 + 78 + 10 \]
A divisão em somas de quadrados aloca a variabilidade nas amostras combinadas em componentes de média, tratamento e residual (erro). Uma análise de variância procede comparando os tamanhos relativos de \(\text{SS}_{\text{trt}}\) e \(\text{SS}_{\text{res}}\). Se \(H_0\) for verdadeira, as variâncias calculadas a partir de \(\text{SS}_{\text{trt}}\) e \(\text{SS}_{\text{res}}\) devem ser aproximadamente iguais.
\[\Diamond\]
A decomposição da soma dos quadrados ilustrada numericamente no Exemplo 6.7 é tão básica que o equivalente algébrico será agora desenvolvido.
Subtraindo \(\bar{x}\) de ambos os lados de (6-34) e elevando ao quadrado, obtemos:
\[ (x_{\ell j} - \bar{x})^2 = ( \bar{x}_{\ell } - \bar{x})^2 + (x_{\ell j} - \bar{x}_{\ell })^2 + 2( \bar{x}_{\ell } - \bar{x})(x_{\ell j} - \bar{x}_{\ell }) \]
Podemos somar ambos os lados sobre \(j\), notando que \(\sum_{i=1}^{n_\ell }{(x_{\ell j} - \bar{x}_{\ell })} = 0\), e obtemos:
\[ \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x})^2 = n_{\ell }(\bar{x}_{\ell } - \bar{x})^2 + \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x}_{\ell })^2 \]
Em seguida, somando ambos os lados sobre $$, obtemos:
\[ \begin{align} \sum_{\ell =1}^{g} \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x})^2 &= \sum_{\ell =1}^{g} n_{\ell }(\bar{x}_{\ell } - \bar{x})^2 + \sum_{\ell =1}^{g} \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x}_{\ell })^2 \\ \text{SS}_{\text{tot}}&=\text{SS}_{\text{trt}}+\text{SS}_{\text{res}} \end{align} \tag{6-35} \]
Em que:
ou
\[ \begin{align} \sum_{\ell =1}^{g} \sum_{j=1}^{n_{\ell }} x_{\ell j}^2 &= \bar{x}^2 \sum_{\ell=1}^{g}{n_\ell} + \sum_{\ell =1}^{g} n_{\ell }(\bar{x}_{\ell } - \bar{x})^2 + \sum_{\ell =1}^{g} \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x}_{\ell })^2 \\ \text{SS}_{\text{obs}}&=\text{SS}_{\text{mean}}+\text{SS}_{\text{trt}}+\text{SS}_{\text{res}} \end{align} \tag{6-36} \]
Em que:
Ao estabelecer (6-36), verificamos que as matrizes representando a média, os efeitos do tratamento e os resíduos são ortogonais. Ou seja, essas matrizes, consideradas como vetores, são perpendiculares independentemente do vetor de observação
\[ \mathbf{y}^{\prime} = [x_{11}, \ldots, x_{1n_1}, x_{21}, \ldots, x_{2n_2},x_{g1}, \ldots, x_{gn_g}] \]
[,1]
[1,] 0
[,1]
[1,] 0
[,1]
[1,] 0
Consequentemente, poderíamos obter \(\text{SS}_{\text{res}}\) por subtração, sem ter que calcular os resíduos individuais, porque
\[ \text{SS}_{\text{res}} = \text{SS}_{\text{obs}} - \text{SS}_{\text{mean}} - \text{SS}_{\text{trt}} \]
No entanto, isso é uma economia falsa porque os gráficos dos resíduos fornecem verificações das suposições do modelo.
As representações vetoriais das matrizes envolvidas na decomposição (6-34) também têm interpretações geométricas que fornecem os graus de liberdade. Para um conjunto arbitrário de observações, o vetor de observação \(\mathbf{y}\) pode estar em qualquer lugar em \(n = \sum_{\ell =1}^{g}{n_\ell }\) dimensões; o vetor médio \(\bar{x}\mathbf{1} = [\bar{x}, \ldots, \bar{x}]^{\prime}\) deve estar ao longo da linha equiangular de \(\mathbf{1}\), e o vetor de efeito de tratamento
reside no hiperplano de combinações lineares dos \(g\) vetores \(\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_g\). Como \(\mathbf{1} = \mathbf{u}_1 + \mathbf{u}_2 + \cdots + \mathbf{u}_g\), o vetor médio também reside neste hiperplano e é sempre perpendicular ao vetor de tratamento. (Consulte o Exercício 6.10). Assim, o vetor médio tem a liberdade de residir em qualquer lugar ao longo da linha equiangular unidimensional, e o vetor de tratamento tem a liberdade de residir em qualquer lugar nas outras \(g - 1\) dimensões. O vetor de resíduo, \(\hat{\mathbf{e}} = \mathbf{y} - \bar{x}\mathbf{1} - ((\bar{x}_i - \bar{x})\mathbf{u}_j + \cdots + (\bar{x}_g - \bar{x})\mathbf{u}_g)\), é perpendicular a ambos, o vetor de média e o vetor de efeito de tratamento, e tem a liberdade de residir em qualquer lugar no subespaço de dimensão \(n - (g -1) - 1 = n - g\), que é perpendicular ao seu hiperplano.
Em resumo, atribuímos 1 grau de liberdade para \(\text{SS}_{\text{mean}}\), \(g - 1\) graus de liberdade para \(\text{SS}_{\text{trt}}\), e \(n - g = \sum_{\ell =1}^{g}{n_\ell } - g\) graus de liberdade para \(\text{SS}_{\text{res}}\). O número total de graus de liberdade é \(\sum_{\ell =1}^{g}{n_\ell }\). Alternativamente, recorrendo à teoria da distribuição univariada, descobrimos que estes são os graus de liberdade para as distribuições qui-quadrado associadas com as respectivas somas dos quadrados. Os cálculos das somas dos quadrados e os graus de liberdade associados são convenientemente resumidos por uma tabela ANOVA.
Fonte de variação | Soma dos Quadrados (SS) | Graus de Liberdade (df) |
---|---|---|
Tratamentos | \(\text{SS}_{\text{trt}} = \sum_{\ell =1}^{g} n_{\ell }(\bar{x}_{\ell } - \bar{x})^2\) | \(g - 1\) |
Resíduo | \(\text{SS}_{\text{res}} = \sum_{\ell =1}^{g} \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x}_{\ell })^2\) | \(\sum_{\ell =1}^{g}{n_\ell } - g\) |
Total | \(\text{SS}_{\text{tot}} = \sum_{\ell =1}^{g} \sum_{j=1}^{n_{\ell }} (x_{\ell j} - \bar{x})^2\) | \(\sum_{\ell =1}^{g}{n_\ell }-1\) |
O teste F usual rejeita \(H_0: \tau_1 = \tau_2 = \cdots = \tau_g = 0\) no nível \(\alpha\) se
\[ F=\dfrac{\dfrac{\text{SS}_{\text{trt}}}{g-1}}{\dfrac{\text{SS}_{\text{res}}}{\sum_{\ell =1}^{g}{n_\ell } - g}} > F_{g-1, \sum_{\ell =1}^{g}{n_\ell } - g}(1-\alpha) \]
Isso é equivalente a rejeitar \(H_0\) para grandes valores de \(\dfrac{\text{SS}_{\text{trt}}}{\text{SS}_{\text{res}}}\) ou para grandes valores de \(1 + \dfrac{\text{SS}_{\text{trt}}}{\text{SS}_{\text{res}}}\).
A estatística apropriada para uma generalização multivariada rejeita \(H_0\) para pequenos valores do recíproco
\[ \dfrac{1}{1 + \dfrac{\text{SS}_{\text{trt}}}{\text{SS}_{\text{res}}}} = \dfrac{\text{SS}_{\text{res}}}{\text{SS}_{\text{res}} + \text{SS}_{\text{trt}}} \tag{6-37} \]
Usando as informações do Exemplo 6.7, temos a seguinte tabela ANOVA:
Fonte de Variação | Soma dos Quadrados | Graus de Liberdade |
---|---|---|
Tratamentos | \(\text{SS}_\text{trt} = 78\) | \(g - 1 = 3 - 1 = 2\) |
Resíduo | \(\text{SS}_\text{res} = 10\) | \(n - g = (3 + 2 + 3) - 3 = 5\) |
Total | \(\text{SS}_\text{tot} = 88\) | \(\sum_{\ell =1}^{g}{n_\ell }-1=7\) |
Consequentemente,
\[ F=\dfrac{\dfrac{\text{SS}_\text{trt}}{g - 1}}{\dfrac{\text{SS}_\text{res}}{n - g}} \sim F_{g-1,n-g} \]
\[ F=\dfrac{\dfrac{\text{SS}_\text{trt}}{g - 1}}{\dfrac{\text{SS}_\text{res}}{n - g}} = \dfrac{\dfrac{78}{2}}{\dfrac{10}{5}}=19.5 \]
Como \(F = 19.5 > F_{2,5}(0.99) = 13.27\), rejeitamos \(H_0: \tau_1 = \tau_2 = \tau_3 = 0\) (sem efeito de tratamento) no nível de significância de 1%.
Por que \(F=(\text{SS}_{\text{trt}}/(g - 1))/(\text{SS}_\text{res}/(n - g)) \sim F_{g-1,n-g}\)?
A resposta é Teorema de Cochran.
O Teorema de Cochran é uma ferramenta fundamental na análise de variância (ANOVA). O teorema fornece a distribuição de somas de quadrados de certos tipos de combinações lineares de variáveis aleatórias normais independentes. Este teorema é frequentemente utilizado em estatísticas para testar hipóteses sobre os parâmetros de populações.
O artigo seminal relacionado ao Teorema de Cochran é:
Outros artigos e referências são:
Kutner, MH et al. (2005) Applied linear statistical model. 5th ed. NY: McGraw-Hill/Irwin, p. 699
Saw, SLC (1992) ANOVA Sums of Squares as Quadratic Forms. The American Statistician 46(4): 288-90.
Tian, Y & Styan, GPH (2006) Cochran’s statistical theorem revisited. Journal of Statistical Planning and Inference 136(8): 2659-67. https://doi.org/10.1016/j.jspi.2004.09.016.
Proof: F-test for main effect in one-way analysis of variance
n1 <- 3
n2 <- 2
n3 <- 3
n <- n1 + n2 + n3
Grupo <- factor(c(rep("1", n1), rep("2", n2), rep("3", n3)))
y <- c(9, 6, 9, 0, 2, 3, 1, 2)
x <- data.frame(Grupo, y)
print(x)
Grupo y
1 1 9
2 1 6
3 1 9
4 2 0
5 2 2
6 3 3
7 3 1
8 3 2
[,1]
[1,] 88
[1] 88
# Soma dos Quadrados do Modelo
x_bar_g <- rep(aggregate(y~Grupo,
FUN="mean",
data=x)$y,
times=c(n1,n2,n3))
SS_trt <- crossprod(x_bar_g - mean(y))
# Soma dos Quadrados dos Resíduos
SS_res <- SS_tot - SS_trt
# Grupos
g <- length(unique(Grupo))
n <- length(y)
# Graus de liberdade
df_trt <- g - 1 # tratamentos
df_res <- n - g # erro
df_tot <- n - 1 # total
# Médias dos Quadrados
MS_trt <- SS_trt / df_trt
MS_res <- SS_res / df_res
# Estatística F
F_val <- MS_trt / MS_res
# Valor p
p_val <- formatC(1 - pf(F_val, df_trt, df_res),
format="e", digits=2)
# Tabela ANOVA
cat("\nAnova Table (Type II tests)\n")
Anova Table (Type II tests)
Response: y
anova_table <- data.frame(Source = c("Grupo", "Residuals", "Total"),
df = c(df_trt, df_res, df_tot),
Sum_Sq = c(SS_trt, SS_res, SS_tot),
Mean_Sq = c(MS_trt, MS_res, ""),
F = c(F_val, "", ""),
"p" = c(p_val, "", ""))
print(anova_table)
Source df Sum_Sq Mean_Sq F p
1 Grupo 2 78 39 19.5 4.35e-03
2 Residuals 5 10 2
3 Total 7 88
# Solução model.matrix, HIJ e OLS
# Kutner et al., 2005, p. 204
# Criar a matriz de design
# One-way ANOVA (offset from reference group): https://en.wikipedia.org/wiki/Design_matrix
X <- model.matrix(~Grupo,
data=x)
X
(Intercept) Grupo2 Grupo3
1 1 0 0
2 1 0 0
3 1 0 0
4 1 1 0
5 1 1 0
6 1 0 1
7 1 0 1
8 1 0 1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$Grupo
[1] "contr.treatment"
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0
[5,] 0 0 0 0 1 0 0 0
[6,] 0 0 0 0 0 1 0 0
[7,] 0 0 0 0 0 0 1 0
[8,] 0 0 0 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 1 1 1 1 1 1 1
[2,] 1 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 1 1 1
[4,] 1 1 1 1 1 1 1 1
[5,] 1 1 1 1 1 1 1 1
[6,] 1 1 1 1 1 1 1 1
[7,] 1 1 1 1 1 1 1 1
[8,] 1 1 1 1 1 1 1 1
1 2 3 4 5 6 7 8
1 0.33 0.33 0.33 0.0 0.0 0.00 0.00 0.00
2 0.33 0.33 0.33 0.0 0.0 0.00 0.00 0.00
3 0.33 0.33 0.33 0.0 0.0 0.00 0.00 0.00
4 0.00 0.00 0.00 0.5 0.5 0.00 0.00 0.00
5 0.00 0.00 0.00 0.5 0.5 0.00 0.00 0.00
6 0.00 0.00 0.00 0.0 0.0 0.33 0.33 0.33
7 0.00 0.00 0.00 0.0 0.0 0.33 0.33 0.33
8 0.00 0.00 0.00 0.0 0.0 0.33 0.33 0.33
SS_trt <- t(y)%*%(H-J/n)%*%y
SS_res <- t(y)%*%(I-H)%*%y
SS_tot <- t(y)%*%(I-J/n)%*%y
df_trt <- sum(diag(H-J/n))
df_res <- sum(diag(I-H))
df_tot <- sum(diag(I-J/n))
# Médias dos Quadrados
MS_trt <- SS_trt / df_trt
MS_res <- SS_res / df_res
# Estatística F
F_val <- MS_trt / MS_res
# Valor p
p_val <- formatC(1 - pf(F_val, df_trt, df_res),
format="e", digits=2)
# Tabela ANOVA
cat("\nAnova Table (Type II tests)\n")
Anova Table (Type II tests)
Response: y
anova_table <- data.frame(Source = c("Grupo", "Residuals", "Total"),
df = c(df_trt, df_res, df_tot),
Sum_Sq = c(SS_trt, SS_res, SS_tot),
Mean_Sq = c(MS_trt, MS_res, ""),
F = c(F_val, "", ""),
"p" = c(p_val, "", ""))
print(anova_table)
Source df Sum_Sq Mean_Sq F p
1 Grupo 2 78 39 19.5 4.35e-03
2 Residuals 5 10 2
3 Total 7 88
# Solução model.matrix e OLS
# https://genomicsclass.github.io/book/pages/expressing_design_formula.html
n1 <- 3
n2 <- 2
n3 <- 3
Grupo <- factor(c(rep("1", n1), rep("2", n2), rep("3", n3)))
y <- c(9, 6, 9, 0, 2, 3, 1, 2)
x <- data.frame(Grupo, y)
# One-way ANOVA (offset from reference group): https://en.wikipedia.org/wiki/Design_matrix
X <- model.matrix(~Grupo,
data=x)
# Computar as somas de quadrados
coeff <- solve(t(X) %*% X) %*% t(X) %*% y # OLS
cat("\nCoefficients:\n")
Coefficients:
Estimate
(Intercept) 8
Grupo2 -7
Grupo3 -6
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8 0.82 9.8 0.00019
Grupo2 -7 1.29 -5.4 0.00289
Grupo3 -6 1.15 -5.2 0.00348
y_hat <- X%*%coeff
SS_trt <- sum((y_hat - mean(y))^2)
SS_res <- sum((y - y_hat)^2)
SS_tot <- sum((y - mean(y))^2)
R2 <- 1-SS_res/SS_tot
cat("R^2 = eta^2 omnibus = ", R2, sep="")
R^2 = eta^2 omnibus = 0.8863636
# Grupos
g <- length(unique(Grupo))
n <- length(y)
# Graus de liberdade
df_trt <- g - 1 # tratamentos
df_res <- n - g # erro
df_tot <- n - 1
# Médias das somas de quadrados
MS_trt <- SS_trt / df_trt
MS_res <- SS_res / df_res
# Estatística F
F_value <- MS_trt / MS_res
# Valor p
p_value <- formatC(1 - pf(F_value, df_trt, df_res),
format="e", digits=2)
# Tabela ANOVA
cat("\nAnova Table (Type II tests)")
Anova Table (Type II tests)
Response: y
anova_table <- data.frame(Source = c("Grupo", "Residuals", "Total"),
df = c(df_trt, df_res, df_tot),
Sum_Sq = c(SS_trt, SS_res, SS_tot),
Mean_Sq = c(MS_trt, MS_res, ""),
F = c(F_value, "", ""),
"p" = c(p_value, "", ""))
print(anova_table)
Source df Sum_Sq Mean_Sq F p
1 Grupo 2 78 39 19.5 4.35e-03
2 Residuals 5 10 2
3 Total 7 88
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
Grupo 2 78 39 19.5 0.004353 **
Residuals 5 10 2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova Table (Type II tests)
Response: y
Sum Sq Df F value Pr(>F)
Grupo 78 2 19.5 0.004353 **
Residuals 10 5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova Table (Type III tests)
Response: y
Sum Sq Df F value Pr(>F)
(Intercept) 192 1 96.0 0.0001885 ***
Grupo 78 2 19.5 0.0043530 **
Residuals 10 5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\[\Diamond\]
\[ \begin{align} \mathbf{X}_{\ell j} &= \boldsymbol{\mu} + \boldsymbol{\tau}_\ell + \mathbf{e}_{\ell j}\\ j &= 1, 2, \ldots, n_\ell \\ \ell &= 1, 2, \ldots, g \end{align} \tag{6-38} \]
em que os \(\mathbf{e}_{\ell j}\) são variáveis independentes \(\mathcal{N}_p(\mathbf{0},\mathbf{\Sigma})\). Aqui, o vetor de parâmetros \(\boldsymbol{\mu}\) é a média geral (nível), e \(\boldsymbol{\tau}_\ell\) representa o efeito de tratamento \(\ell\) com
\[ \sum_{\ell = 1}^{g} n_\ell \boldsymbol{\tau}_\ell = \mathbf{0} \]
De acordo com o modelo em (6-38), cada componente do vetor de observação satisfaz o modelo univariado (6-33). Os erros para os componentes de \(\mathbf{X}_{\ell j}\) estão correlacionados, mas a matriz de covariância \(\mathbf{\Sigma}\) é a mesma para todas as populações.
Um vetor de observações pode ser decomposto conforme sugerido pelo modelo. Assim,
(observação) = (média geral) + (efeito tratamento) + (resíduo)
\[ \begin{align} \mathbf{X}_{\ell j} &= \boldsymbol{\mu} + \boldsymbol{\tau}_\ell + \mathbf{e}_{\ell j}\\ &= \hat{\boldsymbol{\mu}} + \hat{\boldsymbol{\tau}}_\ell + \hat{\mathbf{e}}_{\ell j}\\ &= \hat{\boldsymbol{\mu}} + (\hat{\boldsymbol{\mu}}_\ell - \hat{\boldsymbol{\mu}}) + \hat{\mathbf{e}}_{\ell j}\\ \mathbf{X}_{\ell j} &= \bar{\mathbf{X}} + (\bar{\mathbf{X}}_\ell - \bar{\mathbf{X}}) + (\bar{\mathbf{X}}_{\ell j} - \bar{\mathbf{X}}_\ell) \end{align} \tag{6-39} \]
A decomposição em (6-39) leva ao análogo multivariado da divisão univariada da soma dos quadrados em (6-35). Primeiro, notamos que o produto
\[ (\mathbf{x}_{\ell j} - \bar{\mathbf{x}})(\mathbf{x}_{\ell j} - \bar{\mathbf{x}})^{\prime} \]
pode ser expressa como
\[ \begin{align} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}})(\mathbf{x}_{\ell j} - \bar{\mathbf{x}})^{\prime}&=((\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)+(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}}))((\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)+(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}}))^{\prime}\\ (\mathbf{x}_{\ell j} - \bar{\mathbf{x}})(\mathbf{x}_{\ell j} - \bar{\mathbf{x}})^{\prime}&=(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)^{\prime}+(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})^{\prime}\\ &\quad+(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)^{\prime}+(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})^{\prime} \end{align} \]
A soma sobre \(j\) das duas expressões intermediárias é a matriz zero, porque
\[ \sum_{j=1}^{n_\ell} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)=\mathbf{0} \]
Assim, somando o produto cruzado sobre \(j\) e \(\ell\) obtemos:
Total = Entre(Between) + Dentro(Within)
\[ \begin{align} \sum_{\ell=1}^{g} \sum_{j=1}^{n_\ell} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}})(\mathbf{x}_{\ell j} - \bar{\mathbf{x}})^{\prime} &= \sum_{\ell=1}^{g} n_\ell (\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})^{\prime}+\\ &\sum_{\ell=1}^{g} \sum_{j=1}^{n_\ell} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)^{\prime}\\ \mathbf{T}&=\mathbf{B}+\mathbf{W} \end{align} \tag{6-40} \]
A matriz de somas de quadrados e produtos cruzados dentro (Within) pode ser expressa como:
\[ \begin{align} \mathbf{W} &= \sum_{\ell=1}^{g} \sum_{j=1}^{n_\ell} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)^{\prime}\\ \mathbf{W} &= \sum_{\ell=1}^{g}{(n_\ell-1) \mathbf{S}_\ell} \end{align} \tag{6-41} \]
em que \(\mathbf{S}_\ell\) é a matriz de covariância amostral para a \(\ell\)-ésima condição independente. Esta matriz é uma generalização da matriz agrupada \((n_1 + n_2 - 2)\mathbf{S}_{\text{comb}}\) encontrada no caso de duas condições. Ela desempenha um papel dominante no teste para a presença de efeitos de tratamento.
Análogo ao resultado univariado, a hipótese de não haver efeitos de tratamento
\[ H_0:\;\boldsymbol{\tau}_1=\boldsymbol{\tau}_2=\cdots=\boldsymbol{\tau}_g=\mathbf{0} \]
é testada considerando os tamanhos relativos das somas de quadrados de tratamento e residuais e produtos cruzados. Equivalentemente, podemos considerar os tamanhos relativos das somas de quadrados residuais e totais (corrigidos) e produtos cruzados. Formalmente, resumimos os cálculos que levam à estatística de teste em uma tabela MANOVA.
Tabela MANOVA para Comparação de Vetores de Média Populacionais
Fonte de variação | Matriz de soma de quadrados e produtos cruzados (SSP) | Graus de liberdade (df) |
---|---|---|
Tratamento | \(\mathbf{B} = \sum_{\ell=1}^{g} n_\ell (\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})^{\prime}\) | \(g - 1\) |
Resíduo | \(\mathbf{W} = \sum_{\ell=1}^{g} \sum_{j=1}^{n_\ell} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)(\mathbf{x}_{\ell j} - \bar{\mathbf{x}}_\ell)^{\prime}\) | \(\sum_{\ell=1}^{g}n_\ell - g\) |
Total | \(\mathbf{T}=\sum_{\ell=1}^{g} \sum_{j=1}^{n_\ell} (\mathbf{x}_{\ell j} - \bar{\mathbf{x}})(\mathbf{x}_{\ell j} - \bar{\mathbf{x}})^{\prime}\) | \(\sum_{\ell=1}^{g}n_\ell - 1\) |
Esta tabela tem exatamente a mesma forma, componente por componente, que a tabela ANOVA unifatorial independente de Fisher, exceto que os quadrados de escalares são substituídos por seus contrapartes vetoriais. Por exemplo, \((\bar{\mathbf{x}}_\ell - \bar{x})^2\) torna-se \((\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})(\bar{\mathbf{x}}_{\ell} - \bar{\mathbf{x}})^{\prime}\). Os graus de liberdade correspondem à geometria univariada e também a alguma teoria de distribuição multivariada envolvendo densidades de Wishart. (Veja [1].)
Um teste de \(H_0:\;\boldsymbol{\tau}_1=\boldsymbol{\tau}_2=\cdots=\boldsymbol{\tau}_g=\mathbf{0}\) envolve variâncias generalizadas. Rejeitamos \(H_0\) se a razão das variâncias generalizadas, \(\Lambda^{*}\), for muito pequena:
\[ \Lambda^{*}=\dfrac{|\mathbf{W}|}{|\mathbf{W} + \mathbf{B}|} \tag{6-42} \]
O valor \(\Lambda^{*}\), proposto originalmente por Wilks (veja [25]), corresponde à forma equivalente (6-37) do teste F de \(H_0\): sem efeitos de tratamento no caso univariado. Lambda de Wilks tem a virtude de ser conveniente e relacionado ao critério da razão de verossimilhança. A distribuição exata de \(\Lambda^{*}\) pode ser derivada para os casos especiais listados na Tabela 6.3. Para outros casos e tamanhos de amostras grandes, uma modificação de \(\Lambda^{*}\) devido a Bartlett (veja [4]) pode ser usada para testar \(H_0\).
Lambda de Wilks também pode ser expresso como uma função dos autovalores \(\lambda_1, \lambda_2, \dots, \lambda_r\) de \(\mathbf{W}^{-1}\mathbf{B}\) como:
\[ \Lambda^* = \prod_{i=1}^{r} \dfrac{1}{1 + \lambda_i} \]
em que \(r = \min(p, g - 1)\), sendo \(p\) o posto de \(\mathbf{B}\). Outras estatísticas para verificar a igualdade de várias médias multivariadas, como a estatística de Pillai, a estatística de Lawley-Hotelling e a estatística da maior raiz de Roy também podem ser escritas como funções particulares dos autovalores de \(\mathbf{W}^{-1}\mathbf{B}\). Para amostras grandes, todas essas estatísticas são, essencialmente, equivalentes. (Veja a discussão adicional na página 336).
Tabela 6.3 Distribuição de Lambda de Wilks, \(\Lambda^*\)
No. de variáveis | No. de grupos | Distribuição amostral c/ dados multinormais |
---|---|---|
\(p = 1\) | \(g \ge 2\) | \(\dfrac{\sum{n_\ell}-g}{g-1}\dfrac{1-\Lambda^*}{\Lambda^*}\sim F_{g-1,\sum{n_\ell}-g}\) |
\(p = 2\) | \(g \ge 2\) | \(\dfrac{\sum{n_\ell}-g-1}{g-1}\dfrac{1-\sqrt{\Lambda^*}}{\sqrt{\Lambda^*}}\sim F_{2(g-1),2\left(\sum{n_\ell}-g-1\right)}\) |
\(p \ge 1\) | \(g = 2\) | \(\dfrac{\sum{n_\ell}-p-1}{p}\dfrac{1-\Lambda^*}{\Lambda^*}\sim F_{p,\sum{n_\ell}-p-1}\) |
\(p \ge 1\) | \(g = 3\) | \(\dfrac{\sum{n_\ell}-p-2}{p}\dfrac{1-\sqrt{\Lambda^*}}{\sqrt{\Lambda^*}}\sim F_{2p,2\left(\sum{n_\ell}-p-2\right)}\) |
Bartlett (veja [4]) mostrou que, se \(H_0\) é verdadeira e \(\sum{n_\ell} - p - g \gg 30\),
\[ -\left(n - 1 - \dfrac{p + g}{2}\right) \ln \Lambda^* \underset{a}{\sim} \chi^2_{p(g-1)} \]
Consequentemente, para \(\sum{n_\ell} - p - g \ge 30\) grande, rejeitamos \(H_0\) no nível de significância \(\alpha\) se
\[ -\left(n - 1 - \dfrac{p + g}{2}\right) \ln \Lambda^* > \chi^2_{p(g-1)}(1-\alpha) \]
Suponha que uma variável adicional seja observada juntamente com a variável introduzida no Exemplo 6.7: Os tamanhos de amostra são \(n_1 = 3\), \(n_2 = 2\) e \(n_3 = 3\). Organizando os pares de observações \(\mathbf{x}_{\ell j}\) em linhas, obtemos:
\[ \boldsymbol{x}= \begin{bmatrix} \begin{bmatrix}9\\3\end{bmatrix} &\begin{bmatrix}6\\2\end{bmatrix} &\begin{bmatrix}9\\7\end{bmatrix}\\ \begin{bmatrix}0\\4\end{bmatrix} &\begin{bmatrix}2\\0\end{bmatrix} \\ \begin{bmatrix}3\\8\end{bmatrix} &\begin{bmatrix}1\\9\end{bmatrix} &\begin{bmatrix}2\\7\end{bmatrix} \end{bmatrix} \]
Então:
\[ \bar{\mathbf{x}}_1=\begin{bmatrix}8\\4\end{bmatrix}\quad \bar{\mathbf{x}}_2=\begin{bmatrix}1\\2\end{bmatrix}\quad \bar{\mathbf{x}}_3=\begin{bmatrix}2\\8\end{bmatrix}\\ \bar{\mathbf{x}}=\begin{bmatrix}4\\5\end{bmatrix} \]
[1] "pt_BR.UTF-8"
[1] "LC_COLLATE=pt_BR.UTF-8;LC_CTYPE=pt_BR.UTF-8;LC_MONETARY=pt_BR.UTF-8;LC_NUMERIC=C;LC_TIME=pt_BR.UTF-8"
n1 <- 3
n2 <- 2
n3 <- 3
n <- n1 + n2 + n3
Grupo <- factor(c(rep("1", n1), rep("2", n2), rep("3", n3)))
y <- matrix(c(9, 3, 6, 2, 9, 7,
0, 4, 2, 0,
3, 8, 1, 9, 2, 7),
ncol=2, byrow=TRUE)
x <- data.frame(Grupo, y)
print(x)
Grupo X1 X2
1 1 9 3
2 1 6 2
3 1 9 7
4 2 0 4
5 2 2 0
6 3 3 8
7 3 1 9
8 3 2 7
$`1`
X1 X2
8 4
$`2`
X1 X2
1 2
$`3`
X1 X2
2 8
X1 X2
4 5
B <- n1*(x_bar_i[[1]] - x_bar)%*%t(x_bar_i[[1]] - x_bar) +
n2*(x_bar_i[[2]] - x_bar)%*%t(x_bar_i[[2]] - x_bar) +
n3*(x_bar_i[[3]] - x_bar)%*%t(x_bar_i[[3]] - x_bar)
B
X1 X2
[1,] 78 -12
[2,] -12 48
W <- (y[1,] - x_bar_i[[1]])%*%t(y[1,] - x_bar_i[[1]]) +
(y[2,] - x_bar_i[[1]])%*%t(y[2,] - x_bar_i[[1]]) +
(y[3,] - x_bar_i[[1]])%*%t(y[3,] - x_bar_i[[1]]) +
(y[4,] - x_bar_i[[2]])%*%t(y[4,] - x_bar_i[[2]]) +
(y[5,] - x_bar_i[[2]])%*%t(y[5,] - x_bar_i[[2]]) +
(y[6,] - x_bar_i[[3]])%*%t(y[6,] - x_bar_i[[3]]) +
(y[7,] - x_bar_i[[3]])%*%t(y[7,] - x_bar_i[[3]]) +
(y[8,] - x_bar_i[[3]])%*%t(y[8,] - x_bar_i[[3]])
W
X1 X2
[1,] 10 1
[2,] 1 24
[1] 0.03845535
g <- length(unique(Grupo))
# O livro usou linha 2 da tabela 6.3, mas poderia usar linha 4!
F <- ((1-sqrt(Lambda))/sqrt(Lambda))*(n-g-1)/(g-1)
F
[1] 8.19886
[1] 4
[1] 8
[1] 7.006077
[1] TRUE
pv <- formatC(1 - pf(F, df1, df2),
format="e", digits=2)
cat("\nF(0.99, ", df1, ", ", df2, ") = ",
round(Fcrit,2), ", F = ", round(F,2),
", p = ", pv, "\n", sep="")
F(0.99, 4, 8) = 7.01, F = 8.2, p = 6.23e-03
One-way MANOVA (Bartlett Chi2)
data: x
Wilks' Lambda = 0.038455, Chi2-Value = 14.662, DF = 4.000, p-value =
0.005456
sample estimates:
X1 X2
1 8 4
2 1 2
3 2 8
# rrcov::Wilks.test(x=iris[,1:4],grouping=iris[,5],method="mcd",nrep=1e3)
rrcov::Wilks.test(Grupo~., data=x)
One-way MANOVA (Bartlett Chi2)
data: x
Wilks' Lambda = 0.038455, Chi2-Value = 14.662, DF = 4.000, p-value =
0.005456
sample estimates:
X1 X2
1 8 4
2 1 2
3 2 8
Response X1 :
Call:
lm(formula = X1 ~ Grupo, data = x, na.action = "na.omit")
Residuals:
1 2 3 4 5 6 7 8
1 -2 1 -1 1 1 -1 0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.0000 0.8165 9.798 0.000189 ***
Grupo2 -7.0000 1.2910 -5.422 0.002890 **
Grupo3 -6.0000 1.1547 -5.196 0.003478 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.414 on 5 degrees of freedom
Multiple R-squared: 0.8864, Adjusted R-squared: 0.8409
F-statistic: 19.5 on 2 and 5 DF, p-value: 0.004353
Response X2 :
Call:
lm(formula = X2 ~ Grupo, data = x, na.action = "na.omit")
Residuals:
1 2 3 4 5 6 7
-1.000e+00 -2.000e+00 3.000e+00 2.000e+00 -2.000e+00 5.551e-17 1.000e+00
8
-1.000e+00
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.000 1.265 3.162 0.0250 *
Grupo2 -2.000 2.000 -1.000 0.3632
Grupo3 4.000 1.789 2.236 0.0756 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.191 on 5 degrees of freedom
Multiple R-squared: 0.6667, Adjusted R-squared: 0.5333
F-statistic: 5 on 2 and 5 DF, p-value: 0.06415
Type II MANOVA Tests:
Sum of squares and products for error:
X1 X2
X1 10 1
X2 1 24
------------------------------------------
Term: Grupo
Sum of squares and products for the hypothesis:
X1 X2
X1 78 -12
X2 -12 48
Multivariate Tests: Grupo
Df test stat approx F num Df den Df Pr(>F)
Pillai 2 1.540788 8.388227 4 10 0.0030962 **
Wilks 2 0.038455 8.198860 4 8 0.0062341 **
Hotelling-Lawley 2 9.941423 7.456067 4 6 0.0164318 *
Roy 2 8.076385 20.190963 2 5 0.0040292 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Response X1 :
Df Sum Sq Mean Sq F value Pr(>F)
Grupo 2 78 39 19.5 0.004353 **
Residuals 5 10 2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Response X2 :
Df Sum Sq Mean Sq F value Pr(>F)
Grupo 2 48 24.0 5 0.06415 .
Residuals 5 24 4.8
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type II MANOVA Tests: Wilks test statistic
Df test stat approx F num Df den Df Pr(>F)
Grupo 2 0.0385 8.2 4 8 0.0062 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For one-way between subjects designs, partial eta squared is equivalent
to eta squared. Returning eta squared.
For one-way between subjects designs, partial eta squared is equivalent
to eta squared. Returning eta squared.
# Effect Size for ANOVA (Type I)
Response | Parameter | Eta2 | 95% CI
---------------------------------------------
X1 | Grupo | 0.886 | [0.527, 1.000]
X2 | Grupo | 0.667 | [0.000, 1.000]
Type II Sums of Squares
df X1 X2
Grupo 2 78 48
residuals 5 10 24
F-tests
X1 X2
Grupo 19.5 5
p-values
X1 X2
Grupo 0.004353 0.064150
p-values adjusted (by term) for simultaneous inference by holm method
X1 X2
Grupo 0.0087061 0.0641500
\[\Diamond\]
O Departamento de Saúde e Serviços Sociais de Wisconsin reembolsa os asilos de idosos no estado pelos serviços prestados. O departamento desenvolve um conjunto de fórmulas para taxas para cada instituição, com base em fatores como nível de atendimento, salário médio e salário médio no estado.
Os asilos de idosos podem ser classificados com base na propriedade (particular, organização sem fins lucrativos e governo) e certificação (instalação de enfermagem especializada, instalação de cuidados intermediários ou uma combinação dos dois).
Um dos objetivos de um estudo recente foi investigar os efeitos da propriedade ou certificação (ou ambos) nos custos. Quatro custos, calculados com base no paciente por dia e medidos em horas por dia do paciente, foram selecionados para análise: \(X_1\) = custo da mão-de-obra de enfermagem, \(X_2\) = custo da mão-de-obra dietética, \(X_3\) = custo da mão-de-obra de operação e manutenção da planta, e \(X_4\) = custo da mão-de-obra de limpeza e lavanderia. Um total de \(n = 516\) observações em cada uma das \(p = 4\) variáveis de custo foi inicialmente separado de acordo com a propriedade.
\[ n_1 = 271, \quad n_2 = 138, \quad n_3 = 107 \]
\[ \bar{\mathbf{x}}_1 = \begin{bmatrix} 2.066 \\ 0.480 \\ 0.082 \\ 0.360 \\ \end{bmatrix} \quad \bar{\mathbf{x}}_2 = \begin{bmatrix} 2.167 \\ 0.596 \\ 0.124 \\ 0.418 \\ \end{bmatrix} \quad \bar{\mathbf{x}}_3 = \begin{bmatrix} 2.273 \\ 0.521 \\ 0.125 \\ 0.383 \\ \end{bmatrix} \]
Média geral ponderada:
\[ \bar{\mathbf{x}} = \dfrac{n_1 \cdot \bar{\mathbf{x}}_1 + n_2 \cdot \bar{\mathbf{x}}_2 + n_3 \cdot \bar{\mathbf{x}}_3}{n} \]
O cálculo de \(\mathbf{B}\) necessita apenas dos dados brutos, mas dos tamanhos e médias dos grupos.
A matriz \(B\) é a matriz de somas de quadrados e produtos cruzados entre os grupos (tratamentos) e é calculada como:
\[ \begin{align} \mathbf{B} &= n_1 (\bar{\mathbf{x}}_1 - \bar{\mathbf{x}})(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}})^T + n_2 (\bar{\mathbf{x}}_2 - \bar{\mathbf{x}})(\bar{\mathbf{x}}_2 - \bar{\mathbf{x}})^T + n_3 (\bar{\mathbf{x}}_3 - \bar{\mathbf{x}})(\bar{\mathbf{x}}_3 - \bar{\mathbf{x}})^T\\ \mathbf{B} &= \begin{bmatrix} 3.4688 & 1.0986 & 0.8107 & 0.5860 \\ 1.0986 & 1.2307 & 0.4500 & 0.6157 \\ 0.8107 & 0.4500 & 0.2318 & 0.2311 \\ 0.5860 & 0.6157 & 0.2311 & 0.3086 \\ \end{bmatrix} \end{align} \]
Como as \(\mathbf{s}_\ell\) parecem ser razoavelmente compatíveis, elas foram agrupadas [veja (6-41)] para obter:
\[ \mathbf{W} = (n_1 - 1)\mathbf{s}_1 + (n_2 - 1)\mathbf{s}_2 + (n_3 - 1)\mathbf{s}_3 \]
\[ \begin{align} \mathbf{s}_1 &= \begin{bmatrix} 0.291 & -0.001 & 0.002 & 0.010 \\ -0.001 & 0.011 & 0.000 & 0.003 \\ 0.002 & 0.000 & 0.001 & 0.000 \\ 0.010 & 0.003 & 0.000 & 0.010 \\ \end{bmatrix} \\ \mathbf{s}_2 &= \begin{bmatrix} 0.561 & 0.011 & 0.001 & 0.037 \\ 0.011 & 0.025 & 0.004 & 0.007 \\ 0.001 & 0.004 & 0.005 & 0.002 \\ 0.037 & 0.007 & 0.002 & 0.019 \\ \end{bmatrix} \\ \mathbf{s}_3 &= \begin{bmatrix} 0.261 & 0.030 & 0.003 & 0.018 \\ 0.030 & 0.017 & -0.000 & 0.006 \\ 0.003 & -0.000 & 0.004 & 0.001 \\ 0.018 & 0.006 & 0.001 & 0.013 \\ \end{bmatrix} \end{align} \]
\[ \mathbf{W} = \begin{bmatrix} 183.093 & 4.417 & 0.995 & 9.677 \\ 4.417 & 8.197 & 0.548 & 2.405 \\ 0.995 & 0.548 & 1.379 & 0.380 \\ 9.677 & 2.405 & 0.380 & 6.681 \\ \end{bmatrix} \]
Para testar \(H_0: \boldsymbol{\tau}_1 = \boldsymbol{\tau}_2 = \boldsymbol{\tau}_3\) (sem efeitos de propriedade ou, equivalentemente, sem diferença nos custos médios entre os três tipos de proprietários - privado, sem fins lucrativos e governo), podemos usar o resultado na Tabela 6.3 para \(g = 3\). Então, \(\Lambda^{*}=0.76\).
No entanto, um teste de teoria normal de \(H_0: \mathbf{\Sigma}_1 = \mathbf{\Sigma}_2 = \mathbf{\Sigma}_3\) rejeitaria \(H_0\) em qualquer nível de significância razoável devido aos grandes tamanhos de amostra (veja o Exemplo 6.12).
É informativo comparar os resultados com base no teste “exato” com aqueles obtidos usando o procedimento de grande amostra resumido em (6-43) e (6-44). Para o exemplo atual, \(n - p - g = 516 - 4 - 3 = 509 \ge 30\), e \(H_0\) pode ser testado no nível \(\alpha = 0.01\) comparando
\[ X^2=-\left(n - 1 - \dfrac{p + g}{2}\right)\ln\left(\Lambda^{*}\right) = -511.5 \ln(0.7714) = 138.52 \]
com \(\chi^2_{p(g-1)}(1-\alpha) = \chi^2_8(0.99) = 20.09\). Como \(138.52 > 20.09\), rejeitamos \(H_0\) no nível de 1%.
n1 <- 271
n2 <- 138
n3 <- 107
n <- n1 + n2 + n3
x_bar_1 <- c(2.066, 0.480, 0.082, 0.360)
x_bar_2 <- c(2.167, 0.596, 0.124, 0.418)
x_bar_3 <- c(2.273, 0.521, 0.125, 0.383)
x_bar <- (n1*x_bar_1+n2*x_bar_2+n3*x_bar_3)/n
x_bar
[1] 2.1359360 0.5195252 0.1021492 0.3802810
s1 <- matrix(c(0.291, -0.001, 0.002, 0.010,
-0.001, 0.011, 0.000, 0.003,
0.002, 0.000, 0.001, 0.000,
0.010, 0.003, 0.000, 0.010), nrow=4, byrow=TRUE)
s2 <- matrix(c(0.561, 0.011, 0.001, 0.037,
0.011, 0.025, 0.004, 0.007,
0.001, 0.004, 0.005, 0.002,
0.037, 0.007, 0.002, 0.019), nrow=4, byrow=TRUE)
s3 <- matrix(c(0.261, 0.030, 0.003, 0.018,
0.030, 0.017, -0.000, 0.006,
0.003, -0.000, 0.004, 0.001,
0.018, 0.006, 0.001, 0.013), nrow=4, byrow=TRUE)
B <- n1*(x_bar_1 - x_bar)%*%t(x_bar_1 - x_bar) +
n2*(x_bar_2 - x_bar)%*%t(x_bar_2 - x_bar) +
n3*(x_bar_3 - x_bar)%*%t(x_bar_3 - x_bar)
B
[,1] [,2] [,3] [,4]
[1,] 3.4687989 1.0985713 0.8106779 0.5859503
[2,] 1.0985713 1.2306787 0.4500336 0.6157338
[3,] 0.8106779 0.4500336 0.2317835 0.2311294
[4,] 0.5859503 0.6157338 0.2311294 0.3085943
[,1] [,2] [,3] [,4]
[1,] 183.093 4.417 0.995 9.677
[2,] 4.417 8.197 0.548 2.405
[3,] 0.995 0.548 1.379 0.380
[4,] 9.677 2.405 0.380 6.681
[1] 0.7627582
[1] 18.48786
[1] 8
[1] 1020
[1] 2.528682
[1] TRUE
[1] "0.00e+00"
cat("\nF(0.99, ", df1, ", ", df2, ") = ",
round(Fcrit,2), ", F = ", round(F,2),
", p = ", p_value, "\n", sep="")
F(0.99, 8, 1020) = 2.53, F = 18.49, p = 0.00e+00
[1] 2.511279
cat("\nTCL: F(0.99, ", df1, ", ", df2, ") = ",
round(X2crit,2), ", F = ", round(F,2),
", p = ", p_value, "\n", sep="")
TCL: F(0.99, 8, 1020) = 2.51, F = 18.49, p = 0.00e+00
[1] 138.5215
[1] 20.09024
[1] TRUE
[1] "0.00e+00"
cat("\nTCL (lambda de Wilks): X^2(0.99, ", df, ") = ",
round(X2crit,2), ", X2 = ", round(X2,2),
", p = ", p_value, "\n", sep="")
TCL (lambda de Wilks): X^2(0.99, 8) = 20.09, X2 = 138.52, p = 0.00e+00
\[\Diamond\]
Quando a hipótese de efeitos de tratamento iguais é rejeitada, os efeitos que levaram à rejeição da hipótese são de interesse. Para comparações aos pares, a abordagem de Bonferroni (veja a Seção 5.4) pode ser usada para construir intervalos de confiança simultâneos para os componentes das diferenças \(\boldsymbol{\tau}_k - \boldsymbol{\tau}_\ell\) (ou \(\boldsymbol{\mu}_k - \boldsymbol{\mu}_\ell\)). Esses intervalos são mais curtos do que os obtidos para todos os contrastes e requerem valores críticos apenas para a estatística t univariada.
Seja \(\tau_{ki}\) o i-ésimo componente de \(\boldsymbol{\tau}_k\). Uma vez que \(\boldsymbol{\tau}_k\) é estimado por \(\hat{\boldsymbol{\tau}}_k = \bar{\mathbf{x}}_k - \bar{\mathbf{x}}\),
\[ \hat{\tau}_{ki} = \bar{x}_{ki} - \bar{x}_{i} \]
em que \(\hat{\tau}_{ki} - \hat{\tau}_{\ell i} = \bar{x}_{ki}-\bar{x}_{\ell i}\) é a diferença entre duas médias amostrais independentes. O intervalo de confiança t é válido com um \(\alpha\) apropriadamente modificado. Observe que:
\[ \mathbb{V}\left(\hat{\tau}_{ki} - \hat{\tau}_{\ell i} \right) = \mathbb{V}\left(\bar{X}_{ki} - \bar{X}_{\ell i} \right) = \left(\dfrac{1}{n_k}+\dfrac{1}{n_\ell}\right)\sigma_{ii} \]
Como sugerido por (6-41), \(\mathbb{V}\left(\bar{X}_{ki} - \bar{X}_{\ell i} \right)\) é estimado dividindo-se o elemento correspondente de \(\mathbf{W}\) por seus graus de liberdade. Ou seja,
\[ \widehat{\mathbb{V}}\left(\bar{X}_{ki} - \bar{X}_{\ell i} \right)=\left(\dfrac{1}{n_k}+\dfrac{1}{n_\ell}\right)\dfrac{w_{ii}}{n-g} \]
em que \(w_{ii}\) é o i-ésimo elemento diagonal de \(\mathbf{W}\) e \(n = \sum_{\ell = 1}^{g}{n_\ell}\).
Ainda resta distribuir a probabilidade do erro do tipo I global entre as diversos intervalos de confiança. A relação (5-28) ainda se aplica.
Existem \(p\) variáveis e \(\frac{g(g - 1)}{2}\) diferenças aos pares, então cada intervalo t empregará o valor crítico \(t_{n-g}(\frac{\alpha}{2m})\), em que
\[ m = p\dfrac{g(g - 1)}{2}=pC_{g,2} \tag{6-46} \]
é o número de intervalos de confiança simultâneos.
Resultado 6.5. Seja \(n = \sum_{\ell = 1}^{g}{n_\ell}\). Para o modelo em (6-38), com confiança de pelo menos \((1-\alpha)\),
\[ \text{IC}^{1-\alpha}\left(\tau_{kl} - \tau_{\ell i}\right)=\left[ \bar{x}_{ki} - \bar{x}_{\ell i} \pm t_{n-g}\left(\dfrac{1-\alpha}{pg(g-1)}\right) \sqrt{\left(\dfrac{1}{n_k} + \dfrac{1}{n_\ell}\right)\dfrac{w_{ii}}{n-g}}\right] \]
para todos os componentes \(i = 1, 2, \dots, p\) e todas as diferenças \(\ell < k = 1, 2, \dots, g\). Aqui \(w_{ii}\) é o i-ésimo elemento diagonal de \(\mathbf{W}\).
Ilustraremos a construção de intervalo de confiança simultâneos para as diferenças aos pares nas médias de tratamento usando os dados de asilos de enfermagem introduzidos no Exemplo 6.10.
Vimos no Exemplo 6.10 que os custos médios para asilos de enfermagem diferem, dependendo do tipo de propriedade. Podemos usar o Resultado 6.5 para estimar as magnitudes das diferenças. Uma comparação da variável \(X_3\), custos de operação e mão de obra de manutenção, entre asilos de enfermagem de propriedade privada e asilos de enfermagem de propriedade do governo pode ser feita estimando \(\tau_{13} - \tau_{33}\). Usando (6-39) e as informações no Exemplo 6.10, temos:
\[ \hat{\boldsymbol{\tau}}_1 = \bar{\mathbf{x}}_1 - \bar{\mathbf{x}} = \begin{bmatrix} -0.070 \\ -0.039 \\ -0.020 \\ -0.020 \end{bmatrix} \]
\[ \hat{\boldsymbol{\tau}}_3 = \bar{\mathbf{x}}_3 - \bar{\mathbf{x}} = \begin{bmatrix} 0.137 \\ 0.002 \\ 0.023 \\ 0.003 \end{bmatrix} \]
\[ \mathbf{W} = \begin{bmatrix} 183.093 & 4.417 & 0.995 & 9.677 \\ 4.417 & 8.197 & 0.548 & 2.405 \\ 0.995 & 0.548 & 1.379 & 0.380 \\ 9.677 & 2.405 & 0.380 & 6.681 \\ \end{bmatrix} \]
Consequentemente,
\[ \hat{\tau}_{13} -\hat{\tau}_{33} = -0.020 - 0.023 = -0.043 \]
e \(n = 271 + 138 + 107 = 516\), de forma que
\[ \sqrt{\left(\dfrac{1}{n_1} + \dfrac{1}{n_3}\right)\dfrac{w_{33}}{n - g}} = \sqrt{\left(\dfrac{1}{271} + \dfrac{1}{107}\right)\dfrac{1.484}{516 - 3}} = 0.00614 \]
Uma vez que \(p = 4\) e \(g= 3\), para intervalos de confiança simultâneos de 95%, exigimos que \(t_{513}(0.05/4(3)2) \approx 2.87\). (Consulte o Apêndice, Tabela 1.) O intervalo de confiança simultâneo de 95% é
\[ \begin{align} \text{IC}^{1-0.00208} \left(\tau_{13} - \tau_{33}\right)&= \left[\hat{\tau}_{13} - \hat{\tau}_{33} \pm t_{513}(1-0.00208)\sqrt{\left(\dfrac{1}{n_1} + \dfrac{1}{n_3}\right)\dfrac{n}{n - g}}\right]\\ \text{IC}^{1-0.00208} \left(\tau_{13} - \tau_{33}\right)&= [-0.061, -0.025] \end{align} \]
Concluímos que o custo médio de manutenção e mão de obra para asilos de enfermagem de propriedade do governo é maior em .025 a .061 hora por dia por paciente do que para asilos de enfermagem de propriedade privada. Com a mesma confiança de 95%, podemos dizer que
\[ \text{IC}^{1-0.00208}\left(\tau_{13} - \tau_{23}\right)=[-.058, -.026] \]
e
\[ \text{IC}^{1-0.00208}\left(\tau_{23} - \tau_{33}\right)=[-.021, .019) \]
Assim, existe uma diferença neste custo entre asilos de enfermagem privados e sem fins lucrativos, mas não se observa diferença entre asilos de enfermagem sem fins lucrativos e de propriedade do governo.
# Dados
x_bar_1 <- c(-0.070, -0.039, -0.020, -0.020)
x_bar_3 <- c(0.137, 0.002, 0.023, 0.003)
W <- matrix(c(183.093, 4.417, 0.995, 9.677,
4.417, 8.197, 0.548, 2.405,
0.995, 0.548, 1.379, 0.380,
9.677, 2.405, 0.380, 6.681), nrow=4, byrow=TRUE)
w33 <- W[3,3]
n1 <- 271
n2 <- 138
n3 <- 107
n <- n1 + n2 + n3
g <- 3
p <- length(x_bar_1)
# Cálculos
tau_13 <- x_bar_1[3]
tau_33 <- x_bar_3[3]
tau_13_minus_tau_33 <- tau_13 - tau_33
error_term <- sqrt((1/n1 + 1/n3) * w33 / (n - g))
m <- p*g*(g-1)/2
m
[1] 12
[1] 2.878174
confidence_interval <- c(tau_13_minus_tau_33 - t_value * error_term,
tau_13_minus_tau_33 + t_value * error_term)
# Resultados
cat("IC", round(1-alfa/(2*m),6)*100, "%(tau_13 - tau_33) = ",
confidence_interval)
IC 99.7917 %(tau_13 - tau_33) = -0.06003765 -0.02596235
\[\Diamond\]
Uma das suposições feitas ao comparar dois ou mais vetores médios multivariados é que as matrizes de covariância das potenciais diferentes populações são as mesmas (homocedasticidade multivariada). (Esta suposição aparecerá novamente no Capítulo 11 quando discutirmos discriminação e classificação.) Antes de agrupar a variação entre amostras para formar uma matriz de covariância agrupada ao comparar vetores de média, pode ser útil testar a igualdade das matrizes de covariância da população. Um teste comumente empregado para matrizes de covariância iguais é o teste M de Box ([8], [9]).
Com \(g\) populações, a hipótese nula é
\[ H_0: \mathbf{\Sigma}_1 = \mathbf{\Sigma}_2 = \cdots = \mathbf{\Sigma}_g = \mathbf{\Sigma} \tag{6-47} \]
em que \(\mathbf{\Sigma}_\ell\) é a matriz de covariância para a \(\ell\)-ésima população, \(\ell = 1, 2, \cdots, g\), e \(\mathbf{\Sigma}\) é a matriz de covariância comum presumida. A hipótese alternativa é que pelo menos duas das matrizes de covariância não são iguais.
Assumindo populações multinormais, uma estatística de razão de verossimilhança para testar (6-47) é dada por (veja [1])
\[ \Lambda=\prod_{\ell = 1}^{g}{\left(\dfrac{|\mathbf{S}_\ell|}{|\mathbf{S}_\text{comb}|}\right)^{\frac{n_\ell -1}{2}}} \tag{6-48} \]
Aqui \(n_\ell\) é o tamanho da amostra para o \(\ell\)-ésimo grupo, \(\mathbf{S}_\ell\) é a matriz de covariância amostral do \(\ell\)-ésimo grupo e \(\mathbf{S}_{\text{comb}}\) é a matriz de covariância da amostra combinada dada por
\[ \mathbf{S}_{\text{comb}} = \dfrac{1}{\sum_{\ell =1}^{g}{(n_\ell -1)}} \sum_{\ell =1}^{g}{(n_\ell -1)\mathbf{S}_\ell} \tag{6-49} \]
O teste M de Box é baseado em sua aproximação \(\chi^2\) para a distribuição amostral de \(-2 \ln \Lambda\) (veja Resultado 5.2). Definindo \(-2 \ln \Lambda = M\) (estatística M de Box) temos:
\[ M = \sum_{\ell =1}^{g} {(n_\ell - 1)} \ln |\mathbf{S}_{\text{comb}}| - \sum_{\ell =1}^{g} {(n_\ell - 1) \ln |\mathbf{S}_\ell|} \tag{6-50} \]
Se a hipótese nula for verdadeira, não se espera que as matrizes de covariância amostral individuais difiram muito e, consequentemente, não difiram muito da matriz de covariância combinada. Nesse caso, a razão dos determinantes em (6-48) estará próxima de 1, \(\Lambda\) estará próximo de 1 e a estatística M de Box será pequena. Se a hipótese nula for falsa, as matrizes de covariância amostral podem diferir mais e as diferenças em seus determinantes serão mais pronunciadas. Neste caso, \(\Lambda\) será pequeno e M será relativamente grande. Para ilustrar, observe que o determinante da matriz de covariância combinada, \(|\mathbf{S}_{\text{comb}}|\), estará em algum lugar próximo ao “meio” dos determinantes \(|\mathbf{S}_\ell|\) das matrizes de covariância do grupo individual. À medida que essas quantidades se tornam mais distintas, o produto das razões em (6-44) se aproxima de 0. De fato, à medida que os \(|\mathbf{S}_\ell|\) aumentam em dispersão, \(|\mathbf{S}_{(1)}|/|\mathbf{S}_{\text{comb}}|\) reduz o produto proporcionalmente mais do que \(|\mathbf{S}_{(g)}|/|\mathbf{S}_{\text{comb}}|\) o aumenta, em que \(|\mathbf{S}_{(1)}|\) e \(|\mathbf{S}_{(g)}|\) são os valores determinantes mínimo e máximo, respectivamente.
Seja
\[ u=\left(\sum_{\ell = 1}^{g}{\dfrac{1}{n_\ell - 1}}-\dfrac{1}{\sum_{\ell = 1}^{g}{n_\ell - 1}}\right)\dfrac{2p^2+3p-1}{6(p+1)(g-1)} \tag{6-51} \]
em que \(p\) é o número de variáveis e \(g\) é o número de grupos.
Então
\[ C = (1-u)M \underset{a}{\sim}\chi^2_{p(p + 1)(g - 1)/2} \tag{6-52} \]
No nível de significância \(\alpha\), rejeitar \(H_0\) se \(C > \chi^2_{p(p + 1)(g - 1)/2}(1-\alpha)\).
A aproximação de M de Box funciona bem se cada \(n_\ell - p - g\) exceder 20 e se \(p \le 5\) e \(g \le 5\). Em situações em que essas condições não se mantêm, Box ([7],[8]) forneceu uma aproximação \(F\) mais precisa para a distribuição amostral de M de Box.
Introduzimos os dados de asilos de enfermagem de Wisconsin no Exemplo 6.10. Naquele exemplo, as matrizes de covariância amostral para \(p = 4\) variáveis de custo associadas com \(g = 3\) grupos de asilos de enfermagem são exibidas. Supondo dados multinormais, testamos a hipótese \(H_0: \mathbf{\Sigma}_1 = \mathbf{\Sigma}_2=\mathbf{\Sigma}_3=\mathbf{\Sigma}\).
Usando as informações do Exemplo 6.10, temos \(n_1 = 271\), \(n_2 = 138\), \(n_3 = 107\) e \(|\mathbf{s}_1| = 2.783 \times 10^{-8}\), \(|\mathbf{s}_2| = 89.539 \times 10^{-8}\), \(|\mathbf{s}_3| = 14.579 \times 10^{-8}\), e \(|\mathbf{s}_{\text{pooled}}| = 17.398 \times 10^{-8}\). Tomando os logaritmos naturais dos determinantes, obtemos:
\[ \ln |\mathbf{s}_1| = -17.397, \quad \ln |\mathbf{s}_2| = -13.926, \quad \ln |\mathbf{s}_3| = -15.741\\ \ln |\mathbf{s}_{\text{comb}}| = -15.564 \]
Calculamos
\[ \begin{align} u &= 0.0133\\ M &= (270 + 137 + 106)(-15.564) - (270(-17.397) + 137(-13.926) + 106(-15.741)) \\ M&= 244.15 \end{align} \]
e
\[ C = (1 - 0.0133) \times 244.15 = 240.76 \]
Consultando \(C\) em uma tabela com \(\nu = 4(4 + 1)(3 - 1)/2 = 20\) graus de liberdade, fica claro que \(H_0\) é rejeitado em qualquer nível razoável de significância. Concluímos que as matrizes de covariância das variáveis de custo associadas com as três populações de asilos de enfermagem não são as mesmas (heterocedasticidade multivariada).
# Definindo as matrizes e tamanhos das amostras
n1 <- 271
n2 <- 138
n3 <- 107
n <- n1+n2+n3
s1 <- matrix(c(0.291, -0.001, 0.002, 0.010,
-0.001, 0.011, 0.000, 0.003,
0.002, 0.000, 0.001, 0.000,
0.010, 0.003, 0.000, 0.010), nrow=4, byrow=TRUE)
s2 <- matrix(c(0.561, 0.011, 0.001, 0.037,
0.011, 0.025, 0.004, 0.007,
0.001, 0.004, 0.005, 0.002,
0.037, 0.007, 0.002, 0.019), nrow=4, byrow=TRUE)
s3 <- matrix(c(0.261, 0.030, 0.003, 0.018,
0.030, 0.017, -0.000, 0.006,
0.003, -0.000, 0.004, 0.001,
0.018, 0.006, 0.001, 0.013), nrow=4, byrow=TRUE)
# Definindo p e g
p <- 4
g <- 3
# Calculando a matriz de covariância combinada
s_pooled <- ((n1 - 1)*s1 + (n2 - 1)*s2 + (n3 - 1)*s3) / (n - g)
# Calculando M
u <- (1/(n1-1)+1/(n2-1)+1/(n3-1)-1/(n - g))*
(2*p^2 + 3*p + 1) / (6*(p + 1)*(g - 1))
u
[1] 0.01386571
[1] 244.146
[1] 240.7608
[1] TRUE
pv <- formatC(1 - pchisq(C, df),
format="e", digits=2)
cat("X^2(95%,", df, ") = ", round(X2crit,2),
", X^2 = ", round(C,2), ", p = ", pv, sep="")
X^2(95%,20) = 31.41, X^2 = 240.76, p = 0.00e+00
\[\Diamond\]
O teste M de Box é rotineiramente calculado em muitos pacotes de computação estatística que realizam MANOVA e outros procedimentos que requerem matrizes de covariância iguais. É conhecido que o teste M é sensível a algumas formas de não-normalidade. De forma mais ampla, na presença de não-normalidade, testes de teoria normal em covariâncias são influenciados pela curtose das populações originais (veja [16]). No entanto, com amostras suficientemente grandes, os testes MANOVA de médias ou efeitos de tratamento são bastante robustos à não-normalidade e também funciona o teorema central do limite que permite distribuições multivariadas distintas e eventualmente heterocedásticas nos \(g\) grupos. Assim, o teste M pode rejeitar \(H_0\) em alguns casos de não-normalidade, onde não é prejudicial para os testes MANOVA. Além disso, com tamanhos de amostra iguais (balanceamento perfeito), algumas diferenças nas matrizes de covariância têm pouco efeito nos testes MANOVA. Em resumo, podemos decidir continuar com os testes MANOVA usuais, mesmo que o teste \(M\) leve à rejeição de \(H_0\).
Após nossa abordagem à MANOVA Unifatorial Independente, revisaremos brevemente a análise para um modelo univariado independente de efeitos fixos de dois fatores e, em seguida, generalizaremos simplesmente para o caso multivariado por analogia.
Supomos que as medições são registradas em vários níveis de dois fatores. Em alguns casos, essas condições experimentais representam níveis de um único tratamento organizado em vários blocos. O delineamento experimental particular empregado não nos preocupará neste livro. (Consulte [10] e [17] para discussões sobre delineamento experimental.) No entanto, assumiremos que observações em diferentes combinações de condições experimentais são independentes entre si.
Vamos considerar que os dois conjuntos de condições experimentais sejam os níveis, por exemplo, do fator 1 e do fator 2, respectivamente. Suponha que haja \(g \ge 2\) níveis do fator 1 e \(b \ge 2\) níveis do fator 2, e que \(n \ge 2\) observações independentes possam ser observadas em cada uma das combinações de níveis (tratamentos) \(gb\).
O uso do termo fator para indicar uma condição experimental é conveniente. Os fatores discutidos aqui não devem ser confundidos com os fatores não observáveis considerados no Capítulo 9 no contexto da análise de fatores.
Denotando a r-ésima observação no nível \(\ell\) do fator 1 e nível \(k\) do fator 2 por \(X_{\ell kr}\), especificamos o modelo univariado de dois sentidos como:
\[ X_{\ell kr} = \mu + \tau_\ell + \beta_k + \gamma_{\ell k} + e_{\ell kr} \tag{6-54} \]
Em que:
Sendo que:
\[ \sum_{\ell=1}^{g} \tau_\ell = \sum_{k=1}^{b} \beta_k = \sum_{\ell=1}^{g} \gamma_{\ell k} = \sum_{k=1}^{b} \gamma_{\ell k} = 0 \]
Os \(e_{\ell kr}\) são variáveis aleatórias independentes \(\mathcal{N}(0, \sigma^2)\). Aqui, \(\mu\) representa um nível geral, \(\tau_\ell\) representa o efeito fixo do fator 1, \(\beta_k\) representa o efeito fixo do fator 2, e \(\gamma_{\ell k}\) é a interação entre o fator 1 e o fator 2. A resposta esperada no \(\ell\)-ésimo nível do fator 1 e do \(k\)-ésimo nível do fator 2 é, portanto:
\[ \mathbb{E}(X_{\ell kr}) = \mu + \tau_\ell + \beta_k + \gamma_{\ell k} \]
Em que:
para \(\ell = 1, 2, \ldots, g\) e \(k = 1, 2, \ldots, b\).
A presença da interação, \(\gamma_{\ell k}\), implica que os efeitos dos fatores não são aditivos e a interpretação dos resultados dos resultados torna-se mais complexa.
Figura 6.3 Perfis de médias (a) com interação e (b) sem interação.
As Figuras 6.3(a) e (b) mostram as respostas esperadas em função dos níveis dos fatores com e sem interação, respectivamente. A ausência de interação significa \(\gamma_{\ell k} = 0\) para todos \(\ell\) e \(k\).
De forma análoga a (6-55), cada observação pode ser decomposta como
\[ \begin{align} X_{\ell kr} &= \mu+\tau_{\ell}+\beta_k+\gamma_{\ell k}+e_{\ell kr}\\ &= \mu+\tau_{\ell}+\beta_k+\left(\left(\tau\beta\right)_{\ell k}-\left(\mu+\tau_{\ell}+\beta_k\right)\right)+e_{\ell kr}\\ &= \hat{\mu}+\hat{\tau}_{\ell}+\hat{\beta}_k+\left(\widehat{\left(\tau\beta\right)}_{\ell k}-\left( \hat{\mu}+\hat{\tau}_{\ell}+\hat{\beta}_k\right)\right)+\hat{e}_{\ell kr}\\ &= \bar{X} + (\bar{X}_{\ell} - \bar{X}) + (\bar{X}_{k} - \bar{X}) + (\bar{X}_{\ell k} - (\bar{X}+(\bar{X}_{\ell}- \bar{X}) + (\bar{X}_{k} - \bar{X}))) + (X_{\ell kr} - \bar{X}_{\ell k}) \\ X_{\ell kr} &= \bar{X} + (\bar{X}_{\ell} - \bar{X}) + (\bar{X}_{k} - \bar{X}) + (\bar{X}_{\ell k} - \bar{X}_{\ell} - \bar{X}_{k}+ \bar{X}) + (X_{\ell kr} - \bar{X}_{\ell k}) \end{align} \tag{6-56} \]
em que \(\bar{X}\) é a média geral, \(\bar{X}_{\ell}\) é a média para o nível \(\ell\) do fator 1, \(\bar{X}_k\) é a média para o nível \(k\) do fator 2, e \(X_{\ell k}\) é a média para o nível \(\ell\) do fator 1 e o nível \(k\) do fator 2. Ao quadrar e somar os desvios \((X_{\ell kr} - \bar{X})\), obtemos:
\[ \begin{align} \sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} (X_{\ell kr} - \bar{X})^2 &= b n \sum_{\ell=1}^{g} (\bar{X}_{\ell} - \bar{X})^2 + g n \sum_{k=1}^{b} (\bar{X}_k - \bar{X})^2 +\\ & \quad n \sum_{\ell=1}^{g} \sum_{k=1}^{b} (\bar{X}_{\ell k} - \bar{X}_{\ell} - \bar{X}_k + \bar{X})^2 +\\ & \quad\sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} (X_{\ell kr}- \bar{X}_{\ell k})^2 \end{align} \tag{6-57} \]
Em que:
\[ \begin{align} \bar{X}_{\ell k} &= \dfrac{1}{n} \sum_{r=1}^{n} X_{\ell kr} \\ \bar{X}_{\ell} &= \dfrac{1}{nb} \sum_{k=1}^{b} \sum_{r=1}^{n} X_{\ell kr} \\ \bar{X}_k &= \dfrac{1}{ng} \sum_{\ell=1}^{g} \sum_{r=1}^{n} X_{\ell kr} \\ \bar{X} &= \dfrac{1}{gnb} \sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} X_{\ell kr} \end{align} \]
ou
\[ \text{SS}_{\text{tot}} = \text{SS}_{\text{fat1}} + \text{SS}_{\text{fat2}} + \text{SS}_{\text{int}} + \text{SS}_{\text{res}} \]
Os graus de liberdade correspondentes associados às somas de quadrados na decomposição em (6-57) são:
\[ gbn - 1 = (g - 1) + (b - 1) + (g - 1) (b - 1) + gb(n - 1) \]
Temos que \(\{X_{\ell kr}\}_{\ell=1,k=1,r=1}^{g,b,n}\sim\mathcal{N}\text{IID}(\mu+\tau_{\ell}+\beta_k+\gamma_{\ell k}, \sigma^2)\).
A tabela ANOVA assume a seguinte forma:
Tabela de ANOVA para Modelo Univariado de Dois Fatores com Interação:
Fonte de Variação | Graus de Liberdade (df) | Soma de Quadrados (SS) | Quadrado Médio (MS) | F |
---|---|---|---|---|
Fator 1 | \(g - 1\) | \(\text{SS}_{\text{fat1}}\) | \(\dfrac{\text{SS}_{\text{fat1}}}{g-1}\) | \(\dfrac{\text{MS}_{\text{fat1}}}{\text{MS}_{\text{res}}}\) |
Fator 2 | \(b - 1\) | \(\text{SS}_{\text{fat2}}\) | \(\dfrac{\text{SS}_{\text{fat2}}}{b-1}\) | \(\dfrac{\text{MS}_{\text{fat2}}}{\text{MS}_{\text{res}}}\) |
Interação | \((g - 1)(b - 1)\) | \(\text{SS}_{\text{int}}\) | \(\dfrac{\text{SS}_{\text{int}}} {(g-1)(b-1)}\) | \(\dfrac{\text{MS}_{\text{int}}}{\text{MS}_{\text{res}}}\) |
Resíduo | \(gb(n - 1)\) | \(\text{SS}_{\text{res}}\) | \(\dfrac{\text{SS}_{\text{res}}}{gb(n-1)}\) | - |
Total | \(gbn - 1\) | \(\text{SS}_{\text{tot}}\) | - | - |
A tabela ANOVA necessita que \(n\ge2\), i.e., que o número de réplicas seja pelo menos igual a dois em todas as combinações dos níveis dos dois fatores (tratamentos), pois \(\frac{\text{SS}_{\text{res}}}{gb(n - 1)}\) faz parte do denominador das estatísticas de teste F.
Se \(g=b\), i.e., delineamento perfeitamente balanceado no número de níveis dos dois fatores, e um efeito não é significante, então sua exclusão do modelo não altera os resultados da tabela ANOVA.
Pelo teorema de Cochran, as razões F dos quadrados médios, \(\frac{\text{SS}_{\text{fat1}}}{g - 1}\), \(\frac{\text{SS}_{\text{fat2}}}{b - 1}\), e \(\frac{\text{SS}_{\text{int}}}{(g - 1)(b - 1)}\) em relação ao quadrado médio, \(\frac{\text{SS}_{\text{res}}}{gb(n - 1)}\), podem ser usadas para testar hierarquicamente os efeitos do fator 1, fator 2, e a interação fator 1-fator 2, respectivamente. (Veja [11] para uma discussão sobre análise de variância univariada bifatorial independente).
A hipótese nula omnibus é:
\[ \begin{align} H_{0}^{\text{omnibus}}&: \tau_1=\cdots=\tau_g=\beta_1=\cdots=\beta_b=\gamma_{11}=\cdots=\gamma_{gb}=0\\ H_{1}^{\text{omnibus}}&:\text{algum efeito populacional é diferente de zero} \end{align} \]
Se \(H_{0}^{\text{omnibus}}\) é rejeitada, testar o efeito de interação:
\[ H_{0}^{\text{interação}}: \gamma_{11}=\cdots=\gamma_{gb}=0 \]
Se \(H_{0}^{\text{interação}}\) não é rejeitada, testar os efeitos principais:
\[ H_{0}^{\text{Fator 1}}: \tau_1=\cdots=\tau_g=0\\ H_{0}^{\text{Fator 2}}: \beta_1=\cdots=\beta_b=0 \]
Se \(H_{0}^{\text{interação}}\) é rejeitada, testar os efeitos principais simples.
Figura 1. Árvore de Decisão para conduzir análises post-hoc após uma ANOVA de dois ou três fatores. Howell & Lacroix, 2012
“Primeiramente, numa ANOVA univariada bifatorial independente, teste a presença do efeito de interação. Se o efeito de interação for significante, então o teste de efeitos principais simples deve ser considerado.”
Woodward & Bonett, 1991, p. 256
FIG.1. Estratégia de teste de hipóteses para delineamento de três fatores (Keppel, 1973, p. 293). Woodward & Bonett, 1991
FIG. 2. Estratégia de teste de hipótese para delineamento de três fatores (adaptado de Milliken & Johnson, 1984, p. 198 - veja p. 114 para uma aplicação de análise de estrutura de tratamento de duas vias). Woodward & Bonett, 1991
Tabela 1. Representações univariada e multivariada do GLM (Modelo Linear General). Chartier & Faulkner, 2008
Uma forma de definir com mais precisão o número de variáveis dependentes (DV) é \(pq\). Se \(pq=1\), o modelo é univariado. Se \(pq\ge2\), o modelo é multivariado.
Figura 1. Diagramas de Venn para a) dois preditores intervalares em regressão, b) análise de variância unifatorial, c) ANOVA bifatorial e d) medidas repetidas. Chartier & Faulkner, 2008
porquinho-da-índia
Uma preocupação durante a Segunda Guerra Mundial era o fornecimento de vitamina C aos soldados e, nesse amplo contexto, os efeitos do ácido ascórbico e do suco de laranja foram estudados em animais. Um desses estudos foi:
Crampton EW (1947) O crescimento do odontoblasto dos dentes incisivos como critério de ingestão de vitamina C do porquinho-da-índia. The Journal of Nutrition 33(5): 491–504.
Sessenta porquinhos-da-índia (com 28 dias de idade) receberam um suplemento dietético de vitamina C em uma de três doses (0.5, 1 ou 2 mg/dia) administradas de uma de duas maneiras (como ácido ascórbico ou em suco de laranja). Havia dez porquinhos-da-índia, cinco machos e cinco fêmeas, em cada combinação de dose e método de administração. Após 42 dias na dieta, os porquinhos-da-índia foram sacrificados; os incisivos foram removidos e seccionados para obter medidas do comprimento dos odontoblastos — células que são importantes para o desenvolvimento dos dentes. Poderia haver múltiplas medições (odontoblastos) por animal, então os comprimentos foram médios para fornecer um único valor para cada animal. O resultado de interesse é o comprimento médio dos odontoblastos incisivos; isso foi medido em micrômetros. O interesse é em como a dose e o método de administração influenciam o resultado.
Você pode estar se perguntando o que um estudo de células em porquinhos-da-índia tem a ver com o fornecimento de vitamina C aos humanos. Os pesquisadores haviam estabelecido que o crescimento dos odontoblastos em porquinhos-da-índia era sensível à ingestão de vitamina C e, portanto, medir os odontoblastos poderia fornecer uma maneira de medir de forma confiável a absorção de vitamina C. Havia um debate sobre se as fontes ‘biológicas’ (por exemplo, suco de laranja) de vitamina C eram melhores do que as formas ‘químicas’ (por exemplo, ácido ascórbico), e isso foi explorado em um modelo animal.
ANOVA unifatorial independente de Fisher com o fator tipo de suplemento resulta na não significância deste fator (\(p=0.060\)) com \(\alpha=5\%\). Os resultados a seguir mostram que há o efeito de interação entre tipo de suplemento e dose (\(p=0.022\)) com \(\alpha=5\%\). Há efeitos simples de tipo de suplemento e dose. Então, suco de laranja tem efeito maior do que ácido ascórbico. As três doses têm efeitos progressivamente maiores.
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
alfa <- 0.05
# The Effect of Vitamin C on Tooth Growth in Guinea Pigs
# len: numeric Tooth length
# supp: factor Supplement type (VC or OJ).
# dose: numeric Dose in milligrams/day
Dados <- datasets::ToothGrowth
print.data.frame(Dados)
len supp dose
1 4.2 VC 0.5
2 11.5 VC 0.5
3 7.3 VC 0.5
4 5.8 VC 0.5
5 6.4 VC 0.5
6 10.0 VC 0.5
7 11.2 VC 0.5
8 11.2 VC 0.5
9 5.2 VC 0.5
10 7.0 VC 0.5
11 16.5 VC 1.0
12 16.5 VC 1.0
13 15.2 VC 1.0
14 17.3 VC 1.0
15 22.5 VC 1.0
16 17.3 VC 1.0
17 13.6 VC 1.0
18 14.5 VC 1.0
19 18.8 VC 1.0
20 15.5 VC 1.0
21 23.6 VC 2.0
22 18.5 VC 2.0
23 33.9 VC 2.0
24 25.5 VC 2.0
25 26.4 VC 2.0
26 32.5 VC 2.0
27 26.7 VC 2.0
28 21.5 VC 2.0
29 23.3 VC 2.0
30 29.5 VC 2.0
31 15.2 OJ 0.5
32 21.5 OJ 0.5
33 17.6 OJ 0.5
34 9.7 OJ 0.5
35 14.5 OJ 0.5
36 10.0 OJ 0.5
37 8.2 OJ 0.5
38 9.4 OJ 0.5
39 16.5 OJ 0.5
40 9.7 OJ 0.5
41 19.7 OJ 1.0
42 23.3 OJ 1.0
43 23.6 OJ 1.0
44 26.4 OJ 1.0
45 20.0 OJ 1.0
46 25.2 OJ 1.0
47 25.8 OJ 1.0
48 21.2 OJ 1.0
49 14.5 OJ 1.0
50 27.3 OJ 1.0
51 25.5 OJ 2.0
52 26.4 OJ 2.0
53 22.4 OJ 2.0
54 24.5 OJ 2.0
55 24.8 OJ 2.0
56 30.9 OJ 2.0
57 26.4 OJ 2.0
58 27.3 OJ 2.0
59 29.4 OJ 2.0
60 23.0 OJ 2.0
── Data Summary ────────────────────────
Values
Name Dados
Number of rows 60
Number of columns 3
_______________________
Column type frequency:
factor 1
numeric 2
________________________
Group variables None
── Variable type: factor ───────────────────────────────────────────────────────
skim_variable n_missing complete_rate ordered n_unique top_counts
1 supp 0 1 FALSE 2 OJ: 30, VC: 30
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100
1 len 0 1 18.8 7.65 4.2 13.1 19.2 25.3 33.9
2 dose 0 1 1.17 0.629 0.5 0.5 1 2 2
Dados$supp <- factor(Dados$supp,
levels=c("VC", "OJ"))
Dados$dose <- factor(Dados$dose)
boxplot(len~dose*supp,
data=Dados)
alfaBonf.supp <- alfa/length(unique(Dados$supp))
gplots::plotmeans(len~supp,
connect=FALSE,
col="black",
barcol="black",
p=1-alfaBonf.supp,
main="Porquinho-da-Índia",
data=Dados)
Registered S3 method overwritten by 'gplots':
method from
reorder.factor DescTools
Analysis of Variance Table
Response: len
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.35 3.6683 0.06039 .
Residuals 58 3246.9 55.98
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alfaBonf <- alfa/(length(unique(Dados$supp))*
length((unique(Dados$dose))))
gplots::plotmeans(len~interaction(dose,supp),
connect=FALSE,
col="black",
barcol="black",
p=1-alfaBonf,
main="Porquinho-da-Índia",
data=Dados)
[1] 18.81333
# A tibble: 2 × 2
supp len
<fct> <dbl>
1 VC 17.0
2 OJ 20.7
# A tibble: 3 × 2
dose len
<fct> <dbl>
1 0.5 10.6
2 1 19.7
3 2 26.1
# A tibble: 6 × 3
dose supp len
<fct> <fct> <dbl>
1 0.5 VC 7.98
2 1 VC 16.8
3 2 VC 26.1
4 0.5 OJ 13.2
5 1 OJ 22.7
6 2 OJ 26.1
# Somas de Quadrados
g <- length(unique(Dados$supp))
b <- length(unique(Dados$dose))
n <- nrow(Dados) / (g*b)
p <- 1
SS_total <- sum((Dados$len - X_bar)^2)
SS_A <- b * n * sum((X_bar_l$len - X_bar)^2)
SS_B <- g * n * sum((X_bar_k$len - X_bar)^2)
SS_error <- sum((Dados$len - rep(X_bar_lk$len, times=rep(n,g*b)))^2)
SS_AB <- SS_total - SS_A - SS_B - SS_error
# Graus de Liberdade
df_A <- g - 1
df_B <- b - 1
df_AB <- (g - 1)*(b - 1)
df_error <- g*b*(n - 1)
# Quadrado Médio
MS_A <- SS_A / df_A
MS_B <- SS_B / df_B
MS_AB <- SS_AB / df_AB
MS_error <- SS_error / df_error
# Estatística F
F_A <- MS_A / MS_error
F_B <- MS_B / MS_error
F_AB <- MS_AB / MS_error
# Valor-p (usando a distribuição F)
p_value_A <- formatC(1-pf(F_A, df_A, df_error),
format="e", digits=2)
p_value_B <- formatC(1-pf(F_B, df_B, df_error),
format="e", digits=2)
p_value_AB <- formatC(1-pf(F_AB, df_AB, df_error),
format="e", digits=2)
# Tabela ANOVA
anova_table <- data.frame(
Source = c("Factor A", "Factor B", "Interaction AB", "Error", "Total"),
SS = c(SS_A, SS_B, SS_AB, SS_error, SS_total),
df = c(df_A, df_B, df_AB, df_error, g*b*n - 1),
MS = c(MS_A, MS_B, MS_AB, MS_error, NA),
F = c(F_A, F_B, F_AB, NA, NA),
p_value = c(p_value_A, p_value_B, p_value_AB, NA, NA)
)
print(anova_table, digits=3)
Source SS df MS F p_value
1 Factor A 205 1 205.3 15.57 2.31e-04
2 Factor B 2426 2 1213.2 92.00 0.00e+00
3 Interaction AB 108 2 54.2 4.11 2.19e-02
4 Error 712 54 13.2 NA <NA>
5 Total 3452 59 NA NA <NA>
# Solução model.matrix e OLS
# https://www.maths.usyd.edu.au/u/UG/SM/STAT3022/r/current/Lecture/lecture18_2020JC.html#17
fit <- lm(len~supp*dose,
data=Dados)
print(anova(fit))
Analysis of Variance Table
Response: len
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.35 205.35 15.572 0.0002312 ***
dose 2 2426.43 1213.22 92.000 < 2.2e-16 ***
supp:dose 2 108.32 54.16 4.107 0.0218603 *
Residuals 54 712.11 13.19
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova Table (Type II tests)
Response: len
Sum Sq Df F value Pr(>F)
supp 205.35 1 15.572 0.0002312 ***
dose 2426.43 2 92.000 < 2.2e-16 ***
supp:dose 108.32 2 4.107 0.0218603 *
Residuals 712.11 54
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
y <- Dados$len
X0 <- model.matrix(~ 1, data=Dados)
XA <- model.matrix(~ -1 + supp, data=Dados)
XB <- model.matrix(~ -1 + dose, data=Dados)
XAB <- model.matrix(~ -1 + supp:dose, data=Dados)
PV0 <- X0 %*% solve(t(X0) %*% X0) %*% t(X0)
PWA <- XA %*% solve(t(XA) %*% XA) %*% t(XA) - PV0
PWB <- XB %*% solve(t(XB) %*% XB) %*% t(XB) - PV0
PVAB <- XAB %*% solve(t(XAB) %*% XAB) %*% t(XAB)
PWAB <- PVAB - PWA - PWB - PV0
PRes <- diag(nrow(Dados)) - PVAB
SSA <- t(y) %*% PWA %*% y
SSB <- t(y) %*% PWB %*% y
SSAB <- t(y) %*% PWAB %*% y
SSE <- t(y) %*% PRes %*% y
dfA <- sum(diag(PWA))
dfB <- sum(diag(PWB))
dfAB <- sum(diag(PWAB))
dfE <- sum(diag(PRes))
dfT <- dfA+dfB+dfAB+dfE
SST <- SSA+SSB+SSAB+SSE
R2A <- SSA/SST
R2B <- SSB/SST
R2AB <- SSAB/SST
R2E <- SSE/SST
R2 <- 1 - SSE/SST
eta2partialA <- R2A/(R2A+1-R2)
eta2partialB <- R2B/(R2B+1-R2)
eta2partialAB <- R2AB/(R2AB+1-R2)
FA <- (SSA/dfA)/(SSE/dfE)
pA <- formatC(1-pf(FA, dfA, dfE),
format="e", digits=2)
FB <- (SSB/dfB)/(SSE/dfE)
pB <- formatC(1-pf(FB, dfB, dfE),
format="e", digits=2)
FAB <- (SSAB/dfAB)/(SSE/dfE)
pAB <- formatC(1-pf(FAB, dfAB, dfE),
format="e", digits=2)
R2 <- 1 - SSE/SST
df1 <- dfT-dfE
df2 <- dfE
F_omnibus <- (R2/df1)/((1-R2)/df2) # Pestana & Gageiro, 2005, p. 77
pv <- 1-pf(F_omnibus, df1, df2)
if(pv<2.2e-16) p_omnibus <- "< 2.2e-16"
if(pv>2.2e-16) p_omnibus <- paste0("< ", formatC(pv, format="e", digits=2))
Fcrit <- qf(1-alfa, df1, df2)
cat("Omnibus F(",df1, ",", df2, ") = ", round(Fcrit,4),
", F = ", round(F_omnibus,4),
", p = ", p_omnibus, ", R^2 = eta^2 = ", round(R2,4), sep="")
Omnibus F(5,54) = 2.3861, F = 41.5572, p = < 2.2e-16, R^2 = eta^2 = 0.7937
ANOVA2wtable <- data.frame("Source"=c("supp", "dose", "supp:dose", "Error", "Total"),
# "df"=c(qr(PWA)$rank, qr(PWB)$rank, qr(PWAB)$rank, qr(PRes)$rank),
"df"=c(dfA, dfB, dfAB, dfE, dfT),
"SS"=c(SSA, SSB, SSAB, SSE, SST),
"F"=c(FA, FB, FAB, NA, NA),
"p"=c(pA, pB, pAB, NA, NA),
"R^2"=c(R2A, R2B, R2AB, R2E, NA),
"eta^2 partial"=c(eta2partialA, eta2partialB, eta2partialAB, NA, NA))
print(ANOVA2wtable, digits=2)
Source df SS F p R.2 eta.2.partial
1 supp 1 205 15.6 2.31e-04 0.059 0.22
2 dose 2 2426 92.0 0.00e+00 0.703 0.77
3 supp:dose 2 108 4.1 2.19e-02 0.031 0.13
4 Error 54 712 NA <NA> 0.206 NA
5 Total 59 3452 NA <NA> NA NA
eta2 <- effectsize::eta_squared(anv,
partial=FALSE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=3)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 | 98.3333% CI | interpret
----------------------------------------------
supp | 0.059 | [0.000, 0.253] | small
dose | 0.703 | [0.524, 0.802] | large
supp:dose | 0.031 | [0.000, 0.185] | small
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=3)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 98.3333% CI | interpret
-------------------------------------------------------
supp | 0.224 | [0.034, 0.438] | large
dose | 0.773 | [0.630, 0.850] | large
supp:dose | 0.132 | [0.000, 0.335] | medium
# Regressão Linear Múltipla de two-way ANOVA
fit <- lm(len~supp*dose,
data=Dados)
print(car::Anova(fit))
Anova Table (Type II tests)
Response: len
Sum Sq Df F value Pr(>F)
supp 205.35 1 15.572 0.0002312 ***
dose 2426.43 2 92.000 < 2.2e-16 ***
supp:dose 108.32 2 4.107 0.0218603 *
Residuals 712.11 54
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = len ~ supp * dose, data = Dados)
Residuals:
Min 1Q Median 3Q Max
-8.20 -2.72 -0.27 2.65 8.27
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.980 1.148 6.949 4.98e-09 ***
suppOJ 5.250 1.624 3.233 0.00209 **
dose1 8.790 1.624 5.413 1.46e-06 ***
dose2 18.160 1.624 11.182 1.13e-15 ***
suppOJ:dose1 0.680 2.297 0.296 0.76831
suppOJ:dose2 -5.330 2.297 -2.321 0.02411 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.631 on 54 degrees of freedom
Multiple R-squared: 0.7937, Adjusted R-squared: 0.7746
F-statistic: 41.56 on 5 and 54 DF, p-value: < 2.2e-16
y <- Dados$len
X <- model.matrix(len~supp*dose, data=Dados)
coeff <- solve(t(X) %*% X) %*% t(X) %*% y # OLS
cat("\nCoefficients:\n")
Coefficients:
Estimate
(Intercept) 7.98
suppOJ 5.25
dose1 8.79
dose2 18.16
suppOJ:dose1 0.68
suppOJ:dose2 -5.33
[1] 2740.103
[1] 712.106
[1] 3452.209
R^2 = eta^2 omnibus = 0.7937246
df1 <- dim(X)[2] - 1
N <- dim(X)[1]
df2 <- N-df1-1
F_omnibus <- (R2/df1)/((1-R2)/df2) # Pestana & Gageiro, 2005, p. 77
pv <- 1-pf(F_omnibus, df1, df2)
if(pv<2.2e-16) p_omnibus <- "< 2.2e-16"
if(pv>2.2e-16) p_omnibus <- paste0("< ", formatC(pv, format="e", digits=2))
Fcrit <- qf(1-alfa, df1, df2)
cat("Omnibus F(",df1, ",", df2, ") = ", round(Fcrit,4),
", F = ", round(F_omnibus,4),
", p ", p_omnibus, ", R^2 = eta^2 = ", round(R2,4), sep="")
Omnibus F(5,54) = 2.3861, F = 41.5572, p < 2.2e-16, R^2 = eta^2 = 0.7937
# print(sjPlot::tab_model(fit,
# p.adjust="holm",
# encoding="Windows-1252",
# title="Effect of Vitamin C on Tooth Growth in Guinea Pigs: holm",
# show.reflvl=TRUE,
# show.ngroups=TRUE,
# show.stat=TRUE,
# show.df=TRUE,
# show.se=TRUE,
# digits=4,
# digits.p=4,
# digits.rsq=6,
# digits.re=6))
# Solução GLM
fit <- lm(len~supp,
data=Dados)
print(anova(fit))
Analysis of Variance Table
Response: len
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.35 3.6683 0.06039 .
Residuals 58 3246.9 55.98
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table
Response: len
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.35 205.35 14.017 0.0004293 ***
dose 2 2426.43 1213.22 82.811 < 2.2e-16 ***
Residuals 56 820.43 14.65
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table
Response: len
Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.35 205.35 15.572 0.0002312 ***
dose 2 2426.43 1213.22 92.000 < 2.2e-16 ***
supp:dose 2 108.32 54.16 4.107 0.0218603 *
Residuals 54 712.11 13.19
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova Table (Type II tests)
Response: len
Sum Sq Df F value Pr(>F)
supp 205.35 1 15.572 0.0002312 ***
dose 2426.43 2 92.000 < 2.2e-16 ***
supp:dose 108.32 2 4.107 0.0218603 *
Residuals 712.11 54
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eta2 <- effectsize::eta_squared(anv,
partial=FALSE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=3)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 | 98.3333% CI | interpret
----------------------------------------------
supp | 0.059 | [0.000, 0.253] | small
dose | 0.703 | [0.524, 0.802] | large
supp:dose | 0.031 | [0.000, 0.185] | small
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=3)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 98.3333% CI | interpret
-------------------------------------------------------
supp | 0.224 | [0.034, 0.438] | large
dose | 0.773 | [0.630, 0.850] | large
supp:dose | 0.132 | [0.000, 0.335] | medium
o.par <- par()
alfaBonf_int <- alfa/(length(unique(Dados$supp))*
length((unique(Dados$dose))))
fit.means <- phia::interactionMeans(fit)
plot(fit.means,
errorbar=paste0("ci",
round((1-alfaBonf)*100,4)),
abbrev.levels=FALSE)
par(o.par)
alfaBonf_int <- alfa/(length(unique(Dados$supp))*
length((unique(Dados$dose))))
plot(effects::effect(c("supp", "dose"), fit, confidence.level=1-alfaBonf_int),
multiline=TRUE, ci.style="bars")
NOTE: suppdose is not a high-order term in the model
alfaBonf_fat1 <- alfa/(length(unique(Dados$supp)))
plot(effects::effect(c("supp"), fit, confidence.level=1-alfaBonf_fat1),
ci.style="bars")
NOTE: supp is not a high-order term in the model
alfaBonf_fat2 <- alfa/(length(unique(Dados$dose)))
plot(effects::effect(c("dose"), fit, confidence.level=1-alfaBonf_fat2),
ci.style="bars")
NOTE: dose is not a high-order term in the model
supp simple effect test in the presence of interaction with dose
F Test:
P-value adjustment method: holm
Value SE Df Sum of Sq F Pr(>F)
0.5 -5.25 1.624 1 137.81 10.4505 0.004185 **
1 -5.93 1.624 1 175.82 13.3330 0.001769 **
2 0.08 1.624 1 0.03 0.0024 0.960893
Residuals 54 712.11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dose simple effect test in the presence of interaction with supp
F Test:
P-value adjustment method: holm
dose1 dose2 SE1 SE2 Df Sum of Sq F Pr(>F)
VC -18.16 -9.37 1.624 1.624 2 1649.49 62.541 1.751e-14 ***
OJ -12.83 -3.36 1.624 1.624 2 885.26 33.565 3.363e-10 ***
Residuals 54 712.11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$emmeans
supp dose emmean SE df lower.CL upper.CL
VC 0.5 7.98 1.15 54 5.68 10.3
OJ 0.5 13.23 1.15 54 10.93 15.5
VC 1 16.77 1.15 54 14.47 19.1
OJ 1 22.70 1.15 54 20.40 25.0
VC 2 26.14 1.15 54 23.84 28.4
OJ 2 26.06 1.15 54 23.76 28.4
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
VC dose0.5 - OJ dose0.5 -5.25 1.62 54 -3.233 0.0105
VC dose0.5 - VC dose1 -8.79 1.62 54 -5.413 <.0001
VC dose0.5 - OJ dose1 -14.72 1.62 54 -9.064 <.0001
VC dose0.5 - VC dose2 -18.16 1.62 54 -11.182 <.0001
VC dose0.5 - OJ dose2 -18.08 1.62 54 -11.133 <.0001
OJ dose0.5 - VC dose1 -3.54 1.62 54 -2.180 0.1346
OJ dose0.5 - OJ dose1 -9.47 1.62 54 -5.831 <.0001
OJ dose0.5 - VC dose2 -12.91 1.62 54 -7.949 <.0001
OJ dose0.5 - OJ dose2 -12.83 1.62 54 -7.900 <.0001
VC dose1 - OJ dose1 -5.93 1.62 54 -3.651 0.0035
VC dose1 - VC dose2 -9.37 1.62 54 -5.770 <.0001
VC dose1 - OJ dose2 -9.29 1.62 54 -5.720 <.0001
OJ dose1 - VC dose2 -3.44 1.62 54 -2.118 0.1346
OJ dose1 - OJ dose2 -3.36 1.62 54 -2.069 0.1346
VC dose2 - OJ dose2 0.08 1.62 54 0.049 0.9609
P value adjustment: holm method for 15 tests
supp dose emmean SE df lower.CL upper.CL .group
VC 0.5 7.98 1.15 54 4.83 11.1 a
OJ 0.5 13.23 1.15 54 10.08 16.4 b
VC 1 16.77 1.15 54 13.62 19.9 b
OJ 1 22.70 1.15 54 19.55 25.8 c
OJ 2 26.06 1.15 54 22.91 29.2 c
VC 2 26.14 1.15 54 22.99 29.3 c
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 6 estimates
P value adjustment: holm method for 15 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
NOTE: Results may be misleading due to involvement in interactions
$emmeans
supp emmean SE df lower.CL upper.CL
VC 17.0 0.663 54 15.6 18.3
OJ 20.7 0.663 54 19.3 22.0
Results are averaged over the levels of: dose
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
VC - OJ -3.7 0.938 54 -3.946 0.0002
Results are averaged over the levels of: dose
supp emmean SE df lower.CL upper.CL .group
VC 17.0 0.663 54 15.4 18.5 a
OJ 20.7 0.663 54 19.1 22.2 b
Results are averaged over the levels of: dose
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 2 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
NOTE: Results may be misleading due to involvement in interactions
$emmeans
dose emmean SE df lower.CL upper.CL
0.5 10.6 0.812 54 8.98 12.2
1 19.7 0.812 54 18.11 21.4
2 26.1 0.812 54 24.47 27.7
Results are averaged over the levels of: supp
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
dose0.5 - dose1 -9.13 1.15 54 -7.951 <.0001
dose0.5 - dose2 -15.49 1.15 54 -13.493 <.0001
dose1 - dose2 -6.37 1.15 54 -5.543 <.0001
Results are averaged over the levels of: supp
P value adjustment: holm method for 3 tests
dose emmean SE df lower.CL upper.CL .group
0.5 10.6 0.812 54 8.6 12.6 a
1 19.7 0.812 54 17.7 21.7 b
2 26.1 0.812 54 24.1 28.1 c
Results are averaged over the levels of: supp
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 3 estimates
P value adjustment: holm method for 3 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
\(p\ge2, q=1, g\ge2, b\ge2, n\ge2\)
Carolyn J. Anderson, 2017, Profile Analysis and 2-way MANOVA: https://education.illinois.edu/docs/default-source/carolyn-anderson/edpsy584/lectures/manova_part2_beamer-online.pdf
Procedendo por analogia, especificamos o modelo de efeitos fixos de dois fatores para uma resposta vetorial consistindo em \(p\) componentes [veja (6-54)]
\[ \mathbf{X}_{\ell kr} = \boldsymbol{\mu} + \boldsymbol{\tau}_\ell + \boldsymbol{\beta}_k + \boldsymbol{\gamma}_{\ell k} + \mathbf{e}_{\ell kr} \]
em que \(k = 1,2,..., b\) e \(r = 1,2,...,n\) e
\[ \sum_{\ell=1}^{g} \boldsymbol{\tau}_\ell = \sum_{k=1}^{b} \boldsymbol{\beta}_k = \sum_{\ell=1}^{g} \boldsymbol{\gamma}_{\ell k} = \sum_{k=1}^{b} \boldsymbol{\gamma}_{\ell k} = \mathbf{0} \]
Os vetores são todos da ordem \(p \times 1\), e os \(\mathbf{e}_{\ell kr}\) são vetores aleatórios independentes \(\mathcal{N}_p(\boldsymbol{0}, \boldsymbol{\Sigma})\). Assim, as respostas consistem em \(p\) medições replicadas \(n\) vezes em cada uma das possíveis combinações de níveis dos fatores 1 e 2.
Seguindo (6-56), podemos decompor os vetores de observação \(\mathbf{X}_{\ell kr}\) como:
\[ \begin{align} \mathbf{X}_{\ell kr} &= \mathbf{\bar{X}} + (\mathbf{\bar{X}}_\ell - \mathbf{\bar{X}}) + (\mathbf{\bar{X}}_k - \mathbf{\bar{X}})+ \\ &\quad (\mathbf{\bar{X}}_{\ell k} - \mathbf{\bar{X}}_\ell - \mathbf{\bar{X}}_k + \mathbf{\bar{X}}) + (\mathbf{X}_{\ell kr} - \mathbf{\bar{X}}_{\ell k}) \end{align} \tag{6-60} \]
em que \(\mathbf{\bar{X}}\) é a média geral dos vetores de observação, \(\mathbf{\bar{X}}_\ell\) é a média dos vetores de observação no \(\ell\)-ésimo nível do fator 1, \(\mathbf{\bar{X}}_k\) é a média dos vetores de observação no \(k\)-ésimo nível do fator 2, e \(\mathbf{\bar{X}}_{\ell k}\) é a média dos vetores de observação no \(\ell\)-ésimo nível do fator 1 e no \(k\)-ésimo nível do fator 2.
Generalizações diretas de (6-57) e (6-58) fornecem as decomposições da soma dos quadrados e produtos cruzados e graus de liberdade:
\[ \begin{align} \sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} (\mathbf{X}_{\ell kr} - \mathbf{\bar{X}})(\mathbf{X}_{\ell kr} - \mathbf{\bar{X}})^{\prime} &= b n \sum_{\ell=1}^{g} (\mathbf{\bar{X}}_\ell - \mathbf{\bar{X}})(\mathbf{\bar{X}}_\ell - \mathbf{\bar{X}})^{\prime}+ \\ &\quad g n \sum_{k=1}^{b} (\mathbf{\bar{X}}_k - \mathbf{\bar{X}})(\mathbf{\bar{X}}_k - \mathbf{\bar{X}})^{\prime}+\\ &\quad n \sum_{\ell=1}^{g} \sum_{k=1}^{b} (\mathbf{\bar{X}}_{\ell k} - \mathbf{\bar{X}}_\ell - \mathbf{\bar{X}}_k + \mathbf{\bar{X}})(\mathbf{\bar{X}}_{\ell k} - \mathbf{\bar{X}}_\ell - \mathbf{\bar{X}}_k + \mathbf{\bar{X}})^{\prime}+ \\ &\quad \sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} (\mathbf{X}_{\ell kr} - \mathbf{\bar{X}}_{\ell k})(\mathbf{X}_{\ell kr} - \mathbf{\bar{X}}_{\ell k})^{\prime} \end{align} \tag{6-61} \]
Os graus de liberdade correspondentes são:
\[ gbn - 1 = (g - 1) + (b - 1) + (g - 1)(b - 1) + gb(n - 1) \tag{6-62} \]
Mais uma vez, a generalização da análise univariada para a análise multivariada consiste simplesmente em substituir um escalar, como \((x_\ell - \bar{x})^2\), pela matriz correspondente \((\mathbf{X}_\ell - \mathbf{\bar{X}})(\mathbf{X}_\ell - \mathbf{\bar{X}})^{\prime}\).
Traduzindo e adaptando para a notação solicitada:
Tabela MANOVA
Fonte de Variação | Matriz de Soma de Quadrados e Produtos Cruzados (SSP) | Graus de Liberdade (df) |
---|---|---|
Fator 1 | \(\text{SSP}_{\text{fat1}} = \sum_{\ell=1}^{g} bn (\mathbf{\bar{X}}_{\ell} - \mathbf{\bar{X}})(\mathbf{\bar{X}}_{\ell} - \mathbf{\bar{X}})^{\prime}\) | \(g - 1\) |
Fator 2 | \(\text{SSP}_{\text{fat2}} = \sum_{k=1}^{b} gn (\mathbf{\bar{X}}_k - \mathbf{\bar{X}})(\mathbf{\bar{X}}_k - \mathbf{\bar{X}})^{\prime}\) | \(b - 1\) |
Interação | \(\text{SSP}_{\text{int}} = \sum_{\ell=1}^{g} \sum_{k=1}^{b} n (\mathbf{\bar{X}}_{\ell k} - \mathbf{\bar{X}}_{\ell} - \mathbf{\bar{X}}_k + \mathbf{\bar{X}})(\mathbf{\bar{X}}_{\ell k} - \mathbf{\bar{X}}_{\ell} - \mathbf{\bar{X}}_k + \mathbf{\bar{X}})^{\prime}\) | \((g - 1)(b - 1)\) |
Resíduo | \(\text{SSP}_{\text{res}} = \sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} (\mathbf{X}_{\ell kr} - \mathbf{\bar{X}}_{\ell k})(\mathbf{X}_{\ell kr} - \mathbf{\bar{X}}_{\ell k})^{\prime}\) | \(gb(n - 1)\) |
Total | \(\text{SSP}_{\text{tot}} = \sum_{\ell=1}^{g} \sum_{k=1}^{b} \sum_{r=1}^{n} (\mathbf{X}_{\ell kr} - \mathbf{\bar{X}})(\mathbf{X}_{\ell kr} - \mathbf{\bar{X}})^{\prime}\) | \(gbn - 1\) |
O teste de razão de verossimilhança de (sem efeitos de interação):
\[ \begin{align} H_0&: \boldsymbol{\gamma}_{11} = \boldsymbol{\gamma}_{12} = \cdots = \boldsymbol{\gamma}_{gb} = \mathbf{0} \\ \text{vs.}\\ H_1&: \text{Pelo menos um }\boldsymbol{\gamma}_{\ell k} \ne \mathbf{0} \end{align} \tag{6-63} \]
É conduzido rejeitando \(H_0\) para pequenos valores da razão:
\[ \Lambda^* = \dfrac{|\text{SSP}_\text{res}|}{|\text{SSP}_\text{int} + \text{SSP}_\text{res}|} \tag{6-64} \]
O procedimento de teste de razão de verossimilhança (LR) requer que \(p \leq gb(n - 1)\), para que o \(\text{SSP}_{\text{res}}\) seja definida positiva (com probabilidade 1).
Para amostra grande, \(gbn-p-g-b\ge30\), o lambda de Wilks, \(\Lambda^*\), pode ser referido a um percentil de qui-quadrado.
Usando o multiplicador de Bartlett (veja [6]) para melhorar a aproximação qui-quadrada, rejeitamos \(H_0: \boldsymbol{\gamma}_{11} = \boldsymbol{\gamma}_{12} = \cdots = \boldsymbol{\gamma}_{gb} = \mathbf{0}\) no nível \(\alpha\) se:
\[ X^2 = -\left(gb(n-1)-\dfrac{p+1-(g-1)(b-1)}{2}\right)\ln(\Lambda^*) > \chi^2_{(g-1)(b-1)p}(1-\alpha) \tag{6-65} \]
Normalmente, o teste para interação é realizado antes dos testes para os efeitos principais do fator. Se existirem efeitos de interação, os efeitos do fator não têm uma interpretação clara! Do ponto de vista prático, não é aconselhável prosseguir com os testes adicionais de variáveis. Em vez disso, \(p \ge 2\) análises de variância bifatorial independente univariada (uma para cada resposta) são frequentemente e equivocadamente conduzidas para ver se a interação aparece em algumas respostas, mas não em outras. Aquelas respostas sem interação podem ser interpretadas em termos de efeitos aditivos dos fatores 1 e 2, desde que esses últimos efeitos existam. Em qualquer caso, gráficos de interação semelhantes à Figura 6.3, mas com médias de amostra de tratamento substituindo os valores esperados, esclarecem melhor as magnitudes relativas dos efeitos principais e de interação.
No modelo multivariado, testamos os efeitos principais dos fatores 1 e 2 da seguinte maneira. Primeiro, considere as hipóteses \(H_{\mathbf{0}}: \boldsymbol{\tau}_1 = \boldsymbol{\tau}_2 = \cdots = \boldsymbol{\tau}_g = \mathbf{0}\) e \(H_1:\) pelo menos um \(\boldsymbol{\tau}_\ell \ne\mathbf{0}\) é diferente de 0. Essas hipóteses nula e alternativa especificam que não há efeitos do fator 1 e algum efeito do fator 1, respectivamente. Seja
\[ \Lambda^* = \dfrac{|\text{SSP}_\text{res}|}{|\text{SSP}_\text{fat1} + \text{SSP}_\text{res}|} \tag{6-66} \]
de modo que valores pequenos de \(\Lambda^*\) são consistentes com \(H_1\). Para amostra grande, \(gn-p-g\ge30\), e usando a correção de Bartlett, o teste de razão de verossimilhança é o seguinte:
Rejeite \(H_{\mathbf{0}}: \boldsymbol{\tau}_1 = \boldsymbol{\tau}_2 =\cdots= \boldsymbol{\tau}_g = \mathbf{0}\) (sem efeitos do fator 1) no nível \(\alpha\) se
\[ X^2 = -\left(gb(n-1)-\dfrac{p+1-(g-1)}{2}\right)\ln(\Lambda^*) > \chi^2_{(g-1)p}(1-\alpha) \tag{6-67} \]
De maneira semelhante, os efeitos do fator 2 são testados considerando \(H_{\mathbf{0}}: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_2 = \cdots = \boldsymbol{\beta}_b = \mathbf{0}\) e \(H_1\) onde pelo menos um \(\boldsymbol{\beta}_k \ne\mathbf{0}\). Valores pequenos de
\[ \Lambda^* =\dfrac{|\text{SSP}_\text{res}|}{|\text{SSP}_\text{fat2} + \text{SSP}_\text{res}|} \tag{6-68} \]
são consistentes com \(H_1\). Novamente, para amostra grande, \(bn-p-b\ge30\), e usando a correção de Bartlett:
Rejeite \(H_{\mathbf{0}}: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_2 = \cdots = \boldsymbol{\beta}_b = \mathbf{0}\) (sem efeitos do fator 2) no nível \(\alpha\) se
\[ X^2 = -\left(gb(n-1)-\dfrac{p+1-(b-1)}{2}\right)\ln(\Lambda^*) > \chi^2_{(b-1)p}(1-\alpha) \tag{6-69} \]
Comentário. Consideramos o modelo multivariado de dois fatores com replicações. Ou seja, o modelo permite \(n\ge2\) replicações das respostas em cada combinação de níveis de fatores. Isso nos permite examinar a interação dos fatores. Se apenas um vetor de observação estiver disponível em pelo menos uma combinação de níveis de fatores, o modelo de dois fatores não permite a possibilidade de um termo de interação geral \(\gamma_{\ell k}\). A tabela MANOVA correspondente inclui apenas o fator 1, o fator 2 e as fontes residuais de variação como componentes da variação total. (Veja o Exercício 6.13).
As condições ótimas para extrudir filme plástico foram examinadas usando uma técnica chamada Operação Evolutiva. (Veja [9].) No decorrer do estudo que foi realizado, três respostas - \(X_1\) = resistência ao rasgo, \(X_2\) = brilho, e \(X_3\) = opacidade - foram medidas em dois níveis dos fatores, taxa de extrusão e quantidade de um aditivo. As medições foram repetidas \(n = 5\) vezes em cada combinação dos níveis de fatores. Os dados são exibidos na Tabela 6.4.
porquinho-da-índia
Extrudir é um processo utilizado na produção de objetos com perfil contínuo e homogêneo. Durante a extrusão, um material é empurrado ou puxado através de um molde (conhecido como matriz) que possui um orifício com o perfil desejado, resultando em um produto com a forma desse orifício. Um exemplo comum de extrusão é a produção de massas alimentícias, onde a massa é pressionada através de uma matriz para formar diferentes formatos, como espaguete ou macarrão. Em contextos industriais, a extrusão é amplamente utilizada na produção de itens como tubos, hastes, fios e filmes plásticos. O processo pode ser usado com uma variedade de materiais, incluindo plásticos, metais e cerâmicas.
O valor p omnibus necessita de regressão multivariada múltipla.
porquinho-da-índia
Para testar a interação, calculamos
\[ \Lambda^*=\dfrac{|\text{SSP}_\text{res}|}{|\text{SSP}_\text{int}+\text{SSP}_\text{res}|}=\dfrac{275.7098}{354.7906}=0.7771 \]
Para \((g - l)(b - 1) = 1\),
\[ F=\left(\dfrac{1-\Lambda^*}{\Lambda^*}\right)\dfrac{(gb(n-1)-p+1)/2}{(|(g-1)(b-1)-p|+1)/2} \]
possui uma distribuição F exata com graus de liberdade \(\nu_1 = |(g - l)(b - 1) - p| + 1\) e \(\nu_2=gb(n - 1) - p + 1\). (Veja [1]).
Para o nosso exemplo:
\[ F = \dfrac{1 - 0.7771}{0.7771} \dfrac{(2(2)(4)-3+1)/2}{(|1(1)-3|+1)/2} \]
Uma vez que \(F = 1.34 < F_{3,14}(0.95) = 3.34\), não rejeitamos a hipótese \(H_0: \boldsymbol{\gamma}_{11} = \boldsymbol{\gamma}_{12} = \boldsymbol{\gamma}_{21} = \boldsymbol{\gamma}_{22} = \mathbf{0}\) (sem efeitos de interação).
Observe que a estatística qui-quadrado aproximada para este teste é \(- [2(2)(4) \times (3 + 1 - 2(1))/2] \ln(0.7771) = 3.66\), a partir de (6-65). Uma vez que \(\chi^2(0.95) = 7.81\), chegaríamos à mesma conclusão que foi fornecida pelo teste F exato.
Para testar os efeitos dos fatores 1 e 2 (veja página 317), calculamos:
\[ \Lambda^*_1 = \dfrac{{|\text{SSP}_{\text{res}}|}}{|\text{SSP}_{\text{fac1}} + \text{SSP}_{\text{res}}|} = \dfrac{{275.7098}}{722.0212} = 0.3819 \]
e
\[ \Lambda^*_2 = \dfrac{{|\text{SSP}_{\text{res}}|}}{|\text{SSP}_{\text{fac2}} + \text{SSP}_{\text{res}}|} = \dfrac{{275.7098}}{527.1347} = 0.5230 \]
Para ambos \(g - 1 = 1\) e \(b - 1 = 1\),
\[ F_1 = \left(\dfrac{1 - \Lambda^*_1}{\Lambda^*_1}\right) \dfrac{(gb(n-1) - p + 1)/2}{(|(g-1)(b-1) - p| + 1)/2} \]
e
\[ F_2 = \left(\dfrac{1 - \Lambda^*_2}{\Lambda^*_2}\right) \dfrac{(gb(n-1) - p + 1)/2}{(|(g-1)(b-1) - p| + 1)/2} \]
possuem distribuições F exatas com graus de liberdade \(\nu_1 = |(g - 1)(b - 1) - p| + 1\) e \(\nu_2 = gb(n - 1) - p + 1\), respectivamente. (Veja [1]).
No nosso caso,
\[ F_1 = \dfrac{1 - 0.3819}{0.3819} \dfrac{(2(2)(4) - 3 + 1)/2}{(|(2 - 1)(2 - 1) - 3| + 1)/2} = 7.55 \]
e
\[ F_2 = \dfrac{1 - 0.5230}{0.5230} \dfrac{(2(2)(4) - 3 + 1)/2}{(|(2 - 1)(2 - 1) - 3| + 1)/2} = 4.26 \]
De antes, \(F_{3,14}(0.95) = 3.34\). Temos \(F_1 = 7.55 > F_{3,14}(0.95) = 3.34\), e portanto, rejeitamos \(H_0: \boldsymbol{\tau}_1 = \boldsymbol{\tau}_2 = \mathbf{0}\) (sem efeitos do fator 1) no nível de 5%. Similarmente, \(F_2 = 4.26 > F_{3,14}(0.95) = 3.34\), e rejeitamos \(H_0: \boldsymbol{\beta}_1 = \boldsymbol{\beta}_2 = \mathbf{0}\) (sem efeitos do fator 2) no nível de 5%. Concluímos que tanto a mudança na taxa de extrusão quanto a quantidade de aditivo afetam as respostas, e eles o fazem de maneira aditiva. Outra solução é usar a correção de Bonferroni ou de Kimball (Sidak) para testar os dois efeitos principais.
A natureza dos efeitos dos fatores 1 e 2 nas respostas é explorada no Exercício 6.15. Nesse exercício, intervalos de confiança simultâneos para contrastes nos componentes de \(\boldsymbol{\tau}_\ell\) e \(\boldsymbol{\beta}_k\) são considerados.
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
x <- read.table("JW6Data/T6-4.dat", quote="\"", comment.char="")
alfa <- 0.05
names(x) <- c("F1", "F2", "X1", "X2", "X3")
x$F1 <- factor(x$F1)
x$F2 <- factor(x$F2)
print.data.frame(x)
F1 F2 X1 X2 X3
1 0 0 6.5 9.5 4.4
2 0 0 6.2 9.9 6.4
3 0 0 5.8 9.6 3.0
4 0 0 6.5 9.6 4.1
5 0 0 6.5 9.2 0.8
6 0 1 6.9 9.1 5.7
7 0 1 7.2 10.0 2.0
8 0 1 6.9 9.9 3.9
9 0 1 6.1 9.5 1.9
10 0 1 6.3 9.4 5.7
11 1 0 6.7 9.1 2.8
12 1 0 6.6 9.3 4.1
13 1 0 7.2 8.3 3.8
14 1 0 7.1 8.4 1.6
15 1 0 6.8 8.5 3.4
16 1 1 7.1 9.2 8.4
17 1 1 7.0 8.8 5.2
18 1 1 7.2 9.7 6.9
19 1 1 7.5 10.1 2.7
20 1 1 7.6 9.2 1.9
F1 F2 X1 X2 X3
0:10 0:10 Min. :5.800 Min. : 8.300 Min. :0.800
1:10 1:10 1st Qu.:6.500 1st Qu.: 9.100 1st Qu.:2.525
Median :6.850 Median : 9.350 Median :3.850
Mean :6.785 Mean : 9.315 Mean :3.935
3rd Qu.:7.125 3rd Qu.: 9.625 3rd Qu.:5.325
Max. :7.600 Max. :10.100 Max. :8.400
X1 X2 X3
X1 1.000 -0.169 -0.013
X2 -0.169 1.000 0.098
X3 -0.013 0.098 1.000
Type II MANOVA Tests: Wilks test statistic
Df test stat approx F num Df den Df Pr(>F)
F1 1 0.382 7.55 3 14 0.003 **
F2 1 0.523 4.26 3 14 0.025 *
F1:F2 1 0.777 1.34 3 14 0.302
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
F1 1 0.618 7.55 3 14 0.003 **
F2 1 0.477 4.26 3 14 0.025 *
F1:F2 1 0.223 1.34 3 14 0.302
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa/3,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=3)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 98.3333% CI | interpret
-------------------------------------------------------
F1 | 0.618 | [0.065, 0.813] | large
F2 | 0.477 | [0.000, 0.735] | large
F1:F2 | 0.223 | [0.000, 0.558] | large
# valor p omnibus necessita de regressão multivariada múltipla
# Função para calcular a matriz SSP para o Fator 1
SSP_fat1 <- function(data) {
g <- length(unique(data$F1))
b <- length(unique(data$F2))
n <- nrow(data) / (g * b)
p <- ncol(data[, grep("X",colnames(data))])
ssp_matrix <- matrix(0, p, p)
for (i in unique(data$F1)) {
current_subset <- subset(data, F1 == i)
mean_diff <- colMeans(current_subset[, grep("X",colnames(data))]) -
colMeans(data[, grep("X",colnames(data))])
ssp_matrix <- ssp_matrix + b * n * (matrix(mean_diff) %*% t(matrix(mean_diff)))
}
return(ssp_matrix)
}
# Função para calcular a matriz SSP para o Fator 2
SSP_fat2 <- function(data) {
g <- length(unique(data$F1))
b <- length(unique(data$F2))
n <- nrow(data) / (g * b)
p <- ncol(data[, grep("X",colnames(data))])
ssp_matrix <- matrix(0, p, p)
for (i in unique(data$F2)) {
current_subset <- subset(data, F2 == i)
mean_diff <- colMeans(current_subset[, grep("X",colnames(data))]) -
colMeans(data[, grep("X",colnames(data))])
ssp_matrix <- ssp_matrix + g * n * (matrix(mean_diff) %*% t(matrix(mean_diff)))
}
return(ssp_matrix)
}
# Função para calcular a matriz SSP de Interação
SSP_int <- function(data) {
g <- length(unique(data$F1))
b <- length(unique(data$F2))
n <- nrow(data) / (g * b)
p <- ncol(data[, grep("X",colnames(data))])
ssp_matrix <- matrix(0, p, p)
for (i in unique(data$F1)) {
for (j in unique(data$F2)) {
current_subset <- subset(data, F1 == i & F2 == j)
mean_diff <- colMeans(current_subset[, grep("X",colnames(data))]) -
colMeans(subset(data, F1 == i)[, grep("X",colnames(data))]) -
colMeans(subset(data, F2 == j)[, grep("X",colnames(data))]) +
colMeans(data[, grep("X",colnames(data))])
ssp_matrix <- ssp_matrix + n * (matrix(mean_diff) %*% t(matrix(mean_diff)))
}
}
return(ssp_matrix)
}
# Função para calcular a matriz SSP Residual
SSP_res <- function(data) {
g <- length(unique(data$F1))
b <- length(unique(data$F2))
n <- nrow(data) / (g * b)
p <- ncol(data[, grep("X",colnames(data))])
ssp_matrix <- matrix(0, p, p)
for (i in unique(data$F1)) {
for (j in unique(data$F2)) {
current_subset <- subset(data, F1 == i & F2 == j)
mean_diff <- t(as.matrix(current_subset[, grep("X",colnames(data))] -
matrix(rep(colMeans(current_subset[, grep("X",colnames(data))]), each=nrow(current_subset)), ncol=p)))
ssp_matrix <- ssp_matrix + (mean_diff %*% t(mean_diff))
}
}
return(ssp_matrix)
}
# Agora, vamos calcular cada matriz SSP e reproduzir a tabela
ssp_fat1 <- SSP_fat1(x)
ssp_fat2 <- SSP_fat2(x)
ssp_int <- SSP_int(x)
ssp_res <- SSP_res(x)
ssp_tot <- ssp_fat1 + ssp_fat2 + ssp_int + ssp_res
cat("SSP for Factor 1:\n")
SSP for Factor 1:
[,1] [,2] [,3]
[1,] 1.7405 -1.5045 0.8555
[2,] -1.5045 1.3005 -0.7395
[3,] 0.8555 -0.7395 0.4205
SSP for Factor 2:
[,1] [,2] [,3]
[1,] 0.7605 0.6825 1.9305
[2,] 0.6825 0.6125 1.7325
[3,] 1.9305 1.7325 4.9005
SSP for Interaction:
[,1] [,2] [,3]
[1,] 0.0005 0.0165 0.0445
[2,] 0.0165 0.5445 1.4685
[3,] 0.0445 1.4685 3.9605
SSP for Residual:
X1 X2 X3
X1 1.764 0.020 -3.070
X2 0.020 2.628 -0.552
X3 -3.070 -0.552 64.924
SSP for Total:
X1 X2 X3
X1 4.2655 -0.7855 -0.2395
X2 -0.7855 5.0855 1.9095
X3 -0.2395 1.9095 74.2055
# Teste de razão de verossimilhança por qui-quadrado aproximado
g <- length(unique(x$F1))
b <- length(unique(x$F2))
n <- nrow(x) / (g*b)
p <- 3
df_fat1 <- (g - 1)*p
df_fat2 <- (b - 1)*p
df_int <- (g - 1) * (b - 1)*p
Lambda_star_int <- det(ssp_res) / det(ssp_int + ssp_res)
X2_int <- -(g*b*(n-1) - (p+1-(g-1)*(b-1))/2) * log(Lambda_star_int)
print(Lambda_star_int)
[1] 0.7771058
[1] 3.656593
[1] 7.814728
[1] 0.3010132
Lambda_star_fat1 <- det(ssp_res) / det(ssp_fat1 + ssp_res)
X2_fat1 <- -(g*b*(n-1) - (p+1-(g-1))/2) * log(Lambda_star_fat1)
print(Lambda_star_fat1)
[1] 0.3818584
[1] 13.95923
[1] 7.814728
[1] 0.002961177
Lambda_star_fat2 <- det(ssp_res) / det(ssp_fat2 + ssp_res)
X2_fat2 <- -(g*b*(n-1) - (p+1-(b-1))/2) * log(Lambda_star_fat2)
print(Lambda_star_fat2)
[1] 0.5230349
[1] 9.397553
[1] 7.814728
[1] 0.02444658
# Formatando a saída para se assemelhar à função car::Anova
output <- data.frame(
"Wilks Lambda" = c(Lambda_star_fat1, Lambda_star_fat2, Lambda_star_int),
"Approx X2" = c(X2_fat1, X2_fat2, X2_int),
"df" = c(df_fat1, df_fat2, df_int),
"p" = c(p_fat1, p_fat2, p_int))
rownames(output) <- c("Fator 1", "Fator 2", "Interação")
cat("Teste de razão de verossimilhança por qui-quadrado aproximado\n")
Teste de razão de verossimilhança por qui-quadrado aproximado
Wilks.Lambda Approx.X2 df p
Fator 1 0.382 13.96 3 0.00296
Fator 2 0.523 9.40 3 0.02445
Interação 0.777 3.66 3 0.30101
A análise de perfil refere-se a situações nas quais um conjunto de \(p\) tratamentos (testes, itens etc.) é administrado a dois ou mais grupos de unidades observacionais. Assume-se que as respostas para os diferentes grupos são independentes entre si. Normalmente, poderíamos nos perguntar: os vetores de média da população são os mesmos? Na análise de perfil, a questão da igualdade dos vetores médios é dividida em várias possibilidades específicas.
Figura 6.4 O perfil populacional para p = 4.
Considere as médias populacionais \(\boldsymbol{\mu}_1 = [\mu_{11}\; \mu_{12}\; \mu_{13}\; \mu_{14}]\) representando as respostas médias a quatro tratamentos para o primeiro grupo. Um gráfico dessas médias, conectado por linhas retas, é mostrado na Figura 6.4. Este gráfico de linhas quebradas é o perfil para a população 1. Perfis podem ser construídos para cada população (grupo). Concentraremos em dois grupos.
Seja \(\boldsymbol{\mu}_1 = [\mu_{11}\; \mu_{12}\; \cdots, \;\mu_{1p}]\) e \(\boldsymbol{\mu}_2 = [\mu_{21}\; \mu_{22}\;\cdots\; \mu_{2p}]\) serem as respostas médias para os \(p\) tratamentos para as populações 1 e 2, respectivamente. A hipótese \(H_0: \boldsymbol{\mu}_1 = \boldsymbol{\mu}_2\) implica que os tratamentos têm o mesmo efeito nas duas populações. Em termos dos perfis da população, podemos formular a questão da igualdade de uma forma passo a passo.
Equivalentemente:
\[ H_{01}: \mu_{1i} - \mu_{1(i-1)} = \mu_{2i} - \mu_{2(i-1)}\\ i = 2, 3, \ldots, p \]
É aceitável?
Equivalentemente:
\[ H_{02}: \mu_{1i} = \mu_{2i}\\ i = 1, 2, \ldots, p \]
É aceitável?
A questão “Supondo que os perfis sejam paralelos, os perfis são lineares?” é considerada no Exercício 6.12. A hipótese nula de perfis lineares paralelos pode ser escrita como:
\[ H_{0}: (\mu_{1i} + \mu_{2i}) - \left(\mu_{1(i-1)} + \mu_{2(i-1)}\right)= \left(\mu_{1(i-1)} + \mu_{2(i-1)}\right) - \left(\mu_{1(i-2)} + \mu_{2(i-2)}\right)\\ i = 3,\ldots, p \]
Embora essa hipótese possa ser de interesse em uma situação particular, na prática, a questão de saber se dois perfis paralelos são os mesmos (coincidentes), seja qual for sua natureza, geralmente é de maior interesse.
Equivalentemente:
\[ H_{03}: \mu_{11} = \mu_{12} = \cdots = \mu_{1p} = \mu_{21} = \mu_{22} = \cdots = \mu_{2p} \]
É aceitável?
O teste de paralelismo não exige que os itens intervalares tenha a mesma unidade de medida e que os itens ordinais tenha o mesmo número de níveis. No entanto, os testes de linearidade com ou sem inclinação exigem os itens intervalares tenha a mesma unidade de medida e que os itens ordinais tenha o mesmo número de níveis.
A hipótese nula na etapa 1 pode ser escrita como
\[ H_{01}: \mathbf{C}\boldsymbol{\mu}_1 = \mathbf{C}\boldsymbol{\mu}_2 \]
Em que \(\mathbf{C}\) é a matriz de contraste:
\[ \underset{(p-1) \times p}{\mathbf{C}} = \begin{bmatrix} -1 & 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & -1 & 1 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & -1 & 1 \\ \end{bmatrix} \tag{6-72} \]
Para amostras independentes de tamanhos \(n_1\) e \(n_2\) das duas populações, a hipótese nula pode ser testada construindo as observações transformadas:
\[ \mathbf{C}\mathbf{X}_{1i} \quad i=1,2,\ldots,n_1 \]
e
\[ \mathbf{C}\mathbf{X}_{2j} \quad j=1,2,\ldots,n_2 \]
Estes têm vetores de média de amostra \(\mathbf{C}\bar{\mathbf{X}}_1\) e \(\mathbf{C}\bar{\mathbf{X}}_2\), respectivamente, e matriz de covariância agrupada \(\mathbf{C}\mathbf{S}_{\text{comb}}\mathbf{C}^{\prime}\).
Como os dois conjuntos de observações transformadas têm distribuições \(\mathcal{N}_{p-1}\left(\mathbf{C}\boldsymbol{\mu}_1, \mathbf{C}\boldsymbol{\Sigma}\mathbf{C}^{\prime}\right)\) e \(\mathcal{N}_{p-1}\left(\mathbf{C}\boldsymbol{\mu}_2, \mathbf{C}\boldsymbol{\Sigma}\mathbf{C}^{\prime}\right)\), respectivamente, uma aplicação do Resultado 6.2 fornece um teste para perfis paralelos.
Rejeite \(H_{01}: \mathbf{C}\boldsymbol{\mu}_1 = \mathbf{C}\boldsymbol{\mu}_2\) (perfis paralelos) no nível \(\alpha\) se:
\[ T^2 = \left(\mathbf{C}\left(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\right)\right)^{\prime}\left(\left(\dfrac{1}{{n_1}} + \dfrac{1}{{n_2}}\right)\mathbf{C}\mathbf{s}_{\text{comb}}\mathbf{C}^{\prime}\right)^{-1}\mathbf{C}\left(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\right) >c^2 \tag{6-73} \]
em que
\[ c^2= \dfrac{{(n_1 + n_2 - 2)(p - 1)}}{{n_1 + n_2 - p}}F_{p-1, n_1+n_2-p}(1-\alpha)=T^2_{p-1, n_1+n_2-2}(1-\alpha) \]
Quando os perfis são paralelos, o primeiro está acima do segundo (\(\mu_{1i} > \mu_{2i}\) para todo \(i\)) ou vice-versa. Sob esta condição, os perfis serão coincidentes apenas se as alturas totais \(\mu_{11} + \mu_{12} + \cdots + \mu_{1p} = \mathbf{1}^{\prime}\boldsymbol{\mu}_1\) e \(\mu_{21} + \mu_{22} + \cdots + \mu_{2p} = \mathbf{1}^{\prime}\boldsymbol{\mu}_2\) forem iguais. Portanto, a hipótese nula na etapa 2 pode ser escrita na forma equivalente:
\[ H_{02}: \mathbf{1}^{\prime}\boldsymbol{\mu}_1 = \mathbf{1}^{\prime}\boldsymbol{\mu}_2 \]
Podemos então testar \(H_{02}\) com a usual estatística t com base nas observações univariadas \(\mathbf{1}^{\prime}\mathbf{x}_{1i}\), \(i = 1, 2, \ldots, n_1\) e \(\mathbf{1}^{\prime}\mathbf{x}_{2j}\), \(j = 1, 2, \ldots, n_2\).
Para duas populações normais, rejeite \(H_{02}: \boldsymbol{1}^{\prime}\boldsymbol{\mu}_1 = \boldsymbol{1}^{\prime}\boldsymbol{\mu}_2\) (perfis coincidentes) no nível \(\alpha\) se:
\[ T^2 = \boldsymbol{1}^{\prime}(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)\left(\left(\dfrac{1}{{n_1}} + \dfrac{1}{{n_2}}\right)\boldsymbol{1}^{\prime}\mathbf{s}_{\text{comb}}\boldsymbol{1}\right)\boldsymbol{1}^{\prime}(\bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2)>c^2 \tag{6-74} \]
\[ c^2=F_{1,n_1+n_2-2}(1-\alpha)=T^2_{1, n_1+n_2-2}(1-\alpha) \]
Para perfis coincidentes, \(\mathbf{x}_{11}, \mathbf{x}_{12}, \ldots, \mathbf{x}_{1n_1}\) e \(\mathbf{x}_{21}, \mathbf{x}_{22}, \ldots, \mathbf{x}_{2n_2}\) são todas observações da mesma população normal. O próximo passo é verificar se todas as variáveis têm a mesma média, de modo que o perfil comum seja nivelado.
Quando \(H_{01}\) e \(H_{02}\) são não rejeitáveis, o vetor de média comum é estimado, usando todas as observações \(n_{1} + n_{2}\), por:
\[ \bar{\mathbf{x}} = \dfrac{{n_1 \bar{\mathbf{x}}_1 + n_2 \bar{\mathbf{x}}_2}}{{n_1 + n_2}} \]
Se o perfil comum é nivelado, então \(\mu_1 = \mu_2 = \cdots = \mu_p\), e a hipótese nula na etapa 3 pode ser escrita como:
\[ H_{03}: \mathbf{C}\boldsymbol{\mu} = \mathbf{0} \]
em que \(\mathbf{C}\) é dado por (6-72). Consequentemente, temos o seguinte teste.
Para duas populações normais: Rejeite \(H_{03}: \mathbf{C}\boldsymbol{\mu} = \mathbf{0}\) (perfis nivelados) no nível \(\alpha\) se:
\[ T^2 = (n_1 + n_2)\bar{\mathbf{x}}^{\prime}\mathbf{C}^{\prime}\left(\mathbf{C}\mathbf{s}\mathbf{C}^{\prime}\right)^{-1}\mathbf{C}\bar{\mathbf{x}} > c^2 \tag{6-75} \]
em que \(\mathbf{s}\) é a matriz de covariância amostral baseada em todas as observações \(n_1 + n_2\), e:
\[ c^2 = \dfrac{{(n_1 + n_2 - 2)(p - 1)}}{{n_1 + n_2 - p}}F_{p-1, n_1+n_2-p+1}(1-\alpha)=T^2_{p-1, n_1+n_2-1}(1-\alpha) \]
Como parte de um estudo mais amplo sobre amor e casamento, E. Hatfield, um sociólogo, entrevistou homens e mulheres recém-casados. Eles foram convidados a responder às seguintes perguntas:
Os níveis dos quatro itens Likert de cinco pontos são:
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# profileR::spouse
x <- read.table("JW6Data/T6-14.dat", quote="\"", comment.char="")
alfa <- 0.05
names(x) <- c("X1", "X2", "X3", "X4", "Grupo")
x$Grupo <- factor(x$Grupo)
print.data.frame(x)
X1 X2 X3 X4 Grupo
1 2 3 5 5 Husband
2 5 5 4 4 Husband
3 4 5 5 5 Husband
4 4 3 4 4 Husband
5 3 3 5 5 Husband
6 3 3 4 5 Husband
7 3 4 4 4 Husband
8 4 4 5 5 Husband
9 4 5 5 5 Husband
10 4 4 3 3 Husband
11 4 4 5 5 Husband
12 5 5 4 4 Husband
13 4 4 4 4 Husband
14 4 3 5 5 Husband
15 4 4 5 5 Husband
16 3 3 4 5 Husband
17 4 5 4 4 Husband
18 5 5 5 5 Husband
19 5 5 4 4 Husband
20 4 4 4 4 Husband
21 4 4 4 4 Husband
22 4 4 4 4 Husband
23 3 4 5 5 Husband
24 5 3 5 5 Husband
25 5 5 3 3 Husband
26 3 3 4 4 Husband
27 4 4 4 4 Husband
28 3 3 5 5 Husband
29 4 4 3 3 Husband
30 4 4 5 5 Husband
31 4 4 5 5 Wife
32 4 5 5 5 Wife
33 4 4 5 5 Wife
34 4 5 5 5 Wife
35 4 4 5 5 Wife
36 3 3 4 4 Wife
37 4 3 5 4 Wife
38 3 4 5 5 Wife
39 4 4 5 4 Wife
40 3 4 4 4 Wife
41 4 5 5 5 Wife
42 5 5 5 5 Wife
43 4 4 5 5 Wife
44 4 4 4 4 Wife
45 4 4 5 5 Wife
46 3 4 4 4 Wife
47 5 5 5 5 Wife
48 4 5 4 4 Wife
49 3 4 4 4 Wife
50 5 3 4 4 Wife
51 5 3 4 4 Wife
52 4 5 4 4 Wife
53 2 5 5 5 Wife
54 3 4 5 5 Wife
55 4 3 5 5 Wife
56 4 4 4 4 Wife
57 4 4 5 5 Wife
58 3 4 4 4 Wife
59 4 4 5 4 Wife
60 4 4 5 5 Wife
X1 X2 X3 X4 Grupo
Min. :2.000 Min. :3.000 Min. :3.000 Min. :3.000 Husband:30
1st Qu.:3.000 1st Qu.:4.000 1st Qu.:4.000 1st Qu.:4.000 Wife :30
Median :4.000 Median :4.000 Median :5.000 Median :5.000
Mean :3.867 Mean :4.033 Mean :4.483 Mean :4.467
3rd Qu.:4.000 3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:5.000
Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
$Husband
X1 X2 X3 X4
Min. :2.00 Min. :3.000 Min. :3.000 Min. :3.0
1st Qu.:3.25 1st Qu.:3.000 1st Qu.:4.000 1st Qu.:4.0
Median :4.00 Median :4.000 Median :4.000 Median :4.5
Mean :3.90 Mean :3.967 Mean :4.333 Mean :4.4
3rd Qu.:4.00 3rd Qu.:4.750 3rd Qu.:5.000 3rd Qu.:5.0
Max. :5.00 Max. :5.000 Max. :5.000 Max. :5.0
$Wife
X1 X2 X3 X4
Min. :2.000 Min. :3.00 Min. :4.000 Min. :4.000
1st Qu.:3.250 1st Qu.:4.00 1st Qu.:4.000 1st Qu.:4.000
Median :4.000 Median :4.00 Median :5.000 Median :5.000
Mean :3.833 Mean :4.10 Mean :4.633 Mean :4.533
3rd Qu.:4.000 3rd Qu.:4.75 3rd Qu.:5.000 3rd Qu.:5.000
Max. :5.000 Max. :5.00 Max. :5.000 Max. :5.000
Type II MANOVA Tests: Wilks test statistic
Df test stat approx F num Df den Df Pr(>F)
Grupo 1 0.864 2.16 4 55 0.086 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Grupo 1 0.136 2.16 4 55 0.086 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alfa,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=3)
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-------------------------------------------------------
Grupo | 0.136 | [0.000, 0.274] | medium
Data Summary:
Husband Wife
X1 3.900000 3.833333
X2 3.966667 4.100000
X3 4.333333 4.633333
X4 4.400000 4.533333
Call:
profileR::pbg(data = x[, 1:4], group = x[, 5], original.names = TRUE,
profile.plot = TRUE)
Hypothesis Tests:
$`Ho: Profiles are parallel`
Multivariate.Test Statistic Approx.F num.df den.df p.value
1 Wilks 0.8785726 2.579917 3 56 0.06255945
2 Pillai 0.1214274 2.579917 3 56 0.06255945
3 Hotelling-Lawley 0.1382099 2.579917 3 56 0.06255945
4 Roy 0.1382099 2.579917 3 56 0.06255945
$`Ho: Profiles have equal levels`
Df Sum Sq Mean Sq F value Pr(>F)
group 1 0.234 0.2344 1.533 0.221
Residuals 58 8.869 0.1529
$`Ho: Profiles are flat`
F df1 df2 p-value
1 8.18807 3 56 0.0001310162
Profile Analysis for One Sample with Hotelling's T-Square:
T-Squared F df1 df2
Ho: Ratios of the means over Mu0=1 2391.35697 535.99380 4 26
Ho: All of the ratios are equal to each other 27.30936 8.47532 3 27
p-value
Ho: Ratios of the means over Mu0=1 0.0000000000
Ho: All of the ratios are equal to each other 0.0003953056
Profile Analysis for One Sample with Hotelling's T-Square:
T-Squared F df1 df2
Ho: Ratios of the means over Mu0=1 2939.7205 658.9029 4 26
Ho: All of the ratios are equal to each other 347.6278 107.8845 3 27
p-value
Ho: Ratios of the means over Mu0=1 0.000000e+00
Ho: All of the ratios are equal to each other 3.774758e-15
\[\Diamond\]
Quando o tamanho da amostra é pequeno, uma análise de perfil dependerá da suposição de normalidade. Essa suposição pode ser verificada usando métodos discutidos no Capítulo 4, com as observações originais ou as observações de contraste.
A análise de perfis para várias populações procede de maneira muito semelhante àquela para duas populações. Na verdade, as medidas gerais de comparação são análogas àquelas recém-discutidas. (Veja [13], [18]).
Como mencionamos anteriormente, o termo “medidas repetidas” refere-se a situações em que a mesma característica é observada, em diferentes momentos ou locais, na mesma unidade experimental.
As observações em uma unidade experimental podem corresponder a diferentes tratamentos, como no Exemplo 6.2, em que o tempo entre batimentos cardíacos foi medido sob as combinações de tratamento 2 x 2 aplicadas a cada cão Os tratamentos precisam ser comparados quando as respostas na mesma unidade experimental estão correlacionadas.
Um único tratamento pode ser aplicado a cada unidade experimental e uma única característica observada ao longo de um período de tempo. Por exemplo, poderíamos medir a massa corporal total de um filhote ao nascer e depois uma vez por mês. É a curva traçada por um cão típico que precisa ser modelada. Neste contexto, nos referimos à curva como uma curva de crescimento. Quando algumas unidades experimentais recebem um tratamento e outros outro tratamento, as curvas de crescimento para os tratamentos precisam ser comparadas.
Para ilustrar o modelo de curva de crescimento introduzido por Potthoff e Roy [21], consideramos as medições de cálcio do osso ulna dominante em mulheres idosas. Além de uma leitura inicial, a Tabela 6.5 fornece leituras após um ano, dois anos e três anos para o grupo controle. Leituras obtidas por absorptiometria de fóton do mesmo participante estão correlacionadas, mas aquelas de diferentes participantes devem ser independentes. O modelo assume que a mesma matriz de covariância \(\mathbf{\Sigma}\) se aplica a cada sujeito. Ao contrário das abordagens univariadas, este modelo não requer que as quatro medidas sejam homocedásticas. Um perfil, construído a partir das quatro médias amostrais (\(x_1\), \(x_2\), \(x_3\), \(x_4\)), resume o crescimento que aqui é uma perda de cálcio ao longo do tempo. O padrão de crescimento pode ser adequadamente representado por um polinômio no tempo?
source("summarySEwithin2.R")
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
alfa <- 0.05
xc <- read.table("JW6Data/T6-5.dat", quote="\"", comment.char="")
xc <- data.frame(1:15, xc)
names(xc) <- c("ID", "0", "1", "2", "3")
xc_long <- tidyr::gather(xc,
key="Year",
value="Ca",
-"ID")
xt <- read.table("JW6Data/T6-6.dat", quote="\"", comment.char="")
xt <- data.frame(1:16, xt)
names(xt) <- c("ID", "0", "1", "2", "3")
xt_long <- tidyr::gather(xt,
key="Year",
value="Ca",
-"ID")
Grupo <- rep("c", nrow(xc_long))
xc_lonG <- data.frame(Grupo, xc_long)
xc_lonG$ID <- factor(xc_lonG$ID)
xc_lonG$Grupo <- factor(xc_lonG$Grupo)
xc_lonG$Year <- factor(xc_lonG$Year)
Grupo <- rep("t", nrow(xt_long))
xt_lonG <- data.frame(Grupo, xt_long)
xt_lonG$ID <- factor(xt_lonG$ID)
xt_lonG$Grupo <- factor(xt_lonG$Grupo)
xt_lonG$Year <- factor(xt_lonG$Year)
x_lonG <- rbind(xc_lonG, xt_lonG)
boxplot(Ca ~ Year*Grupo, data=x_lonG)
alfaBonf <- alfa/(length(unique(x_lonG$Year))*
length(unique(x_lonG$Grupo)))
ic <- summarySEwithin2(x_lonG,
measurevar="Ca",
withinvars="Year",
betweenvars="Grupo",
idvar="ID",
na.rm=TRUE,
conf.interval=1-alfaBonf)
print(ic)
Grupo Year N Ca CaNormed sd se ci
1 c 0 15 72.38000 71.97852 4.822324 1.2451186 4.001233
2 c 1 15 73.29333 72.89185 5.004513 1.2921598 4.152402
3 c 2 15 72.47333 72.07185 4.450807 1.1491933 3.692974
4 c 3 15 64.78667 64.38519 3.779052 0.9757470 3.135598
5 t 0 16 69.28750 69.66389 3.278236 0.8195590 2.603825
6 t 1 16 70.65625 71.03264 4.012572 1.0031430 3.187090
7 t 2 16 71.18125 71.55764 4.277150 1.0692875 3.397238
8 t 3 16 64.53125 64.90764 3.743544 0.9358861 2.973408
grf <- ggplot2::ggplot(ic,
ggplot2::aes(x=Year,
y=Ca,
color=Grupo)) +
ggplot2::geom_errorbar(width=.1,
ggplot2::aes(ymin=Ca-ci,
ymax=Ca+ci)) +
ggplot2::geom_point(shape=21,
size=3,
fill="white") +
ggplot2::ylab("Ca") +
ggplot2::ggtitle("Medições de Cálcio na Ulna Dominante\nWithin-subject CI95 Bonferroni") +
ggplot2::theme_bw()
print(grf)
Dados <- MANOVA.RM::EEG
eeg <- tidyr::spread(EEG, feature, resp)
fit <- MANOVA.RM::multRM(cbind(brainrate, complexity) ~ sex * region,
data=eeg,
subject="id",
within="region")
summary(fit)
Call:
cbind(brainrate, complexity) ~ sex * region
A multivariate repeated measures analysis with 1 within-subject factor(s) ( region )and 1 between-subject factor(s).
Descriptive:
sex region n brainrate complexity
1 M central 59 -0.254 -0.302
2 M frontal 59 -0.335 -0.399
3 M temporal 59 -0.294 -0.386
4 W central 101 0.149 0.176
5 W frontal 101 0.195 0.233
6 W temporal 101 0.172 0.226
Wald-Type Statistic (WTS):
Test statistic df p-value
sex "12.45" "2" "0.002"
region "0.192" "4" "0.996"
sex:region "2.79" "4" "0.593"
modified ANOVA-Type Statistic (MATS):
Test statistic
sex 54.203
region 0.048
sex:region 0.703
p-values resampling:
paramBS (WTS) paramBS (MATS)
sex "0.002" "<0.001"
region "0.998" "0.991"
sex:region "0.609" "0.431"
# a multivariate linear model for repeated-measures data
## See ?OBrienKaiser for a description of the data set used in this example.
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
idata
# A tibble: 15 × 2
phase hour
<fct> <ord>
1 pretest 1
2 pretest 2
3 pretest 3
4 pretest 4
5 pretest 5
6 posttest 1
7 posttest 2
8 posttest 3
9 posttest 4
10 posttest 5
11 followup 1
12 followup 2
13 followup 3
14 followup 4
15 followup 5
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5,
post.1, post.2, post.3, post.4, post.5,
fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment*gender,
data=carData::OBrienKaiser)
(av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour))
Type II Repeated Measures MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
(Intercept) 1 0.96954 318.34 1 10 6.532e-09 ***
treatment 2 0.48092 4.63 2 10 0.0376868 *
gender 1 0.20356 2.56 1 10 0.1409735
treatment:gender 2 0.36350 2.86 2 10 0.1044692
phase 1 0.85052 25.61 2 9 0.0001930 ***
treatment:phase 2 0.68518 2.61 4 20 0.0667354 .
gender:phase 1 0.04314 0.20 2 9 0.8199968
treatment:gender:phase 2 0.31060 0.92 4 20 0.4721498
hour 1 0.93468 25.04 4 7 0.0003043 ***
treatment:hour 2 0.30144 0.35 8 16 0.9295212
gender:hour 1 0.29274 0.72 4 7 0.6023742
treatment:gender:hour 2 0.57022 0.80 8 16 0.6131884
phase:hour 1 0.54958 0.46 8 3 0.8324517
treatment:phase:hour 2 0.66367 0.25 16 8 0.9914415
gender:phase:hour 1 0.69505 0.85 8 3 0.6202076
treatment:gender:phase:hour 2 0.79277 0.33 16 8 0.9723693
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
Sum Sq num Df Error SS den Df F value Pr(>F)
(Intercept) 7260.0 1 228.056 10 318.3435 6.532e-09
treatment 211.3 2 228.056 10 4.6323 0.037687
gender 58.3 1 228.056 10 2.5558 0.140974
treatment:gender 130.2 2 228.056 10 2.8555 0.104469
phase 167.5 2 80.278 20 20.8651 1.274e-05
treatment:phase 78.7 4 80.278 20 4.8997 0.006426
gender:phase 1.7 2 80.278 20 0.2078 0.814130
treatment:gender:phase 10.2 4 80.278 20 0.6366 0.642369
hour 106.3 4 62.500 40 17.0067 3.191e-08
treatment:hour 1.2 8 62.500 40 0.0929 0.999257
gender:hour 2.6 4 62.500 40 0.4094 0.800772
treatment:gender:hour 7.8 8 62.500 40 0.6204 0.755484
phase:hour 11.1 8 96.167 80 1.1525 0.338317
treatment:phase:hour 6.3 16 96.167 80 0.3256 0.992814
gender:phase:hour 6.6 8 96.167 80 0.6900 0.699124
treatment:gender:phase:hour 14.2 16 96.167 80 0.7359 0.749562
(Intercept) ***
treatment *
gender
treatment:gender
phase ***
treatment:phase **
gender:phase
treatment:gender:phase
hour ***
treatment:hour
gender:hour
treatment:gender:hour
phase:hour
treatment:phase:hour
gender:phase:hour
treatment:gender:phase:hour
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Mauchly Tests for Sphericity
Test statistic p-value
phase 0.74927 0.27282
treatment:phase 0.74927 0.27282
gender:phase 0.74927 0.27282
treatment:gender:phase 0.74927 0.27282
hour 0.06607 0.00760
treatment:hour 0.06607 0.00760
gender:hour 0.06607 0.00760
treatment:gender:hour 0.06607 0.00760
phase:hour 0.00478 0.44939
treatment:phase:hour 0.00478 0.44939
gender:phase:hour 0.00478 0.44939
treatment:gender:phase:hour 0.00478 0.44939
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
phase 0.79953 7.323e-05 ***
treatment:phase 0.79953 0.01223 *
gender:phase 0.79953 0.76616
treatment:gender:phase 0.79953 0.61162
hour 0.46028 8.741e-05 ***
treatment:hour 0.46028 0.97879
gender:hour 0.46028 0.65346
treatment:gender:hour 0.46028 0.64136
phase:hour 0.44950 0.34573
treatment:phase:hour 0.44950 0.94019
gender:phase:hour 0.44950 0.58903
treatment:gender:phase:hour 0.44950 0.64634
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
HF eps Pr(>F[HF])
phase 0.9278594 2.387543e-05
treatment:phase 0.9278594 8.089765e-03
gender:phase 0.9278594 7.984495e-01
treatment:gender:phase 0.9278594 6.319975e-01
hour 0.5592802 2.014357e-05
treatment:hour 0.5592802 9.887716e-01
gender:hour 0.5592802 6.911521e-01
treatment:gender:hour 0.5592802 6.692976e-01
phase:hour 0.7330608 3.440460e-01
treatment:phase:hour 0.7330608 9.804731e-01
gender:phase:hour 0.7330608 6.552382e-01
treatment:gender:phase:hour 0.7330608 7.080122e-01
## A "doubly multivariate" design with two distinct repeated-measures variables
## (example courtesy of Michael Friendly)
## See ?WeightLoss for a description of the dataset.
imatrix <- matrix(c(
1,0,-1, 1, 0, 0,
1,0, 0,-2, 0, 0,
1,0, 1, 1, 0, 0,
0,1, 0, 0,-1, 1,
0,1, 0, 0, 0,-2,
0,1, 0, 0, 1, 1), 6, 6, byrow=TRUE)
colnames(imatrix) <- c("WL", "SE", "WL.L", "WL.Q", "SE.L", "SE.Q")
rownames(imatrix) <- colnames(WeightLoss)[-1]
(imatrix <- list(measure=imatrix[,1:2], month=imatrix[,3:6]))
$measure
WL SE
wl1 1 0
wl2 1 0
wl3 1 0
se1 0 1
se2 0 1
se3 0 1
$month
WL.L WL.Q SE.L SE.Q
wl1 -1 1 0 0
wl2 0 -2 0 0
wl3 1 1 0 0
se1 0 0 -1 1
se2 0 0 0 -2
se3 0 0 1 1
contrasts(WeightLoss$group) <- matrix(c(-2,1,1, 0,-1,1), ncol=2)
(wl.mod<-lm(cbind(wl1, wl2, wl3, se1, se2, se3)~group, data=WeightLoss))
Call:
lm(formula = cbind(wl1, wl2, wl3, se1, se2, se3) ~ group, data = WeightLoss)
Coefficients:
wl1 wl2 wl3 se1 se2 se3
(Intercept) 5.34444 4.45000 2.17778 14.92778 13.79444 16.28333
group1 0.42222 0.55833 0.04722 0.08889 -0.26944 0.60000
group2 0.43333 1.09167 -0.02500 0.18333 -0.22500 0.71667
Type II Repeated Measures MANOVA Tests: Roy test statistic
Df test stat approx F num Df den Df Pr(>F)
measure 1 86.203 1293.04 2 30 < 2.2e-16 ***
group:measure 2 0.356 5.52 2 31 0.008906 **
month 1 9.407 65.85 4 28 7.807e-14 ***
group:month 2 1.772 12.84 4 29 3.909e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
6.1, 6.5, 6.8, 6.12, 6.13, 6.14, 6.15, 6.18, 6.20, 6.22, 6.28, 6.39, 6.41
Soch, J et al. (2020) The Book of Statistical Proofs: https://statproofbook.github.io/
Batschelet, E (1978) Introdução à matemática para biocientistas. Tradução da 2ª ed. São Paulo: EDUSP e Rio de Janeiro: Interciência.
Batschelet, E (1979) Introduction to mathematics for life scientists. 3rd ed. NY: Springer.
Bickel, PJ & Doksum, KA (2001) Mathematical Statistics: Basic Ideas and Selected Topics, Volume I. USA: CRC.
Bilodeau, M & Brenner, D (1999) Theory of Multivariate Statistics. USA: Springer. T^2 de Hotelling
Denis, DJ (2021) Applied Univariate, Bivariate, and Multivariate Statistics Using Python: A Beginners Guide to Advanced Data Analysis. NJ: Wiley.
Genz, A (1992) Numerical Computation of Multivariate Normal Probabilities. Journal of Computational and Graphical Statistics 1(2): 141–149. https://doi.org/10.2307/1390838
Hardle, W & Hlavka, Z (2007) Multivariate Statistics - Exercises and Solutions. USA: Springer.
KASSAMBARA, A (2017) Practical guide to cluster analysis in R: Unsupervised machine learning. STHDA: http://www.sthda.com.
Kollo. T & von Rosen, D (2005) Advanced Multivariate Statistics with Matrices. USA: Springer.
Kutner, MH et al. (2005) Applied linear statistical model. 5th ed. NY: McGraw-Hill/Irwin.
MAIR, P (2018) Modern psychometrics with R. USA: Springer.
Manly, BFJ & Alberto, JAN (2017) Multivariate Statistical
Methods: A Primer using R. 4th ed. USA: CRC.
Matloff, N (2020) Probability and Statistics for Data Science: Math + R Data_. USA: CRC.
Mirman, D (2014) Growth Curve Analysis and Visualization Using R. USA: CRC.
Pessoa, DGC (2008) Exemplos do livro de Johnson e Wichern.
Rawlings, JO et al. (1998) Applied Regression Analysis: A Research Tool. USA: Springer. T^2 de Hotelling
Schumacker, RE (2016) Using R With Multivariate Statistics. USA: SAGE.
Siqueira, JO (2012) Fundamentos de métodos quantitativos: aplicados em Administração, Economia, Contabilidade e Atuária usando WolframAlpha e SCILAB. São Paulo: Saraiva. Soluções dos exercícios em https://www.researchgate.net/publication/326533772_Fundamentos_de_metodos_quantitativos_-_Solucoes.
Tay, A (2018) OLS using Matrix Algebra. http://www.mysmu.edu/faculty/anthonytay/MFE/OLS_using_Matrix_Algebra.pdf
Teichroew, D (1964) An introduction to management science: deterministic models. NJ: Wiley.
Timm, NH (2002) Applied Multivariate Analysis. USA: Springer. T^2 de Hotelling
Zelterman, D (2015) Applied Multivariate Statistics with R. USA: Springer.
Nguyen, M (2020) A Guide on Data Analysis. Bookdown. https://bookdown.org/mike/data_analysis/
Chaves Neto, A. CE-704 ANÁLISE MULTIVARIADA APLADA À PESQUISA. https://docs.ufpr.br/~soniaisoldi/ce090/CE076AM_2010.pdf
ME 731 - Métodos em Análise Multivariada em R. https://www.ime.unicamp.br/~cnaber/Material_AM_2S_2020.htm
Nogueira, FE (2007) Modelos de regressão multivariada. Dissertação (Mestrado). https://teses.usp.br/teses/disponiveis/45/45133/tde-25062007-163150/publico/dissertacao_4.pdf
Análise Multivariada. https://www.professores.uff.br/samuelcampos/analise-multivariada/
Severn, K (2023) Multivariate Statistics. https://rich-d-wilkinson.github.io/MATH3030/
Powell, J (2019) Multivariate statistics in R. https://www.hiercourse.com/multivariate
Plotting PCA/clustering results using ggplot2 and ggfortify: https://rpubs.com/sinhrks/plot_pca
ggfortify : Extension to ggplot2 to handle some popular packages - R software and data visualization: http://www.sthda.com/english/wiki/ggfortify-extension-to-ggplot2-to-handle-some-popular-packages-r-software-and-data-visualization
A Little Book of Python for Multivariate Analysis: https://python-for-multivariate-analysis.readthedocs.io/index.html