library(daewr)
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
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
data("bread")
bread
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
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"
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
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)
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
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
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)
with(bread, {
std <- tapply(height, time, sd)
weights <- rep(1/std, each = 4)
mod3 <- lm(height ~ time, weights = weights, data = bread)
anova(mod3)
})
teach
modf <- polr(score ~ method, weight = count, data = teach)
modr <- polr(score ~ 1, weight = count, data = teach)
anova(modf, modr)
# 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
sugarbeet
mod4 <- aov(yield ~ treat, data = sugarbeet)
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
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(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
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
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)