Analisis de datos historico de linea de tiempo 1929-2022

knitr::opts_chunk$set(
  echo = FALSE, message = F,
  warning = F, ft.keepnext = F, tab.topcaption = TRUE,
  ft.align = "left", tab.cap.pre = "Tabla ", tab.cap.sep = ": "
)
library(tidyverse)
library(emmeans)
library(agricolae)
library(janitor)
library(data.table)
library(cluster)
library(forcats)
library(gower)
library(officer)
library(ggpubr)
library(flextable)
library(knitr)
library(dumbbell)
library(ggplot2)
library(sf)
#library(rnaturalearth)
#library(rnaturalearthdata)
library(ComplexUpset)
### librarias para modelos
library(gmodels)
library(nnet)
library(emmeans)
library(kableExtra)
source("D:/omar/METRIKA-GROUP/Github/agrograf/count_utils.R")
source("D:/omar/METRIKA-GROUP/Github/agrograf/upset_tools.R")
source("D:/omar/METRIKA-GROUP/Github/agrograf/upset_plot.R")
source("D:/omar/METRIKA-GROUP/Github/agrograf/bar_plot.R")

source("R/utils.R")
source("R/exsity_analysis_genebank.R")
source("R/erosion.R")
set.seed(123)
autonum <- run_autonum(seq_id = "tab", bkm = "TC1", bkm_all = TRUE)
options(knitr.duplicate.label = "allow")

datos_time <- readxl::read_xlsx(path = "D:/omar/METRIKA-GROUP/Github/analisis_datos_cip/TimeLine/Data_Curada/curado_inventario_general/curado_inventario_1927-202_v04.xlsx") %>% as.data.frame()

#datos_time <- datos_time 
## 

1 Descripción del estudio

El Estudio de Linea de Tiempo Historico contiene un total de 12000. Los datos del estudio estan comprendidos entre los años 1927 y 2022. El estudio fue llevado a cabo en Paucartambo, en la región Cusco, en los distritos de Paucartambo, Colquepata, Huancarani, Challabamba, Caicay y Kosñipata.

Debajo una tabla descriptiva con las observaciones estratificadas por año:

res <- datos_time %>% tabyl(ADM3_Name, cu_date_collection)
adm3 <- c(res$ADM3_Name,"Total")
res <- res %>% adorn_totals("row") %>% select(-ADM3_Name) %>% mutate(Total = rowSums(.))
names(res) <- paste0("y", names(res))
res <- res %>% add_column(ADM3_Name=adm3) %>% relocate(ADM3_Name, .before = "y1945")
names(res) <- str_remove_all(names(res),"y")
flextable(res)

ADM3_Name	1945	1947	1951	1953	1963	1970	1971	1972	1975	1982	1986	1988	1990	1991	1993	1994	1997	2001	2005	2008	2009	2013	2016	2019	2022	Total
Challabamba	0	87	0	0	0	0	0	3	12	0	42	0	0	0	0	0	0	0	0	1	2,556	0	0	491	2,935	6,127
Colquepata	0	128	0	0	0	2	0	0	29	0	272	73	0	50	4	110	0	460	0	3	0	26	8	531	2,143	3,839
Paucartambo	50	167	1	1	1	9	70	10	172	1	515	0	21	0	15	109	31	460	1	1	0	8	8	802	2,838	5,291
Total	50	382	1	1	1	11	70	13	213	1	829	73	21	50	19	219	31	920	1	5	2,556	34	16	1,824	7,916	15,257

Numero de observaciones estratificada por año y distrito

ADM3_Name	1941-1960	1961-1980	1981-2000	2001-2020	2021-actualidad	Total
Challabamba	87	15	42	3,048	2,935	6,127
Colquepata	128	31	509	1,028	2,143	3,839
Paucartambo	219	262	692	1,280	2,838	5,291
Total	434	308	1,243	5,356	7,916	15,257

Nota

Se observa entonces que los datos con mayor presencia de observaciones son Paucartambo, Challabamba y Colquepata.

Se observa que el numero de observaciones es completo a partir de 1941-1960, puesto que no existe datos de Challabamba para la decada de 1930-1940.

Teniendo en cuenta las anteriores consideraciones, filtramos los datos con esos criterios

2 Estadística descriptiva

2.1 Frecuencia de variedades unicas por periodo de años

f_period <- count_distinct(datos_time, "period_date_strat", "cu_variety_name",percentage = TRUE)
flextable(f_period)

period_date_strat	n_cu_variety_name	prcnt_cu_variety_name
1941-1960	226	8.9
1961-1980	212	8.4
1981-2000	353	13.9
2001-2020	699	27.6
2021-actualidad	1,047	41.3

2.2 Frecuencia de variedades unicas por cultivar group

f_cgroup <- count_distinct(datos_time, "cultivar_group", "cu_variety_name",percentage = TRUE)
flextable(f_cgroup)

cultivar_group	n_cu_variety_name	prcnt_cu_variety_name
Bitter landrace	138	9.0
Floury landrace	1,368	89.4
Modern variety	24	1.6

2.3 Frecuencia de variedades únicas por ADM3

f_adm3 <- count_distinct(datos_time, "ADM3_Name", "cu_variety_name",percentage = TRUE)
flextable(f_adm3 )

ADM3_Name	n_cu_variety_name	prcnt_cu_variety_name
Challabamba	645	26.6
Colquepata	812	33.4
Paucartambo	972	40.0

2.4 Frecuencia de variedades únicas por cultivar group y periodo de años

f_adm3_period <- count_distinct(datos_time, c("ADM3_Name", "period_date_strat") ,"cu_variety_name") %>% group_by(ADM3_Name) %>% mutate(prcnt_cu_variety_name= pct_fun(n_cu_variety_name))
flextable(f_adm3_period )

ADM3_Name	period_date_strat	n_cu_variety_name	prcnt_cu_variety_name
Challabamba	1941-1960	65	7.1
Challabamba	1961-1980	15	1.6
Challabamba	1981-2000	30	3.3
Challabamba	2001-2020	302	32.9
Challabamba	2021-actualidad	506	55.1
Colquepata	1941-1960	93	7.5
Colquepata	1961-1980	27	2.2
Colquepata	1981-2000	200	16.0
Colquepata	2001-2020	377	30.2
Colquepata	2021-actualidad	551	44.2
Paucartambo	1941-1960	151	9.2
Paucartambo	1961-1980	189	11.5
Paucartambo	1981-2000	260	15.8
Paucartambo	2001-2020	530	32.2
Paucartambo	2021-actualidad	517	31.4

2.5 Frecuencia de variedades unicas por adm3, cultivar group y periodo de años

f_adm3_vari_period <- count_distinct(datos_time, c("ADM3_Name", "period_date_strat","cultivar_group") ,"cu_variety_name") %>% group_by(ADM3_Name,period_date_strat) %>% mutate(pct= pct_fun(n_cu_variety_name)) %>% ungroup()

2.6 Presencia de variedades (únicas) más comunes que estan presentes en diferentes periodos

upstbl <- freq_upset(datos_time, "cu_variety_name", "period_date_strat")
upspam <- pam_upset(upstbl)
graf_ups_period_variety <- upset_plot(upspam, c("1941-1960","1961-1980","1981-2000", "2001-2020", "2021-actualidad"), c("black","green","blue","red","brown"))
graf_ups_period_variety

Consideraciones del gráfico:

1. Hay un grupo de 67 variedades que son persistentes en las 5 periodos: desde 1941 hasta 2021.
2. Existen variedades que estuvieron presentes solamente en alguno de los 5 periodos: 2021-actualidad (527), 2001-2002 (232), 1981-2000 (121), 1961-1980(42), y 1941-1960 (129)
3. Existen variedades que están presentes en 2 periodos, 3 periodos y hasta 4 periodos de forma continua o intermitente.

Dado que hay variabilidad en la presencia-ausencia de variedades en distintos periodos, asi como saltos temporales, es importante estudiar la autocorrelación temporal y su significancia. Con ello identificaremos que tan probable (o no) es la presencia de una variedad dado que estan en algun periódo anterior.

2.7 Modelo Autoregresivo AR1 con respuesta Binomial para presencia-ausencia de variedades

Dado que tenemos datos de variedades que tienen presencia o ausencia en diferentes puntos en el tiempo (discreto), buscamos el grado de autocorrelación temporal en los datos observados.

Para ello hacemos uso del proceso autoregresivo AR(1) con respuesta Binomial (1=Presencia y 0=Ausencia). Utilizamos INLA para modelar este problema:

rowname	mean	sd	0.025quant	0.5quant	0.975quant	mode
Precision for time	0.8971198	0.6804999	0.1277867	0.7238922	2.664238	0.3706502
Rho for time	0.6946634	0.1528557	0.3328123	0.7191617	0.918379	0.7795648
Precision for cu_variety_name	2.9763476	0.5107982	2.1807138	2.9042951	4.174475	2.7286162

De donde obtenemos una autocorrelación temporal (rho) positiva y significativa de 0.69 con un intervalo de credibilidad entre 0.3328123 y 0.918379, el cual indica que en general hay una persistencia temporal en encontrar una especie en el presente dado que estaba en periodo anterior. Hay que tener en cuenta que el intervalo de credibilidad es amplio.

2.8 Probabilidad de estar presente-ausente una variedad a los largo del tiempo

Examinamos la probabilidad de cada variedad a los largo del tiempo

2.9 Análisis de la probabilidad y su significancia de las variedades que estan presentes en los 5 periodos (8 decadas)

Tabla de las probabilidades e intervalos de confianza del grupo de variedades

Veamos la persistencia

2.10 Erosión varietal de los últimos 80 años vs actualidad

2.10.1 Tabla de erosión varietal en los últimos 80 años en Paucartambo

id	variedad	# variedades hace 8 decadas	# variedades 2021-actualidad	GE
1	Maqtillo	251	10	96.02
2	Sunchu	16	1	93.75
3	Compis	59	8	86.44
4	Leqe Chaki	17	3	82.35
5	Waka Wasi	9	2	77.78
6	Allqa Warmi	35	8	77.14
7	Sole	35	8	77.14
8	Soqo	11	3	72.73
9	Waka Chilena	6	2	66.67
10	Churuspiña	32	12	62.50
11	Puka Cusi	8	3	62.50
12	Bole	17	8	52.94
13	Chuqllu	25	12	52.00
14	Qewillo	18	9	50.00
15	Yuraq Miskilla	13	7	46.15
16	Cuchillo Paqui	38	21	44.74
17	Yuraq Chakillu	12	7	41.67
18	Waka Qallun	37	22	40.54
19	Phaspa Sunchu	26	18	30.77
20	Allqa Miskilla	13	9	30.77
21	Amparaeses	13	9	30.77
22	Chimako	34	24	29.41
23	Talaco	18	13	27.78
24	Chilca	11	8	27.27
25	Puka Runtus	26	20	23.08
26	Lluctu Runtus	22	17	22.73
27	Puka Mama	53	41	22.64
28	Lima Wiraqocha	14	11	21.43
29	Puka Bole	124	104	16.13
30	Leqechu	13	12	7.69
31	Muru Qusi	58	54	6.90
32	Waqoto	19	18	5.26
33	Yuraq Ruki	22	21	4.55

2.10.2 Significancia de la erosión varietal

Veamos si hay diferencias significativas usando las proporciones

Se obtuvo la tabla con la significancia de la erosión genética usando las proporciones. Además del

Gráfico de los intervalos de confianza para cada variedad al \(95\%\)

2.11 Erosión varietal de los últimos 60 años vs actualidad

2.11.1 Tabla de erosión varietal en las últimos 60 años en Paucartambo

id	variedad	# variedades hace 6 decadas	# variedades 2021-actualidad	GE
1	Paciencia	12	2	83.33
2	Viruntus	16	4	75.00
3	Puka Viruntus	40	12	70.00
4	Cusi	12	4	66.67
5	Ruki	3	1	66.67
6	Puka Suytu	16	8	50.00
7	Yuraq Rumpus	4	2	50.00
8	Puka Misquilla	5	3	40.00
9	Churillo	24	15	37.50
10	Maralina	3	2	33.33
11	Titeriti	10	8	20.00
12	Yana Ruki	15	13	13.33
13	Azul Ruki	14	13	7.14

2.11.2 Significancia de la erosión varietal

Veamos si hay diferencias significativas usando las proporciones

Se obtuvo la tabla con la significancia de la erosión genética usando las proporciones

Gráfico de los intervalos de confianza para cada variedad al \(95\%\)

2.12 Erosion varietal de variedades de los último 40 años vs actualidad

2.12.1 Tabla de erosión varietal en las últimos 40 años en Paucartambo

id	variedad	# variedades hace 4 decadas	# variedades 2021-actualidad	GE
1	Huayro	200	20	90
2	Yuraq Mestiza	8	1	88
3	Rumpu	10	2	80
4	Yungay	146	38	74
5	Papa Carlos	84	22	74
6	Yana Markilla	10	3	70
7	Chinchero	16	5	69
8	Markilla	16	5	69
9	Wata Kachu	9	3	67
10	Olones	37	15	59
11	Taruka Runtu	7	3	57
12	Tile	14	8	43
13	Muru Viruntus	7	4	43
14	Pusi Bole	44	30	32
15	Dueñas	10	7	30
16	Salamanca	10	7	30
17	Puka Sole	10	8	20
18	Yuraq Waña	43	39	9
19	Puka Chilca	21	20	5

2.12.2 Significancia de la erosión varietal

Veamos si hay diferencias significativas usando las proporciones

Se obtuvo la tabla con la significancia al \(95\%\) de la erosión genética usando las proporciones

Gráfico de los intervalos de confianza para cada variedad al \(95\%\)

2.13 Erosion varietal de variedades de los último 20 años vs actualidad

2.13.1 Tabla de erosión varietal en las últimos 20 años en Paucartambo

id	variedad	# variedades hace 2 decadas	# variedades 2021-actualidad	GE
1	Azul Qeta	7	1	86
2	Canchan	76	11	86
3	Jamachi	6	1	83

2.13.2 Significancia de la erosión varietal

Veamos si hay diferencias significativas usando las proporciones

Se obtuvo la tabla con la significancia de la erosión genética usando las proporciones

Gráfico de los intervalos de confianza para cada variedad al \(95\%\)

2.14 Listado de variedades (conteo único) más comunes en todos los periodos de años

tbl_vunica_period <- upspam %>% filter(`2021-actualidad`==1, `2001-2020`==1, `1981-2000`==1,`1961-1980`==1,`1941-1960`==1) %>% select(cu_variety_name)
vct_vunica_period <- tbl_vunica_period$cu_variety_name %>% as.character()

res <- datos_time %>%  filter(cu_variety_name %in%  vct_vunica_period)

out <- janitor::tabyl(res, period_date_strat ,cu_variety_name) %>% 
        janitor::adorn_totals("row")
top_comunes_period <- out[nrow(out),] %>% t() %>% as.data.frame() %>% rownames_to_column() %>%  slice(-1)
colnames(top_comunes_period) <- c("category","n")
top_comunes_period <- top_comunes_period %>% mutate(n=as.numeric(n)) %>%   mutate(pct=pct_fun(n)) %>% arrange(desc(n)) 

DT::datatable(top_comunes_period)

hbar_plot(top_comunes_period %>% slice(1:10), "category","pct",title = "Top 10 variedades más comunes en 1930-2022")

2.15 Listado de variedades unicas en el periodo 1941-1960

Esta lista contiene un conjunto de variedades que solamente están presentes en el periodo 1941-1960.

tbl_vunica_period_4160 <- upspam %>% filter(`2021-actualidad`==0, `2001-2020`==0, `1981-2000`==0,`1961-1980`==0,`1941-1960`==1) %>% select(cu_variety_name)
vunica_period_4160 <- tbl_vunica_period_4160$cu_variety_name %>% as.character()

res_4160 <- datos_time %>%  filter(cu_variety_name %in%  vunica_period_4160)

El total de variedades unicas durantel periodo 29.

out4160 <- janitor::tabyl(res_4160, period_date_strat ,cu_variety_name) %>% 
        janitor::adorn_totals("row")
top_comunes_period4160 <- out4160[nrow(out4160),] %>% t() %>% as.data.frame() %>% rownames_to_column() %>%  slice(-1)
colnames(top_comunes_period4160) <- c("category","n")
top_comunes_period_4160 <- top_comunes_period4160 %>% mutate(n=as.numeric(n)) %>%   mutate(pct=pct_fun(n)) %>% arrange(desc(n))

Tabla de las variedades y frecuencias de las variedades más comunes de 1941-1960

DT::datatable(top_comunes_period_4160)

2.16 Listado de variedades unicas en el periodo 1961-1980

Esta lista contiene un conjunto de variedades que solamente están presentes en el periodo 1961-1980.

tbl_vunica_period_6180 <- upspam %>% filter(`2021-actualidad`==0, `2001-2020`==0, `1981-2000`==0,`1961-1980`==1,`1941-1960`==0) %>% select(cu_variety_name)
vunica_period_6180 <- tbl_vunica_period_6180$cu_variety_name %>% as.character()

res_6180 <- datos_time %>%  filter(cu_variety_name %in%  vunica_period_6180)

El total de variedades unicas durantel periodo 42.

out6180 <- janitor::tabyl(res_6180, period_date_strat ,cu_variety_name) %>% 
        janitor::adorn_totals("row")
top_comunes_period6180 <- out4160[nrow(out6180),] %>% t() %>% as.data.frame() %>% rownames_to_column() %>%  slice(-1)
colnames(top_comunes_period6180) <- c("category","n")
top_comunes_period_6180 <- top_comunes_period6180 %>% mutate(n=as.numeric(n)) %>%   mutate(pct=pct_fun(n)) %>% arrange(desc(n))

DT::datatable(top_comunes_period_6180)

2.17 Listado de variedades unicas en el periodo 1981-200

Esta lista contiene un conjunto de variedades que solamente están presentes en el periodo 1981-2000.

tbl_vunica_period_8120 <- upspam %>% filter(`2021-actualidad`==0, `2001-2020`==0, `1981-2000`==1,`1961-1980`==0,`1941-1960`==0) %>% select(cu_variety_name)
vunica_period_8120 <- tbl_vunica_period_8120$cu_variety_name %>% as.character()

res_8120 <- datos_time %>%  filter(cu_variety_name %in%  vunica_period_8120)

El total de variedades unicas durantel periodo 121.

out8120 <- janitor::tabyl(res_8120, period_date_strat ,cu_variety_name) %>% 
        janitor::adorn_totals("row")
top_comunes_period8120 <- out8120[nrow(out8120),] %>% t() %>% as.data.frame() %>% rownames_to_column() %>%  slice(-1)
colnames(top_comunes_period8120) <- c("category","n")
top_comunes_period_8120 <- top_comunes_period8120 %>% mutate(n=as.numeric(n)) %>%   mutate(pct=pct_fun(n)) %>% arrange(desc(n))

DT::datatable(top_comunes_period_6180)

2.18 Listado de variedades unicas en el periodo 2000-2020

Esta lista contiene un conjunto de variedades que solamente están presentes en el periodo 2000-2020.

tbl_vunica_period_200020 <- upspam %>% filter(`2021-actualidad`==0, `2001-2020`==1, `1981-2000`==0,`1961-1980`==0,`1941-1960`==0) %>% select(cu_variety_name)
vunica_period_200020 <- tbl_vunica_period_200020$cu_variety_name %>% as.character()

res_200020 <- datos_time %>%  filter(cu_variety_name %in%  vunica_period_200020)

El total de variedades unicas durantel periodo 232.

out200020 <- janitor::tabyl(res_200020, period_date_strat ,cu_variety_name) %>% 
        janitor::adorn_totals("row")
top_comunes_period200020 <- out200020[nrow(out200020),] %>% t() %>% as.data.frame() %>% rownames_to_column() %>%  slice(-1)
colnames(top_comunes_period200020) <- c("category","n")
top_comunes_period_200020 <- top_comunes_period200020 %>% mutate(n=as.numeric(n)) %>%   mutate(pct=pct_fun(n)) %>% arrange(desc(n))

DT::datatable(top_comunes_period_200020)

2.19 Listado de variedades unicas en el periodo 2020-actualidad

Esta lista contiene un conjunto de variedades que solamente están presentes en el periodo 2020-actualidad.

tbl_vunica_period_21actual <- upspam %>% filter(`2021-actualidad`==1, `2001-2020`==0, `1981-2000`==0,`1961-1980`==0,`1941-1960`==0) %>% select(cu_variety_name)
vunica_period_21actual <- tbl_vunica_period_21actual$cu_variety_name %>% as.character()

res_21actual <- datos_time %>%  filter(cu_variety_name %in%  vunica_period_21actual)

El total de variedades unicas durantel periodo 527.

out21actual <- janitor::tabyl(res_21actual, period_date_strat ,cu_variety_name) %>% 
        janitor::adorn_totals("row")
top_comunes_period21actual <- out21actual[nrow(out21actual),] %>% t() %>% as.data.frame() %>% rownames_to_column() %>%  slice(-1)
colnames(top_comunes_period21actual) <- c("category","n")
top_comunes_period_21actual <- top_comunes_period21actual %>% mutate(n=as.numeric(n)) %>%   mutate(pct=pct_fun(n)) %>% arrange(desc(n))

DT::datatable(top_comunes_period_21actual)

2.20 Presencia de variedades (únicas) más comunes que estan presentes en diferentes distritos (Challambamba, Colquepata, Paucartambo)

upstbl <- freq_upset(datos_time, "cu_variety_name", "ADM3_Name")
upspam <- pam_upset(upstbl)
graf_pres_adm3 <- upset_plot(upspam, c("Challabamba", "Colquepata", "Paucartambo"), c("green","blue","red"))
graf_pres_adm3

hbar_freq_plot(top_comunes_adm3 %>% slice(1:10),fill_color = "lightgrey", "category","n",upper_limit = 500 ,title = "Top 10 variedades más comunes en Paucartambo-Challabamba-Colquepata")

2.21 Presencia de variedades más comunes (conteo unico) en los periodos y los distritos

datos_time <- datos_time %>% mutate(perdis=paste(ADM3_Name,period_date_strat,sep="_"))
upstbl <- freq_upset(datos_time, "cu_variety_name", "perdis")
upspam <- pam_upset(upstbl)
graf_ups_period_variety <- upset_plot(upspam, c("Challabamba_1941-1960",
"Challabamba_1961-1980","Challabamba_1981-2000","Challabamba_2001-2020",
"Challabamba_2021-actualidad","Colquepata_1941-1960",
"Colquepata_1961-1980","Colquepata_1981-2000","Colquepata_2001-2020",
"Colquepata_2021-actualidad","Paucartambo_1941-1960",
"Paucartambo_1961-1980","Paucartambo_1981-2000","Paucartambo_2001-2020",
"Paucartambo_2021-actualidad"), colors = rep("black",15))

top10_comunes_dis_period<- upspam %>% 
    adorn_totals(where="col") %>% 
    arrange(desc(Total)) %>% 
    relocate(Total, .before="Challabamba_1941-1960")  %>%  
    slice(1:10)

#top10_comunes_dis_period

flextable(top10_comunes_dis_period)

cu_variety_name	Total	Challabamba_1941-1960	Challabamba_1961-1980	Challabamba_1981-2000	Challabamba_2001-2020	Challabamba_2021-actualidad	Colquepata_1941-1960	Colquepata_1961-1980	Colquepata_1981-2000	Colquepata_2001-2020	Colquepata_2021-actualidad	Paucartambo_1941-1960	Paucartambo_1961-1980	Paucartambo_1981-2000	Paucartambo_2001-2020	Paucartambo_2021-actualidad
Huaman Uma	14	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1
Puka Bole	14	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1
Yana Cusi	14	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1
Chiqchi Puru	13	1	1	0	1	1	1	0	1	1	1	1	1	1	1	1
Muru Qusi	13	1	1	0	1	1	1	0	1	1	1	1	1	1	1	1
Phaspa Sunchu	13	1	1	0	1	1	1	1	1	1	1	0	1	1	1	1
Puka Mama	13	1	1	1	1	1	0	0	1	1	1	1	1	1	1	1
Qequrani	13	1	0	1	1	1	1	0	1	1	1	1	1	1	1	1
Yana Bole	13	1	0	0	1	1	1	1	1	1	1	1	1	1	1	1
Muru Bole	12	1	0	0	1	1	1	0	1	1	1	1	1	1	1	1

hbar_freq_plot(top10_comunes_dis_period , "cu_variety_name","Total",title = "Top 10 variedades más comunes en 1930-2022")

2.22 Presencia de variedades que estan en conservación insitu

dt_exsitu_period <- lapply(X = datos_time$period_date_strat %>% 
                               unique(),function(x) get_freq_genebank(datos_time, 
                                                                "period_date_strat",
                                                                    x, 
                                                                "cu_variety_name", 
                                                                "data_source",
                                                                c("GRIN database",
                                                                 "CIP's GeneBank"),                                                              "cu_variety_name")
                           ) %>% 
                        rbindlist() %>% as.data.frame()

dt_exsitu_cultivar_period <- lapply(1:nrow(dt_exsitu_period), function(x)  
                                get_exsitu_cultivar(datos_time, dt_exsitu_period[x,], 
                                                    "varie", 
                                                dt_exsitu_period[x,]$period_date_strat)
                                ) %>% rbindlist() %>% as.data.frame()

2.23 Indice de Shanon para analisis de diversidad de los 3 distritos y con las variedades más comunes

library(plyr)
tbl_shanon <-  ddply(tbl_freq_varie_adm3,~ADM3_Name,function(x) {
                        data.frame(SHANNON=diversity(x[-1], index="shannon"))
 })

ggpubr::ggdotchart(tbl_shanon, x = "ADM3_Name", y = "SHANNON", color="ADM3_Name", dot.size = 5, y.text.col=TRUE, rotate=TRUE,  ggtheme = theme_pubr(), palette = c("#00AFBB", "#E7B800", "#FC4E07"), # Custom color palette
                   sorting = "descending" ) + theme_cleveland()  + xlab("Indice de Shannon")

2.24 Indice de Pielou para analisis de equidad de los 3 distritos y con las variedades más comunes

library(plyr)
tbl_pileu_adm3 <-  ddply(tbl_freq_varie_adm3,~ADM3_Name,function(x) {
         data.frame(Pielou=diversity(x[-1], index="simpson")/log(sum(x[-1]>0))) })
flextable(tbl_pileu_adm3)

ADM3_Name	Pielou
Challabamba	0.1727401
Colquepata	0.1739207
Paucartambo	0.1737129

ggpubr::ggdotchart(tbl_pileu_adm3, x = "ADM3_Name", y = "Pielou", color="ADM3_Name", dot.size = 5, y.text.col=TRUE, rotate=TRUE,  ggtheme = theme_pubr(), palette = c("#00AFBB", "#E7B800", "#FC4E07"), # Custom color palette
                   sorting = "descending" ) + theme_cleveland()  + xlab("Indice de Shannon")

2.25 Erosión genética

2.26 Tendencia temporal por año y categoria varietal

smry_cgroup_year <- janitor::tabyl(datos_time, cu_date_collection, cultivar_group) %>% 
    pivot_longer(2:4, names_to="cultivar_group") %>% 
    group_by(cu_date_collection) %>%  
    mutate(pct=pct_fun(value,2))
#smry_cgroup_year$cu_date_collection <-  lubridate::ym(sprintf("%d-01",smry_cgroup_year$cu_date_collection))

Gráfico de tendencia en el tiempo segun categoria varietal

library(ggthemes)
library(ggHoriPlot)
gho1 <-  smry_cgroup_year %>%  
  ggplot() +
  ggHoriPlot::geom_horizon(aes(cu_date_collection, 
                   pct,fill = ..Cutpoints..), 
                   origin = 0, horizonscale = c(0,20,40,60,80,100)) +
          scale_fill_hcl(palette = 'RdBu', reverse = T) +
        facet_grid(cultivar_group~.) +
  theme_few() +
  theme(
    panel.spacing.y=unit(0, "lines"),
    strip.text.y = element_text(angle = 0, hjust = 0),
    legend.position = 'top',
    #axis.text.y = element_blank(),
    #axis.title.y = element_blank(),
    #axis.ticks.y = element_blank(),
    panel.border = element_blank()
    )  +
    scale_x_continuous(breaks = smry_cgroup_year$cu_date_collection, labels = smry_cgroup_year$cu_date_collection) + ggpubr::rotate_x_text()
gho1    

2.26.1 Tendencia temporal por estratificacion #2 y categoria varietal

smry_cgroup_period <- janitor::tabyl(datos_time, period_date_strat, cultivar_group) %>%
        pivot_longer(2:4, names_to="cultivar_group") %>% 
        group_by(period_date_strat) %>%  
        mutate(pct=pct_fun(value,2))

library(ggthemes)
library(ggHoriPlot)
gho2 <- smry_cgroup_period %>%  
  ggplot() +
  ggHoriPlot::geom_horizon(aes(period_date_strat, 
                   pct,fill = ..Cutpoints..), 
                   origin = 0, horizonscale = c(0,20,40,60,80,100)) +
  scale_fill_hcl(palette = 'RdBu', reverse = T) +
  facet_grid(cultivar_group~.) +
  theme_few() +
  theme(
    panel.spacing.y=unit(0, "lines"),
    strip.text.y = element_text(angle = 0, hjust = 0),
    legend.position = 'top',
    #axis.text.y = element_blank(),
    #axis.title.y = element_blank(),
    #axis.ticks.y = element_blank(),
    panel.border = element_blank()
    )  +
    ggpubr::rotate_x_text()
gho2

2.27 Relacion entre zonas de altitud vs periodo de años y cultivar_group

a2 <- ggplot(datos_time, aes(x = period_date_strat, y = elevation, color=cultivar_group)) + geom_boxplot() + theme_bw() + stat_summary(fun.y=mean, geom="point", shape=20, size=1, color="yellow", fill="yellow") + facet_wrap(~ADM3_Name)
a2 + ggpubr::rotate_x_text()

lm1 <- lm(elevation~ period_date_strat + cultivar_group, data= datos_time)
aov_lm1 <- aov(lm1)
summary(aov_lm1)

                     Df    Sum Sq  Mean Sq F value   Pr(>F)    
period_date_strat     4 8.028e+07 20068964  231.79  < 2e-16 ***
cultivar_group        2 2.224e+06  1111773   12.84 2.68e-06 ***
Residuals         15250 1.320e+09    86581                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

out <- HSD.test(aov_lm1,"period_date_strat", group=TRUE,console=TRUE)

Study: aov_lm1 ~ "period_date_strat"

HSD Test for elevation 

Mean Square Error:  86580.61 

period_date_strat,  means

                elevation      std    r     Min     Max
1941-1960        3598.765 288.1917  434 2721.68 3917.47
1961-1980        3312.494 671.9209  308 2233.06 4254.93
1981-2000        3649.146 405.7562 1243 2721.68 4425.53
2001-2020        3715.493 201.7777 5356 2721.68 4138.10
2021-actualidad  3766.012 302.7576 7916 2782.43 4518.66

Alpha: 0.05 ; DF Error: 15250 
Critical Value of Studentized Range: 3.858118 

Groups according to probability of means differences and alpha level( 0.05 )

Treatments with the same letter are not significantly different.

                elevation groups
2021-actualidad  3766.012      a
2001-2020        3715.493      b
1981-2000        3649.146      c
1941-1960        3598.765      d
1961-1980        3312.494      e

DT::datatable(out$means)

Study: aov_lm1 ~ "cultivar_group"

HSD Test for elevation 

Mean Square Error:  86580.61 

cultivar_group,  means

                elevation      std     r     Min     Max
Bitter landrace  3752.832 294.4079  1994 2233.06 4518.66
Floury landrace  3721.900 309.0494 12155 2233.06 4518.66
Modern variety   3706.769 246.4412  1108 2782.43 4473.75

Alpha: 0.05 ; DF Error: 15250 
Critical Value of Studentized Range: 3.314818 

Groups according to probability of means differences and alpha level( 0.05 )

Treatments with the same letter are not significantly different.

                elevation groups
Bitter landrace  3752.832      a
Floury landrace  3721.900      b
Modern variety   3706.769      b