小森牧場乳検データ分析レポート

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

📋 基本情報

基本情報

牧場名: haruT

検定日: 2026年6月

検定頭数: 140 頭

---

📊 検定結果統計

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)	35.16	9.65	35.60	10.10	60.50
乳脂率(%)	4.22	0.55	4.23	2.84	6.35
蛋白質率(%)	3.10	0.30	3.07	2.55	4.20
乳糖率(%)	4.60	0.21	4.64	3.55	4.96
MUN(mg/dl)	9.19	1.87	9.05	4.10	15.10
体細胞数(千/ml)	220.41	672.79	70.50	7.00	5727.00
体細胞スコア	2.54	1.82	2.50	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)

繁殖成績統計

項目	平均値	標準偏差	中央値	最小値	最大値
産次	2.41	1.25	2.0	1	6
搾乳日数	200.68	120.00	190.5	8	555
乾乳日数	76.71	35.79	65.0	28	210
分娩間隔	439.34	80.59	431.0	308	715
初産分娩月齢	23.31	2.41	22.0	21	32
空胎日数	116.04	74.08	93.0	7	462
授精回数	1.91	1.08	2.0	1	5

📝 注意事項:

• 乾乳日数：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産	39	32.68	73.59	1.82	27.9
2産	44	32.38	221.50	3.05	31.4
3産以上	57	39.02	320.04	2.65	40.7
合計	140	35.16	220.41	2.54	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産	39	9590.25	13307.12	27.9
2産	44	11944.59	13660.40	31.4
3産以上	57	11425.96	12287.96	40.7
合計	140	11422.92	13003.21	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産	39	195.00	115.31	1.76	27.9
2産	44	234.95	127.75	2.06	31.4
3産以上	57	178.11	107.51	1.90	40.7
合計	140	200.68	116.04	1.91	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	1887	5	72,267
2	1802	6	71,318
3	1820	5	69,879
4	1977	5	64,631
5	2029	5	57,555
6	2221	4	57,230
7	2001	4	56,148
8	2021	4	55,141
9	2047	4	54,744
10	2023	5	54,282
11	2039	4	53,639
12	2102	4	53,610
13	2007	5	50,904
14	2040	4	50,531
15	1975	5	50,361

---

📈 産次別脂肪酸組成

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産	39	22.42	43.94	27.9
2産	44	22.36	43.80	31.4
3産以上	57	22.34	44.70	40.7
合計	140	22.37	44.20	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以下の牛: 12頭**

乳量の高い牛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	2107	4	51	17	60.5
2	1887	5	58	90	56.9
3	2046	5	79	7	53.7
4	2111	3	161	30	53.1
5	2472	3	34	20	52.5
6	2453	3	60	13	51.8
7	2463	2	175	27	51.0
8	1977	5	161	12	50.4
9	2112	4	128	16	50.2
10	2029	5	114	28	49.8

---

🧪 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")
}

## ⚠️ **プロピレングリコール投与推奨: 2頭**

---

🧬 脂肪酸組成分析

デノボ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%未満の牛: 14頭**

プレフォーム(%)

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	2126	4	39	5727	39.1
2	2221	4	307	4623	30.1
3	2685	2	118	2445	35.7
4	752	2	313	1916	10.1
5	2461	3	18	1191	35.7
6	2943	1	217	980	16.1
7	2442	3	199	778	46.1
8	2039	4	260	683	39.7
9	2447	2	435	522	19.8
10	2464	2	423	478	28.2
11	2449	3	166	472	38.1
12	2127	3	364	448	20.7
13	2136	2	395	337	15.8
14	2981	2	321	288	25.6
15	753	2	277	285	22.6
16	2034	4	302	263	31.7
17	2444	3	258	244	25.2
18	2110	4	156	239	46.5
19	2263	2	432	239	19.8
20	2133	3	133	234	40.7

---

📊 体細胞スコア分析

体細胞スコア分布

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-06-26 20:17:30.672897

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