乳検データ分析レポート

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

📋 基本情報

基本情報

牧場名: 小森

検定日: 2026年2月

検定頭数: 144 頭

---

📊 検定結果統計

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)	37.85	9.80	37.10	11.20	64.10
乳脂率(%)	4.00	0.54	4.04	2.53	5.51
蛋白質率(%)	3.32	0.30	3.30	2.62	4.51
乳糖率(%)	4.56	0.20	4.58	3.77	4.93
MUN(mg/dl)	10.46	2.00	10.30	3.80	17.50
体細胞数(千/ml)	202.62	633.34	61.00	5.00	5989.00
体細胞スコア	2.44	1.80	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)

繁殖成績統計

項目	平均値	標準偏差	中央値	最小値	最大値
産次	2.33	1.20	2.0	1	6
搾乳日数	184.17	105.74	178.5	7	497
乾乳日数	72.92	34.04	62.0	28	210
分娩間隔	435.74	83.59	414.0	309	715
初産分娩月齢	23.02	1.57	23.0	21	28
空胎日数	119.48	72.05	96.5	6	348
授精回数	2.10	1.27	2.0	1	7

📝 注意事項:

• 乾乳日数：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産	43	32.75	200.77	1.88	29.9
2産	47	38.66	141.38	2.66	32.6
3産以上	54	41.22	257.41	2.70	37.5
合計	144	37.85	202.62	2.44	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産	43	9850.25	13203.43	29.9
2産	47	12542.70	13710.97	32.6
3産以上	54	12866.50	12493.38	37.5
合計	144	11666.65	13102.82	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産	43	189.81	117.28	1.87	29.9
2産	47	198.96	128.72	2.18	32.6
3産以上	54	166.80	113.19	2.21	37.5
合計	144	184.17	119.48	2.10	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	1802	6	66,056
2	1820	5	64,592
3	1977	5	57,839
4	1976	5	55,300
5	2221	4	53,212
6	2030	4	52,768
7	2001	4	52,099
8	2021	4	51,156
9	2027	4	50,816
10	2023	4	50,081
11	2102	4	49,359
12	2047	4	48,815
13	2039	4	48,077
14	2007	5	46,255
15	2003	4	45,941

---

📈 産次別脂肪酸組成

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産	43	25.25	41.47	29.9
2産	47	25.60	41.21	32.6
3産以上	54	25.79	40.92	37.5
合計	144	25.57	41.18	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以下の牛: 5頭**

乳量の高い牛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	2047	4	57	25	64.1
2	2442	3	72	11	61.2
3	419	4	76	12	60.8
4	2463	2	48	11	60.4
5	2026	4	72	9	57.8
6	2449	3	39	11	57.4
7	1977	5	34	15	57.0
8	2446	3	46	23	55.9
9	2185	2	44	350	55.2
10	2111	3	34	10	54.0

---

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

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

---

🧬 脂肪酸組成分析

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

プレフォーム(%)

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	2221	4	180	5989	39.1
2	2699	1	321	3656	18.6
3	1614	1	46	2877	15.2
4	2127	3	237	898	33.8
5	752	2	186	792	33.4
6	1820	5	245	748	37.9
7	2460	2	335	745	25.7
8	2139	3	8	722	38.7
9	2039	4	133	711	50.9
10	2474	2	256	490	21.2
11	2444	3	131	379	41.5
12	2136	2	268	371	35.9
13	2185	2	44	350	55.2
14	2116	3	412	321	16.7
15	2472	2	269	298	17.8
16	2461	2	342	282	24.2
17	8589	2	411	281	25.5
18	753	2	150	275	33.4
19	2041	4	204	258	30.1
20	1720	4	90	246	22.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-02-26 15:45:45.212854

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