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

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

📋 基本情報

基本情報

牧場名: takaT

検定日: 2026年8月

検定頭数: 154 頭

---

📊 検定結果統計

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)	40.76	10.39	41.85	11.80	62.60
乳脂率(%)	4.18	0.70	4.04	2.90	7.05
蛋白質率(%)	3.40	0.35	3.37	2.63	4.41
乳糖率(%)	4.52	0.16	4.55	3.84	4.85
MUN(mg/dl)	11.11	1.84	11.10	7.60	16.00
体細胞数(千/ml)	115.60	239.78	37.50	4.00	2044.00
体細胞スコア	1.92	1.78	2.00	0.00	7.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.62	1.71	4.0	1	7
搾乳日数	180.06	125.25	168.5	10	815
乾乳日数	69.21	30.66	63.5	26	192
分娩間隔	407.99	65.34	402.0	318	675
初産分娩月齢	23.50	2.71	22.5	21	31
空胎日数	109.75	76.94	85.0	9	636
授精回数	1.65	1.07	1.0	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産	18	30.30	51.56	1.78	11.7
2産	28	38.69	132.57	1.57	18.2
3産以上	108	43.04	121.88	2.03	70.1
合計	154	40.76	115.60	1.92	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産	18	10057.74	12785.22	11.7
2産	28	12826.30	13587.92	18.2
3産以上	108	13298.74	13193.51	70.1
合計	154	12262.14	13220.58	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産	18	252.33	127.17	1.83	11.7
2産	28	151.32	80.32	1.36	18.2
3産以上	108	175.46	114.47	1.67	70.1
合計	154	180.06	109.75	1.65	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	7	92,562
2	1659	6	86,779
3	1551	7	84,985
4	1622	7	84,299
5	1588	7	82,121
6	1906	6	80,766
7	1892	6	78,786
8	1758	6	77,986
9	1955	5	76,506
10	1645	7	71,309
11	1945	6	70,582
12	1918	6	69,025
13	1652	7	68,470
14	1767	6	68,371
15	1790	7	68,311

---

📈 産次別脂肪酸組成

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産	18	26.01	39.26	11.7
2産	28	27.71	36.11	18.2
3産以上	108	28.11	36.42	70.1
合計	154	27.79	36.69	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以下の牛: 4頭**

乳量の高い牛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	2077	4	21	30	62.6
2	2086	4	126	4	61.7
3	1962	6	78	35	59.2
4	1900	5	82	20	58.9
5	2268	4	138	10	58.5
6	1622	7	27	23	57.4
7	2316	3	31	397	57.3
8	1659	6	550	479	56.8
9	2065	5	53	542	55.1
10	1989	5	87	5	54.5

---

🧪 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%未満の牛: 3頭**

プレフォーム(%)

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	2331	2	169	2044	37.8
2	1954	5	399	1307	19.9
3	1652	7	72	1194	21.6
4	1897	6	29	627	46.6
5	1971	6	104	548	43.6
6	2065	5	53	542	55.1
7	1987	5	135	515	43.0
8	2275	4	104	481	48.9
9	1659	6	550	479	56.8
10	2061	5	135	440	37.4
11	2316	3	31	397	57.3
12	2314	3	171	384	53.1
13	2303	3	212	353	40.3
14	2254	4	174	341	35.5
15	2650	2	240	338	11.8
16	1551	7	324	327	18.6
17	2051	5	190	286	28.3
18	1765	7	40	280	54.0
19	1978	5	329	266	23.8
20	1972	6	65	210	51.3

---

📊 体細胞スコア分析

体細胞スコア分布

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-08-30 18:41:17.040482

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