From 130b993f63197131b3e0e75faa4d192b6238f71c Mon Sep 17 00:00:00 2001
From: Walmes Zeviani <walmeszeviani@gmail.com>
Date: Sun, 20 Sep 2015 13:07:50 -0300
Subject: [PATCH] Add dataset of a 3^3 RBD experiment, Pimentel pg 126.

---
 R/legTools.R               |  98 +++++++++++++++++++++++++++++++++++++
 data-raw/sugarcaneYield4.R |  85 ++++++++++++++++++++++++++++++++
 data/sugarcaneYield4.RData | Bin 0 -> 454 bytes
 legTools.bmk               |  38 +++++++++++++-
 man/sugarcaneYield4.Rd     |  98 +++++++++++++++++++++++++++++++++++++
 5 files changed, 318 insertions(+), 1 deletion(-)
 create mode 100644 data-raw/sugarcaneYield4.R
 create mode 100644 data/sugarcaneYield4.RData
 create mode 100644 man/sugarcaneYield4.Rd

diff --git a/R/legTools.R b/R/legTools.R
index a315161..afb1657 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 0000000..3050b59
--- /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
GIT binary patch
literal 454
zcmb2|=3oE==DVT3eun}C>h>sc9@&3Hbh?6IkD;KitHY`Mi@fjFER{%Ry}anyT=oY{
z6%6_hSYvixVy<B9Z`}8-HaV)9lTGZEx7_FNwU=i$^8Nl9W_8ia;Fv|QMS6^2)|Bpk
zj_yr6CV0Pl+WByX+jT3Ot#7_bC7<dlINlN?P;|Hb+_O)Tv8iU}CT%lvV<Wd|Z%MV+
z+AGM*b=@P_WMSLR;$^uzT>kppFg#OtUhvk7_Qi8jSoZ7moAdqa<1CQ6r=^+M8Ngv!
z=GK!YVsroMM{&c0|CAW-)n9D>%UU0u^0~p5GfgT_z?j#(Bl63y;NtE3oOMIyAFT|W
z{W)k=|E2gE*{<T{(-ew5HNVRl-3!cqkl*P4v*f`hXU*5&W#6xO+Be1eL->T{Ts5+n
z{MA3n8eCzIYo6zAv1+%*W78w$?3d;*tWq!+%s-%a_+DGA(|y(}(_QOizrH@O`a`ho
zS#`<t%*#~||IKBXUqAWdDFqvS-oK%W#p|BDoZ<YLBR;kww7j>tJFK*4&hfJIHqC;i
zYr5@p*I8*=X^K9&Q}KUV{HEh`@0E7Vd$5-;M(p)ov1IXtMY&&kLyilIs7-(Hf|>LB
KE0s;`3=9C{;nIfy

literal 0
HcmV?d00001

diff --git a/legTools.bmk b/legTools.bmk
index 93387a0..758a62b 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 0000000..4095578
--- /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}
+
-- 
GitLab