Ejercicio- Analisis Genotipo de área foliar

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)

Gráfico de series de tiempo por tratamiento Temperatura-CO2

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.

Anova de Medidas Repetidas

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.

Post-anova CO2

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.

Postanova Temperatura

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.

Postanova Temp-Genotipo

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.