require(readxl)
## Loading required package: readxl
require(ggplot2)
## Loading required package: ggplot2
require(plotly)
## Loading required package: plotly
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
require(emmeans)
## Loading required package: emmeans
require(multcomp)
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
##
## select
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
require(tidyverse)
## Loading required package: tidyverse
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble 3.1.2 v dplyr 1.0.6
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1
## v purrr 0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks plotly::filter(), stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks MASS::select(), plotly::select()
require(ggpubr)
## Loading required package: ggpubr
require(rstatix)
## Loading required package: rstatix
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:stats':
##
## filter
require(lmerTest)
## Loading required package: lmerTest
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
require(Rmisc)
## Loading required package: Rmisc
## Loading required package: lattice
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:rstatix':
##
## desc, mutate
## The following object is masked from 'package:ggpubr':
##
## mutate
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:purrr':
##
## compact
## The following objects are masked from 'package:plotly':
##
## arrange, mutate, rename, summarise
LeafLength = read_excel("~/2021-II/Bioestadistica/YDRAY-LeafLength.xlsx")
LeafLength$Temp_C02=paste(LeafLength$Temp,"_and_",LeafLength$CO2,sep="")
LeafLength$id=as.factor(paste(LeafLength$Cab,"_",LeafLength$Obs,sep=""))
LeafLength$Genotype=as.factor(LeafLength$Genotype)
gd=summarySE(LeafLength, measurevar="LeafL", groupvars=c("Temp_C02","Time","Genotype"))
g1=ggplot(gd,aes(x=Time,y=LeafL,color=Temp_C02, group =Temp_C02))+
geom_line(data=gd)+geom_point(data=gd, size = 2)+
facet_grid(~Genotype)+theme_bw()+
labs(title = "Leaf Length", x = "Time in days", y = "Length in cm")
g1
Lo arrojado por el gráfico nos permite inferir que las plantas con el genotipo PA107 no altera el área foliar de sus hojas según las condiciones ambientales dominantes de ninguno de los escenarios, mientras que el genotipo SCA6 es si presenta alteraciones teniendo una disminución de area foliar para temperaturas más altas como de 27 a 36° C. Respecto a los niveles de Co2 no parece tener relación alguna con el crecimiento del area foliar.
LeafLength$Time=as.factor(LeafLength$Time)
LeafLength$Temp=as.factor(LeafLength$Temp)
LeafLength$CO2=as.factor(LeafLength$CO2)
anova=aov(LeafL~Temp*CO2*Genotype+Time+Error(id/Genotype/Time),data=LeafLength)
summary(anova)
##
## Error: id
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp 2 6614 3307 17.188 1.06e-05 ***
## CO2 1 214 214 1.113 0.30
## Temp:CO2 2 396 198 1.028 0.37
## Residuals 30 5772 192
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Genotype
## Df Sum Sq Mean Sq F value Pr(>F)
## Genotype 1 871 871 2.016 0.16598
## Temp:Genotype 2 7270 3635 8.409 0.00126 **
## CO2:Genotype 1 46 46 0.108 0.74522
## Temp:CO2:Genotype 2 2673 1337 3.092 0.06013 .
## Residuals 30 12968 432
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: id:Genotype:Time
## Df Sum Sq Mean Sq F value Pr(>F)
## Time 18 17707 983.7 140.2 <2e-16 ***
## Residuals 1278 8964 7.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Se corrobora lo dicho anteriormente pues notamos diferencias significativas para la variable de Temperatura respecto al Genotipo, la variable CO2 no muestra relación con respecto a ninguna de las variables pues esta no tiene diferencias significativas.
bxp = ggboxplot(LeafLength, x = "Time", y = "LeafL",
color = "CO2", palette = "jco")
bxp
pwc = LeafLength %>%
group_by(Time) %>%
pairwise_t_test(LeafL ~ CO2, paired = TRUE,p.adjust.method = "bonferroni")
pwc = pwc %>% add_xy_position(x = "Time")
bxp + stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(subtitle = paste("Anova, ","F(1)", " = ", "1.113, ", "p", " = ", "0.30"),
caption = get_pwc_label(pwc))
Para confirmar la ausencia de relación de CO2 con el crecimiento del area foliar se realizó una correccion de Bonferroni con el fin de que las comparaciones multiples realizadas no interfieran con el resultado y el analisis. Aquí podemos ver que existe solo un par de valores significativos donde el nivel de CO2 de 700 tiene mayor significancia en el crecimiento foliar, sin embargo, estos datos son muy limitados, descartando la posibilidad de implicación de los niveles de CO2 realmente en el crecimiento.
bxp = ggboxplot(LeafLength, x = "Time", y = "LeafL",
color = "Temp", palette = "jco")
bxp
pwc = LeafLength %>%
group_by(Time) %>%
pairwise_t_test(LeafL ~ Temp, paired = TRUE,
p.adjust.method = "bonferroni")
pwc
## # A tibble: 57 x 11
## Time .y. group1 group2 n1 n2 statistic df p p.adj
## * <fct> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 10 LeafL 31_22 33.5_24.5 24 24 -5.14 23 0.0000333 9.99e-5
## 2 10 LeafL 31_22 36_27 24 24 -2.95 23 0.007 2.2 e-2
## 3 10 LeafL 33.5_24.5 36_27 24 24 1.40 23 0.176 5.28e-1
## 4 13 LeafL 31_22 33.5_24.5 24 24 -4.00 23 0.00056 2 e-3
## 5 13 LeafL 31_22 36_27 24 24 -1.08 23 0.291 8.73e-1
## 6 13 LeafL 33.5_24.5 36_27 24 24 1.57 23 0.129 3.87e-1
## 7 16 LeafL 31_22 33.5_24.5 24 24 -1.04 23 0.307 9.21e-1
## 8 16 LeafL 31_22 36_27 24 24 1.45 23 0.161 4.83e-1
## 9 16 LeafL 33.5_24.5 36_27 24 24 3.31 23 0.003 9 e-3
## 10 19 LeafL 31_22 33.5_24.5 24 24 0.421 23 0.678 1 e+0
## # ... with 47 more rows, and 1 more variable: p.adj.signif <chr>
pwc = pwc %>% add_xy_position(x = "Time")
bxp + stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(subtitle = paste("Anova, ","F(2)", " = ", "17.188 , ", "p", " = ", "<0.0001 "),
caption = get_pwc_label(pwc))
El ajuste de Bonferroni para la temperatura confirma las suposisciones realizadas sobre su influencia en el crec del area foliar, siendo el rango de temperatura más optima el de 31_22 y la menos optima la de 36_27.
bxp = ggboxplot(LeafLength, x = "Time", y = "LeafL",
title= "Leaf Length", xlab ="time (days)", ylab= "Length (cm)",
color = "Temp", palette = "lancet", facet.by = "Genotype")
pwc = LeafLength %>%
group_by(Time, Genotype) %>%
pairwise_t_test(LeafL ~ Temp, paired = TRUE,p.adjust.method = "bonferroni")
pwc = pwc %>% add_xy_position(x = "Time")
bxp + stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(subtitle = paste("Anova, ","F(2)", " = ", "8.409 , ", "p", " = ", "0.00126"),
caption = get_pwc_label(pwc))
Para finalizar, se realizó un análisis de como la temperatura influye sobre cada genotipo, se encontró que efectivamente el genotipo PA107 no presenta cambios significativos determinados por las condiciones del ambiente; por su parte, el genotipo SCA6 tiene fluctuaciones en el crecimiento del area de sus hojas de acuerdo a la temperatura a la que sea sometido.