class: split-70 with-border hide-slide-number bg-brand-red background-image: url("images/USydLogo-black.svg") background-size: 200px background-position: 2% 90% .column.white[.content[ <br><br><br> # One-way ANOVA Part II ## .black[STAT3022 Applied Linear Models Lecture 15] <br><br><br> ### .black[2020/02/17] ]] .column.bg-brand-charcoal[.content.white[ ## Today 1. Partitioning ANOVA sums of squares 1. The ANOVA table 1. Treatment sum of squares and its distribution 1. Independent contrasts ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # One-way ANOVA structure ]] .row[.split-90[ .column[.content[ For the one-way ANOVA model `$$\text{M}_1: Y_{ij} = \mu + \alpha_i + \epsilon_{ij}~~~~(i=1,\ldots,t,~~j=1,\ldots,n_i)$$` The ANOVA table seen in the previous lecture with `anova()` is: Source of Variation | Df | Sum Sq | Mean Sq | F value | P value ------------- | ------------- | ------------- | ------------- | ------------- | ------------- Treatment | `\(t-1\)` | `\(\text{RSS}_{H_0}-\text{RSS}_{H_1}\)` | `\(\text{MTSS}=\frac{\text{RSS}_{H_0}-\text{RSS}_{H_1}}{t-1}\)` | `\(f_0=\text{MTSS}/\text{MRSS}\)` | `\(\Pr(f_0>F_{t-1,n-t})\)` Residual | `\(n-t\)` | `\(\text{RSS}_{H_1}\)` | `\(\text{MRSS}=\frac{\text{RSS}_{H_1}}{n-t}\)` | | Total | `\(n-1\)` | `\(\text{RSS}_{H_0}\)` | | | * `\(\text{TSS}=\text{BSS}=\text{RSS}_{H_0}-\text{RSS}_{H_1}=\sum_{i=1}^t\sum_{j=1}^{n_i}(\bar{Y}_{i\bullet} - \bar{Y}_{\bullet\bullet})^2=\sum_{i=1}^tn_i(\bar{Y}_{i\bullet} - \bar{Y}_{\bullet\bullet})^2=\sum_{i=1}^tn_i\bar{Y}_{i\bullet}^2 - n\bar{Y}_{\bullet\bullet}^2\)` * `\(\text{WSS} = \text{RSS}_{H_1}=\sum_{i=1}^t\sum_{j=1}^{n_i}(Y_{ij} - \bar{Y}_{i\bullet})^2=\sum_{i=1}^t\sum_{j=1}^{n_i}Y_{ij}^2 - \sum_{i=1}^tn_i\bar{Y}_{i\bullet}^2\)` * `\(\text{TotSS}=\text{RSS}_{H_0} = \sum_{i=1}^t\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y}_{\bullet\bullet})^2=\sum_{i=1}^t\sum_{j=1}^{n_i}Y_{ij}^2 - n\bar{Y}_{\bullet\bullet}^2\)` * Note also `\(\text{TotSS} = \text{WSS} + \text{BSS} = \text{RSS}_{H_1} + \text{TSS}\)`. * Additionally `\(\text{df}_{\text{Total}} =\text{df}_{\text{Trt}} + \text{df}_{\text{Res}}\)`. ]] .column.bg-brand-gray[.content[ ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # One-way ANOVA structure ]] .row[.split-50[ .column[.content[ ```r caffeine <- read.table(file="data/caffeine.txt", header=TRUE) anova(lm(Taps ~ Dose, data=caffeine)) ``` ``` Analysis of Variance Table Response: Taps Df Sum Sq Mean Sq F value Pr(>F) Dose 2 61.4 30.7000 6.1812 0.006163 ** Residuals 27 134.1 4.9667 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r y=caffeine$Taps; x=caffeine$Dose; m=mean(y) ni=as.vector(table(x)); gm=tapply(y,x,mean) mi=c(rep(gm[3],ni[3]),rep(gm[2],ni[2]),rep(gm[1],ni[1])) ToSS=sum((y-m)^2); RSS=sum((y-mi)^2); TSS=sum((mi-m)^2) c(TSS,RSS,ToSS) ``` ``` [1] 61.4 134.1 195.5 ``` ]] .column.bg-brand-gray[.content[ ```r caffeine %>% #tidyverse; means by treatment group_by(Dose) %>% summarise(TapsSum=sum(Taps), TapsSum2=sum(Taps^2)) ``` ``` # A tibble: 3 x 3 Dose TapsSum TapsSum2 <fct> <int> <dbl> 1 High 2483 616573 2 Low 2464 607168 3 Zero 2448 599322 ``` ```r # CHALLENGE; RSS (yij-barYi)^2 & TSS (barYi-barY)^2 caffeine %>% group_by(Dose) %>% summarise(RSSi=sum((Taps - mean(Taps))^2), TSSi=n() * sum((mean(Taps) - mean(caffeine$Taps))^2)) %>% summarise(RSS=sum(RSSi), TSS=sum(TSSi)) ``` ``` # A tibble: 1 x 2 RSS TSS <dbl> <dbl> 1 134. 61.4 ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # One-way ANOVA in matrix notation ]] .row[.split-40[ .column[.content[ `$$\boldsymbol{Y} = \mathbf{X}\boldsymbol{\beta}+ \boldsymbol{\epsilon}; \quad \boldsymbol{\epsilon}\sim N(\boldsymbol{0}, \sigma^2\mathbf{I}_n)$$` where * `\(\boldsymbol{Y}=(Y_{ij})\)` is a vector of observations ordered by treatment group; * `\(\boldsymbol{\beta}=(\alpha_1, ..., \alpha_t)^\top\)` - note we apply the contraint `\(\mu=0\)` here; and * `\(\boldsymbol{\epsilon}=(\epsilon_{11}, \epsilon_{12}, ..., \epsilon_{t,n_t})\)`. Let `\(\mathbf{P}_{V_T}=\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\)` (the hat matrix). Recall * `\(\mathbf{P}_{V_T}\boldsymbol{Y} = \hat{\boldsymbol{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}}\)`, the fitted values and * `\((\mathbf{I}_n - \mathbf{P}_{V_T})\boldsymbol{Y} = \hat{\boldsymbol{\epsilon}} = \boldsymbol{R}\)`, the residuals. ]] .column.bg-brand-gray[.content[ ```r X <- model.matrix( ~ Dose - 1, data=caffeine); head(X,3) ``` ``` DoseHigh DoseLow DoseZero 1 0 0 1 2 0 0 1 3 0 0 1 ``` `\(\mathbf{P}_{V_0}=\boldsymbol{1}_n(\boldsymbol{1}_n^\top\boldsymbol{1}_n)^{-1}\boldsymbol{1}_n^\top\)` is a `\(n\times n\)` matrix of `\(1/n\)` under `\(H_0\)` with only parameter `\(\mu\)` and `\(\boldsymbol{X}=\boldsymbol{1}_n\)`. Under `\(H_0\)`, `\(\hat{\boldsymbol{Y}}=\mathbf{P}_{V_0}\boldsymbol{Y}=\bar Y_{\bullet \bullet}\)`. ```r n <- nrow(caffeine); y <- caffeine$Taps; mean(y) ``` ``` [1] 246.5 ``` ```r u0 <- rep(1, n) # first column vector of ones in X for H0 P.V0 <- u0 %*% solve(t(u0) %*% u0) %*% t(u0) as.vector(P.V0 %*% y) # fitted under H0 as mean(Y) ``` ``` [1] 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 [13] 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 246.5 [25] 246.5 246.5 246.5 246.5 246.5 246.5 ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Projection Matrices for One-Way ANOVA ]] .row[.split-50[ .column[.content[ `\begin{eqnarray*} \boldsymbol{Y} &=&\underbrace{\mathbf{P}_{V_T}\boldsymbol{Y}}_{\text{fitted values}}+\underbrace{(\mathbf{I} - \mathbf{P}_{V_T})\boldsymbol{Y}}_{\text{residual}} \\ &&\\ \underbrace{(\mathbf{I}_n \color{blue}{- \mathbf{P}_{V_0}})}_{\text{Total}}\boldsymbol{Y}&=& \underbrace{(\mathbf{P}_{V_T} \color{blue}{- \mathbf{P}_{V_0}})}_{\text{Treatment}}\boldsymbol{Y} + \underbrace{(\mathbf{I} - \mathbf{P}_{V_T})}_{\text{Residual}}\boldsymbol{Y} \end{eqnarray*}` Let `\(\mathbf{P}_{W_T} = \mathbf{P}_{V_T} - \mathbf{P}_{V_0}\)` for treatment. Now note that `$$(\mathbf{P}_{W_T}\boldsymbol{Y})^\top\mathbf{P}_{W_T}\boldsymbol{Y} = \boldsymbol{Y}^\top\mathbf{P}_{W_T}^\top\mathbf{P}_{W_T}\boldsymbol{Y}=\boldsymbol{Y}^\top\mathbf{P}_{W_T}\boldsymbol{Y}$$` since the projection matrix `\(\mathbf{P}_{W_T}\)` is idempotent, ie `\(\mathbf{P}_{W_T}\mathbf{P}_{W_T}=\mathbf{P}_{W_T}\)`. Same for `\((\mathbf{I}_n \color{blue}{- \mathbf{P}_{V_0}})\)` and `\((\mathbf{I} - \mathbf{P}_{V_T})\)`. ```r P.VT <- X %*% solve(t(X) %*% X) %*% t(X) P.WT <- P.VT - P.V0 P.Res <- diag(n) - P.VT P.Tot <- diag(n) - P.V0 ``` Note: `diag(n)` = `\(\boldsymbol{I}_n\)`. ]] .column.bg-brand-gray[.content[ <img src="images/Ftest2.png" width="100%" height="100%"> This tests full model against null model in `\(H_0\)`. Hence with respect to regression parameters, the null model is empty. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Checking ranks of projection matrices ]] .row[.split-50[ .column[.content[ ```r data.frame(Source=c("Treatment", "Residual", "Total"), df=c(qr(P.WT)$rank, qr(P.Res)$rank, qr(P.Tot)$rank), SS=c(t(y) %*% P.WT %*% y, t(y) %*% P.Res %*% y, t(y) %*% P.Tot %*% y)) ``` ``` Source df SS 1 Treatment 2 61.4 2 Residual 27 134.1 3 Total 29 195.5 ``` * The function `qr` returns a list which contains `rank` of the matrix. * What do you notice? ]] .column.bg-brand-gray[.content[ ### Check ranks of `\(n\)` by `\(n\)` matrix `\(\mathbf{P}_{V_0}\)` ```r round(P.V0[1:12,1:12],2) #all entities 1/30 ``` ``` [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [2,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [3,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [4,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [5,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [6,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [7,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [8,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [9,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [10,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [11,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 [12,] 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 ``` ```r qr(P.V0)$rank ``` ``` [1] 1 ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Check ranks of `\(n\)` by `\(n\)` matrices `\(\mathbf{P}_{V_T}\)` and `\(\mathbf{P}_{W_T}\)` ]] .row[.split-50[ .column[.content[ ```r round(P.VT[1:12,1:12],2) #3 block diag with 1/10 ``` ``` 1 2 3 4 5 6 7 8 9 10 11 12 1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 3 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 4 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 5 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 6 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 7 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 8 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 10 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.0 0.0 11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.1 12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.1 ``` ```r qr(P.VT)$rank ``` ``` [1] 3 ``` ]] .column.bg-brand-gray[.content[ ```r round(P.WT[1:12,1:12],2) #0.07=1/10-1/30; -0.03=0-1/30 ``` ``` 1 2 3 4 5 6 7 8 9 10 11 12 1 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 2 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 3 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 4 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 5 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 6 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 7 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 8 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 9 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 10 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 -0.03 -0.03 11 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.07 0.07 12 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.07 0.07 ``` Since there are 3 different rows or columns but the sum of all rows or columns is `\(\boldsymbol{0}_{30}\)`, the df is `\(3-1\)`. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Check ranks of `\(\boldsymbol{I}_{n}-\mathbf{P}_{V_T}\)` and `\(\boldsymbol{I}_{n}-\mathbf{P}_{V_0}\)` ]] .row[.split-50[ .column[.content[ ```r round(P.Res[1:12,1:12],2) #1-0.1, or 0-0.1 ``` ``` 1 2 3 4 5 6 7 8 9 10 11 12 1 0.9 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.0 0.0 2 -0.1 0.9 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.0 0.0 3 -0.1 -0.1 0.9 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.0 0.0 4 -0.1 -0.1 -0.1 0.9 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.0 0.0 5 -0.1 -0.1 -0.1 -0.1 0.9 -0.1 -0.1 -0.1 -0.1 -0.1 0.0 0.0 6 -0.1 -0.1 -0.1 -0.1 -0.1 0.9 -0.1 -0.1 -0.1 -0.1 0.0 0.0 7 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.9 -0.1 -0.1 -0.1 0.0 0.0 8 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.9 -0.1 -0.1 0.0 0.0 9 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.9 -0.1 0.0 0.0 10 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.9 0.0 0.0 11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9 -0.1 12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.1 0.9 ``` Since there are 30 different rows or columns but the sum of all rows or columns in each block is `\(\boldsymbol{0}_{30}\)`, the df is `\(30-3\)`. ]] .column.bg-brand-gray[.content[ ```r round(P.Tot[1:12,1:12],2) #0.97=1-1/30; -0.03=0-1/30 ``` ``` [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 [2,] -0.03 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 [3,] -0.03 -0.03 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 [4,] -0.03 -0.03 -0.03 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 [5,] -0.03 -0.03 -0.03 -0.03 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 [6,] -0.03 -0.03 -0.03 -0.03 -0.03 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 [7,] -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.97 -0.03 -0.03 -0.03 -0.03 -0.03 [8,] -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.97 -0.03 -0.03 -0.03 -0.03 [9,] -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.97 -0.03 -0.03 -0.03 [10,] -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.97 -0.03 -0.03 [11,] -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.97 -0.03 [12,] -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 -0.03 0.97 ``` Since there are 30 different rows or columns but the sum of all rows or columns is `\(\boldsymbol{0}_{30}\)`, the df is `\(30-1\)`. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Distribution of treatment sum of squares ]] .row[.split-70[ .column.bg-brand-gray[.content[ Given a constrast `\(c = \sum_{i=1}^t c_i \alpha_i\)`. Using basic properties of `\(\text{E}\)`, `\(\text{Var}\)` and independent normal random variables we can show that `$$\hat c=\sum_{i=1}^t c_i \bar{Y}_{i\bullet} \sim N\Big( c\, ,\, \sigma^2 \sum_{i=1}^t (c_i^2/n_i) \Big)$$` Thus if all `\(\alpha_i\)` are equal, using `\(\sum_{i=1}^t c_i =0\)` `\(\Rightarrow\)` `\(c=\sum_{i=1}^t c_i \alpha_i = 0\)`. So under `\(H_0\)`, `$$\frac{\sum_{i=1}^t c_i \bar{Y}_{i\bullet}}{\sqrt{\sum_{i=1}^t (c_i^2/n_i)}} \sim N\Big( 0\, ,\, \sigma^2 \Big) \;\Rightarrow\;\frac{\left(\sum_{i=1}^t c_i \bar{Y}_{i\bullet}\right)^2}{\sum_{i=1}^t (c_i^2/n_i)} \sim \sigma^2 \chi^2_1.$$` * Each contrast `\(c\)` (say `\(\alpha_1=\alpha_i, \ i \ne 1\)`) determines a 1 df component of the TSS. * The `\(\text{TSS}\sim \sigma^2 \chi^2_{t-1}\)` if `\(H_0\)` is true. * There exist `\((t-1)\)` contrasts to decompose `\(\text{TSS}\)` into independent 1 df components. ]] .column[.content[ .blockquote[ Every contrast `\(\tilde{c}\)` determines a 1 df component of the TSS ] .blockquote[ `$$\frac{\sum_{i=1}^t c_i \bar{Y}_{i\bullet}}{\sqrt{\sum_{i=1}^t (c_i^2/n_i)}} \sim N\Big( 0\, ,\, \sigma^2 \Big)$$` ] .blockquote[ `$$\frac{\left(\sum_{i=1}^t c_i \bar{Y}_{i\bullet}\right)^2}{\sum_{i=1}^t (c_i^2/n_i)} \sim \sigma^2 \chi^2_1.$$` ] ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Good choice of contrast ]] .row[.split-30[ .column[.content[ Choose component-weights so that the contrasts are suitable and meaningful comparisons. * Possibility to determine the source of the differences between the treatments if the null hypothesis is rejected. * Hopefully 1 or 2 components account for most of the TSS. ]] .column.bg-brand-gray[.content[ ### Independent contrasts Given two contrasts * `\(c = \sum_{i=1}^t c_i \alpha_i\)` such that `\(\sum_{i=1}^t c_i = 0,\)` * `\(d = \sum_{i=1}^t d_i \alpha_i\)` such that `\(\sum_{i=1}^t d_i =0,\)` then the corresponding estimates `\(\hat c=\sum_{i=1}^t c_i \bar{Y}_{i\bullet}\;\;\text{ and }\;\; \hat d=\sum_{i=1}^t d_i \bar{Y}_{i\bullet}\)` are independent if `\(\sum_{i=1}^t \frac{c_i d_i}{ n_i} =0\)`. This is true because, `\(\text{Cov} \Big(\sum_{i=1}^t c_i \bar{Y}_{i\bullet}, \sum\limits_{i=1}^t d_i \bar{Y}_{i\bullet}\Big) = \sum\limits_{i=1}^t\sum\limits_{j=1}^t c_i d_j \text{Cov}(\bar{Y}_{i\bullet},\bar{Y}_{j\bullet}) = \sigma^2 \sum\limits_{i=1}^t \frac{c_i d_i}{n_i}\)` as `$$\text{Cov}(\bar{Y}_{i\bullet},\bar{Y}_{j\bullet}) = \begin{cases} 0, & i\neq j, \\ \frac{\sigma^2}{n_i}, & i = j. \end{cases}$$` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Example - Fertilizers ]] .row[.split-50[ .column[.content[ * In a study to compare 4 fertilizers A, B, C and D an experiment was designed and several control areas "0" were incorporated. * Each treatment was applied to 6 plots and the results were: <img src="lecture15_2020JC_files/figure-html/unnamed-chunk-14-1.svg" style="display: block; margin: auto;" /> ]] .column.bg-brand-gray[.content[ Here we have `\(n_i \equiv 6\)` and `\(Y_{\bullet\bullet} = 300\)`. <img src="lecture15_2020JC_files/figure-html/unnamed-chunk-15-1.svg" style="display: block; margin: auto;" /> ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Fertilizers ]] .row[.split-60[ .column[.content[ ```r dat %>% skimr::skim() ``` ``` Skim summary statistics n obs: 30 n variables: 3 -- Variable type:factor ---------------------------------------------------------------------------------------------------- variable missing n n_unique ordered Replicate 0 30 6 FALSE Trt 0 30 5 FALSE -- Variable type:numeric --------------------------------------------------------------------------------------------------- variable missing complete n mean sd p0 p25 p50 p75 p100 hist Y 0 30 30 10 4.81 2 7 9.5 12 22 <U+2583><U+2586><U+2586><U+2587><U+2583><U+2581><U+2581><U+2582> ``` ]] .column.bg-brand-gray[.content[ ```r tapply(dat$Y,dat$Trt,summary) ``` ``` $`0` Min. 1st Qu. Median Mean 3rd Qu. Max. 3.000 4.250 5.000 5.833 6.500 11.000 $A Min. 1st Qu. Median Mean 3rd Qu. Max. 8.00 10.50 12.00 11.50 12.75 14.00 $B Min. 1st Qu. Median Mean 3rd Qu. Max. 8.00 11.00 16.50 15.50 19.75 22.00 $C Min. 1st Qu. Median Mean 3rd Qu. Max. 2.00 5.50 7.50 7.00 8.75 11.00 $D Min. 1st Qu. Median Mean 3rd Qu. Max. 7.00 8.00 10.00 10.17 12.00 14.00 ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Fertilizers: Omnibus F-test ]] .row[.split-50[ .column[.content[ Consider the model `$$Y_{ij} = \mu + \alpha_i +\epsilon_{ij},\;\; \epsilon_{ij}\sim NID(0,\sigma^2)$$` where `\(Y_{ij}\)` denotes the `\(j\)`th observation on the `\(i\)`th treatment, `\(i=1,\ldots,5,\;\; j=1,\ldots,6\)`. .brand-red[Test the claim that all five treatment produce the same mean yield.] ```r anova(lm(Y ~ Trt, dat)) ``` ``` Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) Trt 4 353.33 88.333 6.9299 0.0006743 *** Residuals 25 318.67 12.747 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ]] .column.bg-brand-gray[.content[ The observed test statistic is `\(f_0=6.93\)`. Since `\(F^{-1}_{4,25}(0.999) = 6.49\)`, the p-value `\(p<0.001\)` and so we reject `\(H_0: \alpha_1 = \ldots = \alpha_5 = 0\)`. ### Same as the general F-test ```r Mf <- lm(Y ~ Trt, dat) M0 <- lm(Y ~ 1, dat) *anova(M0, Mf) ``` ``` Analysis of Variance Table Model 1: Y ~ 1 Model 2: Y ~ Trt Res.Df RSS Df Sum of Sq F Pr(>F) 1 29 672.00 2 25 318.67 4 353.33 6.9299 0.0006743 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` `\(672=RSS_{H_0}=TSS_{H_1}+RSS_{H_1}=353.33+318.67\)` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Fertilizers: Partition TSS into independent contrasts ]] .row[.split-50[ .column[.content[ .brand-red[Partition the TSS into 4 independent and one degree of freedom components.] **Solution**: There are several ways of doing this using contrasts. We want to select contrasts that are meaningful. Since all `\(n_i\equiv 6\)` we need to check that `\(\sum_{i=1}^5c_i d_i = 0,\ldots,\sum_{i=1}^5e_i f_i = 0\)`. E.g. Treatments | 0 | A | B | C | D --- | --- | --- | --- | --- A vs D | 0 | 1 | 0 | 0 | -1 | `\(c_i\)` C vs 0 | -1 | 0 | 0 | 1 | 0 | `\(d_i\)` (A,D) vs (C,0) | -1 | 1 | 0 | -1 | 1 | `\(e_i\)` (A,D,C,0) vs B | 1 | 1 | -4 | 1 | 1 | `\(f_i\)` `\(6\bar{Y}_{i\bullet}\)` | 35 | 69 | 93 | 42 | 61 ]] .column.bg-brand-gray[.content[ * `\(\text{Recall: }\frac{\left(\sum_{i=1}^t c_i \bar{Y}_{i\bullet}\right)^2}{\sum_{i=1}^t (c_i^2/n_i)} \sim \sigma^2 \chi^2_1\)` `\(\displaystyle \text{For} \ c: \ \frac{\left( \frac{69}{6} -\frac{61}{6} \right)^2}{\frac{1^2}{6} + \frac{(-1)^2}{6}} = 5.333.\)` `\(\displaystyle \text{For} \ f: \ \frac{\left( \frac{35}{6}+\frac{69}{6} +\frac{42}{6}+\frac{61}{6}-4\frac{93}{6} \right)^2}{4\frac{1^2}{6} + \frac{(-4)^2}{6}} = 226.88.\)` * Similarly for `\(d\)` and `\(e\)`, so TSS `\(= 5.33 + 4.08 + 117.04 + 226.88 = 353.33\)`. ```r m=tapply(dat$Y,dat$Trt,mean); TSSi=6*((m-mean(Y))^2); round(c(TSSi,sum(TSSi)),2) ``` ``` 0 A B C D 104.17 13.50 181.50 54.00 0.17 353.33 ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # TSS Partition in example and basic pairwise ]] .row[.split-60[ .column[.content[ ### TSS Partition in example ```r out <- anova(lm(Y ~ I((Trt=="A")-(Trt=="D"))+ I((Trt=="C")-(Trt=="0"))+ I((Trt=="A")-(Trt=="0")+(Trt=="D")-(Trt=="C"))+ I((Trt=="A")+(Trt=="0")+(Trt=="C")+(Trt=="D")-4*(Trt=="B")), data=dat)) rownames(out) <- NULL; out ``` ``` Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) [1,] 1 5.33 5.333 0.4184 0.5236266 [2,] 1 4.08 4.083 0.3203 0.5764460 [3,] 1 117.04 117.042 9.1821 0.0056148 ** [4,] 1 226.87 226.875 17.7988 0.0002819 *** [5,] 25 318.67 12.747 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ]] .column.bg-brand-gray[.content[ ### TSS Partition in basic pairwise ```r out <- anova(lm(Y ~ I((Trt=="0")-(Trt=="A"))+ I((Trt=="0")-(Trt=="B"))+ I((Trt=="0")-(Trt=="C"))+ I((Trt=="0")-(Trt=="D")), data=dat)) rownames(out) <- NULL; out ``` ``` Analysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) [1,] 1 96.33 96.333 7.5575 0.0109371 * [2,] 1 186.78 186.778 14.6531 0.0007697 *** [3,] 1 70.01 70.014 5.4927 0.0273489 * [4,] 1 0.21 0.208 0.0163 0.8992951 [5,] 25 318.67 12.747 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` There are different ways of partition but TSS is the same. ]] ]]