1. The data

Lets use the dataset ToothGrowth (again). It desribes how much guniea pig teeth grow, based on type of supplement, and the dose of the supplement. We have to make the dosage level a factor (you do not have to worry about this for your data).

the.data = ToothGrowth
the.data$dose = as.factor(the.data$dose)

2. Fitting models.

We now have potentially 5 models of interest. The interaction model, the TFA no-interactions, the factor A effects model, the factor B effects model, and the “empty or null model” (which would be appropriate if there were no factor A or B effects).

To simpify things a bit, I will assume column 1 is numeric, column 2 is factor A, and column 3 is factor B. I will rename my columns for more efficient typing:

names(the.data) = c("Y","A","B")

Fitting the models are simple in R, and can be done with the following commands. Note, I am using AB to denote the interation model, \(A.B\) to denote the non-interation model, and N to denote the empty model.

AB = lm(Y ~ A*B,the.data)
A.B = lm(Y ~ A + B,the.data)
A = lm(Y ~ A,the.data)
B = lm(Y ~ B,the.data)
N = lm(Y ~ 1, the.data)

For your data, you will either have to rename your columns (in the appropriate order to your data), or you will have to replace the values the.data, Y, A, B with your dataset names.

These have fit the “regression version” of ANOVA models, which I may go over in week 7. However, we can use them to find SSE values, and perform the general F test.

3. Finding SSE values

To find a particular SSE value for a particular model, we may use the general command:

sum(the.model$residuals^2)

I have a few commands that will take all the models above and find the SSE for all of them, and label them nicely:

all.models = list(AB,A.B,A,B,N)
SSE = t(as.matrix(sapply(all.models,function(M) sum(M$residuals^2))))
colnames(SSE) = c("AB","(A+B)","A","B","Empty/Null")
rownames(SSE) = "SSE"

You can remove specific models if you so choose.

4. Hypothesis tests

For finding test-statistics and p-values for various tests, the general command is:

anova(smaller.model, larger.model)

This will return all relevant values of SSE, and d.f. For example, if I wanted to test for interactions, I could use the following:

anova(A.B, AB)
## Analysis of Variance Table
## 
## Model 1: Y ~ A + B
## Model 2: Y ~ A * B
##   Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
## 1     56 820.43                             
## 2     54 712.11  2    108.32 4.107 0.02186 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In row 1 we have (in order): \(d.f\{ SSE_R\}\), \(SSE_R\)

In row 2 we have (in order): \(d.f\{ SSE_F\}\), \(SSE_F\), \(d.f\{ SSE_R\} - d.f\{ SSE_F\}\), \(SSE_R - SSE_F\), \(F_S\), and the p-value.

You can access a specific row and column by saving the table, and using [i,j] for the ith row, jth column. For example, if I want the test-statistic and p-value, we can use:

results = anova(A.B,AB)

Then the p-value is 0.0218603, with test-statistic 4.1069911.

If you wanted to do a different hypothesis test, you would simply change what you defined as the “small” model, and what you defined as the “larger” model. I.e., change the names of the models.

5. Finding partial R^2.

I have built a function (surprise!) to find partial R^2 values for you, although it is not strictly necessary. You can find the values of SSE and calculate it yourself, but the function definition is below:

Partial.R2 = function(small.model,big.model){
  SSE1 = sum(small.model$residuals^2)
  SSE2 = sum(big.model$residuals^2)
  PR2 = (SSE1 - SSE2)/SSE1
  return(PR2)
}

Then, after you have copied the definition into R, to find a partial R^2 you can use:

Partial.R2(smaller.model, larger.model)

For example, if I wanted to find \(R^2\{A | \cdot\}\), I could use:

Partial.R2(N,A)
## [1] 0.05948365