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

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

📋 基本情報

基本情報

牧場名: 高橋牧場

検定日: 2025年2月

検定頭数: 147 頭

---

📊 検定結果統計

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)	39.33	9.99	39.90	18.20	59.40
乳脂率(%)	4.51	0.72	4.42	2.72	7.72
蛋白質率(%)	3.72	0.33	3.73	2.95	4.81
乳糖率(%)	4.58	0.18	4.61	4.01	4.97
MUN(mg/dl)	10.87	1.95	10.70	6.00	16.20
体細胞数(千/ml)	165.33	576.76	37.00	4.00	4607.00
体細胞スコア	2.00	1.79	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.30	1.71	3	1	8
搾乳日数	193.35	116.14	179	10	637
乾乳日数	70.26	35.14	64	26	216
分娩間隔	405.94	65.11	398	318	621
初産分娩月齢	23.48	2.46	23	21	31
空胎日数	103.78	73.80	81	9	585
授精回数	1.80	1.03	1	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産	29	31.89	40.72	1.21	19.7
2産	22	38.05	134.64	2.05	15.0
3産以上	96	41.87	210.00	2.23	65.3
合計	147	39.33	165.33	2.00	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産	29	9122.14	12544.84	19.7
2産	22	12440.90	13224.43	15.0
3産以上	96	13437.47	13474.61	65.3
合計	147	12454.82	13264.13	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産	29	215.55	112.38	1.65	19.7
2産	22	182.91	104.68	2.07	15.0
3産以上	96	189.04	100.98	1.80	65.3
合計	147	193.35	103.78	1.80	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	88,319
2	1433	8	87,984
3	1622	6	81,092
4	1659	6	79,755
5	1551	7	78,745
6	1599	7	77,615
7	1588	7	75,690
8	1611	7	73,767
9	1906	5	72,260
10	1892	6	71,167
11	1955	5	71,110
12	1758	6	70,995
13	1652	6	66,466
14	1663	6	65,192
15	1918	5	64,721

---

📈 産次別脂肪酸組成

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産	29	28.22	35.21	19.7
2産	22	28.72	34.55	15.0
3産以上	96	28.18	35.24	65.3
合計	147	28.27	35.13	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	1955	5	123	27	59.4
2	2062	4	164	60	58.6
3	2076	4	69	16	58.2
4	2287	3	61	26	57.3
5	2318	2	51	12	57.3
6	2085	4	38	371	57.1
7	2295	3	101	9	56.9
8	2054	5	50	128	56.2
9	2296	3	64	4	55.5
10	2643	1	412	16	54.7

---

🧪 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	2274	3	340	4607	18.5
2	1978	5	151	4606	41.3
3	2303	3	34	2347	51.9
4	2313	2	82	1494	36.3
5	1659	6	372	670	39.9
6	1652	6	312	455	32.8
7	2316	2	206	404	36.3
8	2286	3	204	399	26.9
9	2051	5	12	374	44.4
10	2085	4	38	371	57.1
11	1551	7	146	342	46.6
12	2074	4	204	306	36.7
13	1611	7	245	272	37.5
14	2077	3	260	266	54.1
15	1765	6	251	257	45.0
16	1973	4	235	253	43.9
17	2677	1	111	253	27.3
18	1960	5	299	230	21.9
19	1947	4	384	199	26.6
20	2649	2	49	187	47.9

---

📊 体細胞スコア分析

体細胞スコア分布

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 16:07:22.651529

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