setwd("/Volumes/GoogleDrive/Mi unidad/Agrosavia/Env_muestra/data")
datos<-read.table("grano3.csv", header=T, sep=',')
datos$curva <- factor(datos$curva, levels = c("1", "2", "3"),
labels = c("P3", "P1", "P2"))
datos$gen<-as.factor(datos$gen)
datos$curva<-as.factor(datos$curva)
datos$id<-as.factor(datos$id)
datos$muestra<-as.factor(datos$muestra)
datos$dia<-as.factor(datos$dia)
library(ggplot2)
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble 3.1.6 ✓ purrr 0.3.4
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.1 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange() masks plyr::arrange()
## x purrr::compact() masks plyr::compact()
## x dplyr::count() masks plyr::count()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::id() masks plyr::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::mutate() masks plyr::mutate()
## x dplyr::rename() masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggpubr)
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
##
## mutate
library(rstatix)
##
## Attaching package: 'rstatix'
## The following objects are masked from 'package:plyr':
##
## desc, mutate
## The following object is masked from 'package:stats':
##
## filter
##Summary statistics
summ<-datos %>%
group_by(curva, gen, dia) %>%
get_summary_stats(cd.grano, type = "mean_sd")
summ %>% as_tibble() %>% print(n=Inf)
## # A tibble: 36 × 7
## curva gen dia variable n mean sd
## <fct> <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 P3 CCN51 0 cd.grano 3 8.68 0.855
## 2 P3 CCN51 2 cd.grano 3 8.40 0.262
## 3 P3 CCN51 5 cd.grano 3 7.49 0.181
## 4 P3 CCN51 6 cd.grano 3 7.60 0.286
## 5 P3 ICS95 0 cd.grano 3 11.5 0.469
## 6 P3 ICS95 2 cd.grano 3 11.3 0.947
## 7 P3 ICS95 5 cd.grano 3 11.3 0.881
## 8 P3 ICS95 6 cd.grano 3 10.3 0.739
## 9 P3 TCS01 0 cd.grano 3 9.03 2.35
## 10 P3 TCS01 2 cd.grano 3 7.98 0.831
## 11 P3 TCS01 5 cd.grano 3 7.75 0.903
## 12 P3 TCS01 6 cd.grano 3 8.28 0.931
## 13 P1 CCN51 0 cd.grano 3 8.50 0.901
## 14 P1 CCN51 2 cd.grano 3 9.09 0.593
## 15 P1 CCN51 5 cd.grano 3 7.99 0.99
## 16 P1 CCN51 6 cd.grano 3 7.95 0.409
## 17 P1 ICS95 0 cd.grano 3 9.17 1.05
## 18 P1 ICS95 2 cd.grano 3 9.36 0.616
## 19 P1 ICS95 5 cd.grano 3 9.46 1.27
## 20 P1 ICS95 6 cd.grano 3 9.49 0.801
## 21 P1 TCS01 0 cd.grano 3 6.86 0.3
## 22 P1 TCS01 2 cd.grano 3 7.42 1.76
## 23 P1 TCS01 5 cd.grano 3 6.68 0.575
## 24 P1 TCS01 6 cd.grano 3 6.78 0.429
## 25 P2 CCN51 0 cd.grano 3 6.52 1.10
## 26 P2 CCN51 2 cd.grano 3 5.71 0.742
## 27 P2 CCN51 5 cd.grano 3 5.75 0.184
## 28 P2 CCN51 6 cd.grano 3 6.26 0.109
## 29 P2 ICS95 0 cd.grano 3 11.0 1.27
## 30 P2 ICS95 2 cd.grano 3 11.2 0.519
## 31 P2 ICS95 5 cd.grano 3 11.1 0.529
## 32 P2 ICS95 6 cd.grano 3 10.4 0.156
## 33 P2 TCS01 0 cd.grano 3 8.23 0.516
## 34 P2 TCS01 2 cd.grano 3 6.69 0.393
## 35 P2 TCS01 5 cd.grano 3 8.11 0.368
## 36 P2 TCS01 6 cd.grano 3 7.28 0.297
##Visualization
bxp <- ggboxplot(
datos, x = "curva", y = "cd.grano",
color = "dia", palette = "jco",
facet.by = "gen"
)
bxp
##Check assumptions
##Outliers
datos %>%
group_by(curva, gen, dia) %>%
identify_outliers(cd.grano)
## [1] curva gen dia muestra id cd.grano
## [7] X X.1 X.2 X.3 X.4 X.5
## [13] X.6 X.7 is.outlier is.extreme
## <0 rows> (or 0-length row.names)
##Normality assumption
##Compute Shapiro-Wilk test for each combinations of factor levels:
norm<-datos %>%
group_by(curva, gen, dia) %>%
shapiro_test(cd.grano)
norm %>% as_tibble() %>% print(n=Inf)
## # A tibble: 36 × 6
## curva gen dia variable statistic p
## <fct> <fct> <fct> <chr> <dbl> <dbl>
## 1 P3 CCN51 0 cd.grano 0.849 0.237
## 2 P3 CCN51 2 cd.grano 0.841 0.216
## 3 P3 CCN51 5 cd.grano 0.803 0.122
## 4 P3 CCN51 6 cd.grano 0.975 0.698
## 5 P3 ICS95 0 cd.grano 0.845 0.227
## 6 P3 ICS95 2 cd.grano 0.999 0.927
## 7 P3 ICS95 5 cd.grano 0.867 0.286
## 8 P3 ICS95 6 cd.grano 0.982 0.740
## 9 P3 TCS01 0 cd.grano 0.988 0.787
## 10 P3 TCS01 2 cd.grano 0.798 0.109
## 11 P3 TCS01 5 cd.grano 0.856 0.256
## 12 P3 TCS01 6 cd.grano 0.992 0.824
## 13 P1 CCN51 0 cd.grano 0.899 0.382
## 14 P1 CCN51 2 cd.grano 0.999 0.940
## 15 P1 CCN51 5 cd.grano 0.906 0.403
## 16 P1 CCN51 6 cd.grano 0.990 0.807
## 17 P1 ICS95 0 cd.grano 0.994 0.852
## 18 P1 ICS95 2 cd.grano 0.962 0.625
## 19 P1 ICS95 5 cd.grano 0.999 0.928
## 20 P1 ICS95 6 cd.grano 0.838 0.209
## 21 P1 TCS01 0 cd.grano 0.909 0.414
## 22 P1 TCS01 2 cd.grano 0.920 0.452
## 23 P1 TCS01 5 cd.grano 1.00 0.984
## 24 P1 TCS01 6 cd.grano 0.815 0.151
## 25 P2 CCN51 0 cd.grano 0.814 0.149
## 26 P2 CCN51 2 cd.grano 0.922 0.460
## 27 P2 CCN51 5 cd.grano 0.975 0.700
## 28 P2 CCN51 6 cd.grano 0.989 0.803
## 29 P2 ICS95 0 cd.grano 0.880 0.325
## 30 P2 ICS95 2 cd.grano 0.971 0.671
## 31 P2 ICS95 5 cd.grano 0.860 0.268
## 32 P2 ICS95 6 cd.grano 0.974 0.690
## 33 P2 TCS01 0 cd.grano 0.992 0.827
## 34 P2 TCS01 2 cd.grano 1.00 0.986
## 35 P2 TCS01 5 cd.grano 0.997 0.892
## 36 P2 TCS01 6 cd.grano 0.988 0.794
##Create QQ plot for each cell of design:
ggqqplot(datos, "cd.grano", ggtheme = theme_bw()) +
facet_grid(dia~ curva*gen, labeller = "label_both")
##Homogneity of variance assumption
##Compute the Levene’s test at each level of the within-subjects factor, here time variable:
lev<-datos %>%
group_by(dia) %>%
levene_test(cd.grano ~ curva*gen)
lev %>% as_tibble() %>% print(n=Inf)
## # A tibble: 4 × 5
## dia df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 0 8 18 0.754 0.645
## 2 2 8 18 0.619 0.751
## 3 5 8 18 0.626 0.745
## 4 6 8 18 0.750 0.649
##Computation
res.aov <- anova_test(
data = datos, dv = cd.grano, wid = id,
within = dia, between = c(curva, gen)
)
get_anova_table(res.aov)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 curva 2.00 18.00 9.467 2.00e-03 * 0.286
## 2 gen 2.00 18.00 94.552 2.83e-10 * 0.800
## 3 dia 2.29 41.25 2.847 6.30e-02 0.089
## 4 curva:gen 4.00 18.00 11.457 8.42e-05 * 0.492
## 5 curva:dia 4.58 41.25 1.581 1.91e-01 0.098
## 6 gen:dia 4.58 41.25 1.027 4.11e-01 0.066
## 7 curva:gen:dia 9.17 41.25 0.799 6.21e-01 0.099
##Splitting dataframe by temperature ramp
## Protocol 3 (P3)
datos.curve1<-filter(datos, curva=="P3")
##Check assumptions
##Outliers
datos.curve1 %>%
group_by(gen, dia) %>%
identify_outliers(cd.grano)
## [1] gen dia curva muestra id cd.grano
## [7] X X.1 X.2 X.3 X.4 X.5
## [13] X.6 X.7 is.outlier is.extreme
## <0 rows> (or 0-length row.names)
##Normality assumption
##Compute Shapiro-Wilk test for each combinations of factor levels:
norm1<-datos.curve1 %>%
group_by(gen, dia) %>%
shapiro_test(cd.grano)
norm1 %>% as_tibble() %>% print(n=Inf)
## # A tibble: 12 × 5
## gen dia variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 CCN51 0 cd.grano 0.849 0.237
## 2 CCN51 2 cd.grano 0.841 0.216
## 3 CCN51 5 cd.grano 0.803 0.122
## 4 CCN51 6 cd.grano 0.975 0.698
## 5 ICS95 0 cd.grano 0.845 0.227
## 6 ICS95 2 cd.grano 0.999 0.927
## 7 ICS95 5 cd.grano 0.867 0.286
## 8 ICS95 6 cd.grano 0.982 0.740
## 9 TCS01 0 cd.grano 0.988 0.787
## 10 TCS01 2 cd.grano 0.798 0.109
## 11 TCS01 5 cd.grano 0.856 0.256
## 12 TCS01 6 cd.grano 0.992 0.824
##Create QQ plot for each cell of design:
ggqqplot(datos.curve1, "cd.grano", ggtheme = theme_bw()) +
facet_grid(dia~ curva*gen, labeller = "label_both")
##Homogneity of variance assumption
##Compute the Levene’s test at each level of the within-subjects factor, here time variable:
lev1<-datos.curve1 %>%
group_by(dia) %>%
levene_test(cd.grano ~ gen)
lev1 %>% as_tibble() %>% print(n=Inf)
## # A tibble: 4 × 5
## dia df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 0 2 6 1.50 0.296
## 2 2 2 6 0.553 0.602
## 3 5 2 6 0.515 0.622
## 4 6 2 6 0.796 0.494
##Computation
res.aov1 <- anova_test(
data = datos.curve1, dv = cd.grano, wid = id,
within = dia, between = gen
)
get_anova_table(res.aov1)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 gen 2 6 28.822 0.000838 * 0.755
## 2 dia 3 18 2.198 0.124000 0.199
## 3 gen:dia 6 18 0.608 0.721000 0.121
#CCN51
datos.ccn<-filter(datos.curve1, gen=="CCN51")
res.aov.ccn1 <- anova_test(
data = datos.ccn, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.ccn1)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 5.222 0.041 * 0.63
#ICS95
datos.ics<-filter(datos.curve1, gen=="ICS95")
res.aov.ics1 <- anova_test(
data = datos.ics, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.ics1)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 2.399 0.167 0.36
#TCS01
datos.tcs<-filter(datos.curve1, gen=="TCS01")
res.aov.tcs1 <- anova_test(
data = datos.tcs, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.tcs1)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 0.478 0.709 0.151
## Protocol 1 (P1)
datos.curve2<-filter(datos, curva=="P1")
##Check assumptions
##Outliers
datos.curve2 %>%
group_by(gen, dia) %>%
identify_outliers(cd.grano)
## [1] gen dia curva muestra id cd.grano
## [7] X X.1 X.2 X.3 X.4 X.5
## [13] X.6 X.7 is.outlier is.extreme
## <0 rows> (or 0-length row.names)
##Normality assumption
##Compute Shapiro-Wilk test for each combinations of factor levels:
norm2<-datos.curve2 %>%
group_by(gen, dia) %>%
shapiro_test(cd.grano)
norm2 %>% as_tibble() %>% print(n=Inf)
## # A tibble: 12 × 5
## gen dia variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 CCN51 0 cd.grano 0.899 0.382
## 2 CCN51 2 cd.grano 0.999 0.940
## 3 CCN51 5 cd.grano 0.906 0.403
## 4 CCN51 6 cd.grano 0.990 0.807
## 5 ICS95 0 cd.grano 0.994 0.852
## 6 ICS95 2 cd.grano 0.962 0.625
## 7 ICS95 5 cd.grano 0.999 0.928
## 8 ICS95 6 cd.grano 0.838 0.209
## 9 TCS01 0 cd.grano 0.909 0.414
## 10 TCS01 2 cd.grano 0.920 0.452
## 11 TCS01 5 cd.grano 1.00 0.984
## 12 TCS01 6 cd.grano 0.815 0.151
##Create QQ plot for each cell of design:
ggqqplot(datos.curve2, "cd.grano", ggtheme = theme_bw()) +
facet_grid(dia~ curva*gen, labeller = "label_both")
##Homogneity of variance assumption
##Compute the Levene’s test at each level of the within-subjects factor, here time variable:
lev2<-datos.curve2 %>%
group_by(dia) %>%
levene_test(cd.grano ~ gen)
lev2 %>% as_tibble() %>% print(n=Inf)
## # A tibble: 4 × 5
## dia df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 0 2 6 0.675 0.544
## 2 2 2 6 0.801 0.491
## 3 5 2 6 0.398 0.688
## 4 6 2 6 0.215 0.813
##Computation
res.aov2 <- anova_test(
data = datos.curve2, dv = cd.grano, wid = id,
within = dia, between = gen
)
get_anova_table(res.aov2)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 gen 2.00 6.00 11.419 0.009 * 0.649
## 2 dia 1.57 9.44 1.170 0.338 0.091
## 3 gen:dia 3.15 9.44 0.540 0.674 0.085
#CCN51
datos.ccn<-filter(datos.curve2, gen=="CCN51")
res.aov.ccn2 <- anova_test(
data = datos.ccn, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.ccn2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 1.687 0.268 0.357
#ICS95
datos.ics<-filter(datos.curve2, gen=="ICS95")
res.aov.ics2 <- anova_test(
data = datos.ics, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.ics2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 0.23 0.872 0.025
#TCS01
datos.tcs<-filter(datos.curve2, gen=="TCS01")
res.aov.tcs2 <- anova_test(
data = datos.tcs, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.tcs2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 0.378 0.772 0.12
## Protocol 2 (P2)
datos.curve3<-filter(datos, curva=="P2")
##Check assumptions
##Outliers
datos.curve3 %>%
group_by(gen, dia) %>%
identify_outliers(cd.grano)
## [1] gen dia curva muestra id cd.grano
## [7] X X.1 X.2 X.3 X.4 X.5
## [13] X.6 X.7 is.outlier is.extreme
## <0 rows> (or 0-length row.names)
##Normality assumption
##Compute Shapiro-Wilk test for each combinations of factor levels:
norm2<-datos.curve3 %>%
group_by(gen, dia) %>%
shapiro_test(cd.grano)
norm2 %>% as_tibble() %>% print(n=Inf)
## # A tibble: 12 × 5
## gen dia variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 CCN51 0 cd.grano 0.814 0.149
## 2 CCN51 2 cd.grano 0.922 0.460
## 3 CCN51 5 cd.grano 0.975 0.700
## 4 CCN51 6 cd.grano 0.989 0.803
## 5 ICS95 0 cd.grano 0.880 0.325
## 6 ICS95 2 cd.grano 0.971 0.671
## 7 ICS95 5 cd.grano 0.860 0.268
## 8 ICS95 6 cd.grano 0.974 0.690
## 9 TCS01 0 cd.grano 0.992 0.827
## 10 TCS01 2 cd.grano 1.00 0.986
## 11 TCS01 5 cd.grano 0.997 0.892
## 12 TCS01 6 cd.grano 0.988 0.794
##Create QQ plot for each cell of design:
ggqqplot(datos.curve3, "cd.grano", ggtheme = theme_bw()) +
facet_grid(dia~ curva*gen, labeller = "label_both")
##Homogneity of variance assumption
##Compute the Levene’s test at each level of the within-subjects factor, here time variable:
lev2<-datos.curve3 %>%
group_by(dia) %>%
levene_test(cd.grano ~ gen)
lev2 %>% as_tibble() %>% print(n=Inf)
## # A tibble: 4 × 5
## dia df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 0 2 6 0.226 0.804
## 2 2 2 6 0.227 0.803
## 3 5 2 6 0.375 0.702
## 4 6 2 6 0.824 0.483
##Computation
res.aov2 <- anova_test(
data = datos.curve3, dv = cd.grano, wid = id,
within = dia, between = gen
)
get_anova_table(res.aov2)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 gen 2 6 158.889 6.36e-06 * 0.941
## 2 dia 3 18 2.841 6.70e-02 0.248
## 3 gen:dia 6 18 2.027 1.15e-01 0.320
#CCN51
datos.ccn<-filter(datos.curve3, gen=="CCN51")
res.aov.ccn2 <- anova_test(
data = datos.ccn, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.ccn2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 1.296 0.359 0.279
#ICS95
datos.ics<-filter(datos.curve3, gen=="ICS95")
res.aov.ics2 <- anova_test(
data = datos.ics, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.ics2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 0.732 0.57 0.221
#TCS01
datos.tcs<-filter(datos.curve3, gen=="TCS01")
res.aov.tcs2 <- anova_test(
data = datos.tcs, dv = cd.grano, wid = id,
within = dia
)
get_anova_table(res.aov.tcs2)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 dia 3 6 10.466 0.008 * 0.788
## Gráficas por réplica y genotipo
datos$dia<-as.numeric(as.character(datos$dia))
##Gráfica por réplica compuesta
pht<- ggplot(datos, aes(x = dia)) +
facet_grid(curva~gen*muestra) +
geom_line(aes(y=cd.grano)) +
geom_point(aes(y=cd.grano)) +
scale_y_continuous(name = expression("Cd (mg*Kg"^"-1)")) + # Etiqueta de la variable continua
scale_x_continuous(name = "dÃa", breaks=seq(0,7,1)) + # Etiqueta de los grupos
theme(axis.line = element_line(colour = "black", # Personalización del tema
size = 0.25)) +
theme(text = element_text(size = 12))
pht
## Gráfica por genotipo
datos2<-summarySE (datos, measurevar = "cd.grano", groupvars = c("curva", "gen","dia"))
write.csv(datos2, "/Volumes/GoogleDrive/Mi unidad/Agrosavia/Env_muestra/data/datos_mean.csv")
pht2<- ggplot(datos2, aes(x = dia)) +
facet_grid(curva~gen) +
geom_errorbar(aes(ymin=cd.grano-ci, ymax=cd.grano+ci), width=.1) +
geom_line(aes(y=cd.grano)) +
geom_point(aes(y=cd.grano)) +
scale_y_continuous(name = expression("Cd (mg*Kg"^"-1)")) + # Etiqueta de la variable continua
scale_x_continuous(name = "dÃa", breaks=seq(0,7,1)) + # Etiqueta de los grupos
theme(axis.line = element_line(colour = "black", # Personalización del tema
size = 0.25)) +
theme(text = element_text(size = 15))
pht2
