-
Model 4 - The Effects Model
- \(Y_{ij}=\mu+\tau_i+\epsilon_{ij}\)
-
where \(\tau_i\) are the deviations of each factor level mean from the overall mean so that \(\sum_{i=1}^T\tau_i=0\).
In the effects model that we discussed in lesson 3, the treatment means were not estimated, but instead, the deviations of treatment means from the overall mean were estimated (the \(\tau_i\)'s). The model must include the overall mean, thus we need to include \(\mu_., \tau_1 , ... , \tau_T\) as elements in the parameter vector, \(\boldsymbol{\beta}\) of the GLM model. Note that, in the case of equal replication at each factor level, the deviations satisfy the following constraint:
\(\sum_{i=1}^{T} \tau_i =0\)
This implies that one of the \(\tau_i\) parameters is not needed since it can be expressed in terms of the other T - 1 parameters and need not be included in the \(\boldsymbol{\beta}\) parameter vector. We shall drop the parameter \(\tau_T\) from the regression equation as it can be expressed in terms of the other T - 1 parameters \(\tau_i\) as follows:
\(\tau_T\ = -\tau_1 - \tau_2 - \cdots \tau_{T-1}\)
Thus, the \(\boldsymbol{\beta}\) vector of the GLM is a \(T \times 1\) vector containing only the parameters \(\mu, \tau_1 , ... , \tau_{T-1}\) for the linear model.
To illustrate how a linear model is developed with this approach, consider again a single-factor study with T = 3 factor levels with \(n_1=n_2=n_3=2\). The \(\mathbf{Y}\), \(\mathbf{X}\), \(\boldsymbol{\beta}\) and \(\boldsymbol{\epsilon}\) matrices for this case are as follows:
\(\mathbf{Y} = \begin{bmatrix}Y_{11}\\ Y_{12}\\ Y_{21}\\ Y_{22}\\ Y_{31}\\ Y_{32}\end{bmatrix}\text{, }\mathbf{X} = \begin{bmatrix}1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & -1 & -1\\ 1 & -1 & -1\end{bmatrix}\text{, }\boldsymbol{\beta} = \begin{bmatrix}\beta_{0}\\ \beta_{1}\\ \beta_{2}\end{bmatrix}\text{, and } \boldsymbol{\epsilon} = \begin{bmatrix}\epsilon_{11}\\ \epsilon_{12}\\ \epsilon_{21}\\ \epsilon_{22}\\ \epsilon_{31}\\ \epsilon_{32}\end{bmatrix}\)
where \(\beta_0, \beta_1 \text{ and } \beta_2\), corresponds to \(\mu, \tau_1 \text{ and } \tau_2\) respectively.
Note that the vector of expected values \(\mathbf{E}\{\mathbf{Y}\} =\mathbf{X}\boldsymbol{\beta}\), yields the following:
\begin{aligned} \mathbf{E}\{\mathbf{Y}\} &= \mathbf{X}\boldsymbol{\beta} \\ \begin{bmatrix}E\{Y_{11}\}\\ E\{Y_{12}\}\\ E\{Y_{21}\}\\ E\{Y_{22}\}\\E\{Y_{31}\}\\E\{Y_{32}\}\end{bmatrix} &= \begin{bmatrix}1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & -1 & -1\\ 1 & -1 & -1\end{bmatrix}\begin{bmatrix}\beta_{0}\\ \beta_{1}\\ \beta_{2}\end{bmatrix} \\ &= \begin{bmatrix}\mu_{.}+\tau_{1}\\ \mu+\tau_{1}\\ \mu+\tau_{2}\\ \mu+\tau_{2}\\ \mu-\tau_{1}-\tau_{2}\\ \mu-\tau_{1}-\tau_{2}\end{bmatrix} \end{aligned}
Since \(\tau_3=-\tau_{1}-\tau_{2}\), see above, we see that \(E\{Y_{31}\} =E\{Y_{32}\}=\mu+\tau_{3}\). Thus, the above \(\mathbf{X}\) matrix and \(\boldsymbol{\beta}\) vector representation provides the appropriate expected values for all factor levels as expressed below:
\(E\{Y_{ij}\} =\mu+\tau_{i}\)
Using Technology Section
Begin by replacing the design matrix in the IML code of Section 4.2 Step 2 with:
/* The Effects Model */
x={
1 1 0,
1 1 0,
1 0 1,
1 0 1,
1 -1 -1,
1 -1 -1};
Here we have another omission of a treatment level, but for a different reason. In the effects model, we have the constraint \(\Sigma \tau_i =0\). As a result, coding for one treatment level can be omitted.
Running the IML program with this design matrix yields:
Regression Coefficients | |
---|---|
Beta_0 | 3.5 |
Beta_1 | -2 |
Beta_2 | 0 |
ANOVA | ||||
---|---|---|---|---|
Treatment | dF | SS | MS | F |
2 | 16 | 8 | 16 | |
Error | 3 | 1.5 | 0.5 | |
Total | 5 | 17.5 |
The regression coefficient Beta_0 is the overall mean and the coefficients \(\beta_1\) and \(\beta_2\) are \(\tau_1\) and \(\tau_2\) respectively. The estimate for \(\tau_3\) is obtained as \(– (\tau_1)-( \tau_2) = 2.0\).
In Minitab, if we change the coding now to be Effect coding (-1,0,+1), which is the default setting, we get the following:
Regression Equation
y = 3.500 - 2.000 trt_A - 0.000 trt_B + 2.000 trt_C
The ANOVA table is the same as for the dummy-variable regression model above. We can also observe that the factor level means and General Linear F Statistics values obtained for all 3 representations (cell means, dummy coded regression and effects coded regression) are identical, confirming that the 3 representations are identical.
The intermediates were:
xprimex | ||
---|---|---|
6 | 0 | 0 |
0 | 4 | 2 |
0 | 2 | 4 |
check | ||
---|---|---|
1 | 0 | 0 |
0 | 1 | 0 |
0 | 0 | 1 |
xprimey |
---|
21 |
-8 |
-4 |
SumY2 |
---|
89.5 |
CF |
---|
73.5 |
xprimexinv | ||
---|---|---|
0.1666667 | 0 | 0 |
0 | 0.3333333 | -0.166667 |
0 | -0.166667 | 0.3333333 |
First, we input the response variable and design matrix.
y = matrix(c(2,1,3,4,6,5), ncol=1)
x = matrix(c(1,1,0,1,1,0,1,0,1,1,0,1,1,-1,-1,1,-1,-1),6,3,byrow=TRUE)
x
[,1] [,2] [,3] [1,] 1 1 0 [2,] 1 1 0 [3,] 1 0 1 [4,] 1 0 1 [5,] 1 -1 -1 [6,] 1 -1 -1
Applying the same code as in the previous sections, we can then calculate the regression coefficients and create the ANOVA table.
beta = solve(t(x)%*%x)%*%(t(x)%*%y)
beta
[,1] [1,] 3.5 [2,] -2.0 [3,] 0.0
n = nrow(y)
p = ncol(x)
J = matrix(1,n,n)
ss_tot = (t(y)%*%y) - (1/n)*(t(y)%*%J)%*%y
ss_trt = t(beta)%*%(t(x)%*%y) - (1/n)*(t(y)%*%J)%*%y
ss_error = ss_tot - ss_trt
total_df = n-1
trt_df = p-1
error_df = n-p
MS_trt = ss_trt/(p-1)
MS_error = ss_error / error_df
Fstat = MS_trt/MS_error
ANOVA <- data.frame(
c ("","Treatment","Error", "Total"),
c("DF", trt_df,error_df,total_df),
c("SS", ss_trt, ss_error, ss_tot),
c("MS", MS_trt, MS_error, ""),
c("F",Fstat,"",""),
stringsAsFactors = FALSE)
names(ANOVA) <- c(" ", " ", " ","","")
print(ANOVA, digits = 3)
1 DF SS MS F 2 Treatment 2 16 8 16 3 Error 3 1.5 0.5 4 Total 5 17.5
The intermediate matrix computations can again be calculated with the following commands (outputs not shown).
xprimex = t(x)%*%x
xprimex
xprimey = t(x)%*%y
xprimey
xprimexinv = solve(t(x)%*%x)
xprimexinv
check = xprimexinv*xprimex
check
SumY2 = t(beta)%*%(t(x)%*%y)
SumY2
CF = (1/n)*(t(y)%*%J)%*%y
CF
To use effect coding in R, we can specify the “sum contrast” option within the lm function as follows.
levels = c(1,1,2,2,3,3)
levels = as.factor(levels)
lm1 = lm(y~levels,contrasts = list(levels = contr.sum))
beta = lm1$coefficients
round(beta,1) # rounding to one decimal point
(Intercept) levels1 levels2 3.5 -2.0 0.0
We can then obtain the Type 3 ANOVA table (again, equivalent to the Type 1 in this case; outputs not shown).
aov3 = car::Anova(lm1,type = 3) # syntax alternative to loading the car package
aov3