Load the libraries

library(daewr)
## Registered S3 method overwritten by 'DoE.base':
##   method           from       
##   factorize.factor conf.design

Simulating data

A simple way of constructing a randomized data collection form, dividing n experimental units into t treatment groups, can be accomplished using base R commands. For example, in the bread rise experiment, if the experimenter wants to examine three different rise times (35 minutes, 40 minutes, and 45 minutes) and test four replicate loaves of bread at each rise time, the following code will create the list.

set.seed(7638)
# 3 treatment factors (t's) each replicated 4 times
f <- factor(rep(c(35,40,45), each = 4))
f
##  [1] 35 35 35 35 40 40 40 40 45 45 45 45
## Levels: 35 40 45
# Randomization
fac <- sample(f, 12)
fac
##  [1] 40 40 35 45 35 45 35 45 40 40 35 45
## Levels: 35 40 45
# Experimental units
eu <- 1:12

plan <- data.frame(loaf = eu, time = fac, height = "")
plan

Let’s use daewr’s bread data

data("bread")
bread

Linear model

mod0 <- lm(height ~ time, data = bread)
summary(mod0)
## 
## Call:
## lm(formula = height ~ time, data = bread)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.812 -1.141  0.000  1.266  2.250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.4375     0.7655   7.104 5.65e-05 ***
## time40        2.8125     1.0825   2.598   0.0288 *  
## time45        2.8750     1.0825   2.656   0.0262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.531 on 9 degrees of freedom
## Multiple R-squared:  0.5056, Adjusted R-squared:  0.3958 
## F-statistic: 4.602 on 2 and 9 DF,  p-value: 0.042

the mod0 is the name of a model previously fit with the R function lm, the string in quotes is the name of the factor in the model whose cell means are compared, and the vector c(1 -1,0) are the contrast coefficients, ci. This produces the result (−2.8125), which is the negative of the second estimate produced

library(gmodels)

fit.contrast(mod0, "time", c(1, -1, 0))
##                   Estimate Std. Error   t value   Pr(>|t|)
## time c=( 1 -1 0 )  -2.8125   1.082532 -2.598076 0.02882905
## attr(,"class")
## [1] "fit_contrast"

In the model for the CRD, the statistical hypothesis of interest is H0 ∶ µ1 = µ2 =… µt or τ1 = τ2 =… = τt versus the alternative Ha ∶ at least two of the τ s differ. If the null hypothesis is true, the model yij = µi+ Eij = µ+τi+Eij simplifies to yij = µ + Eij, which can be represented as a single normal distribution with mean µ and variance σ 2 rather than multiple normal distributions

ANOVA Table

In this table the ssT and msT and the associated degrees of freedom are on the line labeled time, the ssE is on the line labeled Residuals, and the ssT otal can be computed by adding the ssT to the ssE. The F-value is the ratio msT/msE and the last column labeled Pr(>F) is the probability of exceeding the calculated F-value if the null hypothesis is true. This is called the P-value and is illustrated graphically in Figure 2.3. If the experimenter chooses the significance level, α, for his hypothesis test, he would reject the hypothesis if the Pr(>F) value on the aov output is less than the chosen value of α.

For the bread rise experiment there are significant differences among the mean risen dough heights for each rise time at the significance level α = 0.05, since 0.042 < 0.05.

mod1 <- aov(height ~ time, data = bread)
summary(mod1)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## time         2  21.57  10.786   4.602  0.042 *
## Residuals    9  21.09   2.344                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
f_val <- 10.786/2.344
f_val
## [1] 4.601536

Plots of the model

par(mfrow = c(2,2))
plot(mod1, which = 5)
plot(mod1, which = 1)
plot(mod1, which = 2)
plot(residuals(mod1) ~ loaf, main = "Residuals vs. Exp. Units", data = bread)
abline(h = 0, lty = 2)

Box-cox Transformations

Lambda maximizes the log likelihood which is the reciprocal of error sum of squares

library(MASS)
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:daewr':
## 
##     cement, chem
bc <- boxcox(mod1)

lambda <- bc$x[which.max(bc$y)]
lambda
## [1] -0.5050505

Applying transformation to the data and model it

The resulting ANOVA table below shows the P-value for the factor time has decreased to 0.0209 from the 0.042 value shown in the earlier ANOVA of the untransformed data. Therefore the transformation has made the analysis slightly more sensitive.

library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
tbread <- bread %>%
  mutate(theight = height^(lambda))
tbread
mod2 <- aov(theight ~ time, data = tbread)
summary(mod2)
##             Df  Sum Sq  Mean Sq F value Pr(>F)  
## time         2 0.01732 0.008662   6.134 0.0209 *
## Residuals    9 0.01271 0.001412                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Verify assumptions with plots

par(mfrow = c(2,2))
plot(mod2, which = 5)
plot(mod2, which = 1)
plot(mod2, which = 2)
plot(residuals(mod2) ~ loaf, data = tbread, main = "Residuals vs. Exp. Units")
abline(h = 0)

Using weighted least squares for transformation

with(bread, {
  std <- tapply(height, time, sd)
  weights <- rep(1/std, each = 4)
  mod3 <- lm(height ~ time, weights = weights, data = bread)
  anova(mod3)
})

New Data - Teach

teach

The formula in the full model, modf, includes the treatment factor, while the formula in the reduced model, modr, only includes the intercept.

modf <- polr(score ~ method, weight = count, data = teach)
modr <- polr(score ~ 1, weight = count, data = teach)
anova(modf, modr)

Calculating F-power for bread replicates

From this we can see that with r = 5 replicates there would be a 73% chance of detecting a difference in cell means as large as 3.0, and with r = 6 there is a 83% chance. With fewer than five replicates there is at least a 40% chance this difference will be missed. As a rule of thumb, the number of replicates that result in power between 0.80 and 0.90 is usually sufficient for most experimental studies.

# considering a minimum of 2 replicates
rmin <- 2
# considering a maximum of 6 replicates
rmax <- 6
# significance of 0.05
alpha <- rep(0.05, rmax - rmin - 1)
# assume variance of experimental error ~ 2.1
sigma <- sqrt(2.1)
nlev <- 3
nreps <- rmin:rmax
# A delta of 3 is considered significant in response for bread height
Delta = 3
# power calculation
power <- Fpower1(alpha, nlev, nreps, Delta, sigma)
## Warning in cbind(alpha, nlev, nreps, Delta, sigma, power): number of rows
## of result is not a multiple of vector length (arg 1)
power
##      alpha nlev nreps Delta    sigma     power
## [1,]  0.05    3     2     3 1.449138 0.1947995
## [2,]  0.05    3     3     3 1.449138 0.4041857
## [3,]  0.05    3     4     3 1.449138 0.5903406
## [4,]  0.05    3     5     3 1.449138 0.7328895
## [5,]  0.05    3     6     3 1.449138 0.8329923

Rejection of null hypothesis in CRD only tells us that one of the treatment effect means is not equal to the rest

We need to do pre-planned comparisons for further investigation

sugarbeet
mod4 <- aov(yield ~ treat, data = sugarbeet)

The first comparison asks the question: Does a mix of artificial fertilizers change yield? The second comparison asks the question: Is there a difference in yields between plowed and broadcast application of artificial fertilizer? The third comparison asks the question: Does timing of the application change the yield?

con <- matrix(c(1, -1/3, -1/3, -1/3,
                0, 1, -1, 0, 
                0, 0, 1, -1), 4, 3)
con
##            [,1] [,2] [,3]
## [1,]  1.0000000    0    0
## [2,] -0.3333333    1    0
## [3,] -0.3333333   -1    1
## [4,] -0.3333333    0   -1
L <- t(con)

rownames(L) <- c(" - fertilizer effect", " - plowed vs. broadcast", " - Jan vs. April")
L
##                         [,1]       [,2]       [,3]       [,4]
##  - fertilizer effect       1 -0.3333333 -0.3333333 -0.3333333
##  - plowed vs. broadcast    0  1.0000000 -1.0000000  0.0000000
##  - Jan vs. April           0  0.0000000  1.0000000 -1.0000000

Model

mod4
## Call:
##    aov(formula = yield ~ treat, data = sugarbeet)
## 
## Terms:
##                    treat Residuals
## Sum of Squares  291.0050   29.5865
## Deg. of Freedom        3        14
## 
## Residual standard error: 1.453727
## Estimated effects may be unbalanced
fit.contrast(mod4, "treat", coeff = L)
##                              Estimate Std. Error     t value     Pr(>|t|)
## treat - fertilizer effect        -8.8  0.8252024 -10.6640498 4.188782e-08
## treat - plowed vs. broadcast     -3.8  0.9751895  -3.8966785 1.612261e-03
## treat - Jan vs. April             0.1  0.9194175   0.1087645 9.149328e-01
## attr(,"class")
## [1] "fit_contrast"

Contrasts for bread example

contrasts(bread$time) <- contr.poly(3)
contrasts(bread$time)
##               .L         .Q
## 35 -7.071068e-01  0.4082483
## 40 -7.850462e-17 -0.8164966
## 45  7.071068e-01  0.4082483
mod3 <- aov(height ~ time, data = bread)
summary.lm(mod3)
## 
## Call:
## aov(formula = height ~ time, data = bread)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.812 -1.141  0.000  1.266  2.250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.3333     0.4419  16.593 4.68e-08 ***
## time.L        2.0329     0.7655   2.656   0.0262 *  
## time.Q       -1.1227     0.7655  -1.467   0.1765    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.531 on 9 degrees of freedom
## Multiple R-squared:  0.5056, Adjusted R-squared:  0.3958 
## F-statistic: 4.602 on 2 and 9 DF,  p-value: 0.042

Pairwise Comparisons (Tukey’s HSD)

sugarbeet
mod <- aov(yield ~ treat, data = sugarbeet)
mod.tukey <- TukeyHSD(mod, ordered = T)
mod.tukey
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = yield ~ treat, data = sugarbeet)
## 
## $treat
##     diff        lwr       upr     p adj
## B-A  6.3  3.3122236  9.287776 0.0001366
## D-A 10.0  7.1655464 12.834454 0.0000004
## C-A 10.1  7.2655464 12.934454 0.0000003
## D-B  3.7  0.8655464  6.534454 0.0094231
## C-B  3.8  0.9655464  6.634454 0.0077551
## C-D  0.1 -2.5723484  2.772348 0.9995162

Newman-Keuls snk test

library(agricolae)
compare <- SNK.test(mod, "treat", alpha = 0.05)
print(compare)
## $statistics
##    MSerror Df     Mean       CV
##   2.113321 14 45.68333 3.182182
## 
## $parameters
##   test name.t ntr alpha
##    SNK  treat   4  0.05
## 
## $snk
## NULL
## 
## $means
##   yield      std r  Min  Max      Q25   Q50      Q75
## A  38.7 1.645853 4 36.9 40.1 37.50525 38.90 40.09475
## B  45.0 1.334166 4 43.7 46.2 43.92500 45.05 46.12500
## C  48.8 1.102270 5 47.7 50.6 48.20000 48.60 48.90000
## D  48.7 1.677796 5 47.2 51.3 47.30000 48.50 49.20000
## 
## $comparison
## NULL
## 
## $groups
##   yield groups
## C  48.8      a
## D  48.7      a
## B  45.0      b
## A  38.7      c
## 
## attr(,"class")
## [1] "group"

#When using the Dunnett method the glht function (by default) uses the first level of the treatment factor as the control.

library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(mvtnorm)
library(survival)
sugar.dun <- glht(mod4, linfct = mcp(treat = "Dunnett"),
                  alternative = "greater")
summary(sugar.dun)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: aov(formula = yield ~ treat, data = sugarbeet)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value  Pr(>t)    
## B - A <= 0   6.3000     1.0279   6.129 1.3e-05 ***
## C - A <= 0  10.1000     0.9752  10.357 < 1e-05 ***
## D - A <= 0  10.0000     0.9752  10.254 < 1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)