library("tidyverse") library("ggpubr") ### +++ Load data sets +++ #### raw <- read.csv2("MyRepository/data_final.csv") # Remove farms that did not deliver milk raw2 <- filter(raw, no_delivery == 0) # Create subsets of data for increased overview farm <- select(raw2, id_respondant, country, organic_dummy:other_ms_description, own_rearing) suckle <- select(raw2, id_respondant, country, days_suckle, female_suckle_d:female_wean_wk) rear_type <- select(raw2, id_respondant, country, rearing:trt_udder_care) calfsystem <- select(raw2, id_respondant, country, rearing, pct_ge8:calf_no_concentrate) health <- select(raw2, id_respondant, country, cow_fertility:calf_death_rate) drivers <- select(raw2, id_respondant, country, start_year:final_comment) rm(raw) ### +++ General cleanup +++ ##### # Correct variable types farm$id_respondant <- as.factor(farm$id_respondant) farm$country <- as.factor(farm$country) farm$breed1 <- as.factor(farm$breed1) farm$breed2 <- as.factor(farm$breed2) farm$breed3 <- as.factor(farm$breed3) farm$breed4 <- as.factor(farm$breed4) farm$cross <- as.factor(farm$cross) farm$multiple_breeds <- as.factor(farm$multiple_breeds) farm$other_breeds <- as.factor(farm$other_breeds) farm$housing <- as.factor(farm$housing) farm$own_rearing <- as.factor(farm$own_rearing) farm$income_milk <- as.factor(farm$income_milk) suckle$id_respondant <- as.factor(suckle$id_respondant) suckle$country <- as.factor(suckle$country) suckle$female_contact_when <- as.factor(suckle$female_contact_when) rear_type$id_respondant <- as.factor(rear_type$id_respondant) rear_type$country <- as.factor(rear_type$country) rear_type$rearing <- as.factor(rear_type$rearing) calfsystem$id_respondant <- as.factor(calfsystem$id_respondant) calfsystem$country <- as.factor(calfsystem$country) calfsystem$rearing <- as.factor(as.character(calfsystem$rearing)) health$id_respondant <- as.factor(health$id_respondant) health$country <- as.factor(health$country) health$cow_fertility <- as.factor(health$cow_fertility) health$cow_udder_health <- as.factor(health$cow_udder_health) health$calf_health <- as.factor(health$calf_health) health$calf_gain <- as.factor(health$calf_gain) health$calf_diarrhoea <- as.factor(health$calf_diarrhoea) health$calf_cough <- as.factor(health$calf_cough) health$ab_12mo <- as.factor(health$ab_12mo) health$deworm_3y <- as.factor(health$deworm_3y) health$deworm_convent <- as.factor(health$deworm_convent) health$deworm_hpt <- as.factor(health$deworm_hpt) health$deworm_alt <- as.factor(health$deworm_alt) drivers$id_respondant <- as.factor(drivers$id_respondant) drivers$country <- as.factor(drivers$country) # Make income_milk numeric #### farm$num_income <- if_else(farm$income_milk == "0-25", 1, if_else(farm$income_milk == "26-50", 2, if_else(farm$income_milk == "51-75", 3, if_else(farm$income_milk == "76-100", 4, 5)))) # How many milking systems were used per farm? #### farm <- farm %>% mutate(ms_sum = select(., ams_ms:other_ms) %>% rowSums(na.rm = TRUE)) # 3 farms with > 1 milking system; one French farm that had fishbone parlour # both indoors and on pasture, one Italian farm that regularly milked the cows # with mobile milking bucket, but let stressed people hand milk a few of the # cows, and one large Swedish farm that stated both AMS and pipeline milking # without explanation. The Swedish farm had 200 cows and used foster cows that # were not milked for calf rearing - the pipeline cannot have been used for many # cows - maybe calving boxes? # Parlour - bucket - AMS will be used for these 3 farms, as all/most cows were # milked with these systems farm <- farm %>% mutate(ms_type = factor(if_else(ams_ms == 1, "ams", if_else(fishbone_ms == 1 | tandem_ms == 1 | parallell_ms == 1 | butterfly_ms == 1, "parlour", if_else(tie_pipe_ms == 1, "tie", if_else(bucket_ms == 1, "bucket", if_else(hand_ms == 1, "hand", "other"))))))) # Categories CCC allowance #### # Categorize CC-contact as 1) at milking, 2) half-day (between morning and # evening milking or vice verse), 3) permanent (regardless if the suckled cows # were milked or not), or other (combination of different routines, usually # depending on calf age) suckle <- suckle %>% mutate(when2 = factor(if_else(female_contact_when %in% c("after", "before", "both", "during"), "milking", if_else(female_contact_when %in% c("no_milking", "permanent"), "permanent", if_else(female_contact_when == "halfday", "halfday", "other"))))) ### +++ Analyses +++ #### # Q4. Breed #### # Select variables breed <- select(farm, id_respondant, country, breed1, breed2, breed3, breed4, cross, pure_bred, multiple_breeds, other_breeds) # Q6. Herd size #### # Select variables size <- select(farm, id_respondant, country, organic_dummy, cows_tot, ha_crop, ha_perm_pasture, ha_temp_pasture) # Dot-plot herd size vs average herd size per country (ave herd size based on # information from official sources and provided by representatives in each # country) p6 <- ggplot(size, aes(x = country, y = cows_tot)) + geom_jitter(width = 0.15, size = 3, alpha = 0.6) + # labs(x = "Country", y = "Nbr of adult cows") + geom_segment(aes(x = 0.75, y = 18, xend = 1.25, yend = 18), # AT colour = "red", size = 1.5) + geom_segment(aes(x = 1.75, y = 21, xend = 2.25, yend = 21), # CH colour = "red", size = 1.5) + geom_segment(aes(x = 2.75, y = 67, xend = 3.25, yend = 67), # DE colour = "red", size = 1.5) + geom_segment(aes(x = 3.75, y = 62, xend = 4.25, yend = 62), # FR colour = "red", size = 1.5) + geom_segment(aes(x = 4.75, y = 87, xend = 5.25, yend = 87), # IT colour = "red", size = 1.5) + geom_segment(aes(x = 5.75, y = 92, xend = 6.25, yend = 92), # SE colour = "red", size = 1.5) + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 12), legend.title = element_blank(), legend.key.size = unit(1.2, "line"), legend.text = element_text(size = 12), title = element_text(size = 16, face = "bold")) p6 rm(size) # Q8. Use of female calves #### # Note: quite a few confusing answers in this question, with less (or more) # female calves born than specified for different uses. For this reason the # data is summarized on farm level (farm does XX = yes/no) # Select variables female <- select(farm, id_respondant, country, cows_tot, calves_born:heifers_beef, own_rearing) ### Replacement stock ### # Dummy for having raised female calves for recruitment the last 12 months female <- female %>% mutate(own_recruitment = factor(if_else(heifers_kept > 0, "1", "0"))) # Rm farms that did not provide an answer for if they raised replacement stock # on their own farms, and farms that had no recruitment female2 <- filter(female, own_rearing != "no_answer", heifers_kept > 0) # Nbr of farms per country that did, or did not have replacement stock female2$own_rearing <- as.factor(as.character(female2$own_rearing)) own_rear <- female2 %>% group_by(country, own_rearing) %>% summarise(nbr_farms = n()) %>% ungroup() rm(female2, own_rear) ### Fattening ### # Dummy for having fattened female calves for slaughter the last 12 months female <- female %>% mutate(fatten = factor(if_else(heifers_beef > 0, "1", "0"))) # Compute number of farms per country to easily calculate proportion of farms # that do, or do not fatten female calves for slaughter female <- female %>% group_by(country) %>% mutate(tot_farms = n()) %>% ungroup() fatten_c <- female %>% group_by(country, fatten) %>% summarise(nbr_farms = n(), tot_farms = first(tot_farms)) %>% ungroup() # Proportion of farms per country that fatten female calves for slaughter fatten_c <- fatten_c %>% mutate(prp_farms = nbr_farms/tot_farms) rm(female) # Q10. Household income #### # Select variables income <- select(farm, id_respondant, country, cows_tot, income_milk, num_income) # Create median farm size per country, and create dummy farm size over or under # median for country income <- income %>% group_by(country) %>% mutate(median_herd_size_c = median(cows_tot), large_herd_c = if_else(cows_tot >= median_herd_size_c, 1, 0)) %>% ungroup() income <- income %>% mutate(median_herd_size_all = median(cows_tot), large_herd_all = if_else(cows_tot >= median_herd_size_all, 1, 0)) %>% ungroup() # Median herd size per country herd_size <- unique(select(income, country, median_herd_size_c)) # Create separate data set for those that did not answer, so they can be # evaluated for patterns (certain country, small/large farms etc.) no_income <- filter(income, num_income == 5) # Predominantly French farmers - questionnaire issue? # Rm farms that did not answer income_clean <- filter(income, !(num_income == 5)) # Barplot prp of household income from milk vs overall herd size pl_income <- income_clean %>% group_by(income_milk, large_herd_all) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(income_clean$income_milk), unique(income_clean$large_herd_all)) tmp <- tmp %>% rename(income_milk = Var1, large_herd_all = Var2) pl_income2 <- full_join(pl_income, tmp) pl_income2$nbr_farms[is.na(pl_income2$nbr_farms)] <- 0 p10 <- ggplot(data = pl_income2, aes(x = factor(income_milk), y = nbr_farms, fill = factor(large_herd_all))) + geom_col(colour = "black" , width = 0.7, position = position_dodge(0.7)) + scale_y_continuous(expand = c(0, 0)) + scale_x_discrete(breaks = c("0-25", "26-50", "51-75", "76-100"), labels = c("0-25%", "26-50%", "51-75%", "76-100%")) + scale_fill_manual(values = c("grey", "chartreuse4"), name = "Herd size", breaks = c("0", "1"), labels = c("Small", "Large")) + labs(x = "Proportion of household income from milk", y = "Nbr of farms") + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.4, 0.1, 0.4, "cm"), axis.text = element_text(size = 16), legend.key.size = unit(0.5, "cm"), title = element_text(size = 16, face = "bold")) p10 rm(income, income_clean, pl_income, pl_income2) # Q11. Housing system #### housing <- select(farm, id_respondant, country, cows_tot, housing, mixed_description) # Data set summarized by country pl_house <- housing %>% group_by(country, housing) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(housing$country), unique(housing$housing)) tmp <- tmp %>% rename(country = Var1, housing = Var2) pl_house2 <- full_join(pl_house, tmp) pl_house2$nbr_farms[is.na(pl_house2$nbr_farms)] <- 0 p11 <- ggplot(data = pl_house2, aes(x = country, y = nbr_farms, fill = factor(housing, levels = c("mixed", "loosehouse", "freestall", "tiestall", "outside")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("honeydew1", "darkcyan", "lightblue1", "azure3", "burlywood4"), name= "Housing", breaks = c("mixed", "loosehouse", "freestall", "tiestall", "outside"), labels = c("Mixed", "Open pack", "Free stall", "Tie stall", "Outside")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 16), legend.key.size = unit(0.5, "cm"), title = element_text(size = 16, face = "bold")) p11 rm(housing, pl_house, pl_house2) # Q12. Milking system #### ms <- select(farm, id_respondant, country, cows_tot, ms, ams_ms, rotary_ms, fishbone_ms, tandem_ms, parallell_ms, butterfly_ms, tie_pipe_ms, bucket_ms, hand_ms, other_ms, other_ms_description, ms_type) # Relevel factor order ms$ms_type <- factor(ms$ms_type, levels = c("other", "hand", "bucket", "ams", "tie", "parlour")) ### MS ~ country ### pl_ms <- ms %>% group_by(country, ms_type) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(ms$country), unique(ms$ms_type)) tmp <- tmp %>% rename(country = Var1, ms_type = Var2) pl_ms2 <- full_join(pl_ms, tmp) pl_ms2$nbr_farms[is.na(pl_ms2$nbr_farms)] <- 0 p12 <- ggplot(data = pl_ms2, aes(x = country, y = nbr_farms, fill = ms_type)) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("#999999", "mistyrose2", "darkorchid", "#56B4E9", "darkolivegreen4", "#E69F00"), name = "Milking system", breaks = c("other", "hand", "bucket", "ams", "tie", "parlour"), labels = c("Other", "Hand", "Bucket", "AMS", "Pipeline", "Parlour")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 16), legend.key.size = unit(0.5, "cm"), title = element_text(size = 16, face = "bold")) p12 rm(pl_ms, pl_ms2) # Q13. Rearing system #### rear_farm <- inner_join(rear_type, farm, by = c("id_respondant", "country")) # Select variables dam <- select(rear_farm, id_respondant, country, cows_tot, rearing, organic_dummy) # Relevel factor order dam$rearing <- factor(dam$rearing, levels = c("mother", "foster", "mix", "mother_bulk", "bulk_foster")) # Summarize pl_dam <- dam %>% group_by(country, rearing) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back rearing types for countries that lack one or more type for plotting tmp <- expand.grid(unique(dam$country), unique(dam$rearing)) tmp <- tmp %>% rename(country = Var1, rearing = Var2) pl_dam2 <- full_join(pl_dam, tmp) pl_dam2$nbr_farms[is.na(pl_dam2$nbr_farms)] <- 0 p13a <- ggplot(data = pl_dam2, aes(x = country, y = nbr_farms, fill = factor(rearing, levels = c("bulk_foster", "mother_bulk", "mix", "foster", "mother")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("grey40", "goldenrod1", "hotpink2", "darkseagreen4", "lightblue2"), name = "Calf rearing", breaks = c("bulk_foster", "mother_bulk", "mix", "foster", "mother"), labels = c("Manual + Foster", "Dam + Manual", "Mix", "Foster", "Dam")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 16), legend.key.size = unit(0.5, "cm"), title = element_text(size = 18, face = "bold")) p13a # Herd size per rearing type rear_u <- dam %>% group_by(rearing) %>% summarise(n = length(id_respondant), mean_cows = mean(cows_tot), sd_cows = sd(cows_tot), median_cows = median(cows_tot), Q1_cows = quantile(cows_tot, probs = 0.25), Q3_cows = quantile(cows_tot, probs = 0.75), min_cows = min(cows_tot), max_cows = max(cows_tot)) %>% ungroup() ### + Farm size ~ rearing + #### # Scatter/boxplot herd size vs rearing type p13b <- ggplot(data = dam, aes(x = rearing, y = cows_tot, fill = factor(rearing, levels = c("mother", "foster", "mix", "mother_bulk", "bulk_foster")))) + geom_boxplot(outlier.color = NA) + geom_jitter(shape = 16, position = position_jitter(0.25)) + scale_x_discrete(breaks = c("mother", "foster", "mix", "mother_bulk", "bulk_foster"), labels = c("Dam", "Foster", "Mix", "Dam + Manual", "Manual + Foster")) + scale_y_continuous(expand = c(0, 0),limits = c(0, 505), labels = c(0, 100, 200, 300, 400, 500)) + scale_fill_manual(values = c("darkgrey", "goldenrod1", "deeppink2", "darkseagreen4", "darkslategray2"), name = "Calf rearing", breaks = c("mother", "foster", "mix", "mother_bulk", "bulk_foster"), labels = c("Dam", "Foster", "Mix", "Dam + Manual", "Manual + Foster")) + labs(y = "Number of cows") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold"), legend.position = "none") p13b rm(dam, pl_dam, pl_dam2) ### + Days suckling mother_bulk + #### rear_suckle <- inner_join(rear_type, suckle, by = c("id_respondant", "country")) # Select variables suckle_bulk <- select(rear_suckle, id_respondant, country, rearing, female_suckle_d) # Filter for farms with initial dam rearing followed by manual milk feeding suckle_bulk2 <- filter(suckle_bulk, rearing == "mother_bulk") # Remove farms with conflicting answers # 1. Italian farm 157259327 stating calves kept with mother until wk 5, then # kept in group boxes. Allowed to suckle for 20 days, and weaned at 3 weeks # (possibly months, but no way to be sure). # 2. Italian farm 157259332 stating calves kept with cows on pasture for 150 d, # and kept with dam for at least 12 wk (from weekly description), but only # suckle for 60 days? Then are weaned at 25 wks? # 3. Italian farm 157259326 stating 9 days suckling in Q2, then 20 days suckling # for female calves in Q34 but explicitly states that separation takes place at # 9 days of age. Further, the calves have contact with the mother 10 times per # day, but are also housed permanently with the cow, and the female calves are # weaned at 4 weeks of age - the combination of answers seems highly unlikely. suckle_bulk3 <- filter(suckle_bulk2, id_respondant != "157259326", id_respondant != "157259327", id_respondant != "157259332") # Average and SD days suckling mean(suckle_bulk3$female_suckle_d) sd(suckle_bulk3$female_suckle_d) rm(rear_suckle, suckle_bulk, suckle_bulk2, suckle_bulk3) ### + Exploring rearing ~ profit + #### # Dam rearing farms are generally smaller (median 30 cows incl dry), compared to # all other calf rearing strategies (foster median 50 cows, mix median 50 cows, # dam + manual median 40 cows, manual + foster median 210 cows) # Could this be because dam rearing is performed due to ethical considerations, # while the other systems are efficient (either saving labor costs or reducing # the amount of milk lost)? # Here this idea is indirectly explored by looking at the relative frequency of # amount of household income from milk depending on rearing strategy # Select variables profit <- select(rear_farm, id_respondant, country, cows_tot, rearing, income_milk) # Rm farms that did not answer the question about income from milk profit <- filter(profit, income_milk != "unspecified") # Relevel factor order profit$income_milk <- factor(profit$income_milk, levels = c("76-100", "51-75", "26-50", "0-25")) profit$rearing <- factor(profit$rearing, levels = c("mother", "foster", "mix", "mother_bulk", "bulk_foster")) # Summarize pl_profit <- profit %>% group_by(rearing, income_milk) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(profit$rearing), unique(profit$income_milk)) tmp <- tmp %>% rename(rearing = Var1, income_milk = Var2) pl_profit2 <- full_join(pl_profit, tmp) pl_profit2$nbr_farms[is.na(pl_profit2$nbr_farms)] <- 0 p13c <- ggplot(data = pl_profit2, aes(x = rearing, y = nbr_farms, fill = income_milk)) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("black", "forestgreen", "gold1", "orange2", "orangered3"), name = "Share of income", breaks = c("unspecified", "76-100", "51-75", "26-50", "0-25"), labels = c("Unknown", "76-100%", "51-75%", "26-50%", "0-25%")) + labs(x = "Rearing", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p13c rm(profit, pl_profit, pl_profit2) # Q25. Prp with contact - calves #### # Note: Q30 - prp of female/male calves that are allowed to suckle - will be # analyzed in conjunction with Q25. # # Select variables prp_ccc <- select(calfsystem, id_respondant, country, pct_ge8, pct_female_suckle, pct_male_suckle) # Create variable that show difference in prp of calves allowed to suckle, # depending on gender prp_ccc$diff_suckle <- prp_ccc$pct_female_suckle - prp_ccc$pct_male_suckle # Sneak peak at all variables summary(prp_ccc) # Prp of calves with contact more than 7 days p25a <- ggplot(prp_ccc, aes(x = country, y = pct_ge8)) + geom_jitter(size = 4, alpha = 0.5, width = 0.2, height = 0) + labs(x = "Country", y = "Prp of calves w CCC > 7 days (%)") + scale_y_continuous(expand = c(0, 0), limits = c(0,102)) + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p25a # Prp of female calves allowed to suckle p25b <- ggplot(prp_ccc, aes(x = country, y = pct_female_suckle)) + geom_jitter(size = 4, alpha = 0.5, width = 0.2, height = 0) + labs(x = "Country", y = "Female calves allowed to suckle (%)") + scale_y_continuous(expand = c(0, 0), limits = c(0,102)) + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p25b # Prp of male calves allowed to suckle p25c <- ggplot(prp_ccc, aes(x = country, y = pct_male_suckle)) + geom_jitter(size = 4, alpha = 0.5, width = 0.2, height = 0) + labs(x = "Country", y = "Male calves allowed to suckle (%)") + scale_y_continuous(expand = c(0, 0), limits = c(0,102)) + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p25c # Difference in prp of female and male calves that are allowed to suckle per # farm p25d <- ggplot(prp_ccc, aes(x = country, y = diff_suckle)) + geom_hline(yintercept = 0, colour = "red", size = 2) + geom_jitter(size = 4, alpha = 0.5, width = 0.2, height = 0) + labs(x = "Country", y = "(pct_female_suckle - pct_male_suckle") + scale_y_continuous(limits = c(-50, 50)) + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p25d # Farms where < 100% of the calves had contact > 7 days l100_7d <- filter(prp_ccc, pct_ge8 < 100) # Farm where the prp of female and male calves allowed to suckle differed suckle_genderdiff <- filter(prp_ccc, diff_suckle != 0) rm(prp_ccc) # Q26. Calves per foster cow #### # Joining calfsystem with farm so I can compare number of days suckling # depending on certified organic or not rawfoster <- inner_join(calfsystem, farm, by = c("id_respondant", "country")) # Select variables foster <- select(rawfoster, id_respondant, country, organic_dummy, rearing, foster_nbr_calves, foster_start_day) # Keep only farms that use foster cows (n = 45) foster2 <- filter(foster, !is.na(foster_nbr_calves)) # Rm Swiss farm that state mix of mothers and fosters, but put 0 calves to # foster cows (this combination of answers was due to the calves being allowed # to run with the adult cows including the mothers for 160 d, and the farmer # had observed that the calves suckled from more than one cow) foster2 <- filter(foster2, id_respondant != "157259346") foster2$rearing <- as.factor(as.character(foster2$rearing)) # Boxplot to evaluate if there is differences between countries in the number # of calves put to foster cows boxplot(foster_nbr_calves ~ country, ylim = c(0,4), data = foster2) # Summaries per country foster.at <- filter(foster2, country == "AT") summary(foster.at) foster.ch <- filter(foster2, country == "CH") summary(foster.ch) foster.de <- filter(foster2, country == "DE") summary(foster.de) foster.fr <- filter(foster2, country == "FR") summary(foster.fr) foster.it <- filter(foster2, country == "IT") summary(foster.it) foster.se <- filter(foster2, country == "SE") summary(foster.se) # Boxplot to evaluate if there is differences between organic and conventional # farms in the number of calves put to foster cows boxplot(foster_nbr_calves ~ organic_dummy, ylim = c(0,4), data = foster2) # Summaries per farming system foster.org <- filter(foster2, organic_dummy == 1) summary(foster.org) sd(foster.org$foster_nbr_calves) foster.convent <- filter(foster2, organic_dummy == 0) summary(foster.convent) sd(foster.convent$foster_nbr_calves) rm(rawfoster, foster, foster.at, foster.ch, foster.de, foster.fr, foster.it, foster.se, foster.org, foster.convent) # Q27. Days in age to foster cow #### # Rm the 2 farms that did not specify at what age calves were moved to the # foster cows (n = 42) foster3 <- filter(foster2, foster_start_day != 999) # Boxplot to evaluate if there is differences between and withing countries in # when the calves are put to foster cows boxplot(foster_start_day ~ country, data = foster3) # Boxplot to evaluate if there is differences between organic and conventional # farms when the calves are put to foster cows boxplot(foster_start_day ~ organic_dummy, data = foster3) rm(foster2, foster3) # Q28. Colostrum feeding #### # Select variables col <- select(calfsystem, id_respondant, country, rearing, col_suck:col_description) # Which farms has more than one strategy? col <- col %>% mutate(col_sum = select(., col_suck:col_bottle) %>% rowSums(na.rm = TRUE)) # How many farms use what strategy col_u <- col %>% summarise(tot_suck = sum(col_suck), tot_drench = sum(col_drench), tot_bucket = sum(col_bucket), tot_bottle = sum(col_bottle)) # Q32. Forage access - calves #### # Select variables feed <- select(calfsystem, id_respondant, country, rearing, calf_forage_wk:calf_no_concentrate) # Rm farms that did not respond to Q32 roughage <- filter(feed, calf_forage_wk != "unspecified") tmp <- anti_join(feed, roughage) # 21 French farms and 1 Swiss farm did not respond # Make age at roughage access numeric roughage$calf_forage_wk <- as.numeric(roughage$calf_forage_wk) # Overall distribution hist(roughage$calf_forage_wk) # Which farms provided roughage from first wk of life? roughage_early <- filter(roughage, calf_forage_wk < 2) # Which farms were late with forage? roughage_late <- filter(roughage, calf_forage_wk > 1) # Average + SD age per country when first giving forage roughage_age <- roughage %>% group_by(country) %>% summarise(ave_age = mean(calf_forage_wk), sd_age = sd(calf_forage_wk), n = n()) %>% ungroup() rm(roughage, roughage_early, roughage_late) # Q33. Concentrate access - calves #### # Rm farms that did not respond to Q33 concentrate <- filter(feed, calf_concentrate_wk != "unspecified" | is.na(calf_concentrate_wk)) tmp <- anti_join(feed, concentrate) # 21 French farms (19 of which did not answer Q32) did not respond concentrate$calf_concentrate_wk <- as.numeric(concentrate$calf_concentrate_wk) # Nbr of farms that did not provide concentrate per country conc_u <- concentrate %>% group_by(country) %>% summarise(farm_w_answers = length(id_respondant), farms_no_conc = sum(calf_no_concentrate), farms_conc = farm_w_answers - farms_no_conc, prp_farms = farms_conc/farm_w_answers) %>% ungroup() # When is concentrate provided on the farms that do provide it? conc_wk <- filter(concentrate, !is.na(calf_concentrate_wk)) # Average + SD age per country when first giving concentrate conc_age <- conc_wk %>% group_by(country) %>% summarise(ave_age = mean(calf_concentrate_wk), sd_age = sd(calf_concentrate_wk), n = n()) %>% ungroup() rm(concentrate, conc_wk) # Q34. Days suckling #### # 2. Inclusion Q: How many days do calves suckle a cow (mother of foster cow)? # 34. How many days do the female suckler calves suckle? # Create variable female_different, to id how many farms has a different # suckling strategy for female calves specifically (Q2 vs Q34) suckle <- suckle %>% mutate(female_diff = if_else(days_suckle - female_suckle_d != 0, 1, 0)) # --> Weird answers on Q2 - will not be analyzed # Rm farms with conflicting answers: # 1. Italian farm 157259343 stating 90 days suckling in Q2, then 300 d suckling # for female calves in Q34, but 12 weeks for weaning age of female calves. # 2. Italian farm 157259326 stating 9 days suckling in Q2, then 20 days suckling # for female calves in Q34 but explicitly states that separation takes place at # 9 days of age. Further, the calves have contact with the mother 10 times per # day, but are also housed permanently with the cow, and the female calves are # weaned at 4 weeks of age - the combination of answers seems highly unlikely. # 3. Italian farm 157259327 stating calves kept with mother until wk 5, then # kept in group boxes. Allowed to suckle for 20 days, and weaned at 3 weeks # 4. Italian farm 157259332 stating calves kept with cows on pasture for 150 d, # and kept with dam for at least 12 wk (from weekly description), but only # suckle for 60 days? Then are weaned at 25 weeks? suckle2 <- filter(suckle, id_respondant != "157259343", id_respondant != "157259326", id_respondant != "157259327", id_respondant != "157259332") # Joining with farm so I can compare # days suckling depending on organic type suckle_farm <- inner_join(suckle2, farm, by = c("id_respondant", "country")) # Select variables suckle_farm2 <- select(suckle_farm, id_respondant, country, female_suckle_d, organic_dummy) ### X-axis jittered plot days/weeks suckling female calves p34 <- ggplot(suckle_farm2, aes(x = country, y = female_suckle_d/7, colour = as.factor(organic_dummy))) + geom_jitter(size = 5, width = 0.2, height = 0, alpha = 0.4) + stat_summary(fun = median, fun.min = median, fun.max = median, geom = "crossbar", width = 0.6, color = "black", size = 1.1) + scale_y_continuous(expand = c(0, 0), limits = c(0, 50), labels = c(0, 10, 20, 30, 40, 50)) + labs(x = "Country", y = "Suckling period (wk)") + scale_colour_manual(values = c("darkred", "grey55"), breaks = c(0, 1), labels = c("Conventional", "Organic")) + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.3, 0.0, 0.3, "cm"), axis.text = element_text(size = 16), legend.title = element_blank(), legend.key.size = unit(0.3, "cm"), title = element_text(size = 16, face = "bold"), legend.position = "bottom") p34 rm(suckle2, suckle_farm, suckle_farm2) # Q35. Supplemental milk - female calves #### # Select variables extra_milk <- select(suckle, id_respondant, country, female_suckle_d, female_extra) extra_milk2 <- filter(extra_milk, female_extra == 1) rm(extra_milk) # Q36. CCC allowance #### # Select variables contact <- select(suckle, id_respondant, country, when2) # Relevel factor order contact$when2 <- factor(contact$when2, levels = c("permanent", "halfday", "milking", "other")) ### CCC type per country ### # Summarize pl_contact <- contact %>% group_by(country, when2) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(contact$country), unique(contact$when2)) tmp <- tmp %>% rename(country = Var1, when2 = Var2) pl_contact2 <- full_join(pl_contact, tmp) pl_contact2$nbr_farms[is.na(pl_contact2$nbr_farms)] <- 0 p36a <- ggplot(data = pl_contact2, aes(x = country, y = nbr_farms, fill = factor(when2, levels = c("other", "milking", "halfday", "permanent")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("#999999", "gold2", "lightskyblue", "chartreuse4"), name= "Cow-calf \ncontact", breaks = c("other", "milking", "halfday", "permanent"), labels = c("Multiple", "Around milking", "Half day", "Full day")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 16), legend.key.size = unit(0.5, "cm"), title = element_text(size = 16, face = "bold")) p36a rm(pl_contact, pl_contact2) ### + CCC type ~ MS + #### ms_ccc <- full_join(ms, contact) # Remove ms category other, to make graph easier to interpret ms_ccc2 <- filter(ms_ccc, ms_type != "other") ms_ccc2$ms_type <- as.factor(as.character(ms_ccc2$ms_type)) # Relevel factor order ms_ccc2$ms_type <- factor(ms_ccc2$ms_type, levels = c("tie", "parlour", "ams", "bucket", "hand", "other")) # Summarize pl_ms.ccc <- ms_ccc2 %>% group_by(ms_type, when2) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(ms_ccc2$ms_type), unique(ms_ccc2$when2)) tmp <- tmp %>% rename(ms_type = Var1, when2 = Var2) pl_ms.ccc2 <- full_join(pl_ms.ccc, tmp) pl_ms.ccc2$nbr_farms[is.na(pl_ms.ccc2$nbr_farms)] <- 0 p36b <- ggplot(data = pl_ms.ccc2, aes(x = ms_type, y = nbr_farms, fill = factor(when2, levels = c("other", "milking", "halfday", "permanent")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_x_discrete(labels = c("Pipeline", "Parlour", "AMS", "Bucket", "Hand")) + scale_fill_manual(values = c("lightblue2", "gold2", "coral", "chartreuse4"), name= "Cow-calf \ncontact", breaks = c("other", "milking", "halfday", "permanent"), labels = c("Multiple", "At milking", "Half day", "Permanent")) + labs(x = "", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold"), legend.text = element_text(size = 18)) p36b rm(ms, ms_ccc2, pl_ms.ccc, pl_ms.ccc2) ### + CCC-type ~ rearing + #### rear_ccc <- full_join(ms_ccc, rear_type) # Select variables rear_ccc2 <- select(rear_ccc, id_respondant:when2, rearing) # Relevel rearing type rear_ccc2$rearing <- factor(rear_ccc2$rearing, levels = c("mother", "foster", "mix", "mother_bulk", "bulk_foster")) # Summarize pl_rear.ccc <- rear_ccc2 %>% group_by(rearing, when2) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(rear_ccc2$rearing), unique(rear_ccc2$when2)) tmp <- tmp %>% rename(rearing = Var1, when2 = Var2) pl_rear.ccc2 <- full_join(pl_rear.ccc, tmp) pl_rear.ccc2$nbr_farms[is.na(pl_rear.ccc2$nbr_farms)] <- 0 p36c <- ggplot(data = pl_rear.ccc2, aes(x = rearing, y = nbr_farms, # fill = when2)) + fill = factor(when2, levels = c("other", "milking", "halfday", "permanent")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_x_discrete(labels = c("Dam", "Foster", "Mix", "Dam + \nManual", "Manual + \nFoster")) + scale_fill_manual(values = c("lightblue2", "gold2", "coral", "chartreuse4"), name= "Cow-calf \ncontact", breaks = c("other", "milking", "halfday", "permanent"), labels = c("Multiple", "At milking", "Half-day", "Permanent")) + labs(x = "", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold"), legend.text = element_text(size = 18)) p36c rm(rear_ccc, rear_ccc2, pl_rear.ccc, pl_rear.ccc2) # Q37. Nbr of CCC per day #### # Keeping only farms where female calves have restricted access to cows nbr_contact <- filter(suckle, female_contact != "permanent") # Select variables nbr_contact2 <- select(nbr_contact, id_respondant, country, female_suckle_d, female_contact:female_when_description, when2) rm(nbr_contact) # Q41. Weaning #### # Create age in days at weaning to see if suckle period is different than # milk period, allowing for rounding errors since suckle period was registered # in days, while milk period was registered in wks (in some cases one mo was # likely considered 4 wks, i.e. 30 vs. 28 days) suckle$female_wean_d <- suckle$female_wean_wk * 7 suckle$wean_suckle_diff <- suckle$female_wean_d - suckle$female_suckle_d # Rm farms with inconsistent answers: # 1. Italian farm 157259343 stating 90 days suckling in Q2, then 300 d suckling # for female calves in Q34, but 12 weeks for weaning age of female calves. # 2. Swiss farm 157259357 stating 60 days suckling in Q2, then 60 d suckling for # female calves, but 5 weeks for weaning age of female calves - with no weekly # description passed wk 6 # 3. Italian farm 157259326 stating 9 days suckling in Q2, then 20 days suckling # for female calves in Q34 but explicitly states that separation takes place at # 9 days of age. Further, the calves have contact with the mother 10 times per # day, but are also housed permanently with the cow, and the female calves are # weaned at 4 weeks of age - the combination of answers seems highly unlikely. # 4. Italian farm 157259327 stating 30 days suckling in Q2, then 20 d suckling # for female calves in Q34, and 3 weeks for weaning age of female calves - too # early to be biologically possible, likely that weaning was confused with # separation # 5. Swiss farm 157259348 stating 20 days suckling in Q2, then 20 d suckling # for female calves followed by bucket feeding in Q34, but 3 weeks for weaning # age of female calves. As real weaning age is unknown, the farm is rm for this # analysis # 6. Austrian farm 157259390 stating 21 days suckling in Q2, then 21 d suckling # for female calves after which all calves were sold to a rearing farm. As real # weaning age is unknown, the farm is rm for this analysis # 7. Italian farm 157259332 stating 60 days suckling, but kept with dam for # at least 12 weeks based on weekly description, and 150 d on pasture with adult # cows, then weaned at 25 weeks. Unlikely combination of answers wean <- filter(suckle, id_respondant != "157259343", id_respondant != "157259357", id_respondant != "157259326", id_respondant != "157259327", id_respondant != "157259348", id_respondant != "157259390", id_respondant != "157259332") # Note: all farms with negative wean_suckle_diff is due to rounding errors when # using 4 weeks as month in Q41, and 30 days as month in Q34 # Note 2: all farms with long suckling period, but still much shorter than the # milk period are checked and either removed due to inconsistent answers or ok. # Joining with farm so I can compare weaning age depending on organic type wean_farm <- inner_join(wean, farm, by = c("id_respondant", "country")) # Select variables wean_farm2 <- select(wean_farm, id_respondant, country, organic_dummy, female_suckle_d, female_wean_wk, wean_suckle_diff) # Distribution of difference between weaning age and days suckling hist(wean_farm2$wean_suckle_diff, breaks = 14) # Create dummy for number of farms that let female calves suckle from birth # until weaning (cut-off of 10 was determined by checking data row-wise ) wean_farm2$suckle_until_wean <- ifelse(wean_farm2$wean_suckle_diff < 11, 1, 0) ### x-axis jittered plot age (wk) at weaning female calves p41 <- ggplot(wean_farm2, aes(x = country, y = female_wean_wk, colour = as.factor(organic_dummy))) + geom_jitter(size = 5, width = 0.2, height = 0, alpha = 0.4) + stat_summary(fun = median, fun.min = median, fun.max = median, geom = "crossbar", width = 0.6, color = "black", size = 1.1) + scale_y_continuous(expand = c(0, 0), limits = c(0, 52)) + labs(x = "Country", y = "Weaning age (wk)") + scale_colour_manual(values = c("orange2", "grey55"), breaks = c(0, 1), labels = c("Conventional", "Organic")) + theme_classic(base_size = 20) + theme(plot.margin = margin(0.9, 0.3, 0.0, 0.3, "cm"), axis.text = element_text(size = 16), legend.title = element_blank(), legend.key.size = unit(0.3, "cm"), title = element_text(size = 16, face = "bold"), legend.position = "bottom") p41 # Boxplot to evaluate if there is differences between organic and conventional # farms in weaning age for female calves boxplot(female_wean_wk ~ organic_dummy, data = wean_farm2) # Summaries per farming system wean.org <- filter(wean_farm2, organic_dummy == 1) summary(wean.org) sd(wean.org$female_wean_wk) wean.convent <- filter(wean_farm2, organic_dummy == 0) summary(wean.convent) sd(wean.convent$female_wean_wk) # Summaries per country wean.at <- filter(wean_farm2, country == "AT") summary(wean.at) wean.ch <- filter(wean_farm2, country == "CH") summary(wean.ch) wean.de <- filter(wean_farm2, country == "DE") summary(wean.de) wean.fr <- filter(wean_farm2, country == "FR") summary(wean.fr) wean.it <- filter(wean_farm2, country == "IT") summary(wean.it) wean.se <- filter(wean_farm2, country == "SE") summary(wean.se) rm(wean, wean_farm, wean_farm2, wean.org, wean.convent, wean.at, wean.ch, wean.de, wean.fr, wean.it, wean.se) # Q44. Perception cow health #### ### + Fertility + #### # Relevel factor order health$cow_fertility <- factor(health$cow_fertility, levels = c("better", "same", "worse", "dont_know")) # Summarize all data repro_all <- health %>% group_by(cow_fertility) %>% summarise(nbr_farms = length(id_respondant), prp_farms = nbr_farms/104) %>% ungroup() # Summarise per country pl_repro <- health %>% group_by(country, cow_fertility) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(health$country), unique(health$cow_fertility)) tmp <- tmp %>% rename(country = Var1, cow_fertility = Var2) pl_repro2 <- full_join(pl_repro, tmp) pl_repro2$nbr_farms[is.na(pl_repro2$nbr_farms)] <- 0 p44a <- ggplot(data = pl_repro2, aes(x = country, y = nbr_farms, fill = factor(cow_fertility, levels = c("dont_know", "worse", "same", "better")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Fertility", breaks = c("dont_know", "worse", "same", "better"), labels = c("Don't know", "Worse", "Same", "Better")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"), title = element_text(size = 12, face = "bold")) p44a ### + Udder health + #### # Relevel factor order health$cow_udder_health <- factor(health$cow_udder_health, levels = c("better", "same", "worse", "dont_know")) # Summarize all data udder_all <- health %>% group_by(cow_udder_health) %>% summarise(nbr_farms = length(id_respondant), prp_farms = nbr_farms/104) %>% ungroup() # Summarize per country pl_udder <- health %>% group_by(country, cow_udder_health) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(health$country), unique(health$cow_udder_health)) tmp <- tmp %>% rename(country = Var1, cow_udder_health = Var2) pl_udder2 <- full_join(pl_udder, tmp) pl_udder2$nbr_farms[is.na(pl_udder2$nbr_farms)] <- 0 p44b <- ggplot(data=pl_udder2, aes(x = country, y = nbr_farms, fill = factor(cow_udder_health, levels = c("dont_know", "worse", "same", "better")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Udder health", breaks = c("dont_know", "worse", "same", "better"), labels = c("Don't know", "Worse", "Same", "Better")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"), title = element_text(size = 12, face = "bold")) p44b # Q45. Perception calf health #### ### + General health + #### # Relevel factor order health$calf_health <- factor(health$calf_health, levels = c("better", "same", "worse", "dont_know")) # Summarize all data calfhealth_all <- health %>% group_by(calf_health) %>% summarise(nbr_farms = length(id_respondant), prp_farms = nbr_farms/104) %>% ungroup() # Summarize per country pl_calfhealth <- health %>% group_by(country, calf_health) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(health$country), unique(health$calf_health)) tmp <- tmp %>% rename(country = Var1, calf_health = Var2) pl_calfhealth2 <- full_join(pl_calfhealth, tmp) pl_calfhealth2$nbr_farms[is.na(pl_calfhealth2$nbr_farms)] <- 0 p44c <- ggplot(data = pl_calfhealth2, aes(x = country, y = nbr_farms, fill = factor(calf_health, levels = c("dont_know", "worse", "same", "better")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Calf health", breaks = c("dont_know", "worse", "same", "better"), labels = c("Don't know", "Worse", "Same", "Better")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"), title = element_text(size = 12, face = "bold")) p44c ### + Weight gain + #### # Rm farms that did not provide an answer for this question (n = 2) gain <- filter(health, calf_gain != "no_answer") gain$calf_gain <- as.factor(as.character(gain$calf_gain)) # Relevel factor order gain$calf_gain <- factor(gain$calf_gain, levels = c("higher", "same", "lower", "dont_know")) # Summarize all data calfgain_all <- gain %>% group_by(calf_gain) %>% summarise(nbr_farms = length(id_respondant), prp_farms = nbr_farms/102) %>% ungroup() # Summarize per country pl_calfgain <- gain %>% group_by(country, calf_gain) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(gain$country), unique(gain$calf_gain)) tmp <- tmp %>% rename(country = Var1, calf_gain = Var2) pl_calfgain2 <- full_join(pl_calfgain, tmp) pl_calfgain2$nbr_farms[is.na(pl_calfgain2$nbr_farms)] <- 0 p44d <- ggplot(data = pl_calfgain2, aes(x = country, y = nbr_farms, fill = factor(calf_gain, levels = c("dont_know", "lower", "same", "higher")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Weight gain", breaks = c("dont_know", "lower", "same", "higher"), labels = c("Don't know", "Lower", "Same", "Higher")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"), title = element_text(size = 12, face = "bold")) p44d ### + Diarrhea + #### # Note: specifically asked about calf diarrhea events with impaired general # condition # Relevel factor order health$calf_diarrhoea <- factor(health$calf_diarrhoea, levels = c("lower", "same", "higher", "dont_know")) # Summarize all data calfdiarrhoea_all <- health %>% group_by(calf_diarrhoea) %>% summarise(nbr_farms = length(id_respondant), prp_farms = nbr_farms/104) %>% ungroup() # Summarize per country pl_calfdiarrhoea <- health %>% group_by(country, calf_diarrhoea) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(health$country), unique(health$calf_diarrhoea)) tmp <- tmp %>% rename(country = Var1, calf_diarrhoea = Var2) pl_calfdiarrhoea2 <- full_join(pl_calfdiarrhoea, tmp) pl_calfdiarrhoea2$nbr_farms[is.na(pl_calfdiarrhoea2$nbr_farms)] <- 0 p44e <- ggplot(data = pl_calfdiarrhoea2, aes(x = country, y = nbr_farms, fill = factor(calf_diarrhoea, levels = c("dont_know", "higher", "same", "lower")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Diarrhoea", breaks = c("dont_know", "higher", "same", "lower"), labels = c("Don't know", "Higher", "Same", "Lower")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"), title = element_text(size = 12, face = "bold")) p44e ### + Respiratory disease + #### # Rm farms that did not provide an answer for this question (n = 1) cough <- filter(health, calf_cough != "no_answer") cough$calf_cough <- as.factor(as.character(cough$calf_cough)) # Relevel factor order cough$calf_cough <- factor(cough$calf_cough, levels = c("lower", "same", "higher", "dont_know")) # Summarize all data calfcough_all <- cough %>% group_by(calf_cough) %>% summarise(nbr_farms = length(id_respondant), prp_farms = nbr_farms/104) %>% ungroup() # Summarize per country pl_calfcough <- cough %>% group_by(country, calf_cough) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms types for countries that lack one or more type for plotting tmp <- expand.grid(unique(cough$country), unique(cough$calf_cough)) tmp <- tmp %>% rename(country = Var1, calf_cough = Var2) pl_calfcough2 <- full_join(pl_calfcough, tmp) pl_calfcough2$nbr_farms[is.na(pl_calfcough2$nbr_farms)] <- 0 p44f <- ggplot(data = pl_calfcough2, aes(x = country, y = nbr_farms, fill = factor(calf_cough, levels = c("dont_know", "higher", "same", "lower")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Respiratory \ndisease", breaks = c("dont_know", "higher", "same", "lower"), labels = c("Don't know", "Higher", "Same", "Lower")) + labs(x = "Country", y = "Percent of farms (%)") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.4, 0.0, 0.4, "cm"), axis.text = element_text(size = 12), legend.key.size = unit(0.3, "cm"), title = element_text(size = 12, face = "bold")) p44f rm(cough, gain, pl_calfcough, pl_calfcough2, pl_calfdiarrhoea, pl_calfdiarrhoea2, pl_calfgain, pl_calfgain2, pl_calfhealth, pl_calfhealth2, pl_repro, pl_repro2, pl_udder, pl_udder2) # Q48. Mortality 0 - 3 mo #### # Rm the farms that either did not provided information about how many calves # 0-3 mo old that on average dies per year (Q48 = one Austrian and all # French farms), or did not provide information about how many # calves were born the last 12 mo (Q8 = all Swedish farms). # n = 65 mort <- filter(health, calf_death_rate != 999) # Rm Italian farm 157259336, as the large discrepancies in the reported number # of calves born and raised the last 12 mo (fewer calves born in total than # female calves that are kept for recruitment) indicates too poor data to be # able to estimate mortality mort2 <- filter(mort, id_respondant != "157259336") # Select variables mort3 <- select(mort2, id_respondant, country, calf_annual_death, calf_annual_death_description, calf_death_rate) # Histograms all farms with estimated mortality rate p48a <- ggplot(mort3, aes(x = calf_death_rate * 100)) + geom_histogram(colour = "grey22", fill = "grey84", binwidth = 1, size = 1.1) + scale_y_continuous(expand = c(0, 0), limits = c(0, 17), breaks = c(0, 5, 10, 15)) + labs(x = "Calf mortility first 90 days (%)", y = "Number of farms") + theme_classic(base_size = 16) + theme(plot.margin = margin(0.9, 0.3, 0.0, 0.3, "cm"), axis.text = element_text(size = 12), title = element_text(size = 12, face = "bold")) p48a ### + Demographics high mortality farms + #### # Note: I will not compare mortality between organic and conventional farms, as # Italy is the only country with conventional farms that has mortality estimates # > 10% and the CCC-system in Italy was quite different compared to all other # countries mort_farm <- left_join(mort2, farm) # Select variables mort_farm2 <- select(mort_farm, id_respondant:deworm_3y, calf_death_rate, cows_tot:calves_born, own_rearing, income_milk:mixed_description, ms_type) # Create grouping variable annual mortality rate > 10% or not mort_farm2$high_mort <- factor(ifelse(mort_farm2$calf_death_rate >= 0.1, 1, 0)) mort_farm2$acreage <- (mort_farm2$ha_crop + mort_farm2$ha_perm_pasture + mort_farm2$ha_temp_pasture) high_mort <- filter(mort_farm2, high_mort == 1) low_mort <- filter(mort_farm2, high_mort == 0) summary(high_mort) summary(low_mort) # // Herd size #### # Dot-boxplot herd size p48b <- ggplot(mort_farm2, aes(x = high_mort, y = cows_tot)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.15, size = 3, alpha = 0.6) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Herd size") + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p48b # // Acreage #### ### Dot-boxplot crop land p48c <- ggplot(mort_farm2, aes(x = high_mort, y = ha_crop)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.15, size = 3, alpha = 0.6) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Hectar cropland") + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p48c ### Dot-boxplot pastures p48d <- ggplot(mort_farm2, aes(x = high_mort, y = ha_perm_pasture)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.15, size = 3, alpha = 0.6) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Hectar permanent pasture") + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p48d p48e <- ggplot(mort_farm2, aes(x = high_mort, y = ha_temp_pasture)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.15, size = 3, alpha = 0.6) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Hectar temporary pasture") + theme_classic(base_size = 22) + theme(axis.text = element_text(size = 18), title = element_text(size = 20, face = "bold")) p48e mort_farm_suckle <- left_join(mort_farm2, suckle) # Select variables mort_farm_suckle2 <- select(mort_farm_suckle, id_respondant:deworm_3y, calf_death_rate, cows_tot:calves_born, own_rearing, income_milk:mixed_description, ms_type, female_suckle_d) # // household income #### # Relevel factor order mort_farm2$income_milk <- factor(mort_farm2$income_milk, levels = c("unspecified", "76-100", "51-75", "26-50", "0-25")) # Summarize pl_mort_profit <- mort_farm2 %>% group_by(high_mort, income_milk) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms income levels for mortality levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$income_milk)) tmp <- tmp %>% rename(high_mort = Var1, income_milk = Var2) pl_mort_profit2 <- full_join(pl_mort_profit, tmp) pl_mort_profit2$nbr_farms[is.na(pl_mort_profit2$nbr_farms)] <- 0 p48f <- ggplot(data = pl_mort_profit2, aes(x = high_mort, y = nbr_farms, fill = income_milk)) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("black", "forestgreen", "gold1", "orange2", "orangered3"), name = "Share of income", breaks = c("unspecified", "76-100", "51-75", "26-50", "0-25"), labels = c("Unknown", "76-100%", "51-75%", "26-50%", "0-25%")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48f # // ms-type #### # Summarize pl_mort_ms <- mort_farm2 %>% group_by(high_mort, ms_type) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms ms levels for mortality levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$ms_type)) tmp <- tmp %>% rename(high_mort = Var1, ms_type = Var2) pl_mort_ms2 <- full_join(pl_mort_ms, tmp) pl_mort_ms2$nbr_farms[is.na(pl_mort_ms2$nbr_farms)] <- 0 p48g <- ggplot(data = pl_mort_ms2, aes(x = high_mort, y = nbr_farms, fill = ms_type)) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("black", "forestgreen", "gold1", "orange2", "orangered3", "grey"), name = "MS", breaks = c("ams", "bucket", "hand", "other", "parlour", "tie"), labels = c("AMS", "Bucket", "Hand", "Other", "Parlour", "Pipeline")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48g # // housing #### # Summarize pl_mort_house <- mort_farm2 %>% group_by(high_mort, housing) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms housing for mortality levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$housing)) tmp <- tmp %>% rename(high_mort = Var1, housing = Var2) pl_mort_house2 <- full_join(pl_mort_house, tmp) pl_mort_house2$nbr_farms[is.na(pl_mort_house2$nbr_farms)] <- 0 p48h <- ggplot(data = pl_mort_house2, aes(x = high_mort, y = nbr_farms, fill = housing)) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("forestgreen", "gold1", "orange2", "orangered3", "grey"), name = "Cow housing", breaks = c("freestall", "loosehouse", "mixed", "outside", "tiestall"), labels = c("Free stall", "Open pack", "Mixed", "Outside", "Tie stall")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48h # // rearing #### # Select rearing variables rear_for_mort <- select(rear_type, id_respondant:other_rearing_description) mort_rear <- left_join(mort_farm2, rear_for_mort) # Summarize pl_mort_rear <- mort_rear %>% group_by(high_mort, rearing) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms rearing levels for mortality levels that lack levels for plotting tmp <- expand.grid(unique(mort_rear$high_mort), unique(mort_rear$rearing)) tmp <- tmp %>% rename(high_mort = Var1, rearing = Var2) pl_mort_rear2 <- full_join(pl_mort_rear, tmp) pl_mort_rear2$nbr_farms[is.na(pl_mort_rear2$nbr_farms)] <- 0 p48i <- ggplot(data = pl_mort_rear, aes(x = high_mort, y = nbr_farms, fill = rearing)) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("forestgreen", "gold1", "orange2", "orangered3", "grey"), name = "CCC rearing", breaks = c("foster", "mix", "mother", "mother_bulk", "other"), labels = c("Foster", "Mix", "Dam", "Dam + bulk", "Other")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48i # // calf health #### ### General health pl_mort_gen <- mort_farm2 %>% group_by(high_mort, calf_health) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms health levels for mortality levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$calf_health)) tmp <- tmp %>% rename(high_mort = Var1, calf_health = Var2) pl_mort_gen2 <- full_join(pl_mort_gen, tmp) pl_mort_gen2$nbr_farms[is.na(pl_mort_gen2$nbr_farms)] <- 0 p48j <- ggplot(data = pl_mort_gen2, aes(x = high_mort, y = nbr_farms, fill = factor(calf_health, levels = c("dont_know", "worse", "same", "better")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Calf health", breaks = c("dont_know", "worse", "same", "better"), labels = c("Don't know", "Worse", "Same", "Better")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48j ### Weight gain pl_mort_adg <- mort_farm2 %>% group_by(high_mort, calf_gain) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms gain levels for mort levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$calf_gain)) tmp <- tmp %>% rename(high_mort = Var1, calf_gain = Var2) pl_mort_adg2 <- full_join(pl_mort_adg, tmp) pl_mort_adg2$nbr_farms[is.na(pl_mort_adg2$nbr_farms)] <- 0 p48k <- ggplot(data = pl_mort_adg2, aes(x = high_mort, y = nbr_farms, fill = factor(calf_gain, levels = c("no_answer", "dont_know", "lower", "same", "higher")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("black", "gainsboro", "firebrick", "gold", "darkgreen"), name= "Calf gain", breaks = c("no_answer", "dont_know", "lower", "same", "higher"), labels = c("Unknown", "Don't know", "Lower", "Same", "Higher")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48k ### Diarrhea pl_mort_poop <- mort_farm2 %>% group_by(high_mort, calf_diarrhoea) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back diarrhoea levels for mort levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$calf_diarrhoea)) tmp <- tmp %>% rename(high_mort = Var1, calf_diarrhoea = Var2) pl_mort_poop2 <- full_join(pl_mort_poop, tmp) pl_mort_poop2$nbr_farms[is.na(pl_mort_poop2$nbr_farms)] <- 0 p48l <- ggplot(data = pl_mort_poop2, aes(x = high_mort, y = nbr_farms, fill = factor(calf_diarrhoea, levels = c("dont_know", "higher", "same", "lower")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("gainsboro", "firebrick", "gold", "darkgreen"), name= "Diarrhoea", breaks = c("dont_know", "higher", "same", "lower"), labels = c("Don't know", "Higher", "Same", "Lower")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48l ### Cough pl_mort_cough <- mort_farm2 %>% group_by(high_mort, calf_cough) %>% summarise(nbr_farms = length(id_respondant)) %>% ungroup() # Put back ms poop levels for mort levels that lack levels for plotting tmp <- expand.grid(unique(mort_farm2$high_mort), unique(mort_farm2$calf_cough)) tmp <- tmp %>% rename(high_mort = Var1, calf_cough = Var2) pl_mort_cough2 <- full_join(pl_mort_cough, tmp) pl_mort_cough2$nbr_farms[is.na(pl_mort_cough2$nbr_farms)] <- 0 p48m <- ggplot(data = pl_mort_cough2, aes(x = high_mort, y = nbr_farms, fill = factor(calf_cough, levels = c("no_answer", "dont_know", "higher", "same", "lower")))) + geom_bar(colour = "black", position = "fill", stat = "identity", width = 0.65) + scale_y_continuous(expand = c(0, 0), breaks = c(0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1), labels = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) + scale_fill_manual(values = c("black", "gainsboro", "firebrick", "gold", "darkgreen"), name= "Cough", breaks = c("no_answer", "dont_know", "higher", "same", "lower"), labels = c("Unknown", "Don't know", "Higher", "Same", "Lower")) + labs(x = "Low vs high annual mortality 0-3 mo", y = "Percent of farms") + theme_classic(base_size = 18) + theme(axis.text = element_text(size = 18), title = element_text(size = 18, face = "bold")) p48m rm(mort, mort2, mort3, mort_farm, mort_farm2, mort_farm_suckle, rear_for_mort, mort_farm_suckle2, mort_rear, pl_mort_adg, pl_mort_adg2, pl_mort_cough, pl_mort_cough2, pl_mort_gen, pl_mort_gen2, pl_mort_house, pl_mort_house2, pl_mort_ms, pl_mort_ms2, pl_mort_poop, pl_mort_poop2, pl_mort_profit, pl_mort_profit2, pl_mort_rear, pl_mort_rear2) # +++ Effect of contact period +++ #### # Extract Q2 from suckle to evaluate rearing type, CCC allowance, separation # distress, farmer health perception and farmer drivers/barriers for CCC q2 <- select(suckle, id_respondant:days_suckle, when2) # Extract rearing and separation distress rear_stress <- select(rear_type, id_respondant:other_rearing_description, cow_moo:separation_other_description) # Extract health perception health_dur <- select(health, id_respondant:calf_cough, calf_death_rate) # Extract drivers and barriers drivers2 <- select(drivers, id_respondant, country, driver_natural:barrier_other_description) # Combine data sets duration1 <- full_join(q2, rear_stress) duration2 <- full_join(duration1, health_dur) duration3 <- full_join(duration2, drivers2) # Reclassify farms 157259329-IT that in Q2 states 21 days suckling period, but # then specify that female calves suckle 50 days duration3$days_suckle2 <- ifelse(duration3$id_respondant == "157259329", 50, duration3$days_suckle) short_dur <- filter(duration3, days_suckle2 < 29) medium_dur <- filter(duration3, days_suckle2 > 28, days_suckle2 < 91) long_dur <- filter(duration3, days_suckle2 > 90)