diff --git a/data-raw/cowmilkYield.R b/data-raw/cowmilkYield.R
new file mode 100644
index 0000000000000000000000000000000000000000..6df834f87289614b4bc29b79fca20ec49f1a6f29
--- /dev/null
+++ b/data-raw/cowmilkYield.R
@@ -0,0 +1,45 @@
+##----------------------------------------------------------------------
+## Data generation. Pimentel page 269.
+
+cowmilkYield <- expand.grid(casein=c(0, 10, 15, 20),
+                            block=gl(3, 1))
+cowmilkYield$yield <- c(431.4, 687.5, 679.2, 569.7, 485.2, 560.4, 563.3,
+                        502.5, NA, 443, 430.5, 462.4)
+
+
+save(cowmilkYield, file="../data/cowmilkYield.RData")
+
+##----------------------------------------------------------------------
+
+m0 <- lm(yield~block+factor(casein), data=cowmilkYield)
+anova(m0)
+
+## The imputed value is the estimated value by the missing model.
+predict(m0, newdata=cowmilkYield)
+
+cowmilkYield$yield[9] <- 309.8
+m0 <- lm(yield~block+factor(casein), data=cowmilkYield)
+anova(m0)
+
+aggregate(yield~casein, data=cowmilkYield, FUN=mean)
+
+addmargins(with(cowmilkYield,
+                tapply(yield, list(casein, block), FUN=sum)))
+
+##----------------------------------------------------------------------
+## Examples.
+
+library(lattice)
+
+data(cowmilkYield)
+str(cowmilkYield)
+
+xyplot(yield~casein, groups=block,
+       data=cowmilkYield, type="o",
+       ylab=expression(Milk~yield~(kg)),
+       xlab=expression(Casein~(g~day^{-1})))
+
+rm(list=ls())
+load("../data/cowmilkYield.RData")
+ls()
+str(cowmilkYield)
diff --git a/data-raw/cowmilkYield2.R b/data-raw/cowmilkYield2.R
new file mode 100644
index 0000000000000000000000000000000000000000..78e475181ad14208da0d74391ed2a030fee06348
--- /dev/null
+++ b/data-raw/cowmilkYield2.R
@@ -0,0 +1,58 @@
+##----------------------------------------------------------------------
+## Data generation. Pimentel page 272.
+
+cowmilkYield2 <- expand.grid(period=gl(3, 1),
+                             cow=gl(3, 1),
+                             group=gl(4, 1))
+cowmilkYield2$feed <- factor(c(1, 3, 2, 2, 1, 3, 3, 2, 1, 1, 2, 3, 2, 3,
+                               1, 3, 1, 2, 1, 3, 2, 2, 1, 3, 3, 2, 1, 1,
+                               2, 3, 2, 3, 1, 3, 1, 2),
+                             labels=c("A", "B", "C"))
+cowmilkYield2$yield <- c(527L, 883L, 785L, 696L, 635L, 901L, 989L, 899L,
+                         657L, 608L, 715L, 844L, 885L, 1087L, 711L,
+                         940L, 766L, 832L, 472L, 819L, 778L, 734L, 644L,
+                         953L, 897L, 766L, 706L, 586L, 723L, 892L, 635L,
+                         799L, 595L, 805L, 542L, 681L)
+
+names(cowmilkYield2)
+cowmilkYield2 <- cowmilkYield2[, c(3,1,2,4,5)]
+cowmilkYield2 <- cowmilkYield2[
+    with(cowmilkYield2, order(group, period, cow)), ]
+str(cowmilkYield2)
+
+save(cowmilkYield2, file="../data/cowmilkYield2.RData")
+
+##----------------------------------------------------------------------
+
+m0 <- aov(terms(yield~group/(cow+period)+feed+group:feed,
+                keep.order=TRUE),
+          data=cowmilkYield2)
+anova(m0)
+
+m0 <- aov(terms(yield~group/(cow+period)+feed,
+                keep.order=TRUE),
+          data=cowmilkYield2)
+anova(m0)
+
+aggregate(yield~feed, data=cowmilkYield2, FUN=mean)
+
+addmargins(with(cowmilkYield2,
+                tapply(yield, list(feed, group), FUN=sum)))
+
+##----------------------------------------------------------------------
+## Examples.
+
+library(lattice)
+
+data(cowmilkYield2)
+str(cowmilkYield2)
+
+xyplot(yield~feed, groups=group,
+       data=cowmilkYield2, type=c("p", "a"),
+       ylab=expression(Milk~yield~(kg)),
+       xlab="Feed")
+
+rm(list=ls())
+load("../data/cowmilkYield2.RData")
+ls()
+str(cowmilkYield2)
diff --git a/data-raw/wgChickens.R b/data-raw/wgChickens.R
new file mode 100644
index 0000000000000000000000000000000000000000..8bd26bdd43aa840bb357f7bfacbbc1b3d1874dad
--- /dev/null
+++ b/data-raw/wgChickens.R
@@ -0,0 +1,54 @@
+##----------------------------------------------------------------------
+## Data generation. Pimentel page 267.
+
+## wgChickens <- read.table("clipboard", header=FALSE, sep="\t")
+## wgChickens <- wgChickens[with(wgChickens, order(V1, V2)), ]
+## names(wgChickens) <- c("gender", "conc", "n", "twg")
+## dput(wgChickens)
+
+wgChickens <- structure(list(
+    gender = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
+                         2L, 2L, 2L, 2L),
+                       .Label = c("F", "M"),
+                       class = "factor"), 
+    conc = c(0L, 0L, 10L, 10L, 20L, 20L, 30L, 30L, 0L, 0L, 10L, 10L,
+             20L, 20L, 30L, 30L),
+    n = c(12L, 12L, 13L, 12L, 13L, 12L, 13L, 12L, 13L, 13L, 13L, 13L,
+          13L, 13L, 13L, 13L),
+    twg = c(399L, 388L, 503L, 508L, 475L, 437L, 398L, 448L, 548L, 512L,
+            689L, 646L, 543L, 611L, 514L, 537L)),
+            .Names = c("gender", "conc", "n", "twg"),
+            row.names = c(9L, 13L, 10L, 14L, 11L, 15L, 12L, 16L, 1L, 5L,
+                          2L, 6L, 3L, 7L, 4L, 8L),
+            class = "data.frame")
+
+save(wgChickens, file="../data/wgChickens.RData")
+
+##----------------------------------------------------------------------
+
+m0 <- lm(twg~factor(conc)*gender, data=wgChickens)
+anova(m0)
+
+## The number of animals is not considered in the book analysis.
+m0 <- lm((twg/n)~factor(conc)*gender, data=wgChickens, weights=n)
+anova(m0)
+
+##----------------------------------------------------------------------
+## Examples.
+
+library(lattice)
+
+data(wgChickens)
+str(wgChickens)
+
+xyplot(twg/n~conc, groups=gender,
+       data=wgChickens, type=c("p", "a"),
+       auto.key=list(columns=2, corner=c(0.95, 0.95),
+                     title="Gender"),
+       ylab="Mean weight gain (kg)",
+       xlab="Sorghum concentration in the feed (%)")
+
+rm(list=ls())
+load("../data/wgChickens.RData")
+ls()
+str(wgChickens)