function () 
{
    par(mfcol = c(1, 2))
    out.lm = lm(response ~ as.factor(row) + as.factor(treatment), 
        data = fertilizerdata)
    boxplot(response ~ as.factor(treatment), data = fertilizerdata, 
        xlab = "treatment", ylab = expression("nitrogen uptake (ppm" ~ 
            phantom(0) %*% 10^{
                -2
            } ~ ")"))
    boxplot(response ~ as.factor(row), data = fertilizerdata, 
        xlab = "row", ylab = expression("nitrogen uptake (ppm" ~ 
            phantom(0) %*% 10^{
                -2
            } ~ ")"))
    anova(out.lm)
}