Latin Square
I will use Wear experiment as example (p 55-)
First let's enter data:
    > weight<-c(235,251,234,195,236,241,273,270,218,227,274,230,268,229,226,225)
    > position<-as.factor(rep(c(1,2,3,4),c(4,4,4,4)))
    > applications<-as.factor(rep(c(1,2,3,4),4))
    > type<-LETTERS[c(3,1,4,2,4,2,3,1,2,4,1,3,1,3,2,4)]
    > weight
     [1] 235 251 234 195 236 241 273 270 218 227 274 230 268 229 226 225
    > position
     [1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
    > applications
     [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
    > type
     [1] "C" "A" "D" "B" "D" "B" "C" "A" "B" "D" "A" "C" "A" "C" "B" "D"
Let's fit model and get anova:
    >options(contrasts = c("contr.treatment", "contr.poly"))
    > fit<-lm(weight~position+applications+type)
    > summary(fit,cor=F)

    Call: lm(formula = weight ~ position + applications + type)
    Residuals:
       Min    1Q Median    3Q Max
     -9.75 -2.25   0.25 3.687   8

    Coefficients:
                      Value Std. Error   t value  Pr(>|t|)
      (Intercept)  254.7500    6.1872    41.1738    0.0000
        position2   26.2500    5.5340     4.7434    0.0032
        position3    8.5000    5.5340     1.5360    0.1755
        position4    8.2500    5.5340     1.4908    0.1866
    applications2   -2.2500    5.5340    -0.4066    0.6984
    applications3   12.5000    5.5340     2.2588    0.0647
    applications4   -9.2500    5.5340    -1.6715    0.1457
            typeB  -45.7500    5.5340    -8.2671    0.0002
            typeC  -24.0000    5.5340    -4.3368    0.0049
            typeD  -35.2500    5.5340    -6.3697    0.0007

    Residual standard error: 7.826 on 6 degrees of freedom
    Multiple R-Squared: 0.9506
    F-statistic: 12.84 on 9 and 6 degrees of freedom, the p-value is 0.002828
    > anova(fit)
    Analysis of Variance Table

    Response: weight

    Terms added sequentially (first to last)
                 Df Sum of Sq  Mean Sq  F Value      Pr(F)
    position      3    1468.5  489.500  7.99184 0.01616848
    applications  3     986.5  328.833  5.36871 0.03901297
    type          3    4621.5 1540.500 25.15102 0.00084982
    Residuals     6     367.5   61.250
    > model.matrix(fit)
       (Intercept) position2 position3 position4 applications2 applications3
     1           1         0         0         0             0             0
     2           1         0         0         0             1             0
     3           1         0         0         0             0             1
     4           1         0         0         0             0             0
     5           1         1         0         0             0             0
     6           1         1         0         0             1             0
     7           1         1         0         0             0             1
     8           1         1         0         0             0             0
     9           1         0         1         0             0             0
    10           1         0         1         0             1             0
    11           1         0         1         0             0             1
    12           1         0         1         0             0             0
    13           1         0         0         1             0             0
    14           1         0         0         1             1             0
    15           1         0         0         1             0             1
    16           1         0         0         1             0             0
       applications4 typeB typeC typeD
     1             0     0     1     0
     2             0     0     0     0
     3             0     0     0     1
     4             1     1     0     0
     5             0     0     0     1
     6             0     1     0     0
     7             0     0     1     0
     8             1     0     0     0
     9             0     1     0     0
    10             0     0     0     1
    11             0     0     0     0
    12             1     0     1     0
    13             0     0     0     0
    14             0     0     1     0
    15             0     1     0     0
    16             1     0     0     1

Can we fit inteructions? Let's try:
    > fit<-lm(weight~(position+applications+type)^2)
    Error in lm.fit.qr(x, y): computed fit is singular, rank 16
    Dumped
    > fit<-aov(weight~(position+applications+type)^2)
    > summary(fit)
                          Df Sum of Sq  Mean Sq
    position               3    1468.5  489.500
    applications           3     986.5  328.833
    type                   3    4621.5 1540.500
    position:applications  6     367.5   61.250