Cargar datos

library(readxl)
data <- read_excel("C:/Users/sebas/Desktop/Angie Herrera/2021-1/Bioestadistica/data.xlsx")
View(data)

Cargar librerias

require(ggplot2)
require(plotly)
require(emmeans)
require(multcoemp)
require(tidyverse)
require(ggpubr)
require(rstatix)
require(lmerTest)
require(Rmisc)
data$CO2<-as.factor(data$CO2) #Se nombra como un factor por tener dos niveles (400-700)
data$id=as.factor(paste(data$Cab,"_",data$Obvs,sep=""))
names(data)[4]= "Genotype" #Para cambiar el nombre de una de las columnas

Análisis exploratorio

g1<-ggplot(data, aes(x=CO2, y= Flushing, fill= Temp))+ geom_boxplot()+ facet_grid(~Genotype)+ theme_bw()
g1

La gráfica muestra la variabilidad en los dos genotipos en relación de la temperatura y el CO2, logrando evidenciar las diferencias entre temperaturas y no en concentraciones de CO2.

3 way Anova

anova<-aov(Flushing~Temp*CO2*Genotype+Error(id/Genotype),data=data)
summary(anova)
## 
## Error: id
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Temp       2  698.7   349.3  61.892 2.25e-11 ***
## CO2        1    0.5     0.5   0.089    0.768    
## Temp:CO2   2   16.6     8.3   1.469    0.246    
## Residuals 30  169.3     5.6                     
## ---
## 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  88.89   88.89  12.945 0.00114 **
## Temp:Genotype      2  30.53   15.26   2.223 0.12583   
## CO2:Genotype       1   0.50    0.50   0.073 0.78913   
## Temp:CO2:Genotype  2  36.08   18.04   2.627 0.08882 . 
## Residuals         30 206.00    6.87                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se confirma la significancia que tiene la temperatura y genotipo en el modelo.

Post-anova 1

#Diagrama de cajas
bplot<- ggboxplot(data, x = "Temp", y = "Flushing",title= "Flushing",xlab ="Temp", ylab= "Flushing",color = "Temp", palette = "lancet",facet.by = "CO2")
bplot

#Comparaciones
pwc<-data %>% group_by(CO2) %>% pairwise_t_test(Flushing ~ Temp, paired = TRUE,p.adjust.method = "bonferroni")
pwc
## # A tibble: 6 x 11
##   CO2   .y.   group1 group2    n1    n2 statistic    df       p   p.adj
## * <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>   <dbl>
## 1 400   Flus~ 31_22  33.5_~    12    12      5.37    11 2.28e-4 6.84e-4
## 2 400   Flus~ 31_22  36_27     12    12      7.26    11 1.63e-5 4.89e-5
## 3 400   Flus~ 33.5_~ 36_27     12    12      1.52    11 1.56e-1 4.68e-1
## 4 700   Flus~ 31_22  33.5_~    12    12      3.84    11 3.00e-3 8.00e-3
## 5 700   Flus~ 31_22  36_27     12    12      5.44    11 2.04e-4 6.12e-4
## 6 700   Flus~ 33.5_~ 36_27     12    12      3.78    11 3.00e-3 9.00e-3
## # ... with 1 more variable: p.adj.signif <chr>

Se presenta mayor significacia en las variables cuando se encuentran a un nivel 700 de CO2

Post-anova 2

bplot2<- ggboxplot(data, x = "Genotype", y = "Flushing",title= "Flushing",xlab="Genotype", ylab= "Flushing",color = "Genotype", palette = "lancet")
bplot2

pwc2<- data %>% group_by() %>% pairwise_t_test(Flushing ~ Genotype, paired = TRUE,p.adjust.method = "bonferroni")
pwc2
## # A tibble: 1 x 10
##   .y.      group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
## * <chr>    <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
## 1 Flushing PA 107 SCA 06    36    36     -3.38    35 0.002 0.002 **

Existe gran significacia entre los genotipos

Validación de normalidad

data %>%group_by(Temp,CO2) %>%shapiro_test(Flushing)
## # A tibble: 6 x 5
##   CO2   Temp      variable statistic       p
##   <fct> <chr>     <chr>        <dbl>   <dbl>
## 1 400   31_22     Flushing     0.903 0.175  
## 2 700   31_22     Flushing     0.966 0.860  
## 3 400   33.5_24.5 Flushing     0.973 0.936  
## 4 700   33.5_24.5 Flushing     0.797 0.00864
## 5 400   36_27     Flushing     0.916 0.255  
## 6 700   36_27     Flushing     0.903 0.171

Sí se distribuyen normalmente