Piña
David Arango-Londoño
6/3/2021
Diseño Experimento Piña
Se realiza la importación de los datos y se calcula la variable elevación como factor de control en el diseño de experimentos.
##
##Importar Datos y librerias
##
require(readxl)
require(raster)
require(ggplot2)
require(plotly)
require(emmeans)
require(multcomp)
require(tidyverse)
require(ggpubr)
require(rstatix)
require(lmerTest)
require(Rmisc)
require(agricolae)
require(psych)
datos = read_excel("~/Desktop/piña/MATRIZ DE DATOS- DISEÑO DE BCA - JOHN CANACUAN.xlsx", sheet = "R2")Diseño para Variable Respuesta - Rendimiento
Exploratorio
datos$Tratamiento=as.factor(datos$Tratamiento)
g1=ggplot(datos,aes(x=Tratamiento,y=Rendimiento,fill=Tratamiento))+geom_boxplot()+theme_bw()+ scale_fill_grey()+ylab("Rendimiento (ton/ha)")
g1
res=describe.by(x =datos$Rendimiento,group = datos$Tratamiento,mat = TRUE)## Warning: describe.by is deprecated. Please use the describeBy functionres[,c(2,5:7)]## group1 mean sd median
## X11 1 85.13333 4.985352 87.02
## X12 2 91.00667 4.417662 89.94
## X13 3 93.74667 9.496759 94.89
## X14 4 83.74333 6.006166 83.88
## X15 5 92.17667 8.986169 96.56Anova
datos$Tratamiento=as.factor(datos$Tratamiento)
anova=aov(Rendimiento~Tratamiento,data=datos)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 4 237.3 59.33 1.18 0.377
## Residuals 10 502.8 50.28PostAnova - Tratamiento
bxp = ggboxplot(
datos, x = "Tratamiento", y = "Rendimiento",
title= "Postanova",
xlab ="Tratamiento", ylab= "Rendimiento",
color = "Tratamiento", palette = "lancet"
)
pwc = datos %>%
pairwise_t_test(
Rendimiento ~ Tratamiento, paired = FALSE,
p.adjust.method = "bonferroni"
)
pwc = pwc %>% add_xy_position(x = "Tratamiento")
bxp +
stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = paste("Anova, ","F(4)", " = ", "1.179 , ", "p", " = ", "0.377"),
caption = get_pwc_label(pwc)
)
compara=LSD.test(anova,"Tratamiento")
compara$groups## Rendimiento groups
## 3 93.74667 a
## 5 92.17667 a
## 2 91.00667 a
## 1 85.13333 a
## 4 83.74333 apwc[,c(2,3,6)]## # A tibble: 10 x 3
## group1 group2 p
## <chr> <chr> <dbl>
## 1 1 2 0.334
## 2 1 3 0.168
## 3 2 3 0.646
## 4 1 4 0.815
## 5 2 4 0.238
## 6 3 4 0.115
## 7 1 5 0.252
## 8 2 5 0.844
## 9 3 5 0.792
## 10 4 5 0.176Supuesto de Normalidad
datos %>%
group_by(Tratamiento) %>%
shapiro_test(Rendimiento)## # A tibble: 5 x 4
## Tratamiento variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 1 Rendimiento 0.893 0.362
## 2 2 Rendimiento 0.956 0.598
## 3 3 Rendimiento 0.989 0.801
## 4 4 Rendimiento 1.00 0.962
## 5 5 Rendimiento 0.822 0.167Diseño para Variable Respuesta - Grados Brix
Exploratorio
datos$Tratamiento=as.factor(datos$Tratamiento)
g1=ggplot(datos,aes(x=Tratamiento,y=Grados_brix,fill=Tratamiento))+geom_boxplot()+theme_bw()+ scale_fill_grey()+ylab("Grados brix (°Bx)")
g1
res=describe.by(x =datos$Grados_brix,group = datos$Tratamiento,mat = TRUE)## Warning: describe.by is deprecated. Please use the describeBy functionres[,c(2,5:7)]## group1 mean sd median
## X11 1 14.80000 0.6082763 14.5
## X12 2 14.10000 0.1000000 14.1
## X13 3 14.16667 0.5686241 14.0
## X14 4 14.60000 0.2000000 14.6
## X15 5 14.10000 0.1732051 14.0Anova
anova=aov(Grados_brix~Tratamiento,data=datos)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 4 1.271 0.3177 2.054 0.162
## Residuals 10 1.547 0.1547PostAnova - Tratamiento
bxp = ggboxplot(
datos, x = "Tratamiento", y = "Grados_brix",
title= "Postanova",
xlab ="Tratamiento", ylab= "Grados_brix",
color = "Tratamiento", palette = "lancet"
)
pwc = datos %>%
pairwise_t_test(
Grados_brix ~ Tratamiento, paired = TRUE,
p.adjust.method = "bonferroni"
)
pwc = pwc %>% add_xy_position(x = "Tratamiento")
bxp +
stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = paste("Anova, ","F(4)", " = ", "2.054 , ", "p", " = ", "0.162"),
caption = get_pwc_label(pwc)
)
require(agricolae)
compara=LSD.test(anova,"Tratamiento")
compara$group## Grados_brix groups
## 1 14.80000 a
## 4 14.60000 a
## 3 14.16667 a
## 2 14.10000 a
## 5 14.10000 acompara=LSD.test(anova,"Tratamiento")
compara$groups## Grados_brix groups
## 1 14.80000 a
## 4 14.60000 a
## 3 14.16667 a
## 2 14.10000 a
## 5 14.10000 apwc[,c(2,3,8)]## # A tibble: 10 x 3
## group1 group2 p
## <chr> <chr> <dbl>
## 1 1 2 0.225
## 2 1 3 0.348
## 3 1 4 0.51
## 4 1 5 0.109
## 5 2 3 0.866
## 6 2 4 0.102
## 7 2 5 1
## 8 3 4 0.306
## 9 3 5 0.872
## 10 4 5 0.013Supuesto de Normalidad
datos %>%
group_by(Tratamiento) %>%
shapiro_test(Grados_brix)## # A tibble: 5 x 4
## Tratamiento variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 1 Grados_brix 0.818 0.157
## 2 2 Grados_brix 1 1.00
## 3 3 Grados_brix 0.936 0.510
## 4 4 Grados_brix 1 1.00
## 5 5 Grados_brix 0.75 0Diseño para Variable Respuesta - Margen de Ganancia
Exploratorio
datos$Tratamiento=as.factor(datos$Tratamiento)
g1=ggplot(datos,aes(x=Tratamiento,y=Margen_de_ganancia,fill=Tratamiento))+geom_boxplot()+theme_bw()+ scale_fill_grey()+ylab("Margen de ganancia por ha (%)")
g1
res=describe.by(x =datos$Margen_de_ganancia,group = datos$Tratamiento,mat = TRUE)## Warning: describe.by is deprecated. Please use the describeBy functionres[,c(2,5:7)]## group1 mean sd median
## X11 1 45.40000 3.251154 46.7
## X12 2 45.86667 2.596793 45.3
## X13 3 46.66667 5.522982 47.7
## X14 4 43.33333 4.106499 43.6
## X15 5 57.36667 4.398106 59.6Anova
anova=aov(Margen_de_ganancia~Tratamiento,data=datos)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 4 366.7 91.67 5.455 0.0136 *
## Residuals 10 168.0 16.80
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PostAnova - Tratamiento
bxp = ggboxplot(
datos, x = "Tratamiento", y = "Margen_de_ganancia",
title= "Postanova",
xlab ="Tratamiento", ylab= "Margen de Ganancia",
color = "Tratamiento", palette = "lancet"
)
pwc = datos %>%
pairwise_t_test(
Margen_de_ganancia ~ Tratamiento, paired = FALSE,
p.adjust.method = "bonferroni"
)
pwc = pwc %>% add_xy_position(x = "Tratamiento")
bxp +
stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = paste("Anova, ","F(4)", " = ", "6.223 , ", "p", " = ", "0.00883"),
caption = get_pwc_label(pwc)
)
compara=LSD.test(anova,"Tratamiento")
compara$groups## Margen_de_ganancia groups
## 5 57.36667 a
## 3 46.66667 b
## 2 45.86667 b
## 1 45.40000 b
## 4 43.33333 bpwc[,c(2,3,6)]## # A tibble: 10 x 3
## group1 group2 p
## <chr> <chr> <dbl>
## 1 1 2 0.892
## 2 1 3 0.713
## 3 2 3 0.816
## 4 1 4 0.551
## 5 2 4 0.467
## 6 3 4 0.343
## 7 1 5 0.00505
## 8 2 5 0.00638
## 9 3 5 0.00954
## 10 4 5 0.00185Supuesto de Normalidad
datos %>%
group_by(Tratamiento) %>%
shapiro_test(Margen_de_ganancia)## # A tibble: 5 x 4
## Tratamiento variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 1 Margen_de_ganancia 0.880 0.325
## 2 2 Margen_de_ganancia 0.964 0.637
## 3 3 Margen_de_ganancia 0.974 0.689
## 4 4 Margen_de_ganancia 0.997 0.893
## 5 5 Margen_de_ganancia 0.807 0.130Diseño para Variable Respuesta - Costo/Beneficio
Exploratorio
datos$Tratamiento=as.factor(datos$Tratamiento)
g1=ggplot(datos,aes(x=Tratamiento,y=Costo_beneficio,fill=Tratamiento))+geom_boxplot()+theme_bw()+ scale_fill_grey()+ylab("Relación costo/beneficio por planta")
g1
res=describe.by(x =datos$Costo_beneficio,group = datos$Tratamiento,mat = TRUE)## Warning: describe.by is deprecated. Please use the describeBy functionres[,c(2,5:7)]## group1 mean sd median
## X11 1 1.836667 0.11150486 1.88
## X12 2 1.850000 0.09165151 1.83
## X13 3 1.890000 0.19078784 1.91
## X14 4 1.770000 0.13000000 1.77
## X15 5 2.360000 0.22605309 2.47Anova
anova=aov(Costo_beneficio~Tratamiento,data=datos)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 4 0.6797 0.16993 6.784 0.00659 **
## Residuals 10 0.2505 0.02505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PostAnova - Tratamiento
bxp = ggboxplot(
datos, x = "Tratamiento", y = "Costo_beneficio",
title= "Postanova",
xlab ="Tratamiento", ylab= "Costo Beneficio",
color = "Tratamiento", palette = "lancet"
)
pwc = datos %>%
pairwise_t_test(
Costo_beneficio ~ Tratamiento, paired = FALSE,
p.adjust.method = "bonferroni"
)
pwc = pwc %>% add_xy_position(x = "Tratamiento")
bxp +
stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = paste("Anova, ","F(4)", " = ", "6.223 , ", "p", " = ", "0.00883"),
caption = get_pwc_label(pwc)
)
compara=LSD.test(anova,"Tratamiento")
compara$groups## Costo_beneficio groups
## 5 2.360000 a
## 3 1.890000 b
## 2 1.850000 b
## 1 1.836667 b
## 4 1.770000 bpwc[,c(2,3,6)]## # A tibble: 10 x 3
## group1 group2 p
## <chr> <chr> <dbl>
## 1 1 2 0.92
## 2 1 3 0.689
## 3 2 3 0.763
## 4 1 4 0.617
## 5 2 4 0.55
## 6 3 4 0.375
## 7 1 5 0.00232
## 8 2 5 0.00274
## 9 3 5 0.00456
## 10 4 5 0.00103Supuesto de Normalidad
datos %>%
group_by(Tratamiento) %>%
shapiro_test(Costo_beneficio)## # A tibble: 5 x 4
## Tratamiento variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 1 Costo_beneficio 0.887 0.344
## 2 2 Costo_beneficio 0.964 0.637
## 3 3 Costo_beneficio 0.992 0.826
## 4 4 Costo_beneficio 1 1.00
## 5 5 Costo_beneficio 0.822 0.169Diseño para Variable Respuesta - Costo_de_producción_por_kilo
Exploratorio
datos$Tratamiento=as.factor(datos$Tratamiento)
g1=ggplot(datos,aes(x=Tratamiento,y=Costo_de_producción_por_kilo,fill=Tratamiento))+geom_boxplot()+theme_bw()+ scale_fill_grey()+ylab("Costo de producción por kg ($)")
g1
res=describe.by(x =datos$Costo_de_producción_por_kilo,group = datos$Tratamiento,mat = TRUE)## Warning: describe.by is deprecated. Please use the describeBy functionres[,c(2,5:7)]## group1 mean sd median
## X11 1 572.0233 34.49689 558.30
## X12 2 567.2867 27.12129 573.13
## X13 3 558.5067 57.87800 547.94
## X14 4 593.6433 42.79517 590.62
## X15 5 446.7367 46.06216 423.61Anova
anova=aov(Costo_de_producción_por_kilo~Tratamiento,data=datos)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 4 40189 10047 5.444 0.0137 *
## Residuals 10 18457 1846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PostAnova - Tratamiento
bxp = ggboxplot(
datos, x = "Tratamiento", y = "Costo_de_producción_por_kilo",
title= "Postanova",
xlab ="Tratamiento", ylab= "Costo de producción por kilo",
color = "Tratamiento", palette = "lancet"
)
pwc = datos %>%
pairwise_t_test(
Costo_de_producción_por_kilo ~ Tratamiento, paired = FALSE,
p.adjust.method = "bonferroni"
)
pwc = pwc %>% add_xy_position(x = "Tratamiento")
bxp +
stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = paste("Anova, ","F(4)", " = ", "6.223 , ", "p", " = ", "0.00883"),
caption = get_pwc_label(pwc)
)
compara=LSD.test(anova,"Tratamiento")
compara$groups## Costo_de_producción_por_kilo groups
## 4 593.6433 a
## 1 572.0233 a
## 2 567.2867 a
## 3 558.5067 a
## 5 446.7367 bpwc[,c(2,3,6)]## # A tibble: 10 x 3
## group1 group2 p
## <chr> <chr> <dbl>
## 1 1 2 0.895
## 2 1 3 0.708
## 3 2 3 0.807
## 4 1 4 0.551
## 5 2 4 0.47
## 6 3 4 0.34
## 7 1 5 0.00508
## 8 2 5 0.00637
## 9 3 5 0.00971
## 10 4 5 0.00186Supuesto de Normalidad
datos %>%
group_by(Tratamiento) %>%
shapiro_test(Costo_de_producción_por_kilo)## # A tibble: 5 x 4
## Tratamiento variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 1 Costo_de_producción_por_kilo 0.881 0.328
## 2 2 Costo_de_producción_por_kilo 0.965 0.642
## 3 3 Costo_de_producción_por_kilo 0.975 0.697
## 4 4 Costo_de_producción_por_kilo 0.996 0.883
## 5 5 Costo_de_producción_por_kilo 0.811 0.141