4.5 - The Effects Model

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