高橋牧場乳検データ分析レポート

🐄 乳検データ分析レポート

📋 基本情報

基本情報

牧場名: 高橋牧場

検定日: 2026年3月

検定頭数: 151 頭

---

📊 検定結果統計

calc_stats <- function(x) {
  x <- as.numeric(x)
  x <- x[!is.na(x) & is.finite(x)]
  if (length(x) == 0) return(c(mean=NA, sd=NA, median=NA, min=NA, max=NA))
  c(mean   = round(mean(x),   2),
    sd     = round(sd(x),     2),
    median = round(median(x), 2),
    min    = round(min(x),    2),
    max    = round(max(x),    2))
}

items <- c("乳量(kg)", "乳脂率(%)", "蛋白質率(%)", "乳糖率(%)",
           "MUN(mg/dl)", "体細胞数(千/ml)", "体細胞スコア")
cols  <- c("乳量_K", "乳脂率", "蛋白質率", "乳糖率",
           "MUN", "体細胞数", "体細胞スコア")

stats_mat <- t(sapply(cols, function(col) calc_stats(data[[col]])))
stats_df <- data.frame(項目 = items,
                       平均値   = stats_mat[, "mean"],
                       標準偏差 = stats_mat[, "sd"],
                       中央値   = stats_mat[, "median"],
                       最小値   = stats_mat[, "min"],
                       最大値   = stats_mat[, "max"],
                       stringsAsFactors = FALSE)

kable(stats_df, caption = "検定結果統計", row.names = FALSE)

検定結果統計

項目	平均値	標準偏差	中央値	最小値	最大値
乳量(kg)	38.51	10.09	38.00	10.60	61.50
乳脂率(%)	4.36	0.67	4.30	2.46	6.67
蛋白質率(%)	3.67	0.35	3.68	2.63	4.65
乳糖率(%)	4.52	0.19	4.54	3.86	4.88
MUN(mg/dl)	12.29	1.89	12.10	6.60	18.20
体細胞数(千/ml)	223.01	784.25	52.00	6.00	6453.00
体細胞スコア	2.21	1.90	2.00	0.00	9.00

---

🐄 繁殖成績統計

dry_days    <- as.numeric(data$乾乳日数)
dry_days    <- dry_days[!is.na(dry_days) & dry_days != 0]
calving_int <- as.numeric(data$分娩間隔初産月齢[as.numeric(data$産次) != 1])
first_calv  <- as.numeric(data$分娩間隔初産月齢[as.numeric(data$産次) == 1])

repro_items <- c("産次", "搾乳日数", "乾乳日数", "分娩間隔", "初産分娩月齢", "空胎日数", "授精回数")
repro_vals  <- list(as.numeric(data$産次), as.numeric(data$搾乳日数),
                    dry_days, calving_int, first_calv,
                    as.numeric(data$空胎日数), as.numeric(data$授精回数))

repro_mat <- t(sapply(repro_vals, calc_stats))
repro_df <- data.frame(項目 = repro_items,
                       平均値   = repro_mat[, "mean"],
                       標準偏差 = repro_mat[, "sd"],
                       中央値   = repro_mat[, "median"],
                       最小値   = repro_mat[, "min"],
                       最大値   = repro_mat[, "max"],
                       stringsAsFactors = FALSE)

kable(repro_df, caption = "繁殖成績統計", row.names = FALSE)

繁殖成績統計

項目	平均値	標準偏差	中央値	最小値	最大値
産次	3.37	1.74	3	1	8
搾乳日数	190.63	121.49	193	8	667
乾乳日数	70.16	37.36	62	26	216
分娩間隔	403.79	64.27	395	318	621
初産分娩月齢	23.57	2.46	23	21	31
空胎日数	103.44	78.50	75	7	636
授精回数	1.78	1.10	1	1	6

📝 注意事項:

• 乾乳日数：0日のデータは除外して計算
• 分娩間隔：産次1のデータは除外して計算
• 初産分娩月齢：産次1のデータのみを使用して計算

---

📈 産次別乳量と体細胞

data_parity <- data %>%
  mutate(
    産次_num     = as.numeric(産次),
    産次グループ = case_when(
      産次_num == 1 ~ "1産",
      産次_num == 2 ~ "2産",
      産次_num >= 3 ~ "3産以上",
      TRUE          ~ NA_character_
    )
  ) %>%
  filter(!is.na(産次グループ))

parity_summary <- data_parity %>%
  group_by(産次グループ) %>%
  summarise(
    頭数             = n(),
    平均乳量_kg      = round(mean(as.numeric(乳量_K),       na.rm = TRUE), 2),
    平均体細胞数     = round(mean(as.numeric(体細胞数),     na.rm = TRUE), 2),
    平均体細胞スコア = round(mean(as.numeric(体細胞スコア), na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(割合 = round(頭数 / sum(頭数) * 100, 1))

total_row <- data.frame(
  産次グループ     = "合計",
  頭数             = sum(parity_summary$頭数),
  平均乳量_kg      = round(mean(as.numeric(data_parity$乳量_K),       na.rm = TRUE), 2),
  平均体細胞数     = round(mean(as.numeric(data_parity$体細胞数),     na.rm = TRUE), 2),
  平均体細胞スコア = round(mean(as.numeric(data_parity$体細胞スコア), na.rm = TRUE), 2),
  割合             = 100.0,
  stringsAsFactors = FALSE
)

parity_summary <- bind_rows(parity_summary, total_row)
parity_summary$産次グループ <- factor(parity_summary$産次グループ,
                                       levels = c("1産", "2産", "3産以上", "合計"))
parity_summary <- parity_summary[order(parity_summary$産次グループ), ]

kable(parity_summary, caption = "産次別乳量と体細胞", row.names = FALSE)

産次別乳量と体細胞

産次グループ	頭数	平均乳量_kg	平均体細胞数	平均体細胞スコア	割合
1産	28	30.59	52.79	1.50	18.5
2産	22	35.50	313.64	2.05	14.6
3産以上	101	41.37	250.47	2.44	66.9
合計	151	38.51	223.01	2.21	100.0

---

📊 産次別305日乳量と補正乳量

parity_305 <- data_parity %>%
  group_by(産次グループ) %>%
  summarise(
    頭数             = n(),
    平均305日乳量_kg = round(mean(as.numeric(乳量_305日), na.rm = TRUE), 2),
    平均補正乳量_kg  = round(mean(as.numeric(補正乳量),   na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(割合 = round(頭数 / sum(頭数) * 100, 1))

total_305 <- data.frame(
  産次グループ     = "合計",
  頭数             = sum(parity_305$頭数),
  平均305日乳量_kg = round(mean(as.numeric(data_parity$乳量_305日), na.rm = TRUE), 2),
  平均補正乳量_kg  = round(mean(as.numeric(data_parity$補正乳量),   na.rm = TRUE), 2),
  割合             = 100.0,
  stringsAsFactors = FALSE
)

parity_305 <- bind_rows(parity_305, total_305)
parity_305$産次グループ <- factor(parity_305$産次グループ,
                                   levels = c("1産", "2産", "3産以上", "合計"))
parity_305 <- parity_305[order(parity_305$産次グループ), ]

kable(parity_305, caption = "産次別305日乳量と補正乳量", row.names = FALSE)

産次別305日乳量と補正乳量

産次グループ	頭数	平均305日乳量_kg	平均補正乳量_kg	割合
1産	28	9477.48	12692.19	18.5
2産	22	12332.15	13241.48	14.6
3産以上	101	13765.74	13420.63	66.9
合計	151	12288.11	13268.90	100.0

---

🐄 産次別繁殖成績

parity_repro <- data_parity %>%
  group_by(産次グループ) %>%
  summarise(
    頭数         = n(),
    平均搾乳日数 = round(mean(as.numeric(搾乳日数), na.rm = TRUE), 2),
    平均空胎日数 = round(mean(as.numeric(空胎日数), na.rm = TRUE), 2),
    平均授精回数 = round(mean(as.numeric(授精回数), na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(割合 = round(頭数 / sum(頭数) * 100, 1))

total_repro <- data.frame(
  産次グループ = "合計",
  頭数         = sum(parity_repro$頭数),
  平均搾乳日数 = round(mean(as.numeric(data_parity$搾乳日数), na.rm = TRUE), 2),
  平均空胎日数 = round(mean(as.numeric(data_parity$空胎日数), na.rm = TRUE), 2),
  平均授精回数 = round(mean(as.numeric(data_parity$授精回数), na.rm = TRUE), 2),
  割合         = 100.0,
  stringsAsFactors = FALSE
)

parity_repro <- bind_rows(parity_repro, total_repro)
parity_repro$産次グループ <- factor(parity_repro$産次グループ,
                                    levels = c("1産", "2産", "3産以上", "合計"))
parity_repro <- parity_repro[order(parity_repro$産次グループ), ]

kable(parity_repro, caption = "産次別繁殖成績", row.names = FALSE)

産次別繁殖成績

産次グループ	頭数	平均搾乳日数	平均空胎日数	平均授精回数	割合
1産	28	242.14	122.36	1.76	18.5
2産	22	202.73	114.32	2.00	14.6
3産以上	101	173.71	95.82	1.73	66.9
合計	151	190.63	103.44	1.78	100.0

---

🏆 生涯乳量TOP15

top15_lifetime <- data %>%
  filter(!is.na(通算乳量) & !is.na(拡大4桁)) %>%
  mutate(通算乳量_num = as.numeric(通算乳量)) %>%
  arrange(desc(通算乳量_num)) %>%
  head(15) %>%
  mutate(
    順位        = 1:n(),
    通算乳量_kg = format(round(通算乳量_num, 0), big.mark = ",")
  ) %>%
  select(順位, 拡大4桁, 産次, 通算乳量_kg)

kable(top15_lifetime, caption = "生涯乳量TOP15", row.names = FALSE)

生涯乳量TOP15

順位	拡大4桁	産次	通算乳量_kg
1	1765	6	89,564
2	1433	8	89,115
3	1622	6	82,194
4	1659	6	81,059
5	1551	7	80,116
6	1599	7	78,061
7	1588	7	76,936
8	1611	7	74,735
9	1955	5	72,847
10	1892	6	72,815
11	1758	6	72,328
12	1645	7	68,366
13	1663	6	66,493
14	1761	7	66,328
15	1786	5	65,675

---

📈 産次別脂肪酸組成

parity_fa <- data_parity %>%
  group_by(産次グループ) %>%
  summarise(
    頭数             = n(),
    平均デノボFA     = round(mean(as.numeric(デノボFA),     na.rm = TRUE), 2),
    平均プレフォーム = round(mean(as.numeric(プレフォーム), na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(割合 = round(頭数 / sum(頭数) * 100, 1))

total_fa <- data.frame(
  産次グループ     = "合計",
  頭数             = sum(parity_fa$頭数),
  平均デノボFA     = round(mean(as.numeric(data_parity$デノボFA),     na.rm = TRUE), 2),
  平均プレフォーム = round(mean(as.numeric(data_parity$プレフォーム), na.rm = TRUE), 2),
  割合             = 100.0,
  stringsAsFactors = FALSE
)

parity_fa <- bind_rows(parity_fa, total_fa)
parity_fa$産次グループ <- factor(parity_fa$産次グループ,
                                 levels = c("1産", "2産", "3産以上", "合計"))
parity_fa <- parity_fa[order(parity_fa$産次グループ), ]

kable(parity_fa, caption = "産次別脂肪酸組成", row.names = FALSE)

産次別脂肪酸組成

産次グループ	頭数	平均デノボFA	平均プレフォーム	割合
1産	28	28.09	35.56	18.5
2産	22	28.53	35.46	14.6
3産以上	101	27.89	36.13	66.9
合計	151	28.02	35.93	100.0

---

📊 乳量分析

分娩後日数と乳量の関係

milk_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(乳量_K)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    乳量_num     = as.numeric(乳量_K),
    グループ     = ifelse(as.numeric(乳量_K) <= 20, "20kg以下", "21kg以上")
  )

plot_ly(milk_data,
        x = ~搾乳日数_num, y = ~乳量_num,
        color = ~グループ,
        colors = c("20kg以下" = "#F5576C", "21kg以上" = "#667EEA"),
        text = ~paste("拡大4桁:", 拡大4桁, "<br>乳量(kg):", round(乳量_num, 2)),
        hoverinfo = "text", type = "scatter", mode = "markers",
        marker = list(size = 8)) %>%
  layout(title = "分娩後日数と乳量の関係",
         xaxis = list(title = "分娩後日数"),
         yaxis = list(title = "乳量(kg)"),
         hovermode = "closest")

乳量20kg以下の牛のリスト

low_milk <- data %>%
  filter(!is.na(乳量_K) & as.numeric(乳量_K) <= 20) %>%
  mutate(搾乳日数_num = as.numeric(搾乳日数)) %>%
  arrange(搾乳日数_num) %>%
  select(拡大4桁, 産次, 搾乳日数, 乳量_K, 体細胞数) %>%
  head(50)

if (nrow(low_milk) > 0) {
  cat(paste0("⚠️ **乳量20kg以下の牛: ", nrow(low_milk), "頭**\n\n"))
  datatable(low_milk, options = list(pageLength = 10), rownames = FALSE)
} else {
  cat("✅ 乳量20kg以下の牛はいません。\n")
}

## ⚠️ **乳量20kg以下の牛: 6頭**

乳量の高い牛TOP10

top10_milk <- data %>%
  filter(!is.na(乳量_K)) %>%
  mutate(乳量_num = as.numeric(乳量_K)) %>%
  arrange(desc(乳量_num)) %>%
  head(10) %>%
  mutate(順位 = 1:n()) %>%
  select(順位, 拡大4桁, 産次, 搾乳日数, 体細胞数, 乳量_K)

kable(top10_milk, caption = "乳量の高い牛TOP10", row.names = FALSE)

乳量の高い牛TOP10

順位	拡大4桁	産次	搾乳日数	体細胞数	乳量_K
1	2259	4	45	7	61.5
2	1956	6	29	11	59.2
3	1892	6	41	38	59.0
4	1955	5	153	30	56.5
5	2062	4	194	87	56.5
6	2085	4	68	450	56.0
7	2318	2	81	9	55.6
8	1961	5	92	171	55.5
9	2076	4	99	12	55.2
10	2287	3	91	65	54.4

---

🧪 BHBA（ケトン体）分析

分娩後60日以内のBHBA分析

bhba_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(BHB)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    BHB_num      = as.numeric(BHB)
  ) %>%
  filter(搾乳日数_num <= 60) %>%
  mutate(グループ = ifelse(BHB_num >= 0.13, "0.13以上", "問題なし"))

if (nrow(bhba_data) > 0) {
  plot_ly(bhba_data,
          x = ~搾乳日数_num, y = ~BHB_num,
          color = ~グループ,
          colors = c("0.13以上" = "#F5576C", "問題なし" = "#667EEA"),
          text = ~paste("拡大4桁:", 拡大4桁, "<br>BHB(mmol/L):", round(BHB_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8)) %>%
    layout(title = "分娩後60日以内のBHB値",
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "BHB(mmol/L)"),
           hovermode = "closest")
} else {
  cat("データがありません。\n")
}

BHBA高値の牛のリスト（プロピレングリコール投与推奨）

high_bhba <- data %>%
  filter(!is.na(搾乳日数) & !is.na(BHB)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    BHB_num      = as.numeric(BHB)
  ) %>%
  filter(搾乳日数_num <= 60 & BHB_num >= 0.13) %>%
  arrange(desc(BHB_num)) %>%
  select(拡大4桁, 産次, 搾乳日数, 乳量_K, BHB)

if (nrow(high_bhba) > 0) {
  cat(paste0("⚠️ **プロピレングリコール投与推奨: ", nrow(high_bhba), "頭**\n\n"))
  datatable(high_bhba, options = list(pageLength = 10), rownames = FALSE)
} else {
  cat("✅ BHBA高値の牛はいません。\n")
}

## ✅ BHBA高値の牛はいません。

---

🧬 脂肪酸組成分析

デノボMilk(%)

denovo_milk_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(デノボMilk)) %>%
  mutate(
    搾乳日数_num   = as.numeric(搾乳日数),
    デノボMilk_num = as.numeric(デノボMilk)
  )

if (nrow(denovo_milk_data) > 0) {
  all_mean   <- mean(denovo_milk_data$デノボMilk_num, na.rm = TRUE)
  all_median <- median(denovo_milk_data$デノボMilk_num, na.rm = TRUE)
  d60        <- filter(denovo_milk_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$デノボMilk_num, na.rm = TRUE)
  median_60  <- median(d60$デノボMilk_num, na.rm = TRUE)

  plot_ly(denovo_milk_data,
          x = ~搾乳日数_num, y = ~デノボMilk_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>デノボMilk(%):", round(デノボMilk_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = paste0("分娩後日数とデノボMilk(%)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "デノボMilk(%)"),
           hovermode = "closest")
}

デノボFA(%)

denovo_fa_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(デノボFA)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    デノボFA_num  = as.numeric(デノボFA),
    グループ      = ifelse(as.numeric(デノボFA) < 22, "22%未満", "22%以上")
  )

if (nrow(denovo_fa_data) > 0) {
  all_mean   <- mean(denovo_fa_data$デノボFA_num, na.rm = TRUE)
  all_median <- median(denovo_fa_data$デノボFA_num, na.rm = TRUE)
  d60        <- filter(denovo_fa_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$デノボFA_num, na.rm = TRUE)
  median_60  <- median(d60$デノボFA_num, na.rm = TRUE)

  plot_ly(denovo_fa_data,
          x = ~搾乳日数_num, y = ~デノボFA_num,
          color = ~グループ,
          colors = c("22%未満" = "#F5576C", "22%以上" = "#667EEA"),
          text = ~paste("拡大4桁:", 拡大4桁, "<br>デノボFA(%):", round(デノボFA_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8)) %>%
    layout(title = paste0("分娩後日数とデノボFA(%)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "デノボFA(%)"),
           hovermode = "closest")
}

デノボFA 22%未満の牛のリスト（分娩後100日以内）

low_denovo <- data %>%
  filter(!is.na(搾乳日数) & !is.na(デノボFA)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    デノボFA_num  = as.numeric(デノボFA)
  ) %>%
  filter(搾乳日数_num <= 100 & デノボFA_num < 22) %>%
  arrange(デノボFA_num) %>%
  select(拡大4桁, 産次, 搾乳日数, 乳量_K, 体細胞数, デノボFA)

if (nrow(low_denovo) > 0) {
  cat(paste0("⚠️ **分娩後100日以内でデノボFA 22%未満の牛: ", nrow(low_denovo), "頭**\n\n"))
  datatable(low_denovo, options = list(pageLength = 10), rownames = FALSE)
} else {
  cat("✅ デノボFA 22%未満の牛はいません。\n")
}

## ⚠️ **分娩後100日以内でデノボFA 22%未満の牛: 2頭**

プレフォーム(%)

preform_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(プレフォーム)) %>%
  mutate(
    搾乳日数_num    = as.numeric(搾乳日数),
    プレフォーム_num = as.numeric(プレフォーム)
  )

if (nrow(preform_data) > 0) {
  all_mean   <- mean(preform_data$プレフォーム_num, na.rm = TRUE)
  all_median <- median(preform_data$プレフォーム_num, na.rm = TRUE)
  d60        <- filter(preform_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$プレフォーム_num, na.rm = TRUE)
  median_60  <- median(d60$プレフォーム_num, na.rm = TRUE)

  plot_ly(preform_data,
          x = ~搾乳日数_num, y = ~プレフォーム_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>プレフォーム(%):", round(プレフォーム_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = paste0("分娩後日数とプレフォーム(%)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "プレフォーム(%)"),
           hovermode = "closest")
}

---

🧪 乳成分分析

乳脂率（％）

fat_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(乳脂率)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    乳脂率_num   = as.numeric(乳脂率)
  )

if (nrow(fat_data) > 0) {
  all_mean   <- mean(fat_data$乳脂率_num, na.rm = TRUE)
  all_median <- median(fat_data$乳脂率_num, na.rm = TRUE)
  d60        <- filter(fat_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$乳脂率_num, na.rm = TRUE)
  median_60  <- median(d60$乳脂率_num, na.rm = TRUE)

  plot_ly(fat_data,
          x = ~搾乳日数_num, y = ~乳脂率_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>乳脂率(%):", round(乳脂率_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = paste0("分娩後日数と乳脂率(%)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "乳脂率(%)"),
           hovermode = "closest")
}

蛋白質率（％）

protein_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(蛋白質率)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    蛋白質率_num  = as.numeric(蛋白質率)
  )

if (nrow(protein_data) > 0) {
  all_mean   <- mean(protein_data$蛋白質率_num, na.rm = TRUE)
  all_median <- median(protein_data$蛋白質率_num, na.rm = TRUE)
  d60        <- filter(protein_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$蛋白質率_num, na.rm = TRUE)
  median_60  <- median(d60$蛋白質率_num, na.rm = TRUE)

  plot_ly(protein_data,
          x = ~搾乳日数_num, y = ~蛋白質率_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>蛋白質率(%):", round(蛋白質率_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = paste0("分娩後日数と蛋白質率(%)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "蛋白質率(%)"),
           hovermode = "closest")
}

乳糖率（％）

lactose_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(乳糖率)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    乳糖率_num   = as.numeric(乳糖率)
  )

if (nrow(lactose_data) > 0) {
  all_mean   <- mean(lactose_data$乳糖率_num, na.rm = TRUE)
  all_median <- median(lactose_data$乳糖率_num, na.rm = TRUE)
  d60        <- filter(lactose_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$乳糖率_num, na.rm = TRUE)
  median_60  <- median(d60$乳糖率_num, na.rm = TRUE)

  plot_ly(lactose_data,
          x = ~搾乳日数_num, y = ~乳糖率_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>乳糖率(%):", round(乳糖率_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = paste0("分娩後日数と乳糖率(%)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "乳糖率(%)"),
           hovermode = "closest")
}

MUN（mg/dl）

mun_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(MUN)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    MUN_num      = as.numeric(MUN)
  )

if (nrow(mun_data) > 0) {
  all_mean   <- mean(mun_data$MUN_num, na.rm = TRUE)
  all_median <- median(mun_data$MUN_num, na.rm = TRUE)
  d60        <- filter(mun_data, 搾乳日数_num <= 60)
  mean_60    <- mean(d60$MUN_num, na.rm = TRUE)
  median_60  <- median(d60$MUN_num, na.rm = TRUE)

  plot_ly(mun_data,
          x = ~搾乳日数_num, y = ~MUN_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>MUN(mg/dl):", round(MUN_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = paste0("分娩後日数とMUN(mg/dl)<br>全体: 平均 ", round(all_mean, 2),
                          " 中央値 ", round(all_median, 2),
                          " ／ 60日以内: 平均 ", round(mean_60, 2),
                          " 中央値 ", round(median_60, 2)),
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "MUN(mg/dl)"),
           hovermode = "closest")
}

---

🔬 体細胞数分析

分娩後日数と体細胞数

scc_data <- data %>%
  filter(!is.na(搾乳日数) & !is.na(体細胞数)) %>%
  mutate(
    搾乳日数_num = as.numeric(搾乳日数),
    体細胞数_num  = as.numeric(体細胞数)
  )

if (nrow(scc_data) > 0) {
  plot_ly(scc_data,
          x = ~搾乳日数_num, y = ~体細胞数_num,
          text = ~paste("拡大4桁:", 拡大4桁, "<br>体細胞数(千/ml):", round(体細胞数_num, 2)),
          hoverinfo = "text", type = "scatter", mode = "markers",
          marker = list(size = 8, color = "#667EEA")) %>%
    layout(title = "分娩後日数と体細胞数(千/ml)",
           xaxis = list(title = "分娩後日数"),
           yaxis = list(title = "体細胞数(千/ml)"),
           hovermode = "closest")
}

高体細胞数の牛TOP20

top20_scc <- data %>%
  filter(!is.na(体細胞数)) %>%
  mutate(体細胞数_num = as.numeric(体細胞数)) %>%
  arrange(desc(体細胞数_num)) %>%
  head(20) %>%
  mutate(順位 = 1:n()) %>%
  select(順位, 拡大4桁, 産次, 搾乳日数, 体細胞数, 乳量_K)

kable(top20_scc, caption = "高体細胞数の牛TOP20", row.names = FALSE)

高体細胞数の牛TOP20

順位	拡大4桁	産次	搾乳日数	体細胞数	乳量_K
1	2261	4	8	6453	30.4
2	2313	2	112	5651	28.6
3	1978	5	181	3272	31.5
4	2081	4	170	2542	31.3
5	1551	7	176	2151	44.8
6	2074	4	234	1313	35.0
7	2050	4	115	565	50.3
8	1934	6	35	491	45.6
9	1659	6	402	484	46.8
10	1954	5	251	460	28.0
11	2085	4	68	450	56.0
12	1599	7	365	440	10.6
13	2286	3	234	379	23.8
14	2077	3	290	367	52.9
15	1947	4	414	249	25.8
16	2326	2	122	235	33.9
17	2644	1	359	233	30.0
18	1622	6	310	219	31.9
19	1972	5	249	216	44.0
20	1939	5	275	213	33.5

---

📊 体細胞スコア分析

体細胞スコア分布

score_data <- data %>%
  filter(!is.na(体細胞スコア)) %>%
  mutate(体細胞スコア_round = round(as.numeric(体細胞スコア)))

if (nrow(score_data) > 0) {
  score_counts <- score_data %>%
    count(体細胞スコア_round) %>%
    arrange(体細胞スコア_round)

  plot_ly(score_counts,
          x = ~体細胞スコア_round, y = ~n,
          type = "bar",
          marker = list(color = "#667EEA")) %>%
    layout(title = "体細胞スコア分布",
           xaxis = list(title = "体細胞スコア"),
           yaxis = list(title = "頭数"))
}

産次ごとの体細胞スコア割合（100%積み上げ）

score_parity_data <- data %>%
  filter(!is.na(体細胞スコア) & !is.na(産次)) %>%
  mutate(
    体細胞スコア_round = round(as.numeric(体細胞スコア)),
    産次_group = case_when(
      as.numeric(産次) == 1 ~ "1",
      as.numeric(産次) == 2 ~ "2",
      as.numeric(産次) == 3 ~ "3",
      as.numeric(産次) >= 4 ~ "4以上",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(産次_group))

if (nrow(score_parity_data) > 0) {
  score_parity_counts <- score_parity_data %>%
    count(産次_group, 体細胞スコア_round) %>%
    group_by(産次_group) %>%
    mutate(percentage = n / sum(n) * 100) %>%
    ungroup()

  score_parity_counts$産次_group <- factor(score_parity_counts$産次_group,
                                           levels = c("1", "2", "3", "4以上"))

  plot_ly(score_parity_counts,
          x = ~産次_group, y = ~percentage,
          color = ~as.factor(体細胞スコア_round),
          type = "bar",
          text = ~paste0(round(percentage, 1), "%"),
          textposition = "inside") %>%
    layout(title = "産次ごとの体細胞スコア割合（100%積み上げ）",
           xaxis = list(title = "産次"),
           yaxis = list(title = "割合 (%)", range = c(0, 100)),
           barmode = "stack",
           legend = list(title = list(text = "スコア")))
}

---

📝 レポート作成情報

作成日時: 2026-03-27 17:05:41.839355

Rバージョン: R version 4.3.2 (2023-10-31 ucrt)

使用パッケージ: readxl, dplyr, ggplot2, plotly, knitr, DT, tidyr

---

🎯 まとめ

このレポートは以下の分析を含んでいます：

1. ✅ 基本情報 - 牧場名、検定日、検定頭数
2. ✅ 検定結果統計 - 乳量、乳成分、体細胞数の統計
3. ✅ 繁殖成績統計 - 産次、搾乳日数、空胎日数等の統計
4. ✅ 産次別分析 - 乳量、体細胞、305日乳量、繁殖成績、脂肪酸組成
5. ✅ 生涯乳量TOP15 - 通算乳量の上位15頭
6. ✅ 乳量分析 - 散布図、低乳量牛リスト、TOP10
7. ✅ BHBA分析 - ケトン体分析、プロピレングリコール投与推奨リスト
8. ✅ 脂肪酸組成分析 - デノボMilk、デノボFA、プレフォーム
9. ✅ 乳成分分析 - 乳脂率、蛋白質率、乳糖率、MUN
10. ✅ 体細胞分析 - 体細胞数散布図、TOP20、スコア分布

全てのグラフはインタラクティブで、マウスオーバーで詳細情報が表示されます。

---

使用方法:

1. このファイルをRStudioで開く
2. excel_file_path <- "小森20260219.xlsx" を実際のファイルのフルパスに変更
3. Knit ボタンをクリック → HTMLレポートが生成されます

必要なパッケージ: readxl, dplyr, ggplot2, plotly, knitr, DT, tidyr