diff --git a/R/legTools.R b/R/legTools.R
index a3151617e9027b200733a5e61bde5ca383642f78..afb16575aa3d23087ea8de194ac8e75d4675ddf5 100644
--- a/R/legTools.R
+++ b/R/legTools.R
@@ -695,3 +695,101 @@ NULL
#'
#'
NULL
+
+#' @name sugarcaneYield4
+#'
+#' @title Triple factorial NPK fertilization on sugar cane yield
+#'
+#' @description These data are from an \eqn{3^3} factorial experiment
+#' studing the effect of NPK on the yield of sugar cane.
+#'
+#' \itemize{
+#' \item \code{block} a local control factor with 3 levels.
+#' \item \code{rept} factor with 2 levels.
+#' \item \code{N} integer coded nitrogen levels (0, 1, 2).
+#' \item \code{P} integer coded phosphorus levels (0, 1, 2).
+#' \item \code{K} integer coded potassium levels (0, 1, 2).
+#' \item \code{yield} sugar cane yield (ton/ha).
+#' }
+#'
+#' @details There is a missprint in the book for the 9th entry, which
+#' has yield 59.0, that is coded as 202 istead of 220.
+#'
+#' @docType data
+#'
+#' @keywords datasets
+#'
+#' @usage data(sugarcaneYield4)
+#'
+#' @format a \code{data.frame} with 54 records and 6 variables.
+#'
+#' @source Frederico, P. (2009). Curso de Estatística Experimental (15th
+#' ed.). Piracicaba, São Paulo: FEALQ. (page 126)
+#'
+#' @examples
+#'
+#' library(lattice)
+#' library(latticeExtra)
+#' library(multcomp)
+#'
+#' data(sugarcaneYield4)
+#' str(sugarcaneYield4)
+#'
+#' xyplot(yield~N|P, groups=K,
+#' auto.key=list(title="Potassim level", columns=3),
+#' strip=strip.custom(var.name="Phosphorus", strip.names=TRUE,
+#' strip.levels=TRUE, sep=": "),
+#' data=sugarcaneYield4, type=c("p", "a"),
+#' ylab=expression(Yield~(ton~ha^{-1})),
+#' xlab="Nitrogen level level")
+#'
+#' ## Sums in each cell combination.
+#' addmargins(with(sugarcaneYield4, tapply(yield, list(P, N), FUN=sum)))
+#' addmargins(with(sugarcaneYield4, tapply(yield, list(K, N), FUN=sum)))
+#' addmargins(with(sugarcaneYield4, tapply(yield, list(K, P), FUN=sum)))
+#'
+#' sugarcaneYield4 <- transform(sugarcaneYield4,
+#' blockr=interaction(block, rept),
+#' nitro=factor(N),
+#' phosp=factor(P),
+#' potas=factor(K))
+#' str(sugarcaneYield4)
+#'
+#' m0 <- lm(yield~blockr+(nitro+phosp+potas)^3, data=sugarcaneYield4)
+#' par(mfrow=c(2,2)); plot(m0); layout(1)
+#' anova(m0)
+#'
+#' m1 <- update(m0, .~blockr+(nitro+phosp)^2)
+#' par(mfrow=c(2,2)); plot(m1); layout(1)
+#'
+#' anova(m0, m1)
+#' anova(m1)
+#'
+#' m2 <- aov(yield~blockr+nitro/phosp, data=sugarcaneYield4)
+#' anova(m2)
+#'
+#' PinN <- sapply(paste0("nitro", levels(sugarcaneYield4$nitro)),
+#' FUN=grep, x=names(coef(m2))[m2$assign==3L],
+#' simplify=FALSE)
+#'
+#' summary(m2, split=list("nitro:phosp"=PinN))
+#'
+#' X <- model.matrix(m1)
+#' X
+#'
+#' aggregate(X~nitro+phosp, data=sugarcaneYield4, FUN=mean)
+#'
+#' ## It is better use multcomp::LSmatrix().
+#' L <- aggregate(X~nitro+phosp, data=sugarcaneYield4, FUN=mean)
+#' rownames(L) <- with(L, paste0("N", nitro, ":P", phosp))
+#' L <- as.matrix(L[, colnames(X)])
+#' str(L)
+#'
+#' ## Least squares means for N:P combinations.
+#' L%*%coef(m1)
+#'
+#' g1 <- glht(m1, linfct=L)
+#'
+#' confint(g1, calpha=univariate_calpha())
+#'
+NULL
diff --git a/data-raw/sugarcaneYield4.R b/data-raw/sugarcaneYield4.R
new file mode 100644
index 0000000000000000000000000000000000000000..3050b59eef4a0e8f857b92477eb8a8e316df9f2f
--- /dev/null
+++ b/data-raw/sugarcaneYield4.R
@@ -0,0 +1,85 @@
+##----------------------------------------------------------------------
+## Data generation.
+
+sugarcaneYield4 <- read.table("http://www.leg.ufpr.br/~walmes/data/pimentel_npk2.txt",
+ header=TRUE, sep="\t")
+names(sugarcaneYield4)[c(1,6)] <- c("block", "yield")
+sugarcaneYield4 <- transform(sugarcaneYield4, rept=factor(rept))
+str(sugarcaneYield4)
+
+sugarcaneYield4 <- sugarcaneYield4[
+ with(sugarcaneYield4, order(block, rept, N, P, K)),]
+
+save(sugarcaneYield4, file="../data/sugarcaneYield4.RData")
+
+##----------------------------------------------------------------------
+## Examples.
+
+library(lattice)
+library(latticeExtra)
+library(multcomp)
+
+data(sugarcaneYield4)
+str(sugarcaneYield4)
+
+xyplot(yield~N|P, groups=K,
+ auto.key=list(title="Potassim level", columns=3),
+ strip=strip.custom(var.name="Phosphorus", strip.names=TRUE,
+ strip.levels=TRUE, sep=": "),
+ data=sugarcaneYield4, type=c("p", "a"),
+ ylab=expression(Yield~(ton~ha^{-1})),
+ xlab="Nitrogen level level")
+
+## Sums in each cell combination.
+addmargins(with(sugarcaneYield4, tapply(yield, list(P, N), FUN=sum)))
+addmargins(with(sugarcaneYield4, tapply(yield, list(K, N), FUN=sum)))
+addmargins(with(sugarcaneYield4, tapply(yield, list(K, P), FUN=sum)))
+
+sugarcaneYield4 <- transform(sugarcaneYield4,
+ blockr=interaction(block, rept),
+ nitro=factor(N),
+ phosp=factor(P),
+ potas=factor(K))
+str(sugarcaneYield4)
+
+m0 <- lm(yield~blockr+(nitro+phosp+potas)^3, data=sugarcaneYield4)
+par(mfrow=c(2,2)); plot(m0); layout(1)
+anova(m0)
+
+m1 <- update(m0, .~blockr+(nitro+phosp)^2)
+par(mfrow=c(2,2)); plot(m1); layout(1)
+
+anova(m0, m1)
+anova(m1)
+
+m2 <- aov(yield~blockr+nitro/phosp, data=sugarcaneYield4)
+anova(m2)
+
+PinN <- sapply(paste0("nitro", levels(sugarcaneYield4$nitro)),
+ FUN=grep, x=names(coef(m2))[m2$assign==3L],
+ simplify=FALSE)
+
+summary(m2, split=list("nitro:phosp"=PinN))
+
+X <- model.matrix(m1)
+X
+
+aggregate(X~nitro+phosp, data=sugarcaneYield4, FUN=mean)
+
+## It is better use multcomp::LSmatrix().
+L <- aggregate(X~nitro+phosp, data=sugarcaneYield4, FUN=mean)
+rownames(L) <- with(L, paste0("N", nitro, ":P", phosp))
+L <- as.matrix(L[, colnames(X)])
+str(L)
+
+## Least squares means for N:P combinations.
+L%*%coef(m1)
+
+g1 <- glht(m1, linfct=L)
+
+confint(g1, calpha=univariate_calpha())
+
+rm(list=ls())
+load("../data/sugarcaneYield4.RData")
+ls()
+str(sugarcaneYield4)
diff --git a/data/sugarcaneYield4.RData b/data/sugarcaneYield4.RData
new file mode 100644
index 0000000000000000000000000000000000000000..e8821cf112ad01881612e566986a55844359b73c
Binary files /dev/null and b/data/sugarcaneYield4.RData differ
diff --git a/legTools.bmk b/legTools.bmk
index 93387a08a5d6208de2a7f7c6c371ca8a235746ee..758a62b557d41881c488f36e18bf69d696d60edf 100644
--- a/legTools.bmk
+++ b/legTools.bmk
@@ -2,7 +2,43 @@
;;; This format is meant to be slightly human-readable;
;;; nevertheless, you probably don't want to edit it.
;;; -*- End Of Bookmark File Format Version Stamp -*-
-(#1=(#("kornYield" 0 9
+(#1=(#("sugarcaneYield4" 0 15
+ (bmkp-full-record #1#))
+ (filename . "~/GitLab/legTools/R/legTools.R")
+ (buffer-name . "legTools.R")
+ (front-context-string . "@name sugarcaneY")
+ (rear-context-string . "\n#'\n#'\nNULL\n\n#' ")
+ (front-context-region-string)
+ (rear-context-region-string)
+ (visits . 0)
+ (time . #2=(22014 55562 242485 989000))
+ (created . #2#)
+ (position . 21758))
+#1=(#("filterCake" 0 10
+ (bmkp-full-record #1#))
+ (filename . "~/GitLab/legTools/R/legTools.R")
+ (buffer-name . "legTools.R")
+ (front-context-string . "@name filterCake")
+ (rear-context-string . "m1)\n#'\nNULL\n\n#' ")
+ (front-context-region-string)
+ (rear-context-region-string)
+ (visits . 0)
+ (time . #2=(22014 51467 333208 857000))
+ (created . #2#)
+ (position . 20038))
+#1=(#("vinasseFert" 0 11
+ (bmkp-full-record #1#))
+ (filename . "~/GitLab/legTools/R/legTools.R")
+ (buffer-name . "legTools.R")
+ (front-context-string . "@name vinasseFer")
+ (rear-context-string . ")))\n#'\nNULL\n\n#' ")
+ (front-context-region-string)
+ (rear-context-region-string)
+ (visits . 0)
+ (time . #2=(22014 51452 901653 865000))
+ (created . #2#)
+ (position . 18556))
+#1=(#("kornYield" 0 9
(bmkp-full-record #1#))
(filename . "~/GitLab/legTools/R/legTools.R")
(buffer-name . "legTools.R")
diff --git a/man/sugarcaneYield4.Rd b/man/sugarcaneYield4.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..4095578ff1bb766bf613b934cc5d768cfc137b74
--- /dev/null
+++ b/man/sugarcaneYield4.Rd
@@ -0,0 +1,98 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/legTools.R
+\docType{data}
+\name{sugarcaneYield4}
+\alias{sugarcaneYield4}
+\title{Triple factorial NPK fertilization on sugar cane yield}
+\format{a \code{data.frame} with 54 records and 6 variables.}
+\source{
+Frederico, P. (2009). Curso de Estatística Experimental (15th
+ ed.). Piracicaba, São Paulo: FEALQ. (page 126)
+}
+\usage{
+data(sugarcaneYield4)
+}
+\description{
+These data are from an \eqn{3^3} factorial experiment
+ studing the effect of NPK on the yield of sugar cane.
+
+\itemize{
+ \item \code{block} a local control factor with 3 levels.
+ \item \code{rept} factor with 2 levels.
+ \item \code{N} integer coded nitrogen levels (0, 1, 2).
+ \item \code{P} integer coded phosphorus levels (0, 1, 2).
+ \item \code{K} integer coded potassium levels (0, 1, 2).
+ \item \code{yield} sugar cane yield (ton/ha).
+}
+}
+\details{
+There is a missprint in the book for the 9th entry, which
+ has yield 59.0, that is coded as 202 istead of 220.
+}
+\examples{
+library(lattice)
+library(latticeExtra)
+library(multcomp)
+
+data(sugarcaneYield4)
+str(sugarcaneYield4)
+
+xyplot(yield~N|P, groups=K,
+ auto.key=list(title="Potassim level", columns=3),
+ strip=strip.custom(var.name="Phosphorus", strip.names=TRUE,
+ strip.levels=TRUE, sep=": "),
+ data=sugarcaneYield4, type=c("p", "a"),
+ ylab=expression(Yield~(ton~ha^{-1})),
+ xlab="Nitrogen level level")
+
+## Sums in each cell combination.
+addmargins(with(sugarcaneYield4, tapply(yield, list(P, N), FUN=sum)))
+addmargins(with(sugarcaneYield4, tapply(yield, list(K, N), FUN=sum)))
+addmargins(with(sugarcaneYield4, tapply(yield, list(K, P), FUN=sum)))
+
+sugarcaneYield4 <- transform(sugarcaneYield4,
+ blockr=interaction(block, rept),
+ nitro=factor(N),
+ phosp=factor(P),
+ potas=factor(K))
+str(sugarcaneYield4)
+
+m0 <- lm(yield~blockr+(nitro+phosp+potas)^3, data=sugarcaneYield4)
+par(mfrow=c(2,2)); plot(m0); layout(1)
+anova(m0)
+
+m1 <- update(m0, .~blockr+(nitro+phosp)^2)
+par(mfrow=c(2,2)); plot(m1); layout(1)
+
+anova(m0, m1)
+anova(m1)
+
+m2 <- aov(yield~blockr+nitro/phosp, data=sugarcaneYield4)
+anova(m2)
+
+PinN <- sapply(paste0("nitro", levels(sugarcaneYield4$nitro)),
+ FUN=grep, x=names(coef(m2))[m2$assign==3L],
+ simplify=FALSE)
+
+summary(m2, split=list("nitro:phosp"=PinN))
+
+X <- model.matrix(m1)
+X
+
+aggregate(X~nitro+phosp, data=sugarcaneYield4, FUN=mean)
+
+## It is better use multcomp::LSmatrix().
+L <- aggregate(X~nitro+phosp, data=sugarcaneYield4, FUN=mean)
+rownames(L) <- with(L, paste0("N", nitro, ":P", phosp))
+L <- as.matrix(L[, colnames(X)])
+str(L)
+
+## Least squares means for N:P combinations.
+L\%*\%coef(m1)
+
+g1 <- glht(m1, linfct=L)
+
+confint(g1, calpha=univariate_calpha())
+}
+\keyword{datasets}
+