Arithmetic Mean Method: Individual Cow and Herd-Level Means (≥2 Measurements)
Direct Arithmetic Mean Calculation (Manafiazar et al. 2017)
2026-01-13
- 1 OBJECTIVE
- 2 STEP 1: LOAD DATA
- 3 STEP 2: COUNT VISITS AND APPLY MINIMUM VISIT THRESHOLD
- 4 STEP 3: ARITHMETIC MEAN CALCULATION - INDIVIDUAL COWS
- 5 STEP 4: APPLY ≥2 MEASUREMENTS FILTER
- 6 STEP 5: TREATMENT-LEVEL (HERD) MEANS - ARITHMETIC MEAN
- 7 STEP 6: OVERALL HERD MEAN
- 8 STEP 7: WITHIN-ANIMAL VARIATION METRICS
- 9 STEP 8: STATISTICAL COMPARISON
- 10 STEP 9: EFFECT SIZE AND POWER ANALYSIS
- 11 STEP 10: EXPORT RESULTS
- 12 STEP 11: VISUALIZATIONS
- 13 SUMMARY
1 OBJECTIVE
Calculate individual cow means and herd-level (treatment) means using Arithmetic Mean (AM) method per Manafiazar et al. (2017), after removing 20 poor-quality night visits. Report only cows with ≥2 valid measurements for statistical validity. This approach retains more cows than TOD method.
2 STEP 1: LOAD DATA
# Load filtered dataset (3SD outliers removed)
gf_data <- read_excel(
"gf_full_dataset_3SD_filtered.xlsx",
sheet = 1
)
# Display structure
cat("DATA LOADING:\n")## DATA LOADING:cat(strrep("=", 80), "\n")## ================================================================================cat("Original data dimensions:", nrow(gf_data), "rows x", ncol(gf_data), "columns\n")## Original data dimensions: 3355 rows x 21 columnscat("Column names:", colnames(gf_data), "\n\n")## Column names: cow_id grazing_treatment start_time end_time good_data_duration hour_of_day co2_massflow_g_d ch4_massflow_g_d o2_massflow_g_d average_airflow_l_s airflow_cf wind_cf Time_Bin Week_Number Filtered_Status row_id GF_Week Phase Window_14d Visit_Date Cyclehead(gf_data, 10)## # A tibble: 10 × 21
## cow_id grazing_treatment start_time end_time
## <chr> <chr> <dttm> <dttm>
## 1 252H K2_CON 2025-06-27 14:55:30 2025-06-27 14:58:38
## 2 252H K2_CON 2025-06-27 19:47:12 2025-06-27 19:52:15
## 3 252H K2_CON 2025-06-27 22:03:13 2025-06-27 22:05:37
## 4 252H K2_CON 2025-06-28 00:49:31 2025-06-28 00:54:25
## 5 252H K2_CON 2025-06-28 19:20:02 2025-06-28 19:26:41
## 6 252H K2_CON 2025-06-29 01:07:48 2025-06-29 01:10:02
## 7 252H K2_CON 2025-06-29 07:01:34 2025-06-29 07:06:10
## 8 252H K2_CON 2025-06-29 13:06:27 2025-06-29 13:11:07
## 9 252H K2_CON 2025-06-29 20:07:00 2025-06-29 20:10:42
## 10 252H K2_CON 2025-06-30 01:23:07 2025-06-30 01:27:03
## # ℹ 17 more variables: good_data_duration <dttm>, hour_of_day <dbl>,
## # co2_massflow_g_d <dbl>, ch4_massflow_g_d <dbl>, o2_massflow_g_d <dbl>,
## # average_airflow_l_s <dbl>, airflow_cf <dbl>, wind_cf <dbl>, Time_Bin <chr>,
## # Week_Number <dbl>, Filtered_Status <chr>, row_id <dbl>, GF_Week <chr>,
## # Phase <chr>, Window_14d <chr>, Visit_Date <dttm>, Cycle <chr>3 STEP 2: COUNT VISITS AND APPLY MINIMUM VISIT THRESHOLD
# Count visits (unique dates) per cow
visit_counts <- gf_data %>%
group_by(cow_id, grazing_treatment) %>%
summarize(Visit_Count = n_distinct(Visit_Date), .groups = 'drop')
cat("STEP 2: MINIMUM VISIT THRESHOLD (>=20 VISITS)\n")## STEP 2: MINIMUM VISIT THRESHOLD (>=20 VISITS)cat(strrep("=", 80), "\n")## ================================================================================cat("Before filtering: ", nrow(gf_data), " records from ", n_distinct(gf_data$cow_id), " cows\n", sep="")## Before filtering: 3355 records from 37 cowscat("Cows with >=20 visits: ", sum(visit_counts$Visit_Count >= 20), "\n", sep="")## Cows with >=20 visits: 34cat("Cows with <20 visits: ", sum(visit_counts$Visit_Count < 20), "\n\n", sep="")## Cows with <20 visits: 3# Identify and display cows with <20 visits
excluded_cows_visits <- visit_counts %>%
filter(Visit_Count < 20) %>%
arrange(Visit_Count)
if(nrow(excluded_cows_visits) > 0) {
cat("EXCLUDED COWS (<20 visits):\n")
print(excluded_cows_visits)
cat("\nCows to exclude: ", paste(excluded_cows_visits$cow_id, collapse=", "), "\n\n", sep="")
}## EXCLUDED COWS (<20 visits):
## # A tibble: 3 × 3
## cow_id grazing_treatment Visit_Count
## <chr> <chr> <int>
## 1 338D K1_AMP 12
## 2 374E K1_AMP 12
## 3 422J K2_CON 12
##
## Cows to exclude: 338D, 374E, 422J# Filter data to retain only cows with >=20 visits
valid_cow_ids <- visit_counts %>%
filter(Visit_Count >= 20) %>%
pull(cow_id)
gf_data_cleaned <- gf_data %>%
filter(cow_id %in% valid_cow_ids) %>%
filter(!is.na(ch4_massflow_g_d)) %>% # Remove records with missing CH4 values
filter(ch4_massflow_g_d > 0) # Remove zero or negative values
cat("After visit count filtering: ", nrow(gf_data_cleaned), " records from ", n_distinct(gf_data_cleaned$cow_id), " cows\n", sep="")## After visit count filtering: 3306 records from 34 cows4 STEP 3: ARITHMETIC MEAN CALCULATION - INDIVIDUAL COWS
# Calculate arithmetic mean for each cow (direct mean of all valid measurements)
individual_cow_means <- gf_data_cleaned %>%
group_by(cow_id, grazing_treatment) %>%
summarize(
Measurements = n(),
Mean_CH4 = mean(ch4_massflow_g_d, na.rm = TRUE),
Std_Dev = sd(ch4_massflow_g_d, na.rm = TRUE),
CV_percent = (Std_Dev / Mean_CH4) * 100,
Min_CH4 = min(ch4_massflow_g_d, na.rm = TRUE),
Max_CH4 = max(ch4_massflow_g_d, na.rm = TRUE),
Range = Max_CH4 - Min_CH4,
.groups = 'drop'
) %>%
arrange(grazing_treatment, desc(Mean_CH4)) %>%
mutate(across(where(is.numeric), ~round(., 2)))
cat("ARITHMETIC MEAN - INDIVIDUAL COW STATISTICS\n")## ARITHMETIC MEAN - INDIVIDUAL COW STATISTICScat(strrep("=", 80), "\n\n")## ================================================================================# Count cows by number of measurements
measurement_distribution <- individual_cow_means %>%
group_by(Measurements) %>%
summarize(N_Cows = n(), .groups = 'drop')
cat("MEASUREMENT DISTRIBUTION:\n")## MEASUREMENT DISTRIBUTION:print(measurement_distribution)## # A tibble: 31 × 2
## Measurements N_Cows
## <dbl> <int>
## 1 45 1
## 2 60 1
## 3 61 1
## 4 67 1
## 5 72 1
## 6 78 1
## 7 79 1
## 8 80 1
## 9 82 1
## 10 85 2
## # ℹ 21 more rowscat("\n")# Show all cows
print(individual_cow_means, n = Inf)## # A tibble: 34 × 9
## cow_id grazing_treatment Measurements Mean_CH4 Std_Dev CV_percent Min_CH4
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 290G K1_AMP 92 359. 87.1 24.3 184.
## 2 266H K1_AMP 97 351. 89.1 25.4 76.8
## 3 364B K1_AMP 78 347. 99.7 28.7 72.4
## 4 312E K1_AMP 98 330. 79.2 24.0 103.
## 5 322F K1_AMP 103 329. 91.2 27.8 55.9
## 6 454C K1_AMP 72 327. 83.8 25.7 62.5
## 7 322C K1_AMP 61 326. 100. 30.8 119.
## 8 370F K1_AMP 91 323. 84.4 26.1 56.8
## 9 338J K1_AMP 111 319. 69.3 21.7 168.
## 10 314J K1_AMP 132 318. 71.4 22.4 163.
## 11 410C K1_AMP 82 305. 87.3 28.6 145.
## 12 426D K1_AMP 85 304. 86.2 28.3 114.
## 13 370H K1_AMP 102 298. 79.0 26.5 109.
## 14 286J K1_AMP 109 290 88.0 30.3 54.1
## 15 264F K1_AMP 122 290. 69.2 23.9 114.
## 16 436E K1_AMP 89 283. 87.8 31.0 78.6
## 17 258K K1_AMP 121 251. 75.8 30.2 108.
## 18 306K K1_AMP 79 220. 56.8 25.9 47.2
## 19 252H K2_CON 134 382. 83.4 21.8 147.
## 20 270G K2_CON 125 379. 83.3 22.0 139.
## 21 398G K2_CON 60 375. 83.3 22.2 232.
## 22 368J K2_CON 85 368. 91.1 24.7 129.
## 23 296H K2_CON 101 359. 86.7 24.1 49.7
## 24 312F K2_CON 146 344. 77.5 22.5 161.
## 25 332H K2_CON 100 343. 87.3 25.4 66.1
## 26 266D K2_CON 80 337. 98.6 29.2 116.
## 27 404E K2_CON 91 330. 76.7 23.3 199.
## 28 326G K2_CON 100 323. 88.5 27.4 68.6
## 29 400G K2_CON 67 318. 83.9 26.4 162.
## 30 260K K2_CON 126 309. 67.6 21.8 162.
## 31 296G K2_CON 108 308. 76.6 24.9 161.
## 32 442B K2_CON 45 304. 87.2 28.7 49.7
## 33 260C K2_CON 120 300. 83.6 27.9 150.
## 34 446K K2_CON 94 274. 76.8 28.0 71.5
## # ℹ 2 more variables: Max_CH4 <dbl>, Range <dbl>5 STEP 4: APPLY ≥2 MEASUREMENTS FILTER
# Filter to retain only cows with ≥2 measurements
individual_cow_means_filtered <- individual_cow_means %>%
filter(Measurements >= 2)
cat("\nFILTERING: COWS WITH ≥2 MEASUREMENTS\n")##
## FILTERING: COWS WITH ≥2 MEASUREMENTScat(strrep("=", 80), "\n")## ================================================================================cat("Total cows before filtering: ", nrow(individual_cow_means), "\n", sep="")## Total cows before filtering: 34cat("Total cows after filtering: ", nrow(individual_cow_means_filtered), "\n", sep="")## Total cows after filtering: 34# Identify and display excluded cows
excluded_cows <- individual_cow_means %>%
filter(Measurements < 2) %>%
select(cow_id, grazing_treatment, Measurements, Mean_CH4)
if(nrow(excluded_cows) > 0) {
cat("\nEXCLUDED COWS (n=", nrow(excluded_cows), "):\n", sep="")
print(excluded_cows)
cat("\n")
}
# Split by treatment
k1_cows <- individual_cow_means_filtered %>% filter(grazing_treatment == "K1_AMP")
k2_cows <- individual_cow_means_filtered %>% filter(grazing_treatment == "K2_CON")
cat("SAMPLE SIZE BY TREATMENT:\n")## SAMPLE SIZE BY TREATMENT:cat("K1_AMP: ", nrow(k1_cows), " cows\n", sep="")## K1_AMP: 18 cowscat("K2_CON: ", nrow(k2_cows), " cows\n", sep="")## K2_CON: 16 cowscat("Total: ", nrow(individual_cow_means_filtered), " cows\n\n", sep="")## Total: 34 cows6 STEP 5: TREATMENT-LEVEL (HERD) MEANS - ARITHMETIC MEAN
# Method: Calculate treatment mean from individual cow means
# This gives equal weight to each cow regardless of measurement frequency
# Approach 1: Mean of individual cow means (each cow weighted equally)
treatment_summary_cow_level <- individual_cow_means_filtered %>%
group_by(grazing_treatment) %>%
summarize(
Unique_Cows = n_distinct(cow_id),
Mean_CH4_CowLevel = mean(Mean_CH4, na.rm = TRUE),
SD_Between_Cows = sd(Mean_CH4, na.rm = TRUE),
SEM_Cows = sd(Mean_CH4, na.rm = TRUE) / sqrt(n()),
Min_Cow_Mean = min(Mean_CH4, na.rm = TRUE),
Max_Cow_Mean = max(Mean_CH4, na.rm = TRUE),
.groups = 'drop'
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
# Approach 2: Direct mean from all measurements
total_measurements_by_treatment <- gf_data_cleaned %>%
group_by(grazing_treatment) %>%
summarize(
Total_Observations = n(),
Mean_CH4_AllMeas = mean(ch4_massflow_g_d, na.rm = TRUE),
SD_AllMeas = sd(ch4_massflow_g_d, na.rm = TRUE),
SEM_AllMeas = sd(ch4_massflow_g_d, na.rm = TRUE) / sqrt(n()),
.groups = 'drop'
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
# Combine both approaches
treatment_summary <- treatment_summary_cow_level %>%
left_join(total_measurements_by_treatment, by = "grazing_treatment")
cat("TREATMENT-LEVEL MEANS (ARITHMETIC MEAN METHOD)\n")## TREATMENT-LEVEL MEANS (ARITHMETIC MEAN METHOD)cat(strrep("=", 80), "\n\n")## ================================================================================print(treatment_summary)## # A tibble: 2 × 11
## grazing_treatment Unique_Cows Mean_CH4_CowLevel SD_Between_Cows SEM_Cows
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 K1_AMP 18 309. 34.7 8.18
## 2 K2_CON 16 335. 32.1 8.02
## # ℹ 6 more variables: Min_Cow_Mean <dbl>, Max_Cow_Mean <dbl>,
## # Total_Observations <dbl>, Mean_CH4_AllMeas <dbl>, SD_AllMeas <dbl>,
## # SEM_AllMeas <dbl>cat("\nINTERPRETATION:\n")##
## INTERPRETATION:cat("- Mean_CH4_CowLevel: Average of individual cow means (each cow weighted equally)\n")## - Mean_CH4_CowLevel: Average of individual cow means (each cow weighted equally)cat("- Mean_CH4_AllMeas: Direct mean of all measurements (may favor cows with more records)\n")## - Mean_CH4_AllMeas: Direct mean of all measurements (may favor cows with more records)cat("- SD_Between_Cows: Variation in cow means (between-animal variation)\n")## - SD_Between_Cows: Variation in cow means (between-animal variation)cat("- SD_AllMeas: Variation in individual measurements (within + between-animal)\n\n")## - SD_AllMeas: Variation in individual measurements (within + between-animal)# Calculate treatment difference
k1_mean <- treatment_summary$Mean_CH4_CowLevel[treatment_summary$grazing_treatment == "K1_AMP"]
k2_mean <- treatment_summary$Mean_CH4_CowLevel[treatment_summary$grazing_treatment == "K2_CON"]
difference <- k2_mean - k1_mean
percent_diff <- (difference / k1_mean) * 100
cat("TREATMENT COMPARISON (Cow-Level Means):\n")## TREATMENT COMPARISON (Cow-Level Means):cat(strrep("-", 80), "\n")## --------------------------------------------------------------------------------cat(sprintf("K1_AMP Mean: %.2f g/d\n", k1_mean))## K1_AMP Mean: 309.36 g/dcat(sprintf("K2_CON Mean: %.2f g/d\n", k2_mean))## K2_CON Mean: 334.67 g/dcat(sprintf("Difference: %.2f g/d\n", difference))## Difference: 25.31 g/dcat(sprintf("Percent Change: %.2f%%\n", percent_diff))## Percent Change: 8.18%7 STEP 6: OVERALL HERD MEAN
# Overall herd statistics (from all valid measurements)
herd_mean <- mean(gf_data_cleaned$ch4_massflow_g_d, na.rm = TRUE)
herd_sd <- sd(gf_data_cleaned$ch4_massflow_g_d, na.rm = TRUE)
herd_n_meas <- nrow(gf_data_cleaned)
herd_cows <- n_distinct(gf_data_cleaned$cow_id)
cat("OVERALL HERD STATISTICS\n")## OVERALL HERD STATISTICScat(strrep("=", 80), "\n")## ================================================================================cat(sprintf("Total Cows (>=2 measurements): %d\n", herd_cows))## Total Cows (>=2 measurements): 34cat(sprintf("Total Valid Measurements: %d\n", herd_n_meas))## Total Valid Measurements: 3306cat(sprintf("Overall Mean CH4: %.2f g/d\n", herd_mean))## Overall Mean CH4: 321.58 g/dcat(sprintf("Overall Std Dev: %.2f g/d\n", herd_sd))## Overall Std Dev: 89.22 g/dcat(sprintf("Herd Range: %.2f - %.2f g/d\n",
min(gf_data_cleaned$ch4_massflow_g_d, na.rm = TRUE),
max(gf_data_cleaned$ch4_massflow_g_d, na.rm = TRUE)))## Herd Range: 47.15 - 591.19 g/d8 STEP 7: WITHIN-ANIMAL VARIATION METRICS
# Calculate within-animal variation for each cow
within_animal_variation <- individual_cow_means_filtered %>%
group_by(grazing_treatment) %>%
summarize(
N_Cows = n(),
Mean_Std_Dev = mean(Std_Dev, na.rm = TRUE),
SD_of_Std_Dev = sd(Std_Dev, na.rm = TRUE),
Mean_CV_percent = mean(CV_percent, na.rm = TRUE),
SD_of_CV = sd(CV_percent, na.rm = TRUE),
Mean_Range = mean(Range, na.rm = TRUE),
SD_of_Range = sd(Range, na.rm = TRUE),
.groups = 'drop'
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
cat("\nWITHIN-ANIMAL VARIATION (Measurement Variability Within Each Cow)\n")##
## WITHIN-ANIMAL VARIATION (Measurement Variability Within Each Cow)cat(strrep("=", 80), "\n\n")## ================================================================================print(within_animal_variation)## # A tibble: 2 × 8
## grazing_treatment N_Cows Mean_Std_Dev SD_of_Std_Dev Mean_CV_percent SD_of_CV
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 K1_AMP 18 82.6 11.0 26.8 2.87
## 2 K2_CON 16 83.3 7.18 25.0 2.63
## # ℹ 2 more variables: Mean_Range <dbl>, SD_of_Range <dbl>cat("\nINTERPRETATION:\n")##
## INTERPRETATION:cat("- Mean_Std_Dev: Average measurement SD within each cow\n")## - Mean_Std_Dev: Average measurement SD within each cowcat("- Mean_CV_percent: Average coefficient of variation (normalized)\n")## - Mean_CV_percent: Average coefficient of variation (normalized)cat("- Mean_Range: Average difference between max and min measurements per cow\n")## - Mean_Range: Average difference between max and min measurements per cow9 STEP 8: STATISTICAL COMPARISON
# Two-sample t-test on cow-level means
t_test_result <- t.test(
individual_cow_means_filtered$Mean_CH4[individual_cow_means_filtered$grazing_treatment == "K1_AMP"],
individual_cow_means_filtered$Mean_CH4[individual_cow_means_filtered$grazing_treatment == "K2_CON"]
)
# Two-sample t-test on all measurements
t_test_all_meas <- t.test(
gf_data_cleaned$ch4_massflow_g_d[gf_data_cleaned$grazing_treatment == "K1_AMP"],
gf_data_cleaned$ch4_massflow_g_d[gf_data_cleaned$grazing_treatment == "K2_CON"]
)
# Wilcoxon rank-sum test (non-parametric alternative)
wilcox_test_result <- wilcox.test(
individual_cow_means_filtered$Mean_CH4[individual_cow_means_filtered$grazing_treatment == "K1_AMP"],
individual_cow_means_filtered$Mean_CH4[individual_cow_means_filtered$grazing_treatment == "K2_CON"]
)
cat("STATISTICAL TESTS\n")## STATISTICAL TESTScat(strrep("=", 80), "\n\n")## ================================================================================cat("Test 1: Two-Sample t-test (Individual Cow Means)\n")## Test 1: Two-Sample t-test (Individual Cow Means)cat("-" , strrep("-", 76), "\n", sep="")## -----------------------------------------------------------------------------cat("t-statistic: ", round(t_test_result$statistic, 4), "\n")## t-statistic: -2.2094cat("p-value: ", format(t_test_result$p.value, scientific = FALSE), "\n")## p-value: 0.03443767cat("Significant: ", if(t_test_result$p.value < 0.05) "YES (p < 0.05)" else "NO (p ≥ 0.05)", "\n\n")## Significant: YES (p < 0.05)cat("Test 2: t-test (All Measurements)\n")## Test 2: t-test (All Measurements)cat("-" , strrep("-", 76), "\n", sep="")## -----------------------------------------------------------------------------cat("t-statistic: ", round(t_test_all_meas$statistic, 4), "\n")## t-statistic: -9.0076cat("p-value: ", format(t_test_all_meas$p.value, scientific = FALSE), "\n")## p-value: 0.0000000000000000003491114cat("Significant: ", if(t_test_all_meas$p.value < 0.05) "YES (p < 0.05)" else "NO (p ≥ 0.05)", "\n\n")## Significant: YES (p < 0.05)cat("Test 3: Wilcoxon Rank-Sum Test (Non-parametric)\n")## Test 3: Wilcoxon Rank-Sum Test (Non-parametric)cat("-" , strrep("-", 76), "\n", sep="")## -----------------------------------------------------------------------------cat("W-statistic: ", round(wilcox_test_result$statistic, 4), "\n")## W-statistic: 89cat("p-value: ", format(wilcox_test_result$p.value, scientific = FALSE), "\n")## p-value: 0.0593928cat("Significant: ", if(wilcox_test_result$p.value < 0.05) "YES (p < 0.05)" else "NO (p ≥ 0.05)", "\n")## Significant: NO (p ≥ 0.05)10 STEP 9: EFFECT SIZE AND POWER ANALYSIS
# Calculate Cohen's d (effect size)
k1_means <- individual_cow_means_filtered$Mean_CH4[individual_cow_means_filtered$grazing_treatment == "K1_AMP"]
k2_means <- individual_cow_means_filtered$Mean_CH4[individual_cow_means_filtered$grazing_treatment == "K2_CON"]
mean_diff <- mean(k2_means, na.rm = TRUE) - mean(k1_means, na.rm = TRUE)
k1_sd <- sd(k1_means, na.rm = TRUE)
k2_sd <- sd(k2_means, na.rm = TRUE)
# Avoid NaN by checking for valid SD values
if(!is.na(k1_sd) && !is.na(k2_sd) && k1_sd > 0 && k2_sd > 0) {
pooled_sd <- sqrt(((length(k1_means) - 1) * k1_sd^2 +
(length(k2_means) - 1) * k2_sd^2) /
(length(k1_means) + length(k2_means) - 2))
cohens_d <- mean_diff / pooled_sd
cat("\nEFFECT SIZE (Cohen's d)\n")
cat(strrep("=", 80), "\n")
cat(sprintf("Cohen's d: %.4f\n", cohens_d))
interpretation <- if(is.na(cohens_d)) "Unable to calculate"
else if(abs(cohens_d) < 0.2) "Negligible"
else if(abs(cohens_d) < 0.5) "Small"
else if(abs(cohens_d) < 0.8) "Medium"
else "Large"
cat("Interpretation: ", interpretation, "\n")
} else {
cat("\nEFFECT SIZE (Cohen's d)\n")
cat(strrep("=", 80), "\n")
cat("Unable to calculate Cohen's d due to zero or missing standard deviations\n")
cohens_d <- NA
interpretation <- "Unable to calculate"
}##
## EFFECT SIZE (Cohen's d)
## ================================================================================
## Cohen's d: 0.7556
## Interpretation: Medium11 STEP 10: EXPORT RESULTS
# Create comprehensive export workbook
export_list <- list(
"Individual_Cow_Means" = individual_cow_means_filtered,
"Treatment_Summary" = treatment_summary,
"Within_Animal_Variation" = within_animal_variation,
"All_Cows_Including_Excluded" = individual_cow_means,
"Data_Quality_Notes" = data.frame(
Metric = c(
"Original Records",
"After Data Cleaning",
"Records Removed",
"Unique Cows (All)",
"Unique Cows (≥2 Meas)",
"Cows Excluded",
"K1_AMP Sample Size",
"K2_CON Sample Size",
"K1_AMP Mean (g/d)",
"K2_CON Mean (g/d)",
"Difference (g/d)",
"Percent Difference (%)",
"Treatment p-value",
"Cohen's d",
"Effect Size"
),
Value = c(
nrow(gf_data),
herd_n_meas,
nrow(gf_data) - herd_n_meas,
nrow(individual_cow_means),
nrow(individual_cow_means_filtered),
nrow(individual_cow_means) - nrow(individual_cow_means_filtered),
nrow(k1_cows),
nrow(k2_cows),
round(k1_mean, 2),
round(k2_mean, 2),
round(difference, 2),
round(percent_diff, 2),
format(t_test_result$p.value, scientific = FALSE),
if(is.na(cohens_d)) "NA" else round(cohens_d, 4),
interpretation
)
)
)
# Write to Excel
write_xlsx(
export_list,
"AM_Method_Individual_and_Herd_Means_FILTERED_GE2_RESULTS.xlsx"
)
cat("\nRESULTS EXPORTED\n")##
## RESULTS EXPORTEDcat(strrep("=", 80), "\n")## ================================================================================cat("File: AM_Method_Individual_and_Herd_Means_FILTERED_GE2_RESULTS.xlsx\n")## File: AM_Method_Individual_and_Herd_Means_FILTERED_GE2_RESULTS.xlsxcat("Sheets: 5 (Individual_Cow_Means, Treatment_Summary, Within_Animal_Variation,\n")## Sheets: 5 (Individual_Cow_Means, Treatment_Summary, Within_Animal_Variation,cat(" All_Cows_Including_Excluded, Data_Quality_Notes)\n")## All_Cows_Including_Excluded, Data_Quality_Notes)12 STEP 11: VISUALIZATIONS
12.1 Distribution by Treatment
ggplot(gf_data_cleaned, aes(x = grazing_treatment, y = ch4_massflow_g_d, fill = grazing_treatment)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.4, size = 1.5) +
theme_minimal() +
labs(
title = "CH4 Emissions Distribution by Treatment (Arithmetic Mean Method)",
subtitle = paste("n =", nrow(individual_cow_means_filtered), "cows with >=2 measurements"),
x = "Treatment",
y = "CH4 Emissions (g/d)",
fill = "Treatment"
) +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11)
)
12.2 Individual Cow Means Comparison
ggplot(individual_cow_means_filtered, aes(x = reorder(cow_id, Mean_CH4), y = Mean_CH4, fill = grazing_treatment)) +
geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = Mean_CH4 - Std_Dev, ymax = Mean_CH4 + Std_Dev),
width = 0.3, alpha = 0.7) +
facet_wrap(~grazing_treatment, scales = "free_x", ncol = 2) +
theme_minimal() +
labs(
title = "Individual Cow Mean CH4 ± SD (Arithmetic Mean Method, >=2 Measurements)",
x = "Cow ID",
y = "Mean CH4 (g/d)",
fill = "Treatment"
) +
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, size = 9),
text = element_text(size = 11),
plot.title = element_text(face = "bold", size = 14)
)
12.3 Within-Animal Variation (Coefficient of Variation)
ggplot(individual_cow_means_filtered, aes(x = grazing_treatment, y = CV_percent, fill = grazing_treatment)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.6, size = 3) +
theme_minimal() +
labs(
title = "Within-Animal Variation by Treatment (Coefficient of Variation %)",
x = "Treatment",
y = "Coefficient of Variation (%)",
fill = "Treatment"
) +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 14)
)
12.4 Treatment Means Comparison
treatment_means_plot <- data.frame(
Treatment = treatment_summary$grazing_treatment,
Mean = treatment_summary$Mean_CH4_CowLevel,
SE = treatment_summary$SEM_Cows
)
ggplot(treatment_means_plot, aes(x = Treatment, y = Mean, fill = Treatment)) +
geom_bar(stat = "identity", alpha = 0.7) +
geom_errorbar(aes(ymin = Mean - 1.96*SE, ymax = Mean + 1.96*SE),
width = 0.2, size = 1) +
theme_minimal() +
labs(
title = "Treatment Mean CH4 Emissions (Arithmetic Mean Method)",
subtitle = "Error bars = 95% Confidence Interval",
x = "Treatment",
y = "Mean CH4 (g/d)",
fill = "Treatment"
) +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10)
)
13 SUMMARY
Method: Arithmetic Mean (AM) - Direct mean of all valid measurements per Manafiazar et al. (2017)
Data Processing: - Removed poor-quality night visits (data quality issues) - Applied ≥2 measurements filter for statistical validity - Final sample: 34 cows total
Key Results: - K1_AMP: 18 cows, mean = 309.36 g/d - K2_CON: 16 cows, mean = 334.67 g/d - Difference: 25.31 g/d (8.18%) - Statistical Significance: p = 0.03443767 - Effect Size (Cohen’s d): 0.7556 (Medium)
Advantage vs TOD: Arithmetic Mean retains more cows because it doesn’t require 14-day overlapping windows, allowing inclusion of cows with shorter observation periods.