405 lines
16 KiB
R
Executable File
405 lines
16 KiB
R
Executable File
suppressPackageStartupMessages({
|
|
library(tibble)
|
|
library(dplyr)
|
|
library(rlang)
|
|
library(purrr)
|
|
})
|
|
source(here::here('R/common/utils.R'))
|
|
|
|
#' Met en forme le tableau des produits des vaches
|
|
#' Ajoutes les colonnes nécessaires au calcul du rang ecow des vaches
|
|
#' @param produits_vaches dataframe. Tab des produits des vaches du cheptel
|
|
add_data_ecow <- function(liste_produits, czhbc){
|
|
|
|
#Ajout d'une colonne avec le nombre de produits IPG ------------------ TODO PEUT ETRE PLUS UTILE, à comparer avec nb_fin_gestation
|
|
# Croisement avec les données SPIE
|
|
parents_ipg <- get_parents_ipg()
|
|
liste_produits <- base::merge(liste_produits, get_parents_ipg(), by.x='anim', by.y='ANIM', all.x=T, all.y=F)
|
|
|
|
# Calcul REPRO
|
|
liste_produits <- liste_produits %>%
|
|
mutate(
|
|
repro = case_when(
|
|
anim %in% czhbc$ANIM ~ "O",
|
|
!is.na(NBPRODIPG) & NBPRODIPG > 0 ~ "O",
|
|
nbFinGestation > 0 ~ "O", # TODO pas sure que ce soit la bonne variable
|
|
TRUE ~ NA_character_
|
|
)
|
|
) %>%
|
|
# Calcul MORTALITÉ
|
|
mutate(
|
|
age_jours = as.numeric(difftime(dateSortDetenteur, dateNaiss, units = "days")),
|
|
mortnat = if_else(
|
|
!is.na(causeSortDetenteur) & causeSortDetenteur == "M" &
|
|
!is.na(age_jours) & age_jours < 3,
|
|
"O",
|
|
NA_character_
|
|
),
|
|
mortsev = if_else(
|
|
!is.na(causeSortDetenteur) & causeSortDetenteur == "M" &
|
|
!is.na(age_jours) & age_jours >= 3 & age_jours < 211,
|
|
"O",
|
|
NA_character_
|
|
),
|
|
) %>%
|
|
select(-age_jours) # colonne intermédiaire à retirer si inutile
|
|
}
|
|
|
|
#' Ajoute à un tableau de produits les valeurs corrigées par l'effet cheptel
|
|
#' @param produits dataframe. Tab des produits
|
|
apply_effet_chep <- function(produits, effets_chep, rapport_MF){
|
|
effets <- effets_chep %>%
|
|
select(sexe, typeMipg, diff_pn, diff_p120, diff_p210)
|
|
|
|
# Associe les bons effets à chaque produit
|
|
produits <- produits %>%
|
|
dplyr::left_join(effets, by = c("sexe", "typeMipg")) %>%
|
|
mutate(
|
|
nbpp_corr = if_else(sexe == "2", NBPRODIPG * rapport_MF, NBPRODIPG),
|
|
pn_corr = if_else(!is.na(poidsNaiss), poidsNaiss + diff_pn, NA_real_)
|
|
)
|
|
}
|
|
|
|
#' Normalisation des stats cheptel
|
|
#' 1 = max
|
|
norm_chep <- function(x, var, stats_chep) {
|
|
r <- stats_chep %>% dplyr::filter(var == !!var)
|
|
if (nrow(r) == 0) return(rep(NA_real_, length(x)))
|
|
mmin <- r$min[1]; mmax <- r$max[1]
|
|
1 - (abs(mmax - x) / abs(mmax - mmin))
|
|
}
|
|
|
|
#' Calcul des statistiques pour les colonnes d'un tableau
|
|
get_stats_tbl <- function(tab, nom_tab, cols, conditions = NULL, nom_cond = NA) {
|
|
|
|
# Si conditions présentes → filtrer
|
|
if (!is.null(conditions)) {
|
|
tab <- tab %>% filter(!! enquo(conditions))
|
|
}
|
|
|
|
# Pour chaque colonne → calculer les stats
|
|
map_df(cols, function(col) {
|
|
x <- tab[[col]]
|
|
|
|
tibble(
|
|
var = paste0(nom_tab, "$", col),
|
|
cond = nom_cond,
|
|
min = round(min(as.numeric(x), na.rm = TRUE), 1),
|
|
q1 = round(quantile(as.numeric(x), 0.25, na.rm = TRUE), 1),
|
|
med = round(median(as.numeric(x), na.rm = TRUE), 1),
|
|
moy = round(mean(as.numeric(x), na.rm = TRUE), 1),
|
|
q3 = round(quantile(as.numeric(x), 0.75, na.rm = TRUE), 1),
|
|
max = round(max(as.numeric(x), na.rm = TRUE), 1),
|
|
nbval = sum(!is.na(as.numeric(x))),
|
|
nas = sum(is.na(as.numeric(x)))
|
|
)
|
|
})
|
|
}
|
|
|
|
get_synth_prod_parent <- function(vaches, produits, params_ponderation){
|
|
# On récupère les coefficients de pondération pour les pointages adultes
|
|
ppa <- params_ponderation$pointage$adulte
|
|
|
|
# On récupère les coefficients de pondération pour les pointages au sevrage
|
|
pps <- params_ponderation$pointage$sevrage
|
|
|
|
# =========================
|
|
# 1) Synthèse veaux par mère
|
|
# =========================
|
|
veaux_summary <- produits %>%
|
|
group_by(numeroMipg) %>%
|
|
summarise(
|
|
n_veaux = n(),
|
|
|
|
# plage de campagnes velages (max - min + 1), robustes aux NA
|
|
campn_min = {cn <- campagneNaiss[!is.na(campagneNaiss)]; if (length(cn)) min(cn) else NA_integer_},
|
|
campn_max = {cn <- campagneNaiss[!is.na(campagneNaiss)]; if (length(cn)) max(cn) else NA_integer_},
|
|
nbcampvel = ifelse(!is.na(campn_min) & !is.na(campn_max),
|
|
campn_max - campn_min + 1, NA_integer_),
|
|
porteuse = any(embryon == 'O'),
|
|
|
|
# date du dernier anais (dernier événement)
|
|
last_danais = {d <- dateNaiss[!is.na(dateNaiss)]; if (length(d)) max(d) else as.Date(NA)},
|
|
|
|
# moyennes nécessaires pour ptgP et précocité
|
|
mean_devmus = mean(dmSevrage, na.rm = TRUE),
|
|
mean_devsqe = mean(dsSevrage, na.rm = TRUE),
|
|
mean_af = mean(afSevrage, na.rm = TRUE),
|
|
mean_diff_dev = mean(dsSevrage - dmSevrage, na.rm = TRUE),
|
|
|
|
# indicateurs produits
|
|
mort = round((sum(mortsev == "O" | mortnat == "O", na.rm = TRUE)) / n_veaux * 100, 1),
|
|
txrepros = round(sum(repro == "O", na.rm = TRUE) / n_veaux * 100, 1),
|
|
nbpp = sum(NBPRODIPG, na.rm = TRUE),
|
|
nbpp_corr = if ("nbpp_corr" %in% names(.)) sum(nbpp_corr, na.rm = TRUE) else NA_real_,
|
|
txmales = round(sum(sexe == "1", na.rm = TRUE) / n_veaux * 100, 1),
|
|
txvf = round(sum(conditionNaiss %in% c("1","2"), na.rm = TRUE) / n_veaux * 100, 1),
|
|
|
|
# moyennes par sexe
|
|
pn_m = round(mean(poidsNaiss[sexe == "1"], na.rm = TRUE), 1),
|
|
p120_m = round(mean(pat120[sexe == "1"], na.rm = TRUE), 1),
|
|
p210_m = round(mean(pat210[sexe == "1"], na.rm = TRUE), 1),
|
|
pn_f = round(mean(poidsNaiss[sexe == "2"], na.rm = TRUE), 1),
|
|
p120_f = round(mean(pat120[sexe == "2"], na.rm = TRUE), 1),
|
|
p210_f = round(mean(pat210[sexe == "2"], na.rm = TRUE), 1),
|
|
|
|
# versions corrigées
|
|
pn_corr = if ("pn_corr" %in% names(.)) round(mean(pn_corr, na.rm = TRUE), 1) else NA_real_,
|
|
p120_corr = round(mean(pat120Corrige, na.rm = TRUE), 1),
|
|
p210_corr = round(mean(pat210Corrige, na.rm = TRUE), 1),
|
|
.groups = "drop"
|
|
)
|
|
|
|
# =========================
|
|
# 2) Jointure & calculs vaches
|
|
# =========================
|
|
v_ref <- vaches %>%
|
|
dplyr::left_join(veaux_summary, by = c("anim" = "numeroMipg")) %>%
|
|
mutate(
|
|
|
|
# pointage adulte synthétique
|
|
ptgV = ifelse(!is.na(dmcAdulte), round(ppa$dmC * dmcAdulte + ppa$ds * dsAdulte + ppa$af * afAdulte, 1), NA_real_),
|
|
|
|
# précocité & alpha
|
|
precocite = ifelse(!is.na(n_veaux) & n_veaux > 3, round(mean_diff_dev, 2), NA_real_), # TODO à voir pour changer la formule de précocité avec LJ
|
|
alpha = ifelse(is.na(precocite), 1.62, 1.62 - 0.01 * precocite), # TODO A MAJ
|
|
|
|
# estimation du poids adulte (pad)
|
|
pad = dplyr::case_when( # TODO rajouter paramètre cohérence des données, créer max et min
|
|
!is.na(pat24M) ~ round((pat24M - 50 * exp(-720 * alpha * 10^(-3))) / (1 - exp(-720 * alpha * 10^(-3))), 1),
|
|
!is.na(pat18M) ~ round((pat18M - 50 * exp(-540 * alpha * 10^(-3))) / (1 - exp(-540 * alpha * 10^(-3))), 1),
|
|
!is.na(pat12M) ~ round((pat12M - 50 * exp(-360 * alpha * 10^(-3))) / (1 - exp(-360 * alpha * 10^(-3))), 1),
|
|
TRUE ~ NA_real_
|
|
),
|
|
|
|
# temps improductif (e2, e3, e4)
|
|
e2 = dplyr::case_when(
|
|
is.na(ivv1) ~ 0,
|
|
ivv1 < 390 ~ 0,
|
|
TRUE ~ ivv1 - 390
|
|
),
|
|
e3 = dplyr::case_when(
|
|
is.na(ivv2Brut) ~ 0, # TODO voir avec Lauréna si on prend ajust ou brut
|
|
ivv2Brut < 365 ~ 0,
|
|
TRUE ~ ivv2Brut - 365
|
|
),
|
|
days_since_last = as.numeric(difftime(if_else(is.na(dateSortDetenteur), Sys.Date(), dateSortDetenteur), last_danais, units = "days")),
|
|
e4 = dplyr::case_when(
|
|
is.na(days_since_last) ~ 0,
|
|
days_since_last < 365 ~ 0,
|
|
TRUE ~ days_since_last - 365
|
|
),
|
|
|
|
# temps productif (%)
|
|
age_days = time_length( interval( dateNaiss, if_else(is.na(dateSortDetenteur), Sys.Date(), dateSortDetenteur) ), "days" ),
|
|
age_years = round(age_days / 365, 1),
|
|
tempsprod = round( (age_days - (agevel1 * 30.4 + e2 + e3 * (nbcampvel - 2) + e4)) / age_days * 100, 1 ),
|
|
|
|
# prolificité
|
|
prol = round(n_veaux / nbcampvel * 100, 1),
|
|
|
|
# pointage produits
|
|
ptgP = ifelse( is.na(mean_devmus), NA_real_, round(pps$devmus * mean_devmus + pps$devsqe * mean_devsqe + pps$af * mean_af, 1))
|
|
) %>%
|
|
# Conversion des NaN en NA sur certaines moyennes
|
|
mutate(across(
|
|
c(ptgP, pn_m, pn_f, pn_corr, p120_m, p120_f, p120_corr, p210_m, p210_f, p210_corr),
|
|
~ ifelse(is.nan(.), NA_real_, .)
|
|
))
|
|
}
|
|
|
|
get_note_carriere <- function(v_ref, params_ponderation){
|
|
# calcul des stats, valeurs extremes et references pour la normalisation
|
|
stats_chep <- get_stats_tbl(
|
|
tab = v_ref,
|
|
nom_tab = "vaches",
|
|
cols = c("agevel1", "ivv1", "ivv2Brut", "prol", "mort", "txrepros", "nbpp_corr", "txvf", "txmales", "ptgP", "pn_corr", "p120_corr",
|
|
"p210_corr", "pad", "ptgV", "age_years", "tempsprod", "pn_m", "pn_f", "p120_m", "p120_f", "p210_m", "p210_f", "nbpp")
|
|
)
|
|
|
|
# =========================
|
|
# Normalisation
|
|
# =========================
|
|
v_norm <- v_ref %>%
|
|
mutate(
|
|
# age au 1er vêlage normalisé
|
|
agevel1_n = dplyr::case_when(
|
|
is.na(agevel1) ~ NA_real_,
|
|
agevel1 > 48 ~ 0,
|
|
TRUE ~ round(
|
|
-2 * (10^(-6)) * (agevel1 * 30.4)^2 +
|
|
0.0027 * (agevel1 * 30.4) +
|
|
8 * (10^(-15)),
|
|
3
|
|
)
|
|
),
|
|
# ivv1 normalisé
|
|
ivv1_n = dplyr::case_when(
|
|
is.na(ivv1) ~ NA_real_,
|
|
ivv1 > 460 ~ 0,
|
|
ivv1 < 390 ~ 1,
|
|
TRUE ~ round(1 - abs(390 - ivv1) / abs(390 - 460), 3)
|
|
),
|
|
# ivv2+ normalisé
|
|
ivv2p_n = dplyr::case_when(
|
|
is.na(ivv2Brut) | is.nan(ivv2Brut) ~ NA_real_,
|
|
ivv2Brut > 435 ~ 0,
|
|
ivv2Brut < 365 ~ 1,
|
|
TRUE ~ round(1 - abs(365 - ivv2Brut) / abs(365 - 435), 3)
|
|
),
|
|
# normalisations "cheptel 1 = max"
|
|
pad_n = round(norm_chep(pad, "vaches$pad", stats_chep), 3),
|
|
ptgv_n = round(norm_chep(ptgV, "vaches$ptgV", stats_chep), 3),
|
|
txvf_n = round(norm_chep(txvf, "vaches$txvf", stats_chep), 3),
|
|
txm_n = round(norm_chep(txmales, "vaches$txmales", stats_chep), 3),
|
|
txrepros_n = round(norm_chep(txrepros, "vaches$txrepros", stats_chep), 3),
|
|
nbpp_n = round(norm_chep(nbpp_corr, "vaches$nbpp_corr", stats_chep), 3),
|
|
ptgp_n = round(norm_chep(ptgP, "vaches$ptgP", stats_chep), 3),
|
|
p120_n = round(norm_chep(p120_corr, "vaches$p120_corr", stats_chep), 3),
|
|
p210_n = round(norm_chep(p210_corr, "vaches$p210_corr", stats_chep), 3),
|
|
|
|
# prolificité
|
|
prol_n = dplyr::case_when(
|
|
is.na(prol) ~ NA_real_,
|
|
prol >= 100 ~ 1,
|
|
prol < 50 ~ 0,
|
|
TRUE ~ round(1 - (abs(100 - prol) / abs(100 - 50)), 3)
|
|
),
|
|
|
|
# poids naissance corrigé
|
|
pn_n = dplyr::case_when(
|
|
is.na(pn_corr) ~ NA_real_,
|
|
40 < pn_corr & pn_corr < 50 ~ 1,
|
|
22 > pn_corr | pn_corr > 68 ~ 0,
|
|
22 < pn_corr & pn_corr < 40 ~ round(0.056 * (pn_corr - 22), 3),
|
|
TRUE ~ round(1 - 0.056 * (pn_corr - 50), 3)
|
|
),
|
|
|
|
# mortalité
|
|
mort_n = round(1.0 * exp(-0.031 * mort), 3)
|
|
)
|
|
|
|
# =========================
|
|
# Note carrière (pondérée)
|
|
# ========================
|
|
# Récupère les paramètres de pondérations, ATTENTION, il faut que leurs noms soient parfaitement identiques à ceux de v_norm
|
|
weights <- purrr::map_dbl(params_ponderation$carriere, 1)
|
|
|
|
v_final <- v_norm %>%
|
|
rowwise() %>%
|
|
mutate(
|
|
SOMME_tot = {
|
|
x <- c_across(all_of(names(weights)))
|
|
w <- weights
|
|
mask <- !is.na(x) & !is.na(w) &
|
|
is.finite(x) & is.finite(w) &
|
|
w != 0
|
|
|
|
if (sum(mask) == 0) {
|
|
NA_real_
|
|
} else {
|
|
weighted.mean(x[mask], w[mask]) * 10
|
|
}
|
|
},
|
|
ecowcarr = if_else(
|
|
is.na(ptgp_n) & is.na(p120_n) & is.na(p210_n), # règle d'exclusion VA4 : si pas ces trois valeurs ça enlève 1/4 de la note -> pas classable, voir pour créer une alternative pour éleveurs
|
|
NA_real_,
|
|
round(SOMME_tot * 100, 0)
|
|
)
|
|
) %>%
|
|
ungroup()
|
|
|
|
v_final <- v_final %>%
|
|
mutate(
|
|
rg_carr = as.integer(rank(1 / ecowcarr, na.last="keep"))
|
|
)
|
|
|
|
list(synthese = v_final, stats_chep = stats_chep)
|
|
}
|
|
|
|
get_stats_parent <- function(produits, regroupement){
|
|
produits %>%
|
|
group_by({{regroupement}}) %>%
|
|
summarise(
|
|
nb_prod_in_chep = n(),
|
|
utilgen = round(mean(rangVelageMipg == 1, na.rm = TRUE) * 100, 1),
|
|
prol = round(n() / n_distinct(dateNaiss, numeroMipg) * 100, 1),
|
|
mort = round((sum(mortsev == "O" | mortnat == "O", na.rm = TRUE)) / n() * 100, 1),
|
|
|
|
txrepros = round(
|
|
sum(repro == "O", na.rm = TRUE) /
|
|
sum(is.na(mortsev) & is.na(mortnat)) * 100, 1
|
|
),
|
|
|
|
nbpp = sum(NBPRODIPG, na.rm = TRUE),
|
|
txvf = round(mean(conditionNaiss %in% c("1", "2"), na.rm = TRUE) * 100, 1),
|
|
|
|
pnm = round(mean(poidsNaiss[sexe == "1"], na.rm = TRUE), 1),
|
|
pnf = round(mean(poidsNaiss[sexe == "2"], na.rm = TRUE), 1),
|
|
|
|
p120m = round(mean(pat120[sexe == "1"], na.rm = TRUE), 1),
|
|
p120f = round(mean(pat120[sexe == "2"], na.rm = TRUE), 1),
|
|
|
|
p210m = round(mean(pat210[sexe == "1"], na.rm = TRUE), 1),
|
|
p210f = round(mean(pat210[sexe == "2"], na.rm = TRUE), 1),
|
|
|
|
dmsev = round(mean(dmSevrage, na.rm = TRUE), 1),
|
|
dssev = round(mean(dsSevrage, na.rm = TRUE), 1),
|
|
afsev = round(mean(afSevrage, na.rm = TRUE), 1),
|
|
|
|
nb_femelles = sum(sexe == "2", na.rm = TRUE),
|
|
.groups = "drop"
|
|
) %>%
|
|
filter(nb_prod_in_chep >= 5)
|
|
}
|
|
|
|
#' Renvoie les statistiques des produits pour les éléments ayant plus de 3 produits
|
|
#' @param regroupement character. Champs sur lequel on veut faire le regroupement
|
|
#' @param actif booléen. Filtre sur les produits
|
|
#' selon le champs {}
|
|
#' @return Liste d'adhérents
|
|
get_stats_filles <- function(df, prefix, regroupement, actif = FALSE, cheptel) {
|
|
nb <- df %>%
|
|
group_by({{regroupement}}) %>%
|
|
summarise(
|
|
"{prefix}count" := n(),
|
|
.groups = "drop"
|
|
)
|
|
|
|
if (actif) {
|
|
df <- df %>%
|
|
filter(
|
|
is.na(dateSortDetenteur),
|
|
cheptelDetenteur == cheptel # Vérifie pour les descendantes car elles ont pu être vendues
|
|
)
|
|
}
|
|
|
|
res <- df %>% group_by({{regroupement}}) %>%
|
|
summarise(
|
|
"{prefix}nbavecprod" := sum(NBPRODIPG > 0, na.rm = TRUE),
|
|
"{prefix}pctavecprod" := round(sum(NBPRODIPG > 0, na.rm = TRUE) / n() * 100, 1),
|
|
"{prefix}isu" := ifelse(n() >= 3, sum(embryon == "O", na.rm = TRUE), NA),
|
|
"{prefix}age_sort" := ifelse(n() >= 3, round(mean(age_years, na.rm = TRUE), 1), NA),
|
|
"{prefix}agevel1" := ifelse(n() >= 3, round(mean(agevel1, na.rm = TRUE), 1), NA),
|
|
"{prefix}ivv1" := ifelse(n() >= 3, round(mean(ivv1, na.rm = TRUE), 1), NA),
|
|
"{prefix}ivv2p" := ifelse(n() >= 3, round(mean(ivv2Brut, na.rm = TRUE), 1), NA),
|
|
"{prefix}vieprod" := ifelse(n() >= 3, round(mean(tempsprod, na.rm = TRUE), 1), NA),
|
|
|
|
"{prefix}dmad" := ifelse(n() >= 3, round(mean(dmcAdulte, na.rm = TRUE), 1), NA),
|
|
"{prefix}dsad" := ifelse(n() >= 3, round(mean(dsAdulte, na.rm = TRUE), 1), NA),
|
|
"{prefix}afad" := ifelse(n() >= 3, round(mean(afAdulte, na.rm = TRUE), 1), NA),
|
|
|
|
"{prefix}prol" := ifelse(n() >= 3, round(mean(prol, na.rm = TRUE), 1), NA),
|
|
"{prefix}mort" := ifelse(n() >= 3, round(mean(mort, na.rm = TRUE), 1), NA),
|
|
"{prefix}txvf" := ifelse(n() >= 3, round(mean(txvf, na.rm = TRUE), 1), NA),
|
|
|
|
"{prefix}nbprod" := ifelse(n() >= 3, sum(NBPRODIPG, na.rm = TRUE), NA),
|
|
"{prefix}txrepros" := ifelse(n() >= 3, round(mean(txrepros, na.rm = TRUE), 1), NA),
|
|
"{prefix}nbpp" := ifelse(n() >= 3, sum(nbpp, na.rm = TRUE), NA),
|
|
|
|
.groups = "drop"
|
|
)
|
|
|
|
res <- left_join(res, nb, by = rlang::as_name(rlang::ensym(regroupement)))
|
|
}
|