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, noms_parents, czhbc){ #Ajout d'une colonne avec le nombre de produits IPG (croisement base SIG et 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) noms_parents <- as.data.frame(noms_parents) liste_produits$nommere <- noms_parents$V2[match(liste_produits$numeroMipg, noms_parents$V1)] liste_produits$nompere <- noms_parents$V2[match(liste_produits$pereGenetique, noms_parents$V1)] # Calcul REPRO liste_produits <- liste_produits %>% mutate( repro = case_when( anim %in% czhbc$ANIM ~ "O", !is.na(NBPRODIPG) & NBPRODIPG > 0 ~ "O", 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))) ) }) } #' Renvoie la synthèse pour chaque vache des performances de ses produits 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")) ) } get_note_campagne <- function(prod_vaches_corr, params_ponderation){ # Groupement et synthèse des données des produits par vache / date de naissance / rang de vélage de la mère # Ca permet de grouper les jumeaux synth_prod_vache <- prod_vaches_corr %>% group_by(numeroMipg, dateNaiss, rangVelageMipg) %>% summarise( ivv = first(ivv1), # ------------------------------------------------------------------------- TODO vraiment pas sure, à valider pn_c = round(mean(pn_corr, na.rm = TRUE), 1), txvf = round(mean(conditionNaiss %in% c('1','2')) * 100, 1), txm = round(mean(sexe == '1') * 100, 1), ptgp = round(mean(pps$devmus * dmSevrage + pps$devsqe * dsSevrage + pps$af * afSevrage, na.rm = TRUE), 1), p120_c = round(mean(pat120Corrige, na.rm = TRUE), 1), p210_c = round(mean(pat210Corrige, na.rm = TRUE), 1), prol = n() * 100, mort = round(mean(mortsev == "O" | mortnat == "O") * 100, 1), pere = first(pereGenetique), nom_pere = first(nompere), # infos des veaux pour simplifier l'affichage produits = list( pmap( list( nom = nom, anim = anim, sexe = sexe, mortnat = mortnat, mortsev = mortsev, # Nombre de produits tous cheptels confondus NBPRODIPG = NBPRODIPG, embryon = embryon ), list ) ), .groups = "drop" ) # Calcul des stats campagne pour chaque variable stats_camp <- get_stats_tbl( tab = synth_prod_vache, nom_tab = "synth_prod_vache", cols = c( 'ivv', 'prol', 'mort', 'txvf', 'txm', 'ptgp', 'pn_c', 'p120_c', 'p210_c') ) pond_camp_fin <- unlist(params_ponderation$campagne$final) pond_camp_ahp <- unlist(params_ponderation$campagne$AHPtech) # Normalisation des valeurs par rapport aux statistiques de campagnes synth_prod_vache_n <- synth_prod_vache %>% mutate( # pn normalisé pn_n = case_when( is.na(pn_c) ~ NA_real_, pn_c >= 40 & pn_c <= 50 ~ 1, pn_c <= 22 | pn_c >= 68 ~ 0, pn_c > 22 & pn_c < 40 ~ round(0.056 * (pn_c - 22), 3), TRUE ~ round(1 - 0.056 * (pn_c - 50), 3) ), # 5 normalisations linéaires txvf_n = round(1 - abs(stats_camp$max[stats_camp$var=="synth_prod_vache$txvf"] - txvf) / abs(diff(range(stats_camp[stats_camp$var=="synth_prod_vache$txvf",c("min","max")]))), 3), txm_n = round(1 - abs(stats_camp$max[stats_camp$var=="synth_prod_vache$txm"] - txm) / abs(diff(range(stats_camp[stats_camp$var=="synth_prod_vache$txm",c("min","max")]))), 3), ptgp_n = round(1 - abs(stats_camp$max[stats_camp$var=="synth_prod_vache$ptgp"] - ptgp) / abs(diff(range(stats_camp[stats_camp$var=="synth_prod_vache$ptgp",c("min","max")]))), 3), p120_n = round(1 - abs(stats_camp$max[stats_camp$var=="synth_prod_vache$p120_c"] - p120_c) / abs(diff(range(stats_camp[stats_camp$var=="synth_prod_vache$p120_c",c("min","max")]))), 3), p210_n = round(1 - abs(stats_camp$max[stats_camp$var=="synth_prod_vache$p210_c"] - p210_c) / abs(diff(range(stats_camp[stats_camp$var=="synth_prod_vache$p210_c",c("min","max")]))), 3), # prol prol_n = case_when( is.na(prol) ~ NA_real_, prol == 100 ~ 0.8, TRUE ~ 1 ), # mortalité mort_n = round(exp(-0.031 * mort), 3), # IVV ivv_n = case_when( is.na(rangVelageMipg) | rangVelageMipg == 1 | is.na(ivv) ~ NA_real_, # ravelamere == 2 rangVelageMipg == 2 & ivv > 460 ~ 0, rangVelageMipg == 2 & ivv < 390 ~ 1, rangVelageMipg == 2 ~ round(1 - abs(390 - ivv)/abs(390 - 460), 3), # autres ravelamere ivv > 435 ~ 0, ivv < 365 ~ 1, TRUE ~ round(1 - abs(365 - ivv)/abs(365 - 435), 3) ) ) %>% # Calcul du score final ecowcamp rowwise() %>% mutate( perf = list(c_across(c( ivv_n, mort_n, p120_n, p210_n, pn_n, prol_n, ptgp_n, txm_n, txvf_n ))), pond = sum( pond_camp_fin[!(is.na(perf) | is.nan(perf))] ), somme = sum( perf[!(is.na(perf) | is.nan(perf))] * pond_camp_ahp[!(is.na(perf) | is.nan(perf))] ), SOMME_tot = somme / pond * 10, ecowcamp = ifelse( is.na(ptgp_n) & is.na(p120_n) & is.na(p210_n), # Si pas de pointage, on réduit la note NA, round(SOMME_tot * 10, 0) ), produits = toJSON(produits, auto_unbox = TRUE) ) %>% ungroup() } 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 dans le cheptel #' @param regroupement character. Champs sur lequel on veut faire le regroupement #' @param actif booléen. Filtre sur les produits actifs #' @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 %>% filter( !is.na(n_veaux) # Garde seulement les filles ayant produit dans le cheptel ) %>% group_by({{regroupement}}) %>% summarise( "{prefix}nbavecprod" := sum( n_veaux > 0, na.rm = TRUE), "{prefix}pctavecprod" := round(sum(n_veaux > 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(n_veaux, na.rm = TRUE), NA), # On veut uniquement les produits dans le cheptel "{prefix}txrepros" := ifelse(n() >= 3, round(mean(txrepros, na.rm = TRUE), 1), NA), "{prefix}nbpp" := ifelse(n() >= 3, sum(nbpp, na.rm = TRUE), NA), # infos des produits pour l'affichage produits = list( pmap( list( nom = nom, anim = anim, dateNaiss = dateNaiss, rangVelageMipg = rangVelageMipg, numeroMipg = numeroMipg, nommere = nommere, pereGenetique = pereGenetique, nompere = nompere, active = is.na(dateSortDetenteur) & cheptelDetenteur == cheptel ), list ) ), "{prefix}produits" := toJSON(produits, auto_unbox = TRUE), .groups = "drop" ) res <- left_join(res, nb, by = rlang::as_name(rlang::ensym(regroupement))) }