你好！进行玉米的 **NCII（North Carolina Design II，北卡罗来纳II设计...

下面这段 R 脚本可以直接吃你现在这张“宽表”（ph_1 / ph_2 …），自动转成长表，然后对每个性状（ph / ear / mu）分别做：

• NCII 方差分析（ANOVA：Rep + mom + dad + mom:dad）
• 随机效应模型（lmer）估计方差分量：母本GCA、父本GCA、SCA、残差
• 输出 GCA（mom、dad） 和 SCA（mom×dad） 的效应值（用 BLUP：ranef）
• 计算一套常用的 σA（加性）/ σD（显性）/ h²（按组合均值的近似狭义遗传力）
• 把结果写到一个 NCII_results.xlsx（每个性状一组sheet）

你不用改列名；读文件用 file.choose()，会弹窗让你选 Excel。

r
## =========================
## NCII (North Carolina Design II) analysis for your wide-format Excel
## Columns: hybrid, mom, dad, ph_1, ear_1, mu_1, ph_2, ear_2, mu_2
## =========================

## 0) Packages
pkgs <- c("readxl", "dplyr", "tidyr", "stringr", "lme4", "broom", "writexl")
to_install <- pkgs[!pkgs %in% rownames(installed.packages())]
if (length(to_install) > 0) install.packages(to_install, dependencies = TRUE)
invisible(lapply(pkgs, require, character.only = TRUE))

## 1) Read Excel (interactive file picker)
dat0 <- readxl::read_excel(path = file.choose())

## 2) Basic checks
required_cols <- c("hybrid", "mom", "dad")
missing_required <- setdiff(required_cols, names(dat0))
if (length(missing_required) > 0) {
  stop("Missing required columns: ", paste(missing_required, collapse = ", "))
}

## 3) Pivot from wide -> long
##    Any column ending with _<number> will be treated as trait_rep
##    e.g., ph_1, ear_2, mu_2 ...
dat_long <- dat0 %>%
  dplyr::mutate(
    hybrid = as.character(hybrid),
    mom    = as.factor(mom),
    dad    = as.factor(dad)
  ) %>%
  tidyr::pivot_longer(
    cols = dplyr::matches("_.+\\d+$|_\\d+$"),  # tolerant
    names_to = c("trait", "Rep"),
    names_pattern = "^(.+)_([0-9]+)$",
    values_to = "y"
  ) %>%
  dplyr::mutate(
    trait = as.factor(trait),
    Rep   = as.factor(Rep),
    y     = suppressWarnings(as.numeric(y))
  ) %>%
  dplyr::filter(!is.na(y))

if (nrow(dat_long) == 0) stop("No numeric phenotype values found after pivot_longer().")

## Optional: Use sum-to-zero contrasts to make 'effects' more interpretable
options(contrasts = c("contr.sum", "contr.poly"))

## 4) Helper functions
fit_one_trait <- function(df_trait) {
  # Fixed-effects ANOVA (classic factorial)
  aov_fit <- aov(y ~ Rep + mom + dad + mom:dad, data = df_trait)
  aov_tab <- broom::tidy(aov_fit)

  # Random-effects model for variance components
  # Rep kept as fixed (common in field trials with reps)
  lmer_fit <- lme4::lmer(y ~ Rep + (1|mom) + (1|dad) + (1|mom:dad), data = df_trait, REML = TRUE)

  vc <- as.data.frame(lme4::VarCorr(lmer_fit))
  # vc columns: grp, var1, var2, vcov, sdcor

  get_v <- function(grp_name) {
    out <- vc$vcov[vc$grp == grp_name]
    if (length(out) == 0) return(0)
    out[1]
  }

  v_mom    <- get_v("mom")
  v_dad    <- get_v("dad")
  v_momdad <- get_v("mom:dad")
  v_resid  <- get_v("Residual")

  r <- dplyr::n_distinct(df_trait$Rep)

  # Common NCII mapping (diploid, no epistasis assumptions):
  # Var(mom GCA) ~ 1/4 σA ; Var(dad GCA) ~ 1/4 σA  => σA ≈ 2*(v_mom + v_dad)
  # Var(SCA)     ~ 1/4 σD                       => σD ≈ 4*v_momdad
  sigma_A <- 2 * (v_mom + v_dad)
  sigma_D <- 4 * v_momdad

  # Narrow-sense heritability on hybrid MEAN (approx):
  # h2 ≈ σA / (σA + σD + σe/r)
  h2_mean <- ifelse((sigma_A + sigma_D + v_resid / r) > 0,
                    sigma_A / (sigma_A + sigma_D + v_resid / r),
                    NA_real_)

  varcomp <- dplyr::tibble(
    component = c("V_mom(GCA)", "V_dad(GCA)", "V_mom:dad(SCA)", "V_residual",
                  "sigma_A(≈2*(V_mom+V_dad))", "sigma_D(≈4*V_SCA)", "h2_mean(≈sigmaA/(sigmaA+sigmaD+Ve/r))",
                  "nRep"),
    value = c(v_mom, v_dad, v_momdad, v_resid, sigma_A, sigma_D, h2_mean, r)
  )

  # BLUPs for GCA/SCA (random effects)
  re <- lme4::ranef(lmer_fit)

  gca_mom <- re$mom %>%
    tibble::rownames_to_column("mom") %>%
    dplyr::rename(GCA_mom = `(Intercept)`) %>%
    dplyr::arrange(dplyr::desc(GCA_mom))

  gca_dad <- re$dad %>%
    tibble::rownames_to_column("dad") %>%
    dplyr::rename(GCA_dad = `(Intercept)`) %>%
    dplyr::arrange(dplyr::desc(GCA_dad))

  sca <- re[["mom:dad"]] %>%
    tibble::rownames_to_column("mom_dad") %>%
    dplyr::rename(SCA = `(Intercept)`) %>%
    tidyr::separate(mom_dad, into = c("mom", "dad"), sep = ":", remove = TRUE) %>%
    dplyr::arrange(dplyr::desc(SCA))

  # Hybrid means (useful for ranking)
  means <- df_trait %>%
    dplyr::group_by(hybrid, mom, dad) %>%
    dplyr::summarise(mean_y = mean(y, na.rm = TRUE), .groups = "drop") %>%
    dplyr::arrange(dplyr::desc(mean_y))

  list(
    aov_table = aov_tab,
    varcomp   = varcomp,
    gca_mom   = gca_mom,
    gca_dad   = gca_dad,
    sca       = sca,
    means     = means
  )
}

## 5) Run by trait and collect outputs
traits <- levels(dat_long$trait)
if (is.null(traits) || length(traits) == 0) traits <- unique(as.character(dat_long$trait))

out_sheets <- list()

for (tr in traits) {
  df_tr <- dat_long %>% dplyr::filter(trait == tr)

  # skip if too small
  if (nrow(df_tr) < 5) next

  res <- fit_one_trait(df_tr)

  # Put into separate sheets (Excel sheet name limit ~31 chars)
  tr_safe <- stringr::str_replace_all(as.character(tr), "[\$$\$$\\*\\?/\\\\:]", "_")
  tr_safe <- substr(tr_safe, 1, 20)

  out_sheets[[paste0(tr_safe, "_ANOVA")]]    <- res$aov_table
  out_sheets[[paste0(tr_safe, "_VarComp")]]  <- res$varcomp
  out_sheets[[paste0(tr_safe, "_GCA_mom")]]  <- res$gca_mom
  out_sheets[[paste0(tr_safe, "_GCA_dad")]]  <- res$gca_dad
  out_sheets[[paste0(tr_safe, "_SCA")]]      <- res$sca
  out_sheets[[paste0(tr_safe, "_Means")]]    <- res$means
}

if (length(out_sheets) == 0) stop("No trait results generated. Check your trait columns like ph_1, ph_2 ...")

## 6) Write results to Excel
writexl::write_xlsx(out_sheets, path = "NCII_results.xlsx")

message("Done! Results saved to: NCII_results.xlsx")

你跑完以后，重点看：

• *_ANOVA：mom / dad / mom:dad 的显著性（对应 GCA / GCA / SCA）
• *_VarComp：sigma_A、sigma_D、h2_mean
• *_GCA_mom、*_GCA_dad：每个亲本的 GCA（正值通常表示对性状贡献高，取决于性状方向）
• *_SCA：具体组合的 SCA（杂种优势/特配）

如果你希望我帮你把 “ph/ear/mu”重命名成更直观的性状名，或者你有 多地点 Env（需要加 Env 和交互项），把你的表头/示例两三行贴出来，我可以把模型扩展成多环境的 NCII（含 GCA×Env、SCA×Env）。