\[\text{Valeria Riveros Parra y Juan Felipe Paz Paiba}\] \[\text{Analisis de varianza para medidas repetidas de tres vías}\]

library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.2.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.2.3
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.2.3
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(multcomp)
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.2.3
## Loading required package: survival
## Loading required package: TH.data
## Warning: package 'TH.data' was built under R version 4.2.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.2.3
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:rstatix':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(lme4)
## Warning: package 'lme4' was built under R version 4.2.3
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.2.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Datos en formato ancho
datos3 <- read_xlsx ("C:/Users/juanf/Downloads/m_repetidas.xlsx", sheet= "biomasa_3")
datos3
## # A tibble: 30 × 15
##    dosis_fert riego   `0_días` `05_días` `10_días` `15_días` `20_días` `25_días`
##    <chr>      <chr>      <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
##  1 A          Riego_…        0      2.95      3.31      3.12      3.27      3.86
##  2 A          Riego_…        0      2.53      3.02      3.35      3.41      3.73
##  3 A          Riego_…        0      2.68      3.07      3.39      2.72      3.52
##  4 A          Riego_…        0      2.77      2.91      3.11      3.42      2.99
##  5 A          Riego_…        0      2.83      3.15      3.46      3.15      3.08
##  6 A          Riego_…        0      2.88      3.18      3.27      3.75      3.24
##  7 A          Riego_…        0      2.92      3.21      2.67      3.37      3.81
##  8 A          Riego_…        0      2.86      2.96      3.24      2.9       3.45
##  9 A          Riego_…        0      2.99      3.27      3.04      3.51      3.03
## 10 A          Riego_…        0      3.02      3.14      3.19      3.44      3.57
## # ℹ 20 more rows
## # ℹ 7 more variables: `30_días` <dbl>, `35_días` <dbl>, `40_días` <dbl>,
## #   `45_días` <dbl>, `50_días` <dbl>, `55_días` <dbl>, `60_días` <dbl>
# Suponiendo que datos3 es tu dataframe original
datos_larg3 <- datos3 %>%
  gather(key = "tiempo", value = "biomasa", `0_días`, `05_días`, `10_días`, `15_días`, `20_días`, `25_días`, `30_días`, `35_días`, `40_días`, `45_días`, `50_días`, `55_días`, `60_días`) %>%
  mutate(tiempo = as.factor(tiempo))

# Crear la columna id que enumera del 1 al 5 de manera repetitiva
datos_larg3 <- datos_larg3 %>%
  mutate(id = rep(1:5, length.out = n()))

# Verificar el resultado
head(datos_larg3)
## # A tibble: 6 × 5
##   dosis_fert riego             tiempo biomasa    id
##   <chr>      <chr>             <fct>    <dbl> <int>
## 1 A          Riego_cada_2_días 0_días       0     1
## 2 A          Riego_cada_2_días 0_días       0     2
## 3 A          Riego_cada_2_días 0_días       0     3
## 4 A          Riego_cada_2_días 0_días       0     4
## 5 A          Riego_cada_2_días 0_días       0     5
## 6 A          Riego_cada_4_días 0_días       0     1
# Armar la matriz de datos
matriz_datos <- as.matrix(datos_larg3)

# Verificar la matriz de datos
head(matriz_datos)
##      dosis_fert riego               tiempo   biomasa  id 
## [1,] "A"        "Riego_cada_2_días" "0_días" " 0.000" "1"
## [2,] "A"        "Riego_cada_2_días" "0_días" " 0.000" "2"
## [3,] "A"        "Riego_cada_2_días" "0_días" " 0.000" "3"
## [4,] "A"        "Riego_cada_2_días" "0_días" " 0.000" "4"
## [5,] "A"        "Riego_cada_2_días" "0_días" " 0.000" "5"
## [6,] "A"        "Riego_cada_4_días" "0_días" " 0.000" "1"
library(collapsibleTree)
## Warning: package 'collapsibleTree' was built under R version 4.2.3
# Crear un árbol colapsable
tree <- collapsibleTree(datos_larg3, c("tiempo", "dosis_fert","riego", "biomasa"), collapsed = TRUE)

# Visualizar el árbol colapsable
tree
datos_larg3 %>%
  group_by(dosis_fert, riego, tiempo) %>%
  get_summary_stats(biomasa, type = "mean_sd")
## # A tibble: 78 × 7
##    dosis_fert riego             tiempo  variable     n  mean    sd
##    <chr>      <chr>             <fct>   <fct>    <dbl> <dbl> <dbl>
##  1 A          Riego_cada_2_días 0_días  biomasa      5  0    0    
##  2 A          Riego_cada_2_días 05_días biomasa      5  2.75 0.16 
##  3 A          Riego_cada_2_días 10_días biomasa      5  3.09 0.153
##  4 A          Riego_cada_2_días 15_días biomasa      5  3.29 0.159
##  5 A          Riego_cada_2_días 20_días biomasa      5  3.20 0.285
##  6 A          Riego_cada_2_días 25_días biomasa      5  3.44 0.387
##  7 A          Riego_cada_2_días 30_días biomasa      5  3.70 0.278
##  8 A          Riego_cada_2_días 35_días biomasa      5  3.75 0.249
##  9 A          Riego_cada_2_días 40_días biomasa      5  4.25 0.185
## 10 A          Riego_cada_2_días 45_días biomasa      5  4.73 0.293
## # ℹ 68 more rows
bxp <- ggboxplot(
  datos_larg3, x = "dosis_fert", y = "biomasa",
  color = "tiempo", palette = "Zissou1",
  facet.by = "riego", short.panel.labs = FALSE
  )
bxp

#Supuestos
datos_larg3 %>%
    group_by(dosis_fert, riego, tiempo) %>%
  identify_outliers(biomasa)
## # A tibble: 24 × 7
##    dosis_fert riego             tiempo  biomasa    id is.outlier is.extreme
##    <chr>      <chr>             <fct>     <dbl> <int> <lgl>      <lgl>     
##  1 A          Riego_cada_2_días 20_días    2.72     3 TRUE       FALSE     
##  2 A          Riego_cada_2_días 45_días    5.14     2 TRUE       FALSE     
##  3 A          Riego_cada_4_días 10_días    2.96     3 TRUE       FALSE     
##  4 A          Riego_cada_4_días 15_días    2.67     2 TRUE       FALSE     
##  5 A          Riego_cada_4_días 20_días    3.75     1 TRUE       FALSE     
##  6 A          Riego_cada_4_días 20_días    2.9      3 TRUE       TRUE      
##  7 A          Riego_cada_4_días 60_días    5.10     5 TRUE       TRUE      
##  8 B          Riego_cada_2_días 20_días    8.97     1 TRUE       TRUE      
##  9 B          Riego_cada_2_días 20_días    7.51     2 TRUE       TRUE      
## 10 B          Riego_cada_4_días 40_días   10.7      1 TRUE       FALSE     
## # ℹ 14 more rows
#aplicar logaritmos para evitar datos atipicos en la biomasa
datos_larg3$log_biomasa= log(datos_larg3$biomasa)
datos_larg3 %>%
  group_by(dosis_fert, tiempo) %>%
  identify_outliers(log_biomasa)
## # A tibble: 10 × 8
##    dosis_fert tiempo  riego      biomasa    id log_biomasa is.outlier is.extreme
##    <chr>      <fct>   <chr>        <dbl> <int>       <dbl> <lgl>      <lgl>     
##  1 A          05_días Riego_cad…    2.53     2       0.928 TRUE       FALSE     
##  2 A          15_días Riego_cad…    2.67     2       0.982 TRUE       FALSE     
##  3 A          20_días Riego_cad…    2.72     3       1.00  TRUE       FALSE     
##  4 A          60_días Riego_cad…    5.10     5       1.63  TRUE       FALSE     
##  5 B          30_días Riego_cad…    7.41     3       2.00  TRUE       FALSE     
##  6 B          45_días Riego_cad…    8.79     4       2.17  TRUE       FALSE     
##  7 B          50_días Riego_cad…    8.36     1       2.12  TRUE       FALSE     
##  8 C          15_días Riego_cad…    7.86     2       2.06  TRUE       FALSE     
##  9 C          45_días Riego_cad…   24.6      2       3.20  TRUE       FALSE     
## 10 C          45_días Riego_cad…   15.7      4       2.75  TRUE       FALSE
ggqqplot(datos_larg3, "log_biomasa", ggtheme = theme_bw()) +
  facet_grid(dosis_fert + riego ~ tiempo, labeller = "label_both")
## Warning: Removed 30 rows containing non-finite outside the scale range
## (`stat_qq()`).
## Warning: Removed 30 rows containing non-finite outside the scale range
## (`stat_qq_line()`).
## Removed 30 rows containing non-finite outside the scale range
## (`stat_qq_line()`).

res.aov <- anova_test(
    data = datos_larg3, dv = biomasa, wid = id,
  within = c(dosis_fert, riego, tiempo)
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                    Effect DFn DFd        F        p p<.05   ges
## 1              dosis_fert   2   8 7143.187 9.81e-14     * 0.959
## 2                   riego   1   4    2.300 2.04e-01       0.004
## 3                  tiempo  12  48  203.856 7.20e-37     * 0.927
## 4        dosis_fert:riego   2   8    2.551 1.39e-01       0.005
## 5       dosis_fert:tiempo  24  96  127.220 3.44e-62     * 0.873
## 6            riego:tiempo  12  48    0.502 9.03e-01       0.034
## 7 dosis_fert:riego:tiempo  24  96    0.698 8.42e-01       0.036
#como el p valor de la interacción triple es de 0.842 no rechazamos la hipotesis nula, por lo cual, no habría interacción triple
#el p valor de la interacción riego tiempo, es de 0.903, por lo cual no rechazamos la hipotesis nula y no habría interacción doble entre tiempo y riego
#el p valor de la interacción dosis de fertilización y tiempo fue de 0.139, por lo cual no tenemos interacción entre dosis de fertilizante y riego
#Solo tendríamos interacción entre la dosis de fertilización y el tiempo, por lo cual sería necesario ver las medias de la biomasa en cada dosis de fertilización en cada tiempo
#Se puede interpretar el riego puesto que este no presentó ninguna interacción y podemos observar que el riego por si solo tampoco tuvo un efecto en la biomasa
# Convertir las columnas a factores si es necesario
datos_larg3$tiempo <- as.factor(datos_larg3$tiempo)
datos_larg3$dosis_fert <- as.factor(datos_larg3$dosis_fert)
datos_larg3$riego <- as.factor(datos_larg3$riego)

# Modelo Lineal Generalizado Mixto
model <- lmer(biomasa ~ dosis_fert * tiempo + (1 | id), data = datos_larg3)

# Comparaciones por pares con letras de agrupamiento (Tukey)
comp <- glht(model, linfct = mcp(dosis_fert = "Tukey", tiempo = "Tukey"))
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found --
## default contrast might be inappropriate

## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found --
## default contrast might be inappropriate
summary(comp)
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = biomasa ~ dosis_fert * tiempo + (1 | id), data = datos_larg3)
## 
## Linear Hypotheses:
##                                  Estimate Std. Error z value Pr(>|z|)    
## dosis_fert: B - A == 0         -1.951e-13  4.917e-01   0.000   1.0000    
## dosis_fert: C - A == 0         -1.920e-13  4.917e-01   0.000   1.0000    
## dosis_fert: C - B == 0          3.116e-15  4.917e-01   0.000   1.0000    
## tiempo: 05_días - 0_días == 0   2.842e+00  4.917e-01   5.779    <0.01 ***
## tiempo: 10_días - 0_días == 0   3.121e+00  4.917e-01   6.347    <0.01 ***
## tiempo: 15_días - 0_días == 0   3.185e+00  4.917e-01   6.478    <0.01 ***
## tiempo: 20_días - 0_días == 0   3.295e+00  4.917e-01   6.702    <0.01 ***
## tiempo: 25_días - 0_días == 0   3.427e+00  4.917e-01   6.970    <0.01 ***
## tiempo: 30_días - 0_días == 0   3.686e+00  4.917e-01   7.496    <0.01 ***
## tiempo: 35_días - 0_días == 0   3.781e+00  4.917e-01   7.689    <0.01 ***
## tiempo: 40_días - 0_días == 0   4.320e+00  4.917e-01   8.785    <0.01 ***
## tiempo: 45_días - 0_días == 0   4.467e+00  4.917e-01   9.085    <0.01 ***
## tiempo: 50_días - 0_días == 0   4.659e+00  4.917e-01   9.476    <0.01 ***
## tiempo: 55_días - 0_días == 0   5.085e+00  4.917e-01  10.342    <0.01 ***
## tiempo: 60_días - 0_días == 0   5.757e+00  4.917e-01  11.708    <0.01 ***
## tiempo: 10_días - 05_días == 0  2.795e-01  4.917e-01   0.568   1.0000    
## tiempo: 15_días - 05_días == 0  3.437e-01  4.917e-01   0.699   1.0000    
## tiempo: 20_días - 05_días == 0  4.537e-01  4.917e-01   0.923   0.9997    
## tiempo: 25_días - 05_días == 0  5.854e-01  4.917e-01   1.191   0.9959    
## tiempo: 30_días - 05_días == 0  8.444e-01  4.917e-01   1.717   0.9039    
## tiempo: 35_días - 05_días == 0  9.390e-01  4.917e-01   1.910   0.8105    
## tiempo: 40_días - 05_días == 0  1.478e+00  4.917e-01   3.006   0.1258    
## tiempo: 45_días - 05_días == 0  1.626e+00  4.917e-01   3.306   0.0533 .  
## tiempo: 50_días - 05_días == 0  1.818e+00  4.917e-01   3.697   0.0143 *  
## tiempo: 55_días - 05_días == 0  2.244e+00  4.917e-01   4.563    <0.01 ***
## tiempo: 60_días - 05_días == 0  2.916e+00  4.917e-01   5.929    <0.01 ***
## tiempo: 15_días - 10_días == 0  6.420e-02  4.917e-01   0.131   1.0000    
## tiempo: 20_días - 10_días == 0  1.742e-01  4.917e-01   0.354   1.0000    
## tiempo: 25_días - 10_días == 0  3.059e-01  4.917e-01   0.622   1.0000    
## tiempo: 30_días - 10_días == 0  5.649e-01  4.917e-01   1.149   0.9971    
## tiempo: 35_días - 10_días == 0  6.595e-01  4.917e-01   1.341   0.9872    
## tiempo: 40_días - 10_días == 0  1.199e+00  4.917e-01   2.438   0.4344    
## tiempo: 45_días - 10_días == 0  1.346e+00  4.917e-01   2.737   0.2429    
## tiempo: 50_días - 10_días == 0  1.538e+00  4.917e-01   3.128   0.0897 .  
## tiempo: 55_días - 10_días == 0  1.964e+00  4.917e-01   3.995    <0.01 ** 
## tiempo: 60_días - 10_días == 0  2.636e+00  4.917e-01   5.361    <0.01 ***
## tiempo: 20_días - 15_días == 0  1.100e-01  4.917e-01   0.224   1.0000    
## tiempo: 25_días - 15_días == 0  2.417e-01  4.917e-01   0.492   1.0000    
## tiempo: 30_días - 15_días == 0  5.007e-01  4.917e-01   1.018   0.9992    
## tiempo: 35_días - 15_días == 0  5.953e-01  4.917e-01   1.211   0.9952    
## tiempo: 40_días - 15_días == 0  1.134e+00  4.917e-01   2.307   0.5313    
## tiempo: 45_días - 15_días == 0  1.282e+00  4.917e-01   2.607   0.3189    
## tiempo: 50_días - 15_días == 0  1.474e+00  4.917e-01   2.998   0.1278    
## tiempo: 55_días - 15_días == 0  1.900e+00  4.917e-01   3.864    <0.01 ** 
## tiempo: 60_días - 15_días == 0  2.572e+00  4.917e-01   5.230    <0.01 ***
## tiempo: 25_días - 20_días == 0  1.317e-01  4.917e-01   0.268   1.0000    
## tiempo: 30_días - 20_días == 0  3.907e-01  4.917e-01   0.795   1.0000    
## tiempo: 35_días - 20_días == 0  4.853e-01  4.917e-01   0.987   0.9994    
## tiempo: 40_días - 20_días == 0  1.024e+00  4.917e-01   2.084   0.6969    
## tiempo: 45_días - 20_días == 0  1.172e+00  4.917e-01   2.383   0.4747    
## tiempo: 50_días - 20_días == 0  1.364e+00  4.917e-01   2.774   0.2240    
## tiempo: 55_días - 20_días == 0  1.790e+00  4.917e-01   3.641   0.0179 *  
## tiempo: 60_días - 20_días == 0  2.462e+00  4.917e-01   5.007    <0.01 ***
## tiempo: 30_días - 25_días == 0  2.590e-01  4.917e-01   0.527   1.0000    
## tiempo: 35_días - 25_días == 0  3.536e-01  4.917e-01   0.719   1.0000    
## tiempo: 40_días - 25_días == 0  8.928e-01  4.917e-01   1.816   0.8611    
## tiempo: 45_días - 25_días == 0  1.040e+00  4.917e-01   2.115   0.6744    
## tiempo: 50_días - 25_días == 0  1.232e+00  4.917e-01   2.506   0.3865    
## tiempo: 55_días - 25_días == 0  1.659e+00  4.917e-01   3.373   0.0432 *  
## tiempo: 60_días - 25_días == 0  2.330e+00  4.917e-01   4.739    <0.01 ***
## tiempo: 35_días - 30_días == 0  9.460e-02  4.917e-01   0.192   1.0000    
## tiempo: 40_días - 30_días == 0  6.338e-01  4.917e-01   1.289   0.9911    
## tiempo: 45_días - 30_días == 0  7.811e-01  4.917e-01   1.589   0.9459    
## tiempo: 50_días - 30_días == 0  9.734e-01  4.917e-01   1.980   0.7679    
## tiempo: 55_días - 30_días == 0  1.400e+00  4.917e-01   2.846   0.1879    
## tiempo: 60_días - 30_días == 0  2.071e+00  4.917e-01   4.212    <0.01 ** 
## tiempo: 40_días - 35_días == 0  5.392e-01  4.917e-01   1.097   0.9982    
## tiempo: 45_días - 35_días == 0  6.865e-01  4.917e-01   1.396   0.9815    
## tiempo: 50_días - 35_días == 0  8.788e-01  4.917e-01   1.787   0.8742    
## tiempo: 55_días - 35_días == 0  1.305e+00  4.917e-01   2.654   0.2899    
## tiempo: 60_días - 35_días == 0  1.977e+00  4.917e-01   4.020    <0.01 ** 
## tiempo: 45_días - 40_días == 0  1.473e-01  4.917e-01   0.300   1.0000    
## tiempo: 50_días - 40_días == 0  3.396e-01  4.917e-01   0.691   1.0000    
## tiempo: 55_días - 40_días == 0  7.657e-01  4.917e-01   1.557   0.9539    
## tiempo: 60_días - 40_días == 0  1.437e+00  4.917e-01   2.923   0.1555    
## tiempo: 50_días - 45_días == 0  1.923e-01  4.917e-01   0.391   1.0000    
## tiempo: 55_días - 45_días == 0  6.184e-01  4.917e-01   1.258   0.9930    
## tiempo: 60_días - 45_días == 0  1.290e+00  4.917e-01   2.624   0.3082    
## tiempo: 55_días - 50_días == 0  4.261e-01  4.917e-01   0.867   0.9999    
## tiempo: 60_días - 50_días == 0  1.098e+00  4.917e-01   2.233   0.5877    
## tiempo: 60_días - 55_días == 0  6.717e-01  4.917e-01   1.366   0.9848    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
# Convertir las columnas a factores si es necesario
datos_larg3$tiempo <- as.factor(datos_larg3$tiempo)
datos_larg3$riego <- as.factor(datos_larg3$riego)

# Modelo Lineal Generalizado Mixto
model <- lmer(biomasa ~ riego * tiempo + (1 | id), data = datos_larg3)
## boundary (singular) fit: see help('isSingular')
# Comparaciones por pares con letras de agrupamiento (Tukey)
comp <- glht(model, linfct = mcp(riego = "Tukey"))
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found --
## default contrast might be inappropriate
summary(comp)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = biomasa ~ riego * tiempo + (1 | id), data = datos_larg3)
## 
## Linear Hypotheses:
##                                              Estimate Std. Error z value
## Riego_cada_4_días - Riego_cada_2_días == 0 -2.218e-13  2.128e+00       0
##                                            Pr(>|z|)
## Riego_cada_4_días - Riego_cada_2_días == 0        1
## (Adjusted p values reported -- single-step method)
# Obtener las letras de agrupamiento
cld_result <- cld(comp, level = 0.05)
print(cld_result)
## Riego_cada_2_días Riego_cada_4_días 
##               "a"               "a"
bxp <- ggboxplot(
  datos_larg3, x = "dosis_fert", y = "biomasa",
  color = "tiempo", palette = "Zissou1",
  facet.by = "riego", short.panel.labs = FALSE
  )
bxp