Example 1: Split-plot design with one qualitative and one quantitative level factor
Source:R/example1.R
example1.RdGomez & Gomez (1984, p. 143) report a rice experiment with three management practices (minimum, optimum, intensive), five different amounts of nitrogen (N) fertilizer (0, 50, 80, 110, 140 kg/ha), and three varieties (V1, V2, V3). The experiment involved variety and management as qualitative treatment factors and nitrogen fertilizer as a quantitative treatment factor. Overall, there were 45 treatments with three replicates in complete replicate blocks. The fertilizer treatments were applied to main plots, the management practices to split-plots and the varieties to split-split-plots.
Details
Section 1 Section 1 examines treatment effects by fitting qualitative factorial models and the first analysis calculates a full analysis of variance (Table 1) for main plots (nitrogen), split-plots (management) and split-split-plots (variety). Each type of experimental unit (or "stratum") requires a separate error term in the fitted analysis.
The second analysis (Table 2) uses a REML mixed model analysis to find treatment means and SE's for each marginal treatment classification averaged over all the other treatment factors, together with estimates of pairwise contrasts of treatment means and the SE's of the pairwise treatment comparisons. This analysis fits the full set of nitrogen-by-variety interaction effects assuming additive management effects and the fit of the model is tested by a graphical plot of the model residuals. Residual plots provide an important check on model assumptions but many more options for model testing are available and further methods for diagnostic testing are examined in the subsequent examples.
The third analysis (Table 3) shows a mixed model analysis of the full factorial model fitted by REML using the lmer function of the lme4 package. Generally with mixed models, determination of the denominator degrees of freedom for Wald-type F- and t-statistics becomes an issue, and here we use the method proposed by Kenward & Roger (1997).
Section 2 Section 2 examines treatment effects by fitting polynomial models and the first step calculates a full set of four raw polynomials for the 5-levels of N using the poly() function. The N rates are re-scaled by division by 100 to improve numerical stability.
The second step fits a mixed model polynomial analysis of nitrogen effects assuming additive management effects (Table 7). In this analysis, most of the nitrogen treatment effect can be explained by linear and quadratic trend effects. but it is important to note that there is a non-negligible Variety x Cubic N interaction effect. This suggests that not all the varieties responded in a similar way to the N treatments and that some further analysis of the data may be required (see also the N plots of individual varieties and replicates in Fig 1).
The third step fits the required model for the actual fitted model coefficients (Table 8). When estimating model effects, only effects that are significant for the fitted model or that are marginal to those effects (functional marginality) should be included in the model therefore only linear and quadratic nitrogen effects are included in this model. The fitted model for the nitrogen effects fits the actual actual nitrogen levels used in the experiment therefore this model provides the required coefficients for the actual applied nitrogen levels.
Section 3
Section 3 provides checks on some of the assumptions underlying the blocks-by-treatments model.
The first analysis in this section shows a complete partition of the blocks-by-treatments interaction effects into factorial mean square terms where all the terms that contain a replicate:variety interaction effect are estimates of the split-split-plot error variance. If the blocks-by-treatments assumptions are valid, all the estimates of the split-split-plot error variance are expected to have the same error mean square. However, the Replicate:variety effect has a mean square of 1.54 on 4 degrees of freedom whereas the Replicate:management:variety:nitrogen effect has a mean square of 0.26 on 32 degrees of freedom. The ratio of these mean squares is 5.92 with an F-probability of 0.00110 on 4 and 32 degrees of freedom, which means that the Replicate:variety interaction effect is significantly inflated relative to the Replicate:management:variety:nitrogen effect. This shows that the assumptions underlying the blocks-by-treatments analysis of the model are invalid with a high level of probability.
The 4 degrees of freedom in the Replicate:variety interaction effect are the differences between the three varieties differenced between the three replicate blocks. Fig S1 shows graphical plots of variety effects in each replicate block averaged over management effects, and there is clear evidence that the effects of Variety 1 in blocks 1 and 2 were different from the effects of Variety 1 in block 3.
The second analysis in Section 3 shows a complete partition of the blocks-by-treatments interaction effects into factorial mean square terms ignoring Variety 1. This analysis shows a reasonably good fit to the assumed additive block which supports the hypothesis that the non-additivity of the block-and-treatment effects in the full unrestricted analysis is mainly due to Variety 1.
The final analysis in Section 3 shows an analysis of variance of the treatment effects ignoring Variety 1. In this analysis, the management:variety interaction effect becomes significant at the 0.00992 probability level compared with a non-significant management:variety interaction effect in the analysis of the full data set.
Such anomalies are not uncommon in the analysis of real data sets and it is the task of the statistician to identify anomalies as and when they occur. Factorial designs can be very powerful for practical research but, as demonstrated with this data set, the analysis of such designs is complex and anomalies can be easily missed. Unless an anomaly is due to an easily identified cause such as an incorrectly recorded data point, it is likely that the anomaly will need to be investigated by further discussion with the research workers. It is a mistake to suppose that data from a designed experiment can be analysed statistically in isolation from the research workers who conducted the experiment.
References
Gomez, K. A., & Gomez, A. A. (1984). Statistical Procedures for Agricultural Research. John Wiley & Sons.
Kenward, M. G., & Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983–997.
Piepho, H., & Edmondson, R. (2018). A tutorial on the Statistical Analysis of Factorial Experiments with Qualitative and Quantitative treatment factor levels. Journal of Agronomy and Crop Science. (https://onlinelibrary.wiley.com/doi/full/10.1111/jac.12267).
Author
Rodney N. Edmondson (rodney.edmondson@gmail.com)
Hans-Peter Piepho (piepho@uni-hohenheim.de)
Muhammad Yaseen (myaseen208@gmail.com)
Examples
library(broom)
library(broom.mixed)
library(dplyr)
#>
#> Attaching package: ‘dplyr’
#> The following objects are masked from ‘package:stats’:
#>
#> filter, lag
#> The following objects are masked from ‘package:base’:
#>
#> intersect, setdiff, setequal, union
library(emmeans)
#> Welcome to emmeans.
#> Caution: You lose important information if you filter this package's results.
#> See '? untidy'
library(ggfortify)
#> Loading required package: ggplot2
library(ggplot2)
library(lmerTest)
#> Loading required package: lme4
#> Loading required package: Matrix
#>
#> Attaching package: ‘lmerTest’
#> The following object is masked from ‘package:lme4’:
#>
#> lmer
#> The following object is masked from ‘package:stats’:
#>
#> step
library(magrittr)
library(tidyr)
#>
#> Attaching package: ‘tidyr’
#> The following object is masked from ‘package:magrittr’:
#>
#> extract
#> The following objects are masked from ‘package:Matrix’:
#>
#> expand, pack, unpack
options(contrasts = c('contr.treatment', 'contr.poly'))
##----fm1.1----
fm1.1 <- aov(yield ~ Replicate + nitrogen * management * variety +
Error(Replicate/nitrogen/management), rice)
fm1.1.Summary <- broom::tidy(fm1.1)
fm1.1.Summary
#> # A tibble: 11 × 7
#> stratum term df sumsq meansq statistic p.value
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Replicate Repl… 2 0.732 0.366 NA NA
#> 2 Replicate:nitrogen nitr… 4 61.6 15.4 27.7 9.73e- 5
#> 3 Replicate:nitrogen Resi… 8 4.45 0.556 NA NA
#> 4 Replicate:nitrogen:management mana… 2 42.9 21.5 82.0 2.30e-10
#> 5 Replicate:nitrogen:management nitr… 8 1.10 0.138 0.527 8.23e- 1
#> 6 Replicate:nitrogen:management Resi… 20 5.24 0.262 NA NA
#> 7 Within vari… 2 206. 103. 208. 1.06e-27
#> 8 Within nitr… 8 14.1 1.77 3.57 1.92e- 3
#> 9 Within mana… 4 3.85 0.963 1.94 1.15e- 1
#> 10 Within nitr… 16 3.70 0.231 0.467 9.54e- 1
#> 11 Within Resi… 60 29.7 0.496 NA NA
##----fm1.2----
fm1.2 <- lmer(yield ~ Replicate + nitrogen * management * variety +
(1|Replicate:Main) + (1|Replicate:Main:Sub), data = rice)
#> boundary (singular) fit: see help('isSingular')
fm1.2.ANOVA <- anova(fm1.2, ddf = "Kenward-Roger", type = 1)
fm1.2.ANOVA
#> Type I Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> Replicate 0.575 0.288 2 8 0.6578 0.5439096
#> nitrogen 48.424 12.106 4 8 27.6953 9.734e-05 ***
#> management 42.936 21.468 2 20 49.1136 1.919e-08 ***
#> variety 206.013 103.007 2 60 235.6535 < 2.2e-16 ***
#> nitrogen:management 1.103 0.138 8 20 0.3154 0.9508797
#> nitrogen:variety 14.145 1.768 8 60 4.0449 0.0006765 ***
#> management:variety 3.852 0.963 4 60 2.2030 0.0793657 .
#> nitrogen:management:variety 3.699 0.231 16 60 0.5289 0.9210465
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##----fm1.3----
fm1.3 <- lmer(yield ~ Replicate + nitrogen + management + variety + nitrogen:variety +
(1|Replicate:Main) + (1|Replicate:Main:Sub), data = rice)
#> boundary (singular) fit: see help('isSingular')
fm1.3.ANOVA <- anova(fm1.3, ddf = "Kenward-Roger", type = 1)
emmeans::emmeans(fm1.3, ~ nitrogen)
#> NOTE: Results may be misleading due to involvement in interactions
#> nitrogen emmean SE df lower.CL upper.CL
#> 0 5.38 0.144 8 5.05 5.72
#> 50 6.22 0.144 8 5.89 6.55
#> 80 7.00 0.144 8 6.66 7.33
#> 110 6.94 0.144 8 6.61 7.27
#> 140 7.23 0.144 8 6.90 7.56
#>
#> Results are averaged over the levels of: Replicate, management, variety
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
emmeans::emmeans(fm1.3, ~ variety)
#> NOTE: Results may be misleading due to involvement in interactions
#> variety emmean SE df lower.CL upper.CL
#> V1 5.13 0.101 39.7 4.92 5.33
#> V2 6.40 0.101 39.7 6.19 6.60
#> V3 8.14 0.101 39.7 7.94 8.34
#>
#> Results are averaged over the levels of: Replicate, nitrogen, management
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
emmeans::emmeans(fm1.3, ~ nitrogen * variety)
#> nitrogen variety emmean SE df lower.CL upper.CL
#> 0 V1 4.51 0.225 39.7 4.06 4.97
#> 50 V1 4.76 0.225 39.7 4.31 5.22
#> 80 V1 5.83 0.225 39.7 5.38 6.29
#> 110 V1 5.44 0.225 39.7 4.99 5.90
#> 140 V1 5.08 0.225 39.7 4.62 5.53
#> 0 V2 5.16 0.225 39.7 4.71 5.62
#> 50 V2 6.02 0.225 39.7 5.56 6.47
#> 80 V2 6.59 0.225 39.7 6.13 7.04
#> 110 V2 6.92 0.225 39.7 6.47 7.38
#> 140 V2 7.29 0.225 39.7 6.83 7.74
#> 0 V3 6.48 0.225 39.7 6.02 6.93
#> 50 V3 7.88 0.225 39.7 7.43 8.34
#> 80 V3 8.56 0.225 39.7 8.11 9.02
#> 110 V3 8.44 0.225 39.7 7.99 8.90
#> 140 V3 9.34 0.225 39.7 8.88 9.79
#>
#> Results are averaged over the levels of: Replicate, management
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
emmeans::contrast(
emmeans::emmeans(fm1.3, ~ nitrogen|variety)
, alpha = 0.05
, method = "pairwise"
)
#> variety = V1:
#> contrast estimate SE df t.ratio p.value
#> nitrogen0 - nitrogen50 -0.251 0.318 39.7 -0.788 0.9326
#> nitrogen0 - nitrogen80 -1.322 0.318 39.7 -4.158 0.0015
#> nitrogen0 - nitrogen110 -0.932 0.318 39.7 -2.930 0.0420
#> nitrogen0 - nitrogen140 -0.565 0.318 39.7 -1.776 0.4012
#> nitrogen50 - nitrogen80 -1.071 0.318 39.7 -3.370 0.0138
#> nitrogen50 - nitrogen110 -0.681 0.318 39.7 -2.142 0.2230
#> nitrogen50 - nitrogen140 -0.314 0.318 39.7 -0.988 0.8591
#> nitrogen80 - nitrogen110 0.390 0.318 39.7 1.228 0.7356
#> nitrogen80 - nitrogen140 0.757 0.318 39.7 2.382 0.1416
#> nitrogen110 - nitrogen140 0.367 0.318 39.7 1.154 0.7768
#>
#> variety = V2:
#> contrast estimate SE df t.ratio p.value
#> nitrogen0 - nitrogen50 -0.853 0.318 39.7 -2.684 0.0744
#> nitrogen0 - nitrogen80 -1.426 0.318 39.7 -4.485 0.0006
#> nitrogen0 - nitrogen110 -1.762 0.318 39.7 -5.542 <.0001
#> nitrogen0 - nitrogen140 -2.126 0.318 39.7 -6.686 <.0001
#> nitrogen50 - nitrogen80 -0.572 0.318 39.7 -1.800 0.3877
#> nitrogen50 - nitrogen110 -0.908 0.318 39.7 -2.857 0.0500
#> nitrogen50 - nitrogen140 -1.272 0.318 39.7 -4.002 0.0023
#> nitrogen80 - nitrogen110 -0.336 0.318 39.7 -1.057 0.8270
#> nitrogen80 - nitrogen140 -0.700 0.318 39.7 -2.202 0.2002
#> nitrogen110 - nitrogen140 -0.364 0.318 39.7 -1.145 0.7819
#>
#> variety = V3:
#> contrast estimate SE df t.ratio p.value
#> nitrogen0 - nitrogen50 -1.403 0.318 39.7 -4.413 0.0007
#> nitrogen0 - nitrogen80 -2.086 0.318 39.7 -6.561 <.0001
#> nitrogen0 - nitrogen110 -1.965 0.318 39.7 -6.181 <.0001
#> nitrogen0 - nitrogen140 -2.857 0.318 39.7 -8.989 <.0001
#> nitrogen50 - nitrogen80 -0.683 0.318 39.7 -2.147 0.2209
#> nitrogen50 - nitrogen110 -0.562 0.318 39.7 -1.767 0.4064
#> nitrogen50 - nitrogen140 -1.454 0.318 39.7 -4.575 0.0004
#> nitrogen80 - nitrogen110 0.121 0.318 39.7 0.380 0.9954
#> nitrogen80 - nitrogen140 -0.772 0.318 39.7 -2.428 0.1290
#> nitrogen110 - nitrogen140 -0.893 0.318 39.7 -2.808 0.0561
#>
#> Results are averaged over the levels of: Replicate, management
#> Degrees-of-freedom method: kenward-roger
#> P value adjustment: tukey method for comparing a family of 5 estimates
emmeans::contrast(
emmeans::emmeans(fm1.3, ~ variety|nitrogen)
, alpha = 0.05
, method = "pairwise"
)
#> nitrogen = 0:
#> contrast estimate SE df t.ratio p.value
#> V1 - V2 -0.650 0.3 80 -2.169 0.0828
#> V1 - V3 -1.965 0.3 80 -6.559 <.0001
#> V2 - V3 -1.315 0.3 80 -4.390 0.0001
#>
#> nitrogen = 50:
#> contrast estimate SE df t.ratio p.value
#> V1 - V2 -1.253 0.3 80 -4.181 0.0002
#> V1 - V3 -3.117 0.3 80 -10.405 <.0001
#> V2 - V3 -1.865 0.3 80 -6.225 <.0001
#>
#> nitrogen = 80:
#> contrast estimate SE df t.ratio p.value
#> V1 - V2 -0.754 0.3 80 -2.516 0.0366
#> V1 - V3 -2.729 0.3 80 -9.109 <.0001
#> V2 - V3 -1.975 0.3 80 -6.593 <.0001
#>
#> nitrogen = 110:
#> contrast estimate SE df t.ratio p.value
#> V1 - V2 -1.480 0.3 80 -4.940 <.0001
#> V1 - V3 -2.998 0.3 80 -10.007 <.0001
#> V2 - V3 -1.518 0.3 80 -5.068 <.0001
#>
#> nitrogen = 140:
#> contrast estimate SE df t.ratio p.value
#> V1 - V2 -2.211 0.3 80 -7.379 <.0001
#> V1 - V3 -4.258 0.3 80 -14.212 <.0001
#> V2 - V3 -2.047 0.3 80 -6.833 <.0001
#>
#> Results are averaged over the levels of: Replicate, management
#> Degrees-of-freedom method: kenward-roger
#> P value adjustment: tukey method for comparing a family of 3 estimates
##----fm1.3.Plot----
fm1.3.Augment <- broom.mixed::augment(fm1.3)
ggplot(data = fm1.3.Augment, mapping = aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0) +
labs(
x = "Fitted"
, y = "Residuals"
, title = "Full analysis with full nitrogen effects") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
##----fm1.4----
fm1.4 <- lmer(yield ~ Replicate + nitrogen * management * variety + (1|Replicate:Main) +
(1|Replicate:Main:Sub), data = rice)
#> boundary (singular) fit: see help('isSingular')
fm1.4.ANOVA <- anova(fm1.4, ddf = "Kenward-Roger", type = 1)
fm1.4.ANOVA
#> Type I Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> Replicate 0.575 0.288 2 8 0.6578 0.5439096
#> nitrogen 48.424 12.106 4 8 27.6953 9.734e-05 ***
#> management 42.936 21.468 2 20 49.1136 1.919e-08 ***
#> variety 206.013 103.007 2 60 235.6535 < 2.2e-16 ***
#> nitrogen:management 1.103 0.138 8 20 0.3154 0.9508797
#> nitrogen:variety 14.145 1.768 8 60 4.0449 0.0006765 ***
#> management:variety 3.852 0.963 4 60 2.2030 0.0793657 .
#> nitrogen:management:variety 3.699 0.231 16 60 0.5289 0.9210465
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##----fm1.5----
fm1.5 <- lmer(yield ~ Replicate + management + variety * (nrate + I(nrate^2) +
I(nrate^3) + I(nrate^4)) +
(1|Replicate:Main) + (1|Replicate:Main:Sub), data = rice)
#> Warning: Some predictor variables are on very different scales: consider rescaling
#> boundary (singular) fit: see help('isSingular')
#> Warning: Some predictor variables are on very different scales: consider rescaling
fm1.5.ANOVA <- anova(fm1.5, ddf = "Kenward-Roger", type = 1)
fm1.5.ANOVA
#> Type I Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> Replicate 0.531 0.266 2 8 0.6578 0.54391
#> management 42.936 21.468 2 28 53.1500 2.932e-10 ***
#> variety 206.013 103.007 2 80 255.0206 < 2.2e-16 ***
#> nrate 40.624 40.624 1 8 100.5764 8.308e-06 ***
#> I(nrate^2) 2.490 2.490 1 8 6.1642 0.03795 *
#> I(nrate^3) 0.038 0.038 1 8 0.0939 0.76713
#> I(nrate^4) 1.594 1.594 1 8 3.9469 0.08220 .
#> variety:nrate 9.861 4.930 2 80 12.2065 2.363e-05 ***
#> variety:I(nrate^2) 0.804 0.402 2 80 0.9952 0.37417
#> variety:I(nrate^3) 2.783 1.392 2 80 3.4455 0.03669 *
#> variety:I(nrate^4) 0.696 0.348 2 80 0.8620 0.42618
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##----fm1.6----
fm1.6 <- lmer(yield ~ Replicate + management + variety * nrate + I(nrate^2) +
(1|Replicate:Main) + (1|Replicate:Main:Sub), data = rice)
#> Warning: Some predictor variables are on very different scales: consider rescaling
#> boundary (singular) fit: see help('isSingular')
#> Warning: Some predictor variables are on very different scales: consider rescaling
fm1.6.Coef <- summary(fm1.6, ddf = "Kenward-Roger")$coef
fm1.6.Coef
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 3.838742e+00 2.478500e-01 31.46271 15.4881675 2.847131e-16
#> ReplicateR2 1.179333e-01 1.725584e-01 10.00000 0.6834401 5.098535e-01
#> ReplicateR3 -5.922222e-02 1.725584e-01 10.00000 -0.3432010 7.385530e-01
#> managementOptimum 5.857778e-01 1.366639e-01 28.00000 4.2862662 1.941972e-04
#> managementIntensive 1.376333e+00 1.366639e-01 28.00000 10.0709370 8.232324e-11
#> varietyV2 5.374285e-01 2.543624e-01 86.00000 2.1128455 3.751250e-02
#> varietyV3 2.008611e+00 2.543624e-01 86.00000 7.8966494 8.555670e-12
#> nrate 1.612922e-02 5.108418e-03 12.37614 3.1573804 7.981949e-03
#> I(nrate^2) -7.528912e-05 3.327527e-05 10.00000 -2.2626149 4.715858e-02
#> varietyV2:nrate 9.630034e-03 2.822766e-03 86.00000 3.4115589 9.862976e-04
#> varietyV3:nrate 1.322179e-02 2.822766e-03 86.00000 4.6839816 1.043415e-05
##----fm1.6.Coefs----
# fm1.6.Coef[ ,1, drop = FALSE]
# Intercepts
fm1.6.Coef[1, 1] + sum(fm1.6.Coef[2:3, 1])/3 + sum(fm1.6.Coef[4:5, 1])/3
#> [1] 4.512349
fm1.6.Coef[1, 1] + sum(fm1.6.Coef[2:3, 1])/3 + sum(fm1.6.Coef[4:5, 1])/3 +
fm1.6.Coef[6, 1]
#> [1] 5.049778
fm1.6.Coef[1, 1] + sum(fm1.6.Coef[2:3, 1])/3 + sum(fm1.6.Coef[4:5, 1])/3 +
fm1.6.Coef[7, 1]
#> [1] 6.52096
# Linear Slopes
fm1.6.Coef[8, 1]
#> [1] 0.01612922
fm1.6.Coef[8, 1] + fm1.6.Coef[10, 1]
#> [1] 0.02575925
fm1.6.Coef[8, 1] + fm1.6.Coef[11, 1]
#> [1] 0.02935101
# Quadratic Slopes
fm1.6.Coef[9, 1]
#> [1] -7.528912e-05
##----fm1.7----
fm1.7 <- aov(yield ~ Replicate*management * variety * nitrogen, rice)
fm1.7.Summary <- broom::tidy(fm1.7)
fm1.7.Summary
#> # A tibble: 15 × 4
#> term df sumsq meansq
#> <chr> <dbl> <dbl> <dbl>
#> 1 Replicate 2 0.732 0.366
#> 2 management 2 42.9 21.5
#> 3 variety 2 206. 103.
#> 4 nitrogen 4 61.6 15.4
#> 5 Replicate:management 4 0.460 0.115
#> 6 Replicate:variety 4 6.15 1.54
#> 7 management:variety 4 3.85 0.963
#> 8 Replicate:nitrogen 8 4.45 0.556
#> 9 management:nitrogen 8 1.10 0.138
#> 10 variety:nitrogen 8 14.1 1.77
#> 11 Replicate:management:variety 8 2.22 0.278
#> 12 Replicate:management:nitrogen 16 4.78 0.299
#> 13 Replicate:variety:nitrogen 16 13.1 0.820
#> 14 management:variety:nitrogen 16 3.70 0.231
#> 15 Replicate:management:variety:nitrogen 32 8.23 0.257
##----fm1.7.Rice----
Rice1 <-
rice %>%
dplyr::group_by(Replicate, nitrogen, variety) %>%
dplyr::summarise(Yield = mean(yield, na.rm = TRUE))
#> `summarise()` has grouped output by 'Replicate', 'nitrogen'. You can override
#> using the `.groups` argument.
Rice1
#> # A tibble: 45 × 4
#> # Groups: Replicate, nitrogen [15]
#> Replicate nitrogen variety Yield
#> <fct> <fct> <fct> <dbl>
#> 1 R1 0 V1 3.92
#> 2 R1 0 V2 5.92
#> 3 R1 0 V3 6.60
#> 4 R1 50 V1 4.02
#> 5 R1 50 V2 6.32
#> 6 R1 50 V3 7.93
#> 7 R1 80 V1 5.81
#> 8 R1 80 V2 6.26
#> 9 R1 80 V3 8.74
#> 10 R1 110 V1 5.44
#> # ℹ 35 more rows
WideRice1 <-
Rice1 %>%
tidyr::spread(key = nitrogen, value = Yield) %>%
dplyr::ungroup() %>%
dplyr::select(-Replicate, -variety)
##----fm1.7.Plot1----
ggplot(data = Rice1, mapping = aes(x = nitrogen, y = Yield, group = Replicate)) +
geom_line() +
facet_grid(variety ~ Replicate, labeller = label_both) +
labs(
x = "Nitrogen"
, y = "Yield"
, title = "Fig S1. Variety response to nitrogen for individual replicate blocks"
) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))
##----fm1.8----
riceV2V3 <-
rice %>%
dplyr::filter(variety != "V1") %>%
droplevels()
fm1.8 <- aov(yield ~ Replicate*management * variety * nitrogen, riceV2V3)
fm1.8.ANOVA <- broom::tidy(fm1.8)
fm1.8.ANOVA
#> # A tibble: 15 × 4
#> term df sumsq meansq
#> <chr> <dbl> <dbl> <dbl>
#> 1 Replicate 2 3.00 1.50
#> 2 management 2 26.1 13.1
#> 3 variety 1 68.4 68.4
#> 4 nitrogen 4 64.0 16.0
#> 5 Replicate:management 4 0.868 0.217
#> 6 Replicate:variety 2 0.518 0.259
#> 7 management:variety 2 3.68 1.84
#> 8 Replicate:nitrogen 8 6.98 0.872
#> 9 management:nitrogen 8 0.842 0.105
#> 10 variety:nitrogen 4 1.78 0.444
#> 11 Replicate:management:variety 4 1.24 0.311
#> 12 Replicate:management:nitrogen 16 5.29 0.331
#> 13 Replicate:variety:nitrogen 8 4.43 0.554
#> 14 management:variety:nitrogen 8 3.22 0.402
#> 15 Replicate:management:variety:nitrogen 16 4.03 0.252
##----fm1.9----
fm1.9 <- aov(yield ~ Replicate + management * variety * nitrogen +
Error(Replicate/Main/Sub), riceV2V3)
fm1.9.ANOVA <- broom::tidy(fm1.9)
fm1.9.ANOVA
#> # A tibble: 11 × 7
#> stratum term df sumsq meansq statistic p.value
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Replicate Replicate 2 3.00 1.50 NA NA
#> 2 Replicate:Main nitrogen 4 64.0 16.0 18.3 4.30e- 4
#> 3 Replicate:Main Residuals 8 6.98 0.872 NA NA
#> 4 Replicate:Main:Sub management 2 26.1 13.1 42.5 6.34e- 8
#> 5 Replicate:Main:Sub management:nitrog… 8 0.842 0.105 0.342 9.39e- 1
#> 6 Replicate:Main:Sub Residuals 20 6.16 0.308 NA NA
#> 7 Within variety 1 68.4 68.4 201. 7.87e-15
#> 8 Within management:variety 2 3.68 1.84 5.40 9.92e- 3
#> 9 Within variety:nitrogen 4 1.78 0.444 1.30 2.91e- 1
#> 10 Within management:variet… 8 3.22 0.402 1.18 3.43e- 1
#> 11 Within Residuals 30 10.2 0.341 NA NA