Day 2 Workshop
Gasatan, EB
2024-10-10
Loading Libraries
library(readxl)
library(ExpDes)
library(tidyverse)
library(ggstatsplot)
library(ISLR)
library(patchwork)
library(emmeans)
library(multcomp)
library(multcompView)
library(ggpubr)
library(rstatix)
library(Hmisc)
library(GGally)##Analysis of Experiments in RCBD
#Importing the data (Make sure to set the appropriate working directory)
agrondata <- read_excel("agrondata.xlsx")
head(agrondata) #displays the 1st 6 rows of the data frame## # A tibble: 6 × 10
## treatment block DTE HT15 HT30 HT45 HT60 DTF DTM HERBAGE
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 4 56.2 102. 152. 210. 44 60 158000
## 2 1 2 4 53.9 92.9 140 194. 44 60 192000
## 3 1 3 4 53.4 90.4 134. 189. 44 60 180000
## 4 1 4 4 52.2 89 129. 183. 44 60 155000
## 5 1 5 4 54.6 102. 150. 207. 44 60 185000
## 6 2 1 7 47.3 73.5 131. 173. 76 90 255000#ANOVA for RCBD experiment with post hoc
with(agrondata, rbd(treat = treatment,
block = block,
resp = HT15,
quali = TRUE,
mcomp = "tukey"))## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Treatament 2 2416.67 1208.33 232.588 0.000000
## Block 4 74.27 18.57 3.574 0.059085
## Residuals 8 41.56 5.20
## Total 14 2532.50
## ------------------------------------------------------------------------
## CV = 5.7 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.1688107
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## ------------------------------------------------------------------------
## Homogeneity of variances test
## p-value: 0.07919401
## According to the test of oneillmathews at 5% of significance, the variances can be considered homocedastic.
## ------------------------------------------------------------------------
##
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 1 54.06
## b 2 42.6
## c 3 23.3
## ------------------------------------------------------------------------#Presenting the results visually
#Run ANOVA the traditional way
agrondata$Trt <- as.factor(agrondata$treatment)
agrondata$Blk <- as.factor(agrondata$block)
head(agrondata)## # A tibble: 6 × 12
## treatment block DTE HT15 HT30 HT45 HT60 DTF DTM HERBAGE Trt Blk
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct>
## 1 1 1 4 56.2 102. 152. 210. 44 60 158000 1 1
## 2 1 2 4 53.9 92.9 140 194. 44 60 192000 1 2
## 3 1 3 4 53.4 90.4 134. 189. 44 60 180000 1 3
## 4 1 4 4 52.2 89 129. 183. 44 60 155000 1 4
## 5 1 5 4 54.6 102. 150. 207. 44 60 185000 1 5
## 6 2 1 7 47.3 73.5 131. 173. 76 90 255000 2 1out1 <- aov(HT15 ~ Trt + Blk,
data = agrondata)
anova(out1)## Analysis of Variance Table
##
## Response: HT15
## Df Sum Sq Mean Sq F value Pr(>F)
## Trt 2 2416.67 1208.33 232.588 8.171e-08 ***
## Blk 4 74.27 18.57 3.574 0.05908 .
## Residuals 8 41.56 5.20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TRT.means <- agrondata |>
group_by(treatment) |>
dplyr::summarize(Mean=mean(HT15),
SD = sd(HT15),
n = length(HT15)) |>
dplyr::mutate(SE = SD/sqrt(n))
TRT.means## # A tibble: 3 × 5
## treatment Mean SD n SE
## <dbl> <dbl> <dbl> <int> <dbl>
## 1 1 54.1 1.48 5 0.663
## 2 2 42.6 4.90 5 2.19
## 3 3 23.3 1.67 5 0.746#Compute summary statistics by Stocks
trt.means <- emmeans(out1, specs="Trt")
trt.means## Trt emmean SE df lower.CL upper.CL
## 1 54.1 1.02 8 51.7 56.4
## 2 42.6 1.02 8 40.2 45.0
## 3 23.3 1.02 8 20.9 25.7
##
## Results are averaged over the levels of: Blk
## Confidence level used: 0.95#Add the letters
trt.means.cld <- cld(trt.means,
Letters=letters,
decreasing = TRUE)
trt.means.cld <- trt.means.cld |>
arrange(Trt) |>
mutate(SE1=TRT.means$SE)
trt.means.cld## Trt emmean SE df lower.CL upper.CL .group SE1
## 1 1 54.06 1.01933 8 51.70942 56.41058 a 0.6630234
## 2 2 42.60 1.01933 8 40.24942 44.95058 b 2.1897488
## 3 3 23.30 1.01933 8 20.94942 25.65058 c 0.7463243#Create graphs with SE and letters
ggplot(trt.means.cld, aes(Trt, emmean,label = .group)) +
geom_col(fill = "lightblue", col="black") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.2,
size = 0.7) +
geom_text(vjust=-1.5, size = 5) +
geom_text(vjust = 3, label = round(trt.means.cld$emmean,2)) +
labs(x= "Treatment", y= "Mean Plant Height (15 days)") +
lims(y = c(0,60)) +
theme_classic()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
##Analysis of Factorial Experiments in CRD
#Importing the data
fact_crd <- read_excel("Factorial.xlsx",
sheet = "CRD")
#Display the first 6 rows of the data frame
head(fact_crd)## # A tibble: 6 × 4
## A B Rep Height
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 46.5
## 2 1 1 2 55.9
## 3 1 1 3 78.7
## 4 1 2 1 49.5
## 5 1 2 2 59.5
## 6 1 2 3 78.7#ANOVA Table
with(fact_crd,fat2.crd(factor1 = A,
factor2 = B,
resp = Height,
quali = c(TRUE,TRUE),
mcomp = "tukey",
fac.names = c("Spacing", "Age"),
))## ------------------------------------------------------------------------
## Legend:
## FACTOR 1: Spacing
## FACTOR 2: Age
## ------------------------------------------------------------------------
##
##
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Spacing 1 409.0 3 1.5207 0.241120
## Age 2 12846.3 5 23.8835 0.000066
## Spacing*Age 2 996.6 4 1.8529 0.198942
## Residuals 12 3227.2 2
## Total 17 17479.1 1
## ------------------------------------------------------------------------
## CV = 20.36 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.4529232
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## No significant interaction: analyzing the simple effect
## ------------------------------------------------------------------------
## Spacing
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 85.30000
## 2 2 75.76667
## ------------------------------------------------------------------------
## Age
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 118.0833
## b 2 65.36667
## b 1 58.15
## ------------------------------------------------------------------------#Graphical presentation
#Create the ANOVA model
fact_crd$Spacing = as.factor(fact_crd$A)
fact_crd$Age = as.factor(fact_crd$B)
head(fact_crd)## # A tibble: 6 × 6
## A B Rep Height Spacing Age
## <dbl> <dbl> <dbl> <dbl> <fct> <fct>
## 1 1 1 1 46.5 1 1
## 2 1 1 2 55.9 1 1
## 3 1 1 3 78.7 1 1
## 4 1 2 1 49.5 1 2
## 5 1 2 2 59.5 1 2
## 6 1 2 3 78.7 1 2#Create the ANOVA model
fact_mod1 <- aov(Height ~ Spacing + Age + Spacing:Age,
data = fact_crd)
#Display the ANOVA table
summary(fact_mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## Spacing 1 409 409 1.521 0.241
## Age 2 12846 6423 23.883 6.55e-05 ***
## Spacing:Age 2 997 498 1.853 0.199
## Residuals 12 3227 269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fact_crd <- fact_crd |>
group_by(Age) |>
dplyr::summarize(Mean = mean(Height),
SD = sd(Height),
n = length(Height)) |>
dplyr::mutate(SE = SD/sqrt(n))
fact_crd## # A tibble: 3 × 5
## Age Mean SD n SE
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 58.2 12.0 6 4.89
## 2 2 65.4 10.4 6 4.24
## 3 3 118. 26.0 6 10.6#Generate the summary statistics for the levels of Factor B
B.means <- emmeans(fact_mod1, specs="Age")## NOTE: Results may be misleading due to involvement in interactionsB.means#Print the table of means## Age emmean SE df lower.CL upper.CL
## 1 58.1 6.69 12 43.6 72.7
## 2 65.4 6.69 12 50.8 80.0
## 3 118.1 6.69 12 103.5 132.7
##
## Results are averaged over the levels of: Spacing
## Confidence level used: 0.95#Assign the letter designations to the means of the Factor B
B.means.cld <- cld(B.means,
Letters=letters) |>
arrange(Age) |>
mutate(SE1 = fact_crd$SE)
B.means.cld #Print the table of means## Age emmean SE df lower.CL upper.CL .group SE1
## 1 1 58.15000 6.694975 12 43.56290 72.73710 a 4.888882
## 2 2 65.36667 6.694975 12 50.77957 79.95376 a 4.237269
## 3 3 118.08333 6.694975 12 103.49624 132.67043 b 10.610008#Generate the two-way plot
ggplot(B.means.cld, aes(Age, emmean,label = .group)) +
geom_col(fill="lightblue") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.2, size = 0.7) +
geom_text(vjust=-3.2, size = 5) +
geom_text(label=round(B.means.cld$emmean,2),vjust=5, color = "black") +
lims(y = c(0,160)) +
labs(x= "Age", y= "Mean Height") +
theme_classic()
##Analysis of Factorial Experiments in RCBD
#Importing the data into R
fact_rbd <- read_excel("Factorial.xlsx",
sheet = "RCBD")
head(fact_rbd)## # A tibble: 6 × 4
## Block P N Yield
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 0.9
## 2 1 1 2 1.2
## 3 1 1 3 1.3
## 4 1 1 4 1.8
## 5 1 1 5 1.1
## 6 1 2 1 0.9#Factorial ANOVA (RCBD)
with(fact_rbd,fat2.rbd(factor1 = P,
factor2 = N,
block = Block,
resp = Yield,
quali = c(TRUE,TRUE),
fac.names = c("Phosporous","Nitrogen"),
mcomp="tukey"))## ------------------------------------------------------------------------
## Legend:
## FACTOR 1: Phosporous
## FACTOR 2: Nitrogen
## ------------------------------------------------------------------------
##
##
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Block 2 0.0640 3 1.631 0.213772
## Phosporous 2 0.1693 4 4.316 0.023244
## Nitrogen 4 2.4902 6 31.732 0.000000
## Phosporous*Nitrogen 8 1.0151 5 6.468 0.000090
## Residuals 28 0.5493 2
## Total 44 4.2880 1
## ------------------------------------------------------------------------
## CV = 11.12 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.1377132
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
##
##
## Significant interaction: analyzing the interaction
## ------------------------------------------------------------------------
##
## Analyzing Phosporous inside of each level of Nitrogen
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr.Fc
## Block 2 0.06400 0.032 1.6311 0.2138
## Nitrogen 4 2.49022 0.62256 31.7322 0
## Phosporous:Nitrogen 1 2 0.01556 0.00778 0.3964 0.6764
## Phosporous:Nitrogen 2 2 0.12667 0.06333 3.2282 0.0548
## Phosporous:Nitrogen 3 2 0.01556 0.00778 0.3964 0.6764
## Phosporous:Nitrogen 4 2 0.62000 0.31 15.801 0
## Phosporous:Nitrogen 5 2 0.40667 0.20333 10.3641 4e-04
## Residuals 28 0.54933 0.01962
## Total 44 4.28800
## ------------------------------------------------------------------------
##
##
##
## Phosporous inside of the level 1 of Nitrogen
##
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 0.9333333
## 2 2 0.8333333
## 3 3 0.8666667
## ------------------------------------------------------------------------
##
##
## Phosporous inside of the level 2 of Nitrogen
##
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 1.2333333
## 2 2 0.9666667
## 3 3 1.2000000
## ------------------------------------------------------------------------
##
##
## Phosporous inside of the level 3 of Nitrogen
##
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 1.400000
## 2 2 1.300000
## 3 3 1.366667
## ------------------------------------------------------------------------
##
##
## Phosporous inside of the level 4 of Nitrogen
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 1 1.933333
## b 3 1.433333
## b 2 1.333333
## ------------------------------------------------------------------------
##
##
## Phosporous inside of the level 5 of Nitrogen
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 1.666667
## b 1 1.233333
## b 3 1.2
## ------------------------------------------------------------------------
##
##
##
## Analyzing Nitrogen inside of each level of Phosporous
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr.Fc
## Block 2 0.06400 0.032 1.6311 0.2138
## Phosporous 2 0.16933 0.08467 4.3155 0.0232
## Nitrogen:Phosporous 1 4 1.63067 0.40767 20.7791 0
## Nitrogen:Phosporous 2 4 1.29733 0.32433 16.5316 0
## Nitrogen:Phosporous 3 4 0.57733 0.14433 7.3568 4e-04
## Residuals 28 0.54933 0.01962
## Total 44 4.28800
## ------------------------------------------------------------------------
##
##
##
## Nitrogen inside of the level 1 of Phosporous
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 4 1.933333
## b 3 1.4
## bc 2 1.233333
## bc 5 1.233333
## c 1 0.9333333
## ------------------------------------------------------------------------
##
##
## Nitrogen inside of the level 2 of Phosporous
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 5 1.666667
## b 4 1.333333
## b 3 1.3
## c 2 0.9666667
## c 1 0.8333333
## ------------------------------------------------------------------------
##
##
## Nitrogen inside of the level 3 of Phosporous
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 4 1.433333
## a 3 1.366667
## a 2 1.2
## a 5 1.2
## b 1 0.8666667
## ------------------------------------------------------------------------#Visualizing the results
#Converting the integer codes to factors
fact_rbd$Block = factor(fact_rbd$Block)
fact_rbd$P = factor(fact_rbd$P)
fact_rbd$N = factor(fact_rbd$N)
head(fact_rbd)## # A tibble: 6 × 4
## Block P N Yield
## <fct> <fct> <fct> <dbl>
## 1 1 1 1 0.9
## 2 1 1 2 1.2
## 3 1 1 3 1.3
## 4 1 1 4 1.8
## 5 1 1 5 1.1
## 6 1 2 1 0.9#Create the ANOVA model
fact_mod2 <- aov(Yield ~ Block + P + N + P:N,
data = fact_rbd)
#Display the ANOVA table
anova(fact_mod2)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 0.06400 0.03200 1.6311 0.21377
## P 2 0.16933 0.08467 4.3155 0.02324 *
## N 4 2.49022 0.62256 31.7322 4.946e-10 ***
## P:N 8 1.01511 0.12689 6.4676 8.979e-05 ***
## Residuals 28 0.54933 0.01962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pn.means <- fact_rbd |>
group_by(P,N) |>
dplyr::summarise(Mean = mean(Yield),
SD = sd(Yield),
n = length(Yield)) |>
dplyr::mutate(SE = SD/sqrt(n))## `summarise()` has grouped output by 'P'. You can override using the `.groups`
## argument.pn.means## # A tibble: 15 × 6
## # Groups: P [3]
## P N Mean SD n SE
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 1 1 0.933 0.0577 3 0.0333
## 2 1 2 1.23 0.0577 3 0.0333
## 3 1 3 1.4 0.1 3 0.0577
## 4 1 4 1.93 0.153 3 0.0882
## 5 1 5 1.23 0.153 3 0.0882
## 6 2 1 0.833 0.0577 3 0.0333
## 7 2 2 0.967 0.115 3 0.0667
## 8 2 3 1.3 0.2 3 0.115
## 9 2 4 1.33 0.252 3 0.145
## 10 2 5 1.67 0.208 3 0.120
## 11 3 1 0.867 0.153 3 0.0882
## 12 3 2 1.2 0.2 3 0.115
## 13 3 3 1.37 0.0577 3 0.0333
## 14 3 4 1.43 0.0577 3 0.0333
## 15 3 5 1.2 0.1 3 0.0577#Generate summary statistics for all treatment combinations
PN.means <- emmeans(fact_mod2, specs=c("P", "N"))
PN.means#Print the table of means## P N emmean SE df lower.CL upper.CL
## 1 1 0.933 0.0809 28 0.768 1.099
## 2 1 0.833 0.0809 28 0.668 0.999
## 3 1 0.867 0.0809 28 0.701 1.032
## 1 2 1.233 0.0809 28 1.068 1.399
## 2 2 0.967 0.0809 28 0.801 1.132
## 3 2 1.200 0.0809 28 1.034 1.366
## 1 3 1.400 0.0809 28 1.234 1.566
## 2 3 1.300 0.0809 28 1.134 1.466
## 3 3 1.367 0.0809 28 1.201 1.532
## 1 4 1.933 0.0809 28 1.768 2.099
## 2 4 1.333 0.0809 28 1.168 1.499
## 3 4 1.433 0.0809 28 1.268 1.599
## 1 5 1.233 0.0809 28 1.068 1.399
## 2 5 1.667 0.0809 28 1.501 1.832
## 3 5 1.200 0.0809 28 1.034 1.366
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95#Assign the letters to the treatment means
PN.means.cld <- cld(PN.means,
Letters=letters,
decreasing = TRUE)
PN.means.cld## P N emmean SE df lower.CL upper.CL .group
## 1 4 1.933 0.0809 28 1.768 2.099 a
## 2 5 1.667 0.0809 28 1.501 1.832 ab
## 3 4 1.433 0.0809 28 1.268 1.599 bc
## 1 3 1.400 0.0809 28 1.234 1.566 bc
## 3 3 1.367 0.0809 28 1.201 1.532 bcd
## 2 4 1.333 0.0809 28 1.168 1.499 bcde
## 2 3 1.300 0.0809 28 1.134 1.466 bcde
## 1 5 1.233 0.0809 28 1.068 1.399 cdef
## 1 2 1.233 0.0809 28 1.068 1.399 cdef
## 3 5 1.200 0.0809 28 1.034 1.366 cdef
## 3 2 1.200 0.0809 28 1.034 1.366 cdef
## 2 2 0.967 0.0809 28 0.801 1.132 def
## 1 1 0.933 0.0809 28 0.768 1.099 ef
## 3 1 0.867 0.0809 28 0.701 1.032 f
## 2 1 0.833 0.0809 28 0.668 0.999 f
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 15 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.PN.means.cld <- PN.means.cld |>
dplyr::arrange(P,N) |>
dplyr::mutate(SE1 = pn.means$SE)
PN.means.cld## P N emmean SE df lower.CL upper.CL .group SE1
## 1 1 1 0.9333333 0.0808683 28 0.7676821 1.0989845 ef 0.03333333
## 2 1 2 1.2333333 0.0808683 28 1.0676821 1.3989845 cdef 0.03333333
## 3 1 3 1.4000000 0.0808683 28 1.2343488 1.5656512 bc 0.05773503
## 4 1 4 1.9333333 0.0808683 28 1.7676821 2.0989845 a 0.08819171
## 5 1 5 1.2333333 0.0808683 28 1.0676821 1.3989845 cdef 0.08819171
## 6 2 1 0.8333333 0.0808683 28 0.6676821 0.9989845 f 0.03333333
## 7 2 2 0.9666667 0.0808683 28 0.8010155 1.1323179 def 0.06666667
## 8 2 3 1.3000000 0.0808683 28 1.1343488 1.4656512 bcde 0.11547005
## 9 2 4 1.3333333 0.0808683 28 1.1676821 1.4989845 bcde 0.14529663
## 10 2 5 1.6666667 0.0808683 28 1.5010155 1.8323179 ab 0.12018504
## 11 3 1 0.8666667 0.0808683 28 0.7010155 1.0323179 f 0.08819171
## 12 3 2 1.2000000 0.0808683 28 1.0343488 1.3656512 cdef 0.11547005
## 13 3 3 1.3666667 0.0808683 28 1.2010155 1.5323179 bcd 0.03333333
## 14 3 4 1.4333333 0.0808683 28 1.2676821 1.5989845 bc 0.03333333
## 15 3 5 1.2000000 0.0808683 28 1.0343488 1.3656512 cdef 0.05773503#Comparing the levels of N for each level of P
ggplot(PN.means.cld, aes(x = P,
y = emmean,
fill = N,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.3,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-4,
position = position_dodge(0.9),
size = 3) +
labs(x= "P", y= "Mean Yield") +
lims(y=c(0,2.5))+
scale_fill_brewer(palette = "Blues") +
theme_classic()
##Analysis of Split-Plot Experiments in RCBD
#Importing the data
spplot <- read_excel("splitplot.xlsx")
head(spplot)## # A tibble: 6 × 4
## Block Variety Concentration Yield
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 42.9
## 2 1 1 4 44.4
## 3 1 1 3 49.5
## 4 1 1 2 53.8
## 5 1 2 3 59.8
## 6 1 2 1 53.3#ANOVA for split plot experiments
with(spplot, split2.rbd(factor1 = Variety,
factor2 = Concentration,
block = Block,
resp = Yield,
quali = c(TRUE, TRUE),
fac.names = c("Variety", "Concentration"),
mcomp = "tukey",
unfold = "1"))## ------------------------------------------------------------------------
## Legend:
## FACTOR 1 (plot): Variety
## FACTOR 2 (split-plot): Concentration
## ------------------------------------------------------------------------
##
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr(>Fc)
## Variety 2 3631.5 1815.77 60.775 1.5e-05 ***
## Block 4 2931.3 732.82 24.528 0.000152 ***
## Error a 8 239.0 29.88
## Concentration 3 350.2 116.73 8.316 0.000248 ***
## Variety*Concentration 6 368.3 61.38 4.372 0.002049 **
## Error b 36 505.4 14.04
## Total 59 8025.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ------------------------------------------------------------------------
## CV 1 = 10.6146 %
## CV 2 = 7.275885 %
##
## No significant interaction: analyzing the simple effects
## ------------------------------------------------------------------------
## Variety
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 60.255
## b 2 52.88
## c 1 41.35
## ------------------------------------------------------------------------
##
## Concentration
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 55.45333
## b 3 51.22
## b 4 50.29333
## b 1 49.01333
## ------------------------------------------------------------------------
##
##
##
##
## Significant interaction: analyzing the interaction
## ------------------------------------------------------------------------
##
## Analyzing Variety inside of each level of Concentration
## ------------------------------------------------------------------------
## DF SS MS Fc p.value
## Variety : Concentration 1 2.00000 1573.2013 786.60067 43.70568 0.000000
## Variety : Concentration 2 2.00000 495.6253 247.81267 13.76915 0.000048
## Variety : Concentration 3 2.00000 422.6680 211.33400 11.74229 0.000148
## Variety : Concentration 4 2.00000 1508.3253 754.16267 41.90333 0.000000
## Pooled Error 32.22141 579.9106 17.99768 NA NA
## ------------------------------------------------------------------------
##
##
## Variety inside of Concentration 1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 60.84
## b 2 50.34
## c 1 35.86
## ------------------------------------------------------------------------
##
## Variety inside of Concentration 2
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 62.78
## b 2 54.84
## b 1 48.74
## ------------------------------------------------------------------------
##
## Variety inside of Concentration 3
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 56.96
## a 2 52.54
## b 1 44.16
## ------------------------------------------------------------------------
##
## Variety inside of Concentration 4
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 60.44
## b 2 53.8
## c 1 36.64
## ------------------------------------------------------------------------
##
##
## Analyzing Concentration inside of each level of Variety
## ------------------------------------------------------------------------
## DF SS MS Fc p.value
## Concentration : Variety 1 3 574.1620 191.38733 13.633626 0.000004
## Concentration : Variety 2 3 56.2760 18.75867 1.336288 0.277818
## Concentration : Variety 3 3 88.0455 29.34850 2.090663 0.118618
## Error b 36 505.3640 14.03789 NA NA
## ------------------------------------------------------------------------
##
##
## Concentration inside of Variety 1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 48.74
## a 3 44.16
## b 4 36.64
## b 1 35.86
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
##
##
## Concentration inside of Variety 2
## ------------------------------------------------------------------------
## According to F test, the means of this factor are not different.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 50.34
## 2 2 54.84
## 3 3 52.54
## 4 4 53.80
## ------------------------------------------------------------------------
##
## Concentration inside of Variety 3
## ------------------------------------------------------------------------
## According to F test, the means of this factor are not different.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 60.84
## 2 2 62.78
## 3 3 56.96
## 4 4 60.44
## ------------------------------------------------------------------------#Another way to generate ANOVA for split plot experiments
#Converts integer codes to nominal codes
spplot <- spplot |>
dplyr::mutate(Blk = as.factor(Block),
Var = as.factor(Variety),
Conc = as.factor(Concentration))
#Generates ANOVA
spplot.out <- aov(Yield ~ Blk + Var + Blk:Var + Conc + Var:Conc,
data = spplot)
anova(spplot.out)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Blk 4 2931.3 732.82 52.2028 1.691e-14 ***
## Var 2 3631.5 1815.77 129.3477 < 2.2e-16 ***
## Conc 3 350.2 116.73 8.3156 0.0002483 ***
## Blk:Var 8 239.0 29.88 2.1283 0.0583461 .
## Var:Conc 6 368.3 61.38 4.3725 0.0020493 **
## Residuals 36 505.4 14.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Extract the Mean Square for *Blk:Var* and use this as *MSE(a)*
MSE_a <- anova(spplot.out)["Mean Sq"][4,1]
MS_Blk <- anova(spplot.out)["Mean Sq"][1,1]
MS_Var <- anova(spplot.out)["Mean Sq"][2,1]
print(c(MSE_a,MS_Blk, MS_Var))## [1] 29.87704 732.81692 1815.76850#Generate the visualization of the interaction effect
VC.means <- spplot |>
group_by(Var, Conc) |>
dplyr::summarize(Mean = mean(Yield),
SD = sd(Yield),
n = length(Yield)) |>
dplyr::mutate(SE = SD/sqrt(n))## `summarise()` has grouped output by 'Var'. You can override using the `.groups`
## argument.VC.means## # A tibble: 12 × 6
## # Groups: Var [3]
## Var Conc Mean SD n SE
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 1 1 35.9 6.27 5 2.80
## 2 1 2 48.7 7.20 5 3.22
## 3 1 3 44.2 7.10 5 3.18
## 4 1 4 36.6 6.46 5 2.89
## 5 2 1 50.3 12.7 5 5.66
## 6 2 2 54.8 9.92 5 4.44
## 7 2 3 52.5 10.1 5 4.52
## 8 2 4 53.8 7.45 5 3.33
## 9 3 1 60.8 9.57 5 4.28
## 10 3 2 62.8 5.67 5 2.54
## 11 3 3 57.0 9.78 5 4.37
## 12 3 4 60.4 10.0 5 4.47#Generate summary statistics for all treatment combinations
VC.means1 <- emmeans(spplot.out, specs=c("Var", "Conc"))
VC.means1#Print the table of means## Var Conc emmean SE df lower.CL upper.CL
## 1 1 35.9 1.68 36 32.5 39.3
## 2 1 50.3 1.68 36 46.9 53.7
## 3 1 60.8 1.68 36 57.4 64.2
## 1 2 48.7 1.68 36 45.3 52.1
## 2 2 54.8 1.68 36 51.4 58.2
## 3 2 62.8 1.68 36 59.4 66.2
## 1 3 44.2 1.68 36 40.8 47.6
## 2 3 52.5 1.68 36 49.1 55.9
## 3 3 57.0 1.68 36 53.6 60.4
## 1 4 36.6 1.68 36 33.2 40.0
## 2 4 53.8 1.68 36 50.4 57.2
## 3 4 60.4 1.68 36 57.0 63.8
##
## Results are averaged over the levels of: Blk
## Confidence level used: 0.95#Assign the letters to the treatment means
VC.means.cld <- cld(VC.means1,
Letters=letters,
decreasing = TRUE) |>
dplyr::arrange(Var, Conc) |>
dplyr::mutate(SE1 = VC.means$SE)
VC.means.cld#Print the table of means## Var Conc emmean SE df lower.CL upper.CL .group SE1
## 1 1 1 35.86 1.675583 36 32.46176 39.25824 g 2.802963
## 2 1 2 48.74 1.675583 36 45.34176 52.13824 de 3.217856
## 3 1 3 44.16 1.675583 36 40.76176 47.55824 ef 3.176256
## 4 1 4 36.64 1.675583 36 33.24176 40.03824 fg 2.890778
## 5 2 1 50.34 1.675583 36 46.94176 53.73824 de 5.659205
## 6 2 2 54.84 1.675583 36 51.44176 58.23824 abcd 4.435380
## 7 2 3 52.54 1.675583 36 49.14176 55.93824 cd 4.520133
## 8 2 4 53.80 1.675583 36 50.40176 57.19824 bcd 3.331516
## 9 3 1 60.84 1.675583 36 57.44176 64.23824 ab 4.280958
## 10 3 2 62.78 1.675583 36 59.38176 66.17824 a 2.537597
## 11 3 3 56.96 1.675583 36 53.56176 60.35824 abcd 4.373168
## 12 3 4 60.44 1.675583 36 57.04176 63.83824 abc 4.472762#Comparing the levels of Concentration for each level of Variety
ggplot(VC.means.cld, aes(x = Var,
y = emmean,
fill = Conc,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-3.2,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "Variety", y= "Mean Yield") +
lims(y = c(0, 75)) +
scale_fill_brewer(palette = "Blues") +
theme_classic()
#Comparing the levels of Variety for each level of Concentration
ggplot(VC.means.cld, aes(x = Conc,
y = emmean,
fill = Var,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-3.5,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "Concentration", y= "Mean Yield") +
lims(y = c(0, 75)) +
scale_fill_brewer(palette = "Blues") +
theme_classic()
##Analysis of Experiments with Repeated Measurements
head(agrondata)## # A tibble: 6 × 12
## treatment block DTE HT15 HT30 HT45 HT60 DTF DTM HERBAGE Trt Blk
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct>
## 1 1 1 4 56.2 102. 152. 210. 44 60 158000 1 1
## 2 1 2 4 53.9 92.9 140 194. 44 60 192000 1 2
## 3 1 3 4 53.4 90.4 134. 189. 44 60 180000 1 3
## 4 1 4 4 52.2 89 129. 183. 44 60 155000 1 4
## 5 1 5 4 54.6 102. 150. 207. 44 60 185000 1 5
## 6 2 1 7 47.3 73.5 131. 173. 76 90 255000 2 1#Subset the agrondata to include only height measurements
HT <- agrondata |>
select(treatment, block, HT15, HT30, HT45, HT60)
HT## # A tibble: 15 × 6
## treatment block HT15 HT30 HT45 HT60
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 56.2 102. 152. 210.
## 2 1 2 53.9 92.9 140 194.
## 3 1 3 53.4 90.4 134. 189.
## 4 1 4 52.2 89 129. 183.
## 5 1 5 54.6 102. 150. 207.
## 6 2 1 47.3 73.5 131. 173.
## 7 2 2 44.4 72.4 119. 159.
## 8 2 3 36.6 70.1 118. 158.
## 9 2 4 38.2 68.8 119. 148.
## 10 2 5 46.5 74.8 133. 174.
## 11 3 1 24.4 38.7 50.4 74.2
## 12 3 2 21.9 36.6 45.9 69.5
## 13 3 3 22.2 37.3 46.7 70
## 14 3 4 22.3 36 45.6 71.8
## 15 3 5 25.7 39.1 51.4 77.4HT_Long <- HT |>
gather(key = "Time", value = "HT", HT15, HT30, HT45, HT60) |>
convert_as_factor(treatment, block, Time) |>
dplyr::rename(Treatment = treatment) |>
dplyr::mutate(Days = recode(Time,
"HT15" = "15",
"HT30" = "30",
"HT45" = "45",
"HT60" = "60"))
HT_Long## # A tibble: 60 × 5
## Treatment block Time HT Days
## <fct> <fct> <fct> <dbl> <fct>
## 1 1 1 HT15 56.2 15
## 2 1 2 HT15 53.9 15
## 3 1 3 HT15 53.4 15
## 4 1 4 HT15 52.2 15
## 5 1 5 HT15 54.6 15
## 6 2 1 HT15 47.3 15
## 7 2 2 HT15 44.4 15
## 8 2 3 HT15 36.6 15
## 9 2 4 HT15 38.2 15
## 10 2 5 HT15 46.5 15
## # ℹ 50 more rows#Testing for presence of outliers
HT_Long |>
group_by(Treatment, Days) |>
identify_outliers(HT)## [1] Treatment Days block Time HT is.outlier is.extreme
## <0 rows> (or 0-length row.names)HT_Long |>
group_by(Treatment, Days) |>
shapiro_test(HT)## # A tibble: 12 × 5
## Treatment Days variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 1 15 HT 0.989 0.975
## 2 1 30 HT 0.837 0.157
## 3 1 45 HT 0.930 0.595
## 4 1 60 HT 0.911 0.476
## 5 2 15 HT 0.866 0.252
## 6 2 30 HT 0.961 0.813
## 7 2 45 HT 0.762 0.0384
## 8 2 60 HT 0.912 0.478
## 9 3 15 HT 0.839 0.161
## 10 3 30 HT 0.929 0.592
## 11 3 45 HT 0.837 0.156
## 12 3 60 HT 0.921 0.536rep.aov <- anova_test(
data = HT_Long,
dv = HT,
wid = block,
within = c(Treatment, Days))
#To display the ANOVA table
get_anova_table(rep.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Treatment 1.04 4.15 915.983 4.91e-06 * 0.969
## 2 Days 1.02 4.09 1289.421 2.80e-06 * 0.979
## 3 Treatment:Days 6.00 24.00 297.420 2.65e-21 * 0.877#To display the results of sphericity test
rep.aov$`Mauchly's Test for Sphericity` ## Effect W p p<.05
## 1 Treatment 0.071000 0.019 *
## 2 Days 0.000758 0.003 *head(HT_Long)## # A tibble: 6 × 5
## Treatment block Time HT Days
## <fct> <fct> <fct> <dbl> <fct>
## 1 1 1 HT15 56.2 15
## 2 1 2 HT15 53.9 15
## 3 1 3 HT15 53.4 15
## 4 1 4 HT15 52.2 15
## 5 1 5 HT15 54.6 15
## 6 2 1 HT15 47.3 15aov.rep <- aov(HT ~ Treatment + Treatment:block + Days + Days:Treatment,
data = HT_Long)
anova(aov.rep)## Analysis of Variance Table
##
## Response: HT
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 2 61993 30996.3 2347.9997 < 2.2e-16 ***
## Days 3 91141 30380.3 2301.3391 < 2.2e-16 ***
## Treatment:block 12 1524 127.0 9.6226 5.516e-08 ***
## Treatment:Days 6 14314 2385.6 180.7149 < 2.2e-16 ***
## Residuals 36 475 13.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Extract the Mean Square for Blk:Var and use this as MSE(a)
MSE_a <- anova(aov.rep)["Mean Sq"][3,1]
MS_Trt <- anova(aov.rep)["Mean Sq"][1,1]
print(c(MSE_a, MS_Trt))## [1] 127.0289 30996.2702#Visualization
#Generate summary statistics for all treatment combinations
DT.means <- emmeans(aov.rep, specs=c("Days", "Treatment"))## NOTE: A nesting structure was detected in the fitted model:
## block %in% TreatmentDT.means#Print the table of means## Days Treatment emmean SE df lower.CL upper.CL
## 15 1 54.1 1.62 36 50.8 57.4
## 30 1 95.2 1.62 36 91.9 98.5
## 45 1 141.1 1.62 36 137.8 144.4
## 60 1 196.3 1.62 36 193.0 199.6
## 15 2 42.6 1.62 36 39.3 45.9
## 30 2 71.9 1.62 36 68.6 75.2
## 45 2 123.9 1.62 36 120.6 127.2
## 60 2 162.6 1.62 36 159.3 165.9
## 15 3 23.3 1.62 36 20.0 26.6
## 30 3 37.5 1.62 36 34.2 40.8
## 45 3 48.0 1.62 36 44.7 51.3
## 60 3 72.6 1.62 36 69.3 75.9
##
## Results are averaged over the levels of: block
## Confidence level used: 0.95#Assign the letters to the treatment means
DT.means.cld <- cld(DT.means,
Letters=letters,
decreasing = TRUE) |>
dplyr::arrange(Treatment, Days)
DT.means.cld#Print the table of means## Days Treatment emmean SE df lower.CL upper.CL .group
## 15 1 54.1 1.62 36 50.8 57.4 g
## 30 1 95.2 1.62 36 91.9 98.5 e
## 45 1 141.1 1.62 36 137.8 144.4 c
## 60 1 196.3 1.62 36 193.0 199.6 a
## 15 2 42.6 1.62 36 39.3 45.9 hi
## 30 2 71.9 1.62 36 68.6 75.2 f
## 45 2 123.9 1.62 36 120.6 127.2 d
## 60 2 162.6 1.62 36 159.3 165.9 b
## 15 3 23.3 1.62 36 20.0 26.6 j
## 30 3 37.5 1.62 36 34.2 40.8 i
## 45 3 48.0 1.62 36 44.7 51.3 gh
## 60 3 72.6 1.62 36 69.3 75.9 f
##
## Results are averaged over the levels of: block
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 12 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.HT_Long |>
group_by(Days, Treatment) |>
dplyr::summarize(Mean = mean(HT),
SD = sd(HT),
n = length(HT)) |>
dplyr::mutate(SE = SD/sqrt(n)) |>
ggplot(aes(x = Days, y = Mean, group=Treatment)) +
geom_line(aes(color = Treatment)) +
geom_point() +
geom_errorbar(aes(ymin = Mean - SE,
ymax = Mean + SE,
colour = Treatment),
width = 0.2,
size = 1) +
labs(y = "Mean Height") +
theme_classic() +
theme(legend.position="bottom") ## `summarise()` has grouped output by 'Days'. You can override using the
## `.groups` argument.
##Correlation analysis
#Loading a data set into R
plant_data <- read_excel("ACIAR Data2.xlsx", sheet = "Plant")
head(plant_data)## # A tibble: 6 × 7
## Treatment Replication Variety ODW TotalN TotalP TotalK
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 47.1 2.74 5.16 2.84
## 2 1 2 1 15.6 2.26 4.09 3.00
## 3 1 3 1 42.1 2.41 3.74 2.81
## 4 1 4 1 27.7 2.56 5.22 2.97
## 5 2 1 1 62.8 2.39 4.19 2.66
## 6 2 2 1 60.3 2.56 4.51 3.04#Correlation coefficient
#Correlation between two variables
plant_data |>
select(ODW, TotalN) |>
cor()## ODW TotalN
## ODW 1.0000000 0.3285493
## TotalN 0.3285493 1.0000000#Test of significance of the correlation
plant_data |>
select(ODW, TotalN) |>
cor_test()## # A tibble: 1 × 8
## var1 var2 cor statistic p conf.low conf.high method
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 ODW TotalN 0.33 2.56 0.0134 0.0719 0.544 Pearson#Correlation between several variables
plant_data |>
select(-c(Treatment, Variety, Replication)) |>
cor()## ODW TotalN TotalP TotalK
## ODW 1.0000000 0.3285493 -0.5482409 0.4746828
## TotalN 0.3285493 1.0000000 -0.6361370 0.7781900
## TotalP -0.5482409 -0.6361370 1.0000000 -0.7443045
## TotalK 0.4746828 0.7781900 -0.7443045 1.0000000#Correlation between several variables
plant_data |>
select(-c(Treatment, Variety, Replication)) |>
cor_test()## # A tibble: 16 × 8
## var1 var2 cor statistic p conf.low conf.high method
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 ODW ODW 1 493147422. 0 1 1 Pearson
## 2 ODW TotalN 0.33 2.56 1.34e- 2 0.0719 0.544 Pearson
## 3 ODW TotalP -0.55 -4.82 1.22e- 5 -0.709 -0.333 Pearson
## 4 ODW TotalK 0.47 3.96 2.19e- 4 0.242 0.656 Pearson
## 5 TotalN ODW 0.33 2.56 1.34e- 2 0.0719 0.544 Pearson
## 6 TotalN TotalN 1 Inf 0 1 1 Pearson
## 7 TotalN TotalP -0.64 -6.06 1.37e- 7 -0.770 -0.448 Pearson
## 8 TotalN TotalK 0.78 9.11 1.69e-12 0.648 0.864 Pearson
## 9 TotalP ODW -0.55 -4.82 1.22e- 5 -0.709 -0.333 Pearson
## 10 TotalP TotalN -0.64 -6.06 1.37e- 7 -0.770 -0.448 Pearson
## 11 TotalP TotalP 1 348707886. 0 1 1 Pearson
## 12 TotalP TotalK -0.74 -8.19 4.88e-11 -0.842 -0.599 Pearson
## 13 TotalK ODW 0.47 3.96 2.19e- 4 0.242 0.656 Pearson
## 14 TotalK TotalN 0.78 9.11 1.69e-12 0.648 0.864 Pearson
## 15 TotalK TotalP -0.74 -8.19 4.88e-11 -0.842 -0.599 Pearson
## 16 TotalK TotalK 1 348707886. 0 1 1 Pearson#Correlation between several variables
plant_data |>
select(-c(Treatment, Variety, Replication)) |>
correlation::correlation(include_factors = TRUE, method = "Pearson")## # Correlation Matrix (Pearson-method)
##
## Parameter1 | Parameter2 | r | 95% CI | t(54) | p
## --------------------------------------------------------------------
## ODW | TotalN | 0.33 | [ 0.07, 0.54] | 2.56 | 0.013*
## ODW | TotalP | -0.55 | [-0.71, -0.33] | -4.82 | < .001***
## ODW | TotalK | 0.47 | [ 0.24, 0.66] | 3.96 | < .001***
## TotalN | TotalP | -0.64 | [-0.77, -0.45] | -6.06 | < .001***
## TotalN | TotalK | 0.78 | [ 0.65, 0.86] | 9.11 | < .001***
## TotalP | TotalK | -0.74 | [-0.84, -0.60] | -8.19 | < .001***
##
## p-value adjustment method: Holm (1979)
## Observations: 56#Visualization #Scatter Plot
#from the ggstatplot package
ggscatterstats(
data = plant_data,
x = ODW,
y = TotalN,
bf.message = FALSE,
marginal = FALSE # remove histograms
)
plant_data |>
select(-c(Treatment, Variety, Replication)) |>
ggpairs() #from the ggally package
#Correlogram
plant_data |>
select(-c(Treatment, Variety, Replication)) |>
ggcorr() #from the ggally package
plant_data |>
mutate(TRT = as.factor(Treatment),
VAR = as.factor(Variety)) |>
ggpairs(columns = 4:7, ggplot2::aes(colour=VAR))