4
$\begingroup$

I tried to reproduce an example from John Lawson's book "Design and Analysis of Experiments with R" and got a different result. The analysis is rather simple:

library(daewr)
data(bioequiv)
lm.out = lm(y ~ Subject + Period + Treat + Carry, data=bioequiv,
            contrasts = list(Subject=contr.sum, Period=contr.sum, 
                             Treat=contr.sum, Carry=contr.sum))
car::Anova(lm.out, type="III", singular.ok=TRUE)

The result reads:

Anova Table (Type III tests)

Response: y
          Sum Sq Df F values Pr(>F)    
Subject   403586 35  21.9194 <2e-16 ***
Period        38  1   0.0724 0.7888    
Treat       2209  1   4.1993 0.0443 *  
Carry       1051  1   1.9970 0.1622    
Residuals  35772 68

There are two main differences:

  1. The degrees of freedom in Period. As Period is a factor with three levels, there should be two degrees of freedom. This can also be cross-checked by using the standard anova (type I):
anova(lm.out)
  1. The result does not contain an intercept.
Improve this question
$\endgroup$
2
  • $\begingroup$ Questions about code are off topic here, but it looks like R is treating period as if it was continuous. $\endgroup$ Commented yesterday
  • $\begingroup$ That was exactly my first guess, so I explicitly defined Period as factor. In addition, I replaced Period = {1,2,3} by Period={'A','B','C'}. The results does not change. $\endgroup$ Commented yesterday

1 Answer 1

Reset to default
8
$\begingroup$

The problem here is twofold:

  1. Carry == none is completely aliased with Period == 1: this makes your model of less than full rank and is likely what prompted you to include singular.ok=TRUE in the ANOVA call. Really only one of the two remaining Carry values contributes any information because they are also aliased with the remaining periods.
  2. When provided with such model car::Anova changes its default behaviour, as indicated by the message that comes out of your call:

Note: model has aliased coefficients
sums of squares computed by model comparison

In practice this message means that they swap to car::Anova.glm which in this case calls the internal car:::Anova_III_F_glm method. That in turn performs inference on individual terms via stats::drop1 (see here).

Now, the problem is that when Period is dropped from the model, Carry is no longer aliased with it. An extra coefficient just became estimable. Because drop1 calculates the degrees of freedom by comparing the ranks of the full and reduced model, and this last one suddenly is of full rank, it is off by one. You can see this behaviour using just drop1 or by inspecting the models:

## Full model (car::Anova will convert to glm if not of full rank)
glm <- glm(y ~ Subject + Treat + Period + Carry, data=bioequiv)

coef(glm)[37:41]
#>  TreatB Period2 Period3  CarryB Carrynone 
#>  9.5940 -7.7646 -6.3104  7.6397        NA   ## Aliased

drop1(glm, "Period")
#>        Df Deviance    AIC
#> Period  1    35810 1013.3   ## Df=1?

## Reduced model
glm0 <- glm(y ~ Subject + Treat + Carry, data=bioequiv)

coef(glm0)[37:39]
#>  TreatB  CarryB Carrynone 
#>  9.5940  7.6397    7.0375   ## No longer aliased!

glm$rank - glm0$rank
#> 1  ## Oops... expected 2

Is this a bug? Maybe? I don't see anything in drop1's documentation about this behaviour, but on the other hand it's not exactly good form to be working with such less-than-full-rank models... The real solution is probably to parametrize the Carry effect as a single dummy for any of the levels other than none. This will give consistent results, for example:

contrasts <- list(Subject=contr.sum, Period=contr.sum, Treat=contr.sum)

df <- within(bioequiv, {CB <- 0+(Carry=="B")})
lm <- lm(y~Subject + Period + Treat + CB, data=df, contrasts=contrasts)

car::Anova(lm, type="III")
#>             Sum Sq Df   F value Pr(>F)    
#> (Intercept) 633029  1 1203.3274 <2e-16 ***
#> Subject     403586 35   21.9194 <2e-16 ***
#> Period         930  2    0.8835 0.4180 
#> ...

## n.b. car did not use stats::drop1() in the above Anova call
drop1(lm, "Period")
#>        Df Sum of Sq   RSS    AIC
#> Period  2    929.54 36702 705.47   ## Matches
Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ I am impressed! Thank you! $\endgroup$ Commented yesterday
  • 3
    $\begingroup$ +1. A good object lesson in the critical importance of actually reading and reporting error and warning messages issued by the code! $\endgroup$ Commented yesterday

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.