##Just Run This Chunk 
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
## Warning: package 'multcomp' was built under R version 4.1.1
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.1.1
## Loading required package: survival
## Loading required package: TH.data
## Warning: package 'TH.data' was built under R version 4.1.1
## 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
## Warning: package 'tidyverse' was built under R version 4.1.1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.1.3     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1
## v purrr   0.3.4
## Warning: package 'readr' was built under R version 4.1.1
## Warning: package 'forcats' was built under R version 4.1.1
## -- 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
## Warning: package 'ggpubr' was built under R version 4.1.1
require(rstatix)
## Loading required package: rstatix
## Warning: package 'rstatix' was built under R version 4.1.1
## 
## 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
## Warning: package 'lmerTest' was built under R version 4.1.1
## 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
## Warning: package 'Rmisc' was built under R version 4.1.1
## Loading required package: lattice
## Loading required package: plyr
## Warning: package 'plyr' was built under R version 4.1.1
## ------------------------------------------------------------------------------
## 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("C:/Users/Usuario/Desktop/Universidad/semestre 5/Bioestadistica/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)

##Exploratoring data

##Time Series Plot by Temp-CO2 Treatment

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

En la gráfica anterior se están comparando los comportamientos de dos genotipos con respecto a la cantidad de area foliar que producen a diferentes temperaturas y cantidades de CO2 en el transcurso de un periodo de tiempo. Para el primer genotipo PA107 no se ven diferencias significativas en la cantidad de area foliar que se obtiene al varias la temperatura y el CO2, mientras que para el segundo genotipo SCA6 se tienen cambios relevantes en la cantidad de area foliar al variar las temperaturas y concentraciones de CO2.

Para ambas gráficas se puede notar que antes de la semana 20 se tiene un crecimiento elevado del area foliar, pero que luego de esta semana se estabiliza el crecimiento y se empiezan a evidenciar o no las diferencias en la cantidad de area foliar según las temperaturas y concentraciones de CO2.

##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

Del análisis con el anova se puede ver que hay un efecto significativo de la temperatura, quiere decir que hay diferencia en la cantidad de area foliar entre las distintas temperaturas. Seguidamente, se analiza el genotipo en relación con las distintas variables; temperatura y CO2. Con respecto a la temperatura se evidenció que el genotipo tiende a variar según los cambios de esta, haciendo que responda de manera diferente y significativa. Mientras que con la variable CO2 no se registraron cambios significativos en la cantidad de area foliar de los genotipos. Por último, también se evaluó el tiempo, para el cual se obtuvo cambios significativos en la cantidad de area foliar para los genotipos.

##postanova - CO2

##Postanova


bxp = ggboxplot(
  LeafLength, x = "Time", y = "LeafL",
  color = "CO2", palette = "jco"
)


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)
  )

Según la gráfica de postanova para CO2 se ve que la diferencia en la concentración de este no estimula cambios significativos en la longitud del area foliar en ningun momento del tiempo. Para el CO2 no tiene sentido hacer el postanova, por la falta de cambios significativos, pero en este caso puede ser de utilidad para ver como da el postanova de la temperatura.

##postanova - temp

  ##POstanova Temp
  
  
bxp = ggboxplot(
  LeafLength, x = "Time", y = "LeafL",
  color = "Temp", palette = "jco"
)

pwc = LeafLength %>%
  group_by(Time) %>%
  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)", " = ", "17.188 , ", "p", " = ", "<0.0001 "),
    caption = get_pwc_label(pwc)
  )

En está gráfica del postanova para la temperatura se ven tres niveles de temperatura; baja en color azul, media en color amarillo y alta en color gris. Hasta la semana 19 se ve un crecimiento ascendente, que después de esta semana se estabiliza y desde ahí se evidencian las diferencias en la cantidad de area foliar según la temperatura. Por la gráfica se concluye que a temperaturas altas entre 27 a 36, la planta se ve afectada bajando a la menor cantidad de area foliar, mientras que la mayor cantidad de area foliar se registra a temperaturas bajas entre 22 a 31. Estas variaciones representan unas diferencias significativas en la interacción entre el area foliar y la temperatura.

##postanova - interaccion temp vs genotipo

  ##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)
  )

En esta gráfica se está considerando la interacción entre temperatura y genotipo sobre la cantidad de area foliar resultante. Para el genotipo 1 se ve que no hay variaciones significativas en la cantidad de area foliar según la temperatura a la que se exponga, por este motivo este genotipo será más tolerante ante los cambios de temperatura pudiendo mantener sus area foliar. Por otro lado, para el genotipo 2 se ve que sí hay respuesta significativa ante los cambios de temperatura, lo que se puede entender como que es un genotipo menos tolerante a esta variable, haciendo que sea suceptible a tener menos area foliar al aumento de la temperatura y a tener más area foliar a la disminución de la temperatura.