FIP606-Proyecto: Análisis Estadistico

Published

July 3, 2024

Preparación Base de Datos

Cargar paquetes de análisis

Codigo

library(tidyverse)
library(gsheet)
library(cowplot)
library(patchwork)
library(ggthemes)
library(viridis)
library(epifitter)
library(ggplot2)
library("writexl")
library(nlme)
library(lme4)
library(DHARMa)
library(performance)
library(report)
library(emmeans)
library(multcompView)
library(multcomp)
library(corrplot)
library(see)
library(lubridate)
library(agridat)
library(cowplot)
library(agricolae)
library(sf)
library(lme4)
library(broom)
library(lattice)
library(car)
library(scales)
library(readxl)
library(dplyr)
library(knitr)
library(kableExtra)
library(easyanova)
library(tidyr)
library(PerformanceAnalytics)
library(magrittr)
library(car)
library(ggpubr)
library(rstatix)
library(reshape)
library(kableExtra)
library(formattable)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggh4x)
library(gridExtra)
library(stringr)
library(epiR)
library(ggridges)
library(RColorBrewer)
library(DT)
library(gsheet)

Preparación de la Base de datos

Codigo

dat<-read.csv2("DB_PAT104022.csv")

dat |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   Buttons = c('excel', "csv"))) |> 
                        formatRound(c('Inc_promedio','Def_calculada','Sev_total','Severidad','Sev_condicional'), 2)

Uso de la funcion select para convertir una tabla resumida

Codigo

CLR_ <-dat |> 
      select(Evaluacion,Time_,Time_E,Eva_E,Eva_,Parcela,Surco,Arbol_1,Arbol,Rep,Htotal,Inc_promedio,Def_calculada,Severidad,Sev_condicional) |> 
      filter((Eva_ >= 4 & Eva_ <= 7) | (Eva_ >= 9 & Eva_ <= 12))
  

CLR_ |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   Buttons = c('excel', "csv")))

ANOVA + Prueba de Supuestos

Incidencia

Codigo

# Crear el modelo ANOVA
aov_inc <- aov(Inc_promedio~ Parcela, data = CLR_)

# Resumen del modelo
summary(aov_inc)

              Df  Sum Sq Mean Sq F value Pr(>F)    
Parcela        2 1045239  522620    1015 <2e-16 ***
Residuals   5392 2776486     515                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
365 observations deleted due to missingness

homocedasticidad

Codigo

check_heteroscedasticity(aov_inc)

Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

Normalidad

Codigo

check_normality(aov_inc)

Warning: Non-normality of residuals detected (p < .001).

check

Codigo

plot(simulateResiduals(aov_inc))

Codigo

check_model(aov_inc)

No cumplen con los supuestos de normalidad y homocedasticidad, hay varias alternativas que puedes considerar para realizar análisis estadísticos apropiados.

1. Transformaciones de Datos
1. Pruebas No Paramétricas El análisis no paramétrico no asume una distribución específica para los datos y es más flexible en términos de los supuestos que hace.
1. Modelos Generalizados
1. Resampling Methods
1. Análisis Robustos
1. Análisis de Varianza no Paramétrico
Área Bajo la Curva del Progreso de la Enfermedad (AUDPC)

Area Bajo la curva Incidencia

Codigo

dat_1<-CLR_ |> 
  group_by(Evaluacion,Time_E,Parcela,Arbol) |> 
  summarize(Incidencia = mean(Inc_promedio,na.rm=TRUE))


dat_1 |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound('Incidencia', 2)

Codigo

dat_audps <- dat_1 |> 
          group_by(Evaluacion,Parcela,Arbol) |> 
          summarize(audps = AUDPS(Time_E,Incidencia,))


dat_audps |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound('audps', 2)

Codigo

t<-90

dat1_audpc <- dat_audps |> 
  mutate(audpc2 = audps / t)


dat1_audpc  %>%
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound(c ('audps','audpc2'), 2)

Codigo

dat1_audpc_f<-  dat1_audpc  |> 
                  group_by(Evaluacion,Parcela,Arbol) |> 
                  summarize(audpc = sum(audps,na.rm=TRUE),TPD=mean(audpc2,na.rm=TRUE))

dat1_audpc_f  |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound(c('audpc','TPD'), 2)

Codigo

dat1_audpc_f_<-  dat1_audpc  |> 
                  group_by(Parcela) |> 
                  summarize(audpc = mean(audps,na.rm=TRUE),TPD=mean(audpc2,na.rm=TRUE))

dat1_audpc_f_  |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound(c('audpc','TPD'), 2)

Codigo

dat1_audpc_f |>
            ggplot(aes(Parcela, audpc, fill = Parcela)) +
            geom_boxplot() +
            facet_wrap(~ Evaluacion)+
            theme_clean()+
              theme(axis.text = element_text(size = 8),
                  axis.text.x = element_text( size =10),
                  axis.title = element_text(size = 10),
                  strip.text.x=element_text(face="bold",size =10,margin=margin(1,0,1,0)),
                  plot.background = element_rect(colour = "white"),
                  legend.background = element_rect(colour = "white"),
                  legend.title = element_text( face="bold",size = 10),
                  legend.text = element_text( size = 10),
                  plot.title = element_text(size = 15,face="bold"),
                  plot.subtitle = element_text(size = 15),
                  legend.position = "none" ) +
              scale_fill_manual(values = c( "#669933","#FFCC66","#990000"))+ 
              labs(
              y = "% Defoliaci",
              title = "A. Boxplot ABCPE Incidencia",subtitle =""
            )

Codigo

# Test de Kruskal-Wallis
kruskal.test(audpc ~ Parcela, data = dat1_audpc_f)


    Kruskal-Wallis rank sum test

data:  audpc by Parcela
Kruskal-Wallis chi-squared = 220.4, df = 2, p-value < 2.2e-16

Interpretación del Test de Kruskal-Wallis: Este test no requiere la asunción de normalidad. Si el p-valor es menor que 0.05, hay diferencias significativas entre los grupos.

Codigo

kruskal(dat1_audpc_f$audpc, dat1_audpc_f$Parcela, console = TRUE)


Study: dat1_audpc_f$audpc ~ dat1_audpc_f$Parcela
Kruskal-Wallis test's
Ties or no Ties

Critical Value: 220.4003
Degrees of freedom: 2
Pvalue Chisq  : 0 

dat1_audpc_f$Parcela,  means of the ranks

               dat1_audpc_f.audpc   r
Tto Fungicida            109.8937 160
Tto Nuevo                284.4688 160
Tto Testigo SA           327.1375 160

Post Hoc Analysis

t-Student: 1.96495
Alpha    : 0.05
Minimum Significant Difference: 22.43448 

Treatments with the same letter are not significantly different.

               dat1_audpc_f$audpc groups
Tto Testigo SA           327.1375      a
Tto Nuevo                284.4688      b
Tto Fungicida            109.8937      c

Codigo

aov_inc2 <- aov(audpc~ Parcela, data = dat1_audpc_f)

aov_inc2

Call:
   aov(formula = audpc ~ Parcela, data = dat1_audpc_f)

Terms:
                   Parcela  Residuals
Sum of Squares  3490469797 3788596830
Deg. of Freedom          2        477

Residual standard error: 2818.253
Estimated effects may be unbalanced

Codigo

aov_inc_means<- cld(emmeans(aov_inc2, ~ Parcela),alpha = 0.05, Letters = LETTERS,reverse=F)

aov_inc_means|> 
      DT::datatable(
        extensions = 'Buttons', 
        options = list(dom = 'Bfrtip', 
                       buttons = c('excel', "csv"))) |> 
                            formatRound(c('emmean','SE','lower.CL','upper.CL'), 2)

Defoliación

Codigo

# Crear el modelo ANOVA
aov_def <- aov(Def_calculada~ Parcela, data = CLR_)

# Resumen del modelo
summary(aov_def)

              Df  Sum Sq Mean Sq F value Pr(>F)    
Parcela        2  252680  126340   445.5 <2e-16 ***
Residuals   5757 1632470     284                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

homocedasticidad

Codigo

check_heteroscedasticity(aov_def)

OK: Error variance appears to be homoscedastic (p = 0.316).

Normalidad

Codigo

check_normality(aov_def)

Warning: Non-normality of residuals detected (p = 0.006).

check

Codigo

plot(simulateResiduals(aov_def))

Codigo

check_model(aov_def)

Se cumplen con el supuestos de normalidad, pero no con la homocedasticidad, hay varias alternativas que puedes considerar para realizar análisis estadísticos apropiados.

1. Transformaciones de Datos
1. Pruebas No Paramétricas
1. Modelos Generalizados
1. Resampling Methods
1. Análisis Robustos
1. Análisis de Varianza no Paramétrico
1. Transformaciones de Datos
1. Pruebas No Paramétricas El análisis no paramétrico no asume una distribución específica para los datos y es más flexible en términos de los supuestos que hace.
1. Modelos Generalizados
1. Resampling Methods
1. Análisis Robustos
1. Análisis de Varianza no Paramétrico
Área Bajo la Curva del Progreso de la Enfermedad (AUDPC)

Area Bajo la curva Defoliación

Codigo

dat_1<-CLR_ |> 
  group_by(Evaluacion,Time_E,Parcela,Arbol) |> 
  summarize(Defoliacion = mean(Def_calculada,na.rm=TRUE))


dat_1 |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound('Defoliacion', 2)

Codigo

dat_audps <- dat_1 |> 
          group_by(Evaluacion,Parcela,Arbol) |> 
          summarize(audps = AUDPS(Time_E,Defoliacion))


dat_audps |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound('audps', 2)

Codigo

t<-90

dat1_audpc <- dat_audps |> 
  mutate(audpc2 = audps / t)


dat1_audpc  %>%
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound(c ('audps','audpc2'), 2)

Codigo

dat1_audpc_f<-  dat1_audpc  |> 
                  group_by(Evaluacion,Parcela,Arbol) |> 
                  summarize(audpc = sum(audps,na.rm=TRUE),TPD=mean(audpc2,na.rm=TRUE))

dat1_audpc_f  |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound(c('audpc','TPD'), 2)

Codigo

dat1_audpc_f_<-  dat1_audpc  |> 
                  group_by(Parcela) |> 
                  summarize(audpc = mean(audps,na.rm=TRUE),TPD=mean(audpc2,na.rm=TRUE))

dat1_audpc_f_  |> 
  DT::datatable(
    extensions = 'Buttons', 
    options = list(dom = 'Bfrtip', 
                   buttons = c('excel', "csv"))) |> 
                        formatRound(c('audpc','TPD'), 2)

Codigo

dat1_audpc_f |>
            ggplot(aes(Parcela, audpc, fill = Parcela)) +
            geom_boxplot() +
            facet_wrap(~ Evaluacion)+
            theme_clean()+
              theme(axis.text = element_text(size = 8),
                  axis.text.x = element_text( size =10),
                  axis.title = element_text(size = 10),
                  strip.text.x=element_text(face="bold",size =10,margin=margin(1,0,1,0)),
                  plot.background = element_rect(colour = "white"),
                  legend.background = element_rect(colour = "white"),
                  legend.title = element_text( face="bold",size = 10),
                  legend.text = element_text( size = 10),
                  plot.title = element_text(size = 15,face="bold"),
                  plot.subtitle = element_text(size = 15),
                  legend.position = "none" ) +
              scale_fill_manual(values = c( "#669933","#FFCC66","#990000"))+ 
              labs(
              y = "% Defoliaci",
              title = "B. Boxplot ABCPE Defoliación",subtitle =""
            )

Codigo

# Test de Kruskal-Wallis
kruskal.test(audpc ~ Parcela, data = dat1_audpc_f)


    Kruskal-Wallis rank sum test

data:  audpc by Parcela
Kruskal-Wallis chi-squared = 220.4, df = 2, p-value < 2.2e-16

Interpretación del Test de Kruskal-Wallis: Este test no requiere la asunción de normalidad. Si el p-valor es menor que 0.05, hay diferencias significativas entre los grupos.

Codigo

kruskal(dat1_audpc_f$audpc, dat1_audpc_f$Parcela, console = TRUE)


Study: dat1_audpc_f$audpc ~ dat1_audpc_f$Parcela
Kruskal-Wallis test's
Ties or no Ties

Critical Value: 220.4003
Degrees of freedom: 2
Pvalue Chisq  : 0 

dat1_audpc_f$Parcela,  means of the ranks

               dat1_audpc_f.audpc   r
Tto Fungicida            109.8937 160
Tto Nuevo                284.4688 160
Tto Testigo SA           327.1375 160

Post Hoc Analysis

t-Student: 1.96495
Alpha    : 0.05
Minimum Significant Difference: 22.43448 

Treatments with the same letter are not significantly different.

               dat1_audpc_f$audpc groups
Tto Testigo SA           327.1375      a
Tto Nuevo                284.4688      b
Tto Fungicida            109.8937      c

Codigo

aov_inc2 <- aov(audpc~ Parcela, data = dat1_audpc_f)

aov_inc2

Call:
   aov(formula = audpc ~ Parcela, data = dat1_audpc_f)

Terms:
                   Parcela  Residuals
Sum of Squares  3490469797 3788596830
Deg. of Freedom          2        477

Residual standard error: 2818.253
Estimated effects may be unbalanced

Codigo

aov_inc_means<- cld(emmeans(aov_inc2, ~ Parcela),alpha = 0.05, Letters = LETTERS,reverse=F)

aov_inc_means|> 
      DT::datatable(
        extensions = 'Buttons', 
        options = list(dom = 'Bfrtip', 
                       buttons = c('excel', "csv"))) |> 
                            formatRound(c('emmean','SE','lower.CL','upper.CL'), 2)

いいですか、私たちの神は主おひとりです。

Preparación Base de Datos

Cargar paquetes de análisis

Preparación de la Base de datos

ANOVA + Prueba de Supuestos

Incidencia

homocedasticidad

Normalidad

check

Área Bajo la Curva del Progreso de la Enfermedad (AUDPC)

Area Bajo la curva Incidencia

Defoliación

homocedasticidad

Normalidad

check

Área Bajo la Curva del Progreso de la Enfermedad (AUDPC)

Area Bajo la curva Defoliación