Chapter 36 Modeling Interaction

36.1 One cell interaction

Recall the carbon wire data.

> library(cfcdae)
> wire <- read.table("carbonwire.txt",header=TRUE)
> wire <- within(wire, {degas<-factor(degas);temp<-factor(temp);diff<-factor(diff)})
> wire.fit <- lm(resist~degas*temp*diff,data=wire)
> anova(wire.fit)
Analysis of Variance Table

Response: resist
                Df Sum Sq Mean Sq  F value    Pr(>F)    
degas            1   0.47   0.467   1.5918    0.2125    
temp             2 410.69 205.346 699.6024 < 2.2e-16 ***
diff             2  80.54  40.270 137.1989 < 2.2e-16 ***
degas:temp       2   0.31   0.157   0.5342    0.5892    
degas:diff       2   0.70   0.351   1.1957    0.3104    
temp:diff        4  14.81   3.704  12.6177 2.566e-07 ***
degas:temp:diff  4   0.87   0.219   0.7450    0.5656    
Residuals       54  15.85   0.294                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The temp:diff interaction is very significant. When we look at the plot, it seems like the diff=1,temp=3 factor combination is lower than it should be.

> interactplot(diff,temp,resist,data=wire)
Interaction plot of diff:temp for carbon wire data

Figure 36.1: Interaction plot of diff:temp for carbon wire data

One way to model a one-cell interaction is to create a dummy variable (mostly zeroes with ones in a few special places) that indicates the combination of interest.

> cell <- 1*with(wire,diff==1 & temp==3)
> cell
 [1] 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0
[49] 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0

What we do is fit the dummy variable before the corresponding interaction. The test for the interaction then considers whether there is evidence for anything beyond the one cell.

The data are balanced, but the cell term is one degree of freedom within the temp:diff interaction. We don’t want a Type II sum of squares for cell in this case, because cell would be adjusted for temp:diff, and there wouldn’t be any SS or df left for cell after temp:diff. A Type II SS for temp:diff adjusted for cell works fine.

What we see is that temp:diff is not significant once we have cell in the model. Thus the one-cell interaction is a good model.

> wire.fit.1 <- lm(resist~degas+temp+diff+cell+degas*temp*diff,data=wire)
> anova(wire.fit.1)
Analysis of Variance Table

Response: resist
                Df Sum Sq Mean Sq  F value    Pr(>F)    
degas            1   0.47   0.467   1.5918    0.2125    
temp             2 410.69 205.346 699.6024 < 2.2e-16 ***
diff             2  80.54  40.270 137.1989 < 2.2e-16 ***
cell             1  13.01  13.005  44.3073 1.488e-08 ***
degas:temp       2   0.31   0.157   0.5342    0.5892    
degas:diff       2   0.70   0.351   1.1957    0.3104    
temp:diff        3   1.81   0.603   2.0546    0.1171    
degas:temp:diff  4   0.87   0.219   0.7450    0.5656    
Residuals       54  15.85   0.294                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> car::Anova(wire.fit.1,type=2)
Note: model has aliased coefficients
      sums of squares computed by model comparison
Anova Table (Type II tests)

Response: resist
                Sum Sq Df  F value    Pr(>F)    
degas             0.47  1   1.5918    0.2125    
temp            372.65  2 634.7992 < 2.2e-16 ***
diff             28.91  2  49.2484 6.708e-13 ***
cell                    0                       
degas:temp        0.31  2   0.5342    0.5892    
degas:diff        0.70  2   1.1957    0.3104    
temp:diff         1.81  3   2.0546    0.1171    
degas:temp:diff   0.87  4   0.7450    0.5656    
Residuals        15.85 54                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actually, we need to take a little more care when talking about signficance. The one degree of freedom for the one cell interaction is, in effect, a contrast suggested by the data. Thus a Scheffe adjustment is appropriate. There are 4 df in the temp:diff interaction, so we must divide the F for cell by 4 and treat it as having 4 and 54 df. It is still very significant after this adjustment.

> pf(44.30/4,4,54,lower.tail = FALSE)
[1] 1.245092e-06

If you want to look at the coefficients with the cell indicator, you need to refit without insignificant terms, because everything is unbalanced now.

The interacting cell is about 1.9 units below what the pattern of the rest of the data would predict.

> wire.fit.1.reduced <- lm(resist~degas+temp+diff+cell,data=wire)
> summary(wire.fit.1.reduced)

Call:
lm(formula = resist ~ degas + temp + diff + cell, data = wire)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26806 -0.31840 -0.07639  0.39028  1.13194 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 22.62083    0.07226 313.047  < 2e-16 ***
degas1       0.08056    0.06463   1.246  0.21710    
temp1       -3.15000    0.09695 -32.492  < 2e-16 ***
temp2       -0.18750    0.09695  -1.934  0.05746 .  
diff1       -1.04583    0.11195  -9.342 1.24e-13 ***
diff2        0.28750    0.09695   2.966  0.00422 ** 
cell        -1.91250    0.29084  -6.576 9.68e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5484 on 65 degrees of freedom
Multiple R-squared:  0.9627,    Adjusted R-squared:  0.9593 
F-statistic: 279.7 on 6 and 65 DF,  p-value: < 2.2e-16

36.2 Polynomial models

When factors are quantitative, we know that we can model a factor effect as a polynomial. We can do the same thing in factorial models. Our motivation in doing this is basically the same as in single factor models: we would like to be able to predict at unobserved values of the quantitative factors, and we would generally like to use as few degrees of freedom as possible in our models for the mean.

Recall the amylase activity data. Analysis temperature is quantitative, so we can consider using polynomials instead of simply 8 groups. Here we use orthogonal polynomials. I have selected as high an order of polynomial as possible, but 7 is an unrealistically high order in standard practice.

Growth temperature is also quantitative, but it only has two levels. In such a case, the “linear” degree of freedom and the “between groups” degree of freedom are just describing the same variability with different parameterizations. We’ll look at both.

Recall that we analyzed these data on the log scale.

> library(cfcdae)
> data(AmylaseActivity)
> amylase.quantquant <- lm(log(amylase)~variety*poly(gTemp.z,1)*poly(aTemp.z,7),
+                          data=AmylaseActivity)
> amylase.qualquant <- lm(log(amylase)~variety*gTemp*poly(aTemp.z,7),
+                          data=AmylaseActivity)

We start by looking at the model with gTemp as qualitative. Just looking at the ANOVA table, we would say there is no evidence of a three way interaction, no evidence of a variety by analysis temperature interaction, and equivocal evidence of a growth temperature by analysis temperature interaction.

> anova(amylase.qualquant)
Analysis of Variance Table

Response: log(amylase)
                               Df  Sum Sq Mean Sq  F value    Pr(>F)    
variety                         1 0.58957 0.58957 107.9085 2.305e-15 ***
gTemp                           1 0.00438 0.00438   0.8016 0.3739757    
poly(aTemp.z, 7)                7 3.01613 0.43088  78.8628 < 2.2e-16 ***
variety:gTemp                   1 0.08599 0.08599  15.7392 0.0001863 ***
variety:poly(aTemp.z, 7)        7 0.02758 0.00394   0.7212 0.6543993    
gTemp:poly(aTemp.z, 7)          7 0.08106 0.01158   2.1195 0.0539203 .  
variety:gTemp:poly(aTemp.z, 7)  7 0.04764 0.00681   1.2457 0.2916176    
Residuals                      64 0.34967 0.00546                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

But now look a little closer at the coefficients. A significant linear or quadratic effect can be hidden by being combined with multiple nonsignificant higher order terms. Nothing higher than cubic in analysis temperature is anywhere near significant.

> summary(amylase.qualquant)

Call:
lm(formula = log(amylase) ~ variety * gTemp * poly(aTemp.z, 7), 
    data = AmylaseActivity)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.13079 -0.04861  0.00117  0.05163  0.12640 

Coefficients:
                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        5.805966   0.007544 769.610  < 2e-16 ***
variety1                           0.078367   0.007544  10.388 2.31e-15 ***
gTemp1                            -0.006754   0.007544  -0.895 0.373976    
poly(aTemp.z, 7)1                  0.935612   0.073916  12.658  < 2e-16 ***
poly(aTemp.z, 7)2                 -1.445585   0.073916 -19.557  < 2e-16 ***
poly(aTemp.z, 7)3                 -0.204923   0.073916  -2.772 0.007280 ** 
poly(aTemp.z, 7)4                  0.053281   0.073916   0.721 0.473641    
poly(aTemp.z, 7)5                 -0.001156   0.073916  -0.016 0.987566    
poly(aTemp.z, 7)6                  0.058510   0.073916   0.792 0.431536    
poly(aTemp.z, 7)7                 -0.052764   0.073916  -0.714 0.477924    
variety1:gTemp1                   -0.029929   0.007544  -3.967 0.000186 ***
variety1:poly(aTemp.z, 7)1        -0.033144   0.073916  -0.448 0.655382    
variety1:poly(aTemp.z, 7)2         0.137003   0.073916   1.853 0.068422 .  
variety1:poly(aTemp.z, 7)3        -0.046678   0.073916  -0.632 0.529962    
variety1:poly(aTemp.z, 7)4         0.029442   0.073916   0.398 0.691721    
variety1:poly(aTemp.z, 7)5        -0.061100   0.073916  -0.827 0.411531    
variety1:poly(aTemp.z, 7)6         0.002074   0.073916   0.028 0.977707    
variety1:poly(aTemp.z, 7)7         0.030507   0.073916   0.413 0.681185    
gTemp1:poly(aTemp.z, 7)1          -0.188227   0.073916  -2.546 0.013298 *  
gTemp1:poly(aTemp.z, 7)2          -0.009436   0.073916  -0.128 0.898820    
gTemp1:poly(aTemp.z, 7)3          -0.170622   0.073916  -2.308 0.024224 *  
gTemp1:poly(aTemp.z, 7)4          -0.078812   0.073916  -1.066 0.290325    
gTemp1:poly(aTemp.z, 7)5           0.082983   0.073916   1.123 0.265775    
gTemp1:poly(aTemp.z, 7)6          -0.031378   0.073916  -0.425 0.672615    
gTemp1:poly(aTemp.z, 7)7          -0.048450   0.073916  -0.655 0.514517    
variety1:gTemp1:poly(aTemp.z, 7)1  0.001186   0.073916   0.016 0.987250    
variety1:gTemp1:poly(aTemp.z, 7)2 -0.016851   0.073916  -0.228 0.820396    
variety1:gTemp1:poly(aTemp.z, 7)3  0.201712   0.073916   2.729 0.008195 ** 
variety1:gTemp1:poly(aTemp.z, 7)4  0.005656   0.073916   0.077 0.939247    
variety1:gTemp1:poly(aTemp.z, 7)5 -0.076204   0.073916  -1.031 0.306445    
variety1:gTemp1:poly(aTemp.z, 7)6 -0.001836   0.073916  -0.025 0.980261    
variety1:gTemp1:poly(aTemp.z, 7)7  0.028729   0.073916   0.389 0.698809    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07392 on 64 degrees of freedom
Multiple R-squared:  0.9168,    Adjusted R-squared:  0.8765 
F-statistic: 22.74 on 31 and 64 DF,  p-value: < 2.2e-16

Fitting a model with just up to cubic terms indicates that we didn’t need those higher order terms.

> amylase.qualquant3 <- lm(log(amylase)~variety*gTemp*poly(aTemp.z,3),
+                          data=AmylaseActivity)
> anova(amylase.qualquant3,amylase.qualquant)
Analysis of Variance Table

Model 1: log(amylase) ~ variety * gTemp * poly(aTemp.z, 3)
Model 2: log(amylase) ~ variety * gTemp * poly(aTemp.z, 7)
  Res.Df     RSS Df Sum of Sq     F Pr(>F)
1     80 0.38735                          
2     64 0.34967 16   0.03768 0.431 0.9683

Looking back at the summary, the cubic term in the three way looked significant, but the linear and quadratic did not. Keeping cubic means we have to keep linear and quadratic. Is cubic three way sufficiently important if viewed with all three df? Fit without the three way and compare.

The three-way can reasonably be removed.

> amylase.qualquant3b <- lm(log(amylase)~(variety+gTemp+poly(aTemp.z,3))^2,
+                          data=AmylaseActivity)
> anova(amylase.qualquant3b,amylase.qualquant3,amylase.qualquant)
Analysis of Variance Table

Model 1: log(amylase) ~ (variety + gTemp + poly(aTemp.z, 3))^2
Model 2: log(amylase) ~ variety * gTemp * poly(aTemp.z, 3)
Model 3: log(amylase) ~ variety * gTemp * poly(aTemp.z, 7)
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1     83 0.42832                              
2     80 0.38735  3  0.040973 2.4997 0.06739 .
3     64 0.34967 16  0.037680 0.4310 0.96825  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now think about two way interactions. The only order even vaguely significant in variety by analysis temperature was quadratic (p-value .07). We should be able to safely remove that interaction.

> amylase.qualquant3c <- lm(log(amylase)~variety*gTemp+gTemp*poly(aTemp.z,3)^2,
+                          data=AmylaseActivity)
> anova(amylase.qualquant3c,amylase.qualquant3b,amylase.qualquant3,amylase.qualquant)
Analysis of Variance Table

Model 1: log(amylase) ~ variety * gTemp + gTemp * poly(aTemp.z, 3)^2
Model 2: log(amylase) ~ (variety + gTemp + poly(aTemp.z, 3))^2
Model 3: log(amylase) ~ variety * gTemp * poly(aTemp.z, 3)
Model 4: log(amylase) ~ variety * gTemp * poly(aTemp.z, 7)
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1     86 0.45037                              
2     83 0.42832  3  0.022047 1.3451 0.26764  
3     80 0.38735  3  0.040973 2.4997 0.06739 .
4     64 0.34967 16  0.037680 0.4310 0.96825  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The growth temperature by analysis temperature cubic model is probably needed (two of the three had small p-values), but lets check for completeness.

> amylase.qualquant3d <- lm(log(amylase)~variety*gTemp+poly(aTemp.z,3)^2,
+                          data=AmylaseActivity)
> anova(amylase.qualquant3d,amylase.qualquant3c,amylase.qualquant3b,amylase.qualquant3,amylase.qualquant)
Analysis of Variance Table

Model 1: log(amylase) ~ variety * gTemp + poly(aTemp.z, 3)^2
Model 2: log(amylase) ~ variety * gTemp + gTemp * poly(aTemp.z, 3)^2
Model 3: log(amylase) ~ (variety + gTemp + poly(aTemp.z, 3))^2
Model 4: log(amylase) ~ variety * gTemp * poly(aTemp.z, 3)
Model 5: log(amylase) ~ variety * gTemp * poly(aTemp.z, 7)
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1     89 0.51500                              
2     86 0.45037  3  0.064630 3.9431 0.01207 *
3     83 0.42832  3  0.022047 1.3451 0.26764  
4     80 0.38735  3  0.040973 2.4997 0.06739 .
5     64 0.34967 16  0.037680 0.4310 0.96825  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As expected, that last bit was fairly significant, so we will retain it and our final model has all main effects plus the interation of growth temperature and variety, plus the interaction of growth temperature and analysis temperature, but only up to cubic in analysis temperature.

Now that we know the terms we need, it is often easier to think about the model somewhat differently.

  1. Go back to regular polynomials for better understanding.
  2. The poly(aTemp.z,3) + gTemp:poly(aTemp.z,3) part of the model says that there is an overall cubic coefficient shared by all levels of gTemp plus the cubic coefficient gets adjusted up and down for each different level of gTemp.
  3. If we use instead gTemp:poly(aTemp.z,3) (without the main effect of aTemp), then we just get two separate cubics, one for each level of gTemp.
  4. The gTemp*variety contribution can be thought of as an overall intercept, an adjustment to the intercept depending on gTemp, and adjustment to the intercept depending on variety, and an adjustment to the intercept depending on both gTemp and variety.
  5. If we use gTemp:variety-1, then we just get four intercepts, one for each of the gTemp by variety combinations.

Putting it all together, we have a cubic model in analysis temperature, with different coefficients depending on growth temperature, and the intercepts are determined by the growth temperature by variety combinations.

One interesting thing is that the cubic term only appears to be needed for growth temperature 13.

> amylase.final <- lm(log(amylase)~gTemp:variety-1+gTemp:(aTemp.z +I(aTemp.z^2)+I(aTemp.z^3)),data=AmylaseActivity)
> summary(amylase.final)

Call:
lm(formula = log(amylase) ~ gTemp:variety - 1 + gTemp:(aTemp.z + 
    I(aTemp.z^2) + I(aTemp.z^3)), data = AmylaseActivity)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.144875 -0.049499 -0.001341  0.049588  0.172718 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
gTemp13:varietyB73    5.426e+00  1.957e-01  27.725  < 2e-16 ***
gTemp25:varietyB73    4.813e+00  1.957e-01  24.595  < 2e-16 ***
gTemp13:varietyOh43   5.329e+00  1.957e-01  27.230  < 2e-16 ***
gTemp25:varietyOh43   4.596e+00  1.957e-01  23.489  < 2e-16 ***
gTemp13:aTemp.z      -2.864e-03  2.793e-02  -0.103 0.918577    
gTemp25:aTemp.z       9.043e-02  2.793e-02   3.238 0.001712 ** 
gTemp13:I(aTemp.z^2)  2.566e-03  1.199e-03   2.140 0.035156 *  
gTemp25:I(aTemp.z^2) -1.386e-03  1.199e-03  -1.156 0.251040    
gTemp13:I(aTemp.z^3) -5.816e-05  1.585e-05  -3.670 0.000421 ***
gTemp25:I(aTemp.z^3) -5.312e-06  1.585e-05  -0.335 0.738320    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07237 on 86 degrees of freedom
Multiple R-squared:  0.9999,    Adjusted R-squared:  0.9998 
F-statistic: 6.187e+04 on 10 and 86 DF,  p-value: < 2.2e-16

Let’s look at the results. We will plot the predictions for analysis temperature from 10 to 40, with a different line for each growth temperature by variety combination.

Lines 1 and 2 (growth temperature 13) are non-symmetric, rising more slowly to their max and then falling off more quickly (that’s the cubic contribution). Conversely, lines 3 and 4 (growth temperature 25) are pretty symmetric. The gap between lines 1 and 2 (variety effect) is smaller than the gap between lines 3 and 4; this difference is showing the interaction between growth temperature and variety.

> # set up the prediction frame
> aTemp.z <- rep(10:40,4)
> gTemp <- rep(c(13,25),each=62)
> gTemp <- factor(gTemp)
> variety <- rep(c("B73","Oh43"),each=31,length=124)
> variety <- factor(variety)
> pred.frame <- data.frame(gTemp,aTemp.z,variety)
> # get the predictions and rearrage by gtemp/variety
> preds <- predict(amylase.final,pred.frame)
> preds <- matrix(preds,ncol=4)
> # make a blank plot
> plot(c(10,40),range(preds),xlab="Analysis Temperature",ylab="log(amylase)", type='n')
> # add lines
> for(i in 1:4) lines(10:40,preds[,i],lty=i)
> # add a legend
> legend(25,5.75,c("13 B73","13 Oh43","25 B73","25 Oh43"),lty=1:4)
Predictions using final amylase model

Figure 36.2: Predictions using final amylase model

36.3 Tukey One Degree of Freedom

Let’s look at the dye adsorption data from problem 8.14 in the text. Assume that we only had the subset of the data where temperature was 60 and the oxidant was 355. This gives us a 5x3 single replication in time and initial concentration.

> data(DyeRemoval)
> head(DyeRemoval)
  time.z time oxidant.z oxidant initial.z initial temp.z temp  CRE
1      5    5       355     355        10      10     60   60 14.6
2      5    5       355     355        10      10     70   70 33.5
3      5    5       355     355        15      15     60   60 12.4
4      5    5       355     355        15      15     70   70 20.4
5      5    5       355     355        20      20     60   60 10.6
6      5    5       355     355        20      20     70   70 19.5
> dye.subset <- with(DyeRemoval,DyeRemoval[temp==60 & oxidant==355,])
> with(dye.subset,tapply(CRE,list(time,initial),mean))
     10   15   20
5  14.6 12.4 10.6
10 26.1 18.8 17.3
15 34.3 25.9 21.4
20 39.6 27.6 23.8
25 44.0 30.4 26.1

With a single replication, we fit the additive model. Both factors are highly significant. Residuals have the Tukey pattern (positives in opposite corners, negatives in opposite corners), and a dramatic flopping fish.

> dye.fit <- lm(CRE~time+initial,data=dye.subset)
> anova(dye.fit)
Analysis of Variance Table

Response: CRE
          Df Sum Sq Mean Sq F value    Pr(>F)    
time       4 837.18 209.294  23.746 0.0001708 ***
initial    2 378.23 189.114  21.456 0.0006096 ***
Residuals  8  70.51   8.814                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> with(dye.subset,tapply(residuals(dye.fit),list(time,initial),mean))
          10          15        20
5  -4.793333  1.70666667  3.086667
10 -1.493333 -0.09333333  1.586667
15  0.240000  0.54000000 -0.780000
20  2.406667 -0.89333333 -1.513333
25  3.640000 -1.26000000 -2.380000
> plot(dye.fit,which=1)
Residual plot of dye adsorption (subset)

Figure 36.3: Residual plot of dye adsorption (subset)

So get the rescaled squared predicted values, add to the model, and test.

Tukey is highly significant. That one degree of freedom (out of 8) accounted for 97% of the sum of squares.

The coefficient for rspv is 1.41, so the estimated transformation is power -.5 (roughly). Box-Cox will give a very similar suggestion.

Keep in mind that Box-Cox requires positive responses, but the Tukey model does not (it does if you want to transform, but not just to fit the model).

> rspv <- predict(dye.fit)^2/2/coef(dye.fit)[1]
> dye.fitTukey <- lm(CRE~time+initial+rspv,data=dye.subset)
> anova(dye.fitTukey)
Analysis of Variance Table

Response: CRE
          Df Sum Sq Mean Sq F value    Pr(>F)    
time       4 837.18 209.294  641.41 4.719e-09 ***
initial    2 378.23 189.114  579.56 1.676e-08 ***
rspv       1  68.23  68.228  209.09 1.803e-06 ***
Residuals  7   2.28   0.326                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> summary(dye.fitTukey)

Call:
lm(formula = CRE ~ time + initial + rspv, data = dye.subset)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.67260 -0.33849  0.04872  0.27123  0.78478 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.98911    1.38208   3.610  0.00863 ** 
time1        2.36223    1.05779   2.233  0.06069 .  
time2        2.80808    0.56304   4.987  0.00159 ** 
time3        0.46366    0.32226   1.439  0.19339    
time4       -1.52734    0.56693  -2.694  0.03090 *  
initial1    -3.45636    0.74330  -4.650  0.00234 ** 
initial2     1.38104    0.30517   4.526  0.00271 ** 
rspv         1.41333    0.09774  14.460  1.8e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5712 on 7 degrees of freedom
Multiple R-squared:  0.9982,    Adjusted R-squared:  0.9964 
F-statistic:   562 on 7 and 7 DF,  p-value: 4.384e-09

Note that both time and initial concentration are quantitative factors. We could also have approached this interaction as a linear by linear interaction.