setwd("G:/Mi unidad/Agrosavia/SOCODEVI/Informe final")
sensorial <- read.table("senso.csv", header=T, sep=";")
sensorial$Vereda<-as.factor(sensorial$Vereda)
attach(sensorial)
# Packages
library(vegan)
## Warning: package 'vegan' was built under R version 4.4.3
## Cargando paquete requerido: permute
## Warning: package 'permute' was built under R version 4.4.3
## Cargando paquete requerido: lattice
library(ggplot2)
library(concaveman)
## Warning: package 'concaveman' was built under R version 4.4.3
library(ggforce)
## Warning: package 'ggforce' was built under R version 4.4.3
# Generando la matriz general de perfil sensorial - extrayendo columnas con la información de sabor
sens = sensorial[,4:ncol(sensorial)]
# Convirtiendo las columnas extraidas en una matriz
m_sens=as.matrix(sens)
# Obteniendo las distancias Bray-Curtis de la matriz
set.seed(123)
nmds = metaMDS(m_sens, distance = "bray")
## Run 0 stress 0.04670144
## Run 1 stress 0.06446481
## Run 2 stress 0.06446472
## Run 3 stress 0.04670143
## ... New best solution
## ... Procrustes: rmse 8.259902e-05 max resid 0.000233247
## ... Similar to previous best
## Run 4 stress 0.04670149
## ... Procrustes: rmse 6.741703e-05 max resid 0.0001896569
## ... Similar to previous best
## Run 5 stress 0.04389918
## ... New best solution
## ... Procrustes: rmse 0.05534243 max resid 0.1312288
## Run 6 stress 0.04389919
## ... Procrustes: rmse 4.054449e-05 max resid 6.058739e-05
## ... Similar to previous best
## Run 7 stress 0.04670143
## Run 8 stress 0.04389909
## ... New best solution
## ... Procrustes: rmse 0.0002490384 max resid 0.000451389
## ... Similar to previous best
## Run 9 stress 0.06446471
## Run 10 stress 0.0644647
## Run 11 stress 0.04670145
## Run 12 stress 0.04670143
## Run 13 stress 0.06446477
## Run 14 stress 0.3071451
## Run 15 stress 0.2937874
## Run 16 stress 0.06446472
## Run 17 stress 0.04389921
## ... Procrustes: rmse 0.0002275235 max resid 0.0004090331
## ... Similar to previous best
## Run 18 stress 0.04670147
## Run 19 stress 0.04670145
## Run 20 stress 0.04389919
## ... Procrustes: rmse 0.0002429038 max resid 0.0004401424
## ... Similar to previous best
## *** Best solution repeated 3 times
nmds
##
## Call:
## metaMDS(comm = m_sens, distance = "bray")
##
## global Multidimensional Scaling using monoMDS
##
## Data: m_sens
## Distance: bray
##
## Dimensions: 2
## Stress: 0.04389909
## Stress type 1, weak ties
## Best solution was repeated 3 times in 20 tries
## The best solution was from try 8 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'm_sens'
# Generando el gráfico Non-Metric multidimensional scaling
plot(nmds)
# extrayendo los puntajes NMDS (coordenadas x and y)
data.scores = as.data.frame(scores(nmds))
# adicionando columnas al data frame
data.scores$Vereda = sensorial$Vereda
head(data.scores)
## sites.NMDS1 sites.NMDS2 species.NMDS1 species.NMDS2 Vereda
## 1 -0.07274958 -0.03993710 0.11549862 0.009949239 Caño Baul
## 2 -0.10384857 -0.08367856 0.14308485 0.014095816 Caño Baul
## 3 -0.15249141 -0.11021291 0.07435605 -0.008475395 Caño Baul
## 4 0.02373042 -0.03095687 -0.07442190 0.030805700 Caño Baul
## 5 -0.02466335 0.10800730 -0.16444339 -0.025591914 Mataredonda
## 6 -0.18570148 0.01671731 -0.04304244 -0.045680196 Mataredonda
# Creando escala de color para graficar
cbp1 <- c("#999999", "#E69F00", "#56B4E9", "#009E73","#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# Gráfica del análisis
xx = ggplot(data.scores, aes(x = sites.NMDS1, y = sites.NMDS2)) +
geom_point(size = 4, aes(shape = Vereda)) +
geom_mark_hull(concavity = 5,expand=0,radius=0,aes(fill=Vereda))+
theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"),
axis.text.x = element_text(colour = "black", face = "bold", size = 12),
legend.text = element_text(size = 12, face ="bold", colour ="black"),
legend.position = "right", axis.title.y = element_text(face = "bold", size = 14),
axis.title.x = element_text(face = "bold", size = 14, colour = "black"),
legend.title = element_text(size = 14, colour = "black", face = "bold"),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
legend.key=element_blank()) +
labs(x = "NMDS1", y = "NMDS2", shape ="Vereda", fill="Vereda") +
scale_colour_manual(values = cbp1)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
xx + scale_fill_discrete(name ="Vereda") + scale_shape_discrete(name="Vereda")+
geom_text(x=0.275, y=-0.07, label="Stress = 0.04") + xlim(-0.3, 0.4)
## Probando si el perfil del sabor es estadísticamente distinto basado en los niveles de Vereda
# Análisis de similaridad identificando las diferencias en el perfil de sabor para Vereda
anoVereda <- anosim(m_sens, Vereda, distance = "bray", permutations = 9999)
anoVereda
##
## Call:
## anosim(x = m_sens, grouping = Vereda, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.3056
## Significance: 0.0363
##
## Permutation: free
## Number of permutations: 9999
# PostHoc análisis de las diferencias cuantitativas del perfil de sabor entre los niveles de Vereda
library(pairwiseAdonis)
## Cargando paquete requerido: cluster
pairwise.adonis2(m_sens ~ Vereda, data = sensorial)
## $parent_call
## [1] "m_sens ~ Vereda , strata = Null , permutations 999"
##
## $`Caño Baul_vs_Mataredonda`
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.015581 0.20601 1.5568 0.249
## Residual 6 0.060051 0.79399
## Total 7 0.075632 1.00000
##
## $`Caño Baul_vs_Mira Valle`
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.036926 0.48305 5.6065 0.056 .
## Residual 6 0.039517 0.51695
## Total 7 0.076443 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Mataredonda_vs_Mira Valle`
## Df SumOfSqs R2 F Pr(>F)
## Model 1 0.030657 0.29605 2.5233 0.157
## Residual 6 0.072898 0.70395
## Total 7 0.103556 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"
## Análisis de sabor indicador: este análisis permite identificar el sabor que es encontrado más a menudo
## Se comparan los niveles de Vereda
library(indicspecies)
## Warning: package 'indicspecies' was built under R version 4.4.3
## Análisis de sabor indicador
invgen <- multipatt(m_sens, Vereda, func = "r.g", control = how(nperm=9999), duleg = T)
summary(invgen, alpha = 1)
##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 1
##
## Total number of species: 12
## Selected number of species: 12
## Number of species associated to 1 group: 12
## Number of species associated to 2 groups: 0
##
## List of species associated to each combination:
##
## Group Caño Baul #sps. 6
## stat p.value
## Madera 0.760 0.0257 *
## Frutal 0.625 0.0893 .
## Floral 0.568 0.1576
## Vegetal 0.490 0.3343
## Dulce 0.369 0.5944
## Ácido 0.263 0.8549
##
## Group Mataredonda #sps. 4
## stat p.value
## Lácteo 0.700 0.066 .
## Especiado 0.632 0.264
## Frutos.secos 0.378 0.652
## Cacao 0.300 0.627
##
## Group Mira Valle #sps. 2
## stat p.value
## Astringente 0.547 0.206
## Amargo 0.309 0.698
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1