linetester <- function(replications, lines, testers, y) {

  replications <- as.factor(replications)
  lines        <- as.factor(lines)
  testers      <- as.factor(testers)

  datos <- data.frame(
    Replications = replications,
    Lines        = lines,
    Testers      = testers,
    Y            = y
  )

  n_r <- nlevels(datos$Replications)
  n_l <- nlevels(datos$Lines)
  n_t <- nlevels(datos$Testers)

  # ---- crosses (complete rows) and parents (line or tester absent) ----
  datos2 <- na.omit(datos)                                        # crosses
  datos3 <- subset(datos, is.na(datos$Lines) | is.na(datos$Testers))  # parents

  # ---- ANOVA with parents: Replications, Treatments, Error ----
  Treatments <- as.factor(paste(datos$Lines, datos$Testers))
  modelo1 <- aov(Y ~ Replications + Treatments, data = datos)
  matriz1 <- as.matrix(anova(modelo1))

  # ---- Lines, Testers, Lines x Testers: fitted on the CROSSES only ----
  modelo4 <- aov(Y ~ Lines * Testers, data = datos2)
  matriz4 <- as.matrix(anova(modelo4))

  # ---- Parents ----
  Parents <- as.factor(paste(datos3$Lines, datos3$Testers))
  modelo2 <- aov(Y ~ Parents, data = datos3)
  matriz2 <- as.matrix(anova(modelo2))

  # ---- Crosses ----
  Crosses <- as.factor(paste(datos2$Lines, datos2$Testers))
  modelo3 <- aov(Y ~ Crosses, data = datos2)
  matriz3 <- as.matrix(anova(modelo3))

  # ---- cell means of the crosses ----
  mm_cross <- tapply(datos2$Y, datos2[, c("Lines", "Testers")], mean, na.rm = TRUE)
  gm_cell  <- mean(mm_cross, na.rm = TRUE)

  # ---- SCA: cell - line marginal - tester marginal + grand ----
  SCA <- round(
    mm_cross -
      outer(rowMeans(mm_cross, na.rm = TRUE),
            colMeans(mm_cross, na.rm = TRUE), "+") +
      gm_cell,
    3
  )

  # ---- GCA of lines and testers, from the cell means ----
  GCA.lines   <- round(rowMeans(mm_cross, na.rm = TRUE) - gm_cell, 3)
  GCA.testers <- round(colMeans(mm_cross, na.rm = TRUE) - gm_cell, 3)

  # ---- assemble the ten-line table ----
  # Parents vs. Crosses by subtraction: Treatments - Parents - Crosses.
  # Only df and SS are meaningful here; columns 3-5 are recomputed below.
  matriz5 <- matriz1[2, ] - matriz2[1, ] - matriz3[1, ]

  matriz <- rbind(
    matriz1[1:2, ],   # Replications, Treatments
    matriz2[1, ],     # Parents
    matriz5,          # Parents vs. Crosses
    matriz3[1, ],     # Crosses
    matriz4[1:3, ],   # Lines, Testers, Lines x Testers
    matriz1[3, ]      # Error
  )
  matriz <- rbind(matriz, c(sum(matriz1[, 1]), sum(matriz1[, 2]), NA, NA, NA))

  rownames(matriz) <- c(
    "Replications", "Treatments", "Parents", "Parents vs. Crosses",
    "Crosses", "Lines", "Testers", "Lines x Testers", "Error", "Total"
  )

  # ---- F-tests: Lines (6) and Testers (7) against Lines x Testers (8);
  #      every other source against Error (9) ----
  for (i in 1:9) {
    matriz[i, 3] <- matriz[i, 2] / matriz[i, 1]
    matriz[i, 4] <- matriz[i, 3] / matriz[9, 3]
    matriz[i, 5] <- 1 - pf(matriz[i, 4], matriz[i, 1], matriz[9, 1])

    if (i == 6 || i == 7) {
      matriz[i, 4] <- matriz[i, 3] / matriz[8, 3]
      matriz[i, 5] <- 1 - pf(matriz[i, 4], matriz[i, 1], matriz[8, 1])
    }
  }

  # ---- standard errors, with Me the Error mean square of the full table ----
  Me <- matriz[9, 3]
  SE <- list(
    gca_line     = sqrt(Me / (n_r * n_t)),
    gca_tester   = sqrt(Me / (n_r * n_l)),
    sca          = sqrt(Me / n_r),
    gi_gj_line   = sqrt(2 * Me / (n_r * n_t)),
    gi_gj_tester = sqrt(2 * Me / (n_r * n_l)),
    sij_skl      = sqrt(2 * Me / n_r)
  )

  # ---- genetic components ----
  MS_l  <- matriz[6, 3]
  MS_t  <- matriz[7, 3]
  MS_lt <- matriz[8, 3]

  cov1 <- (MS_l - MS_lt) / (n_r * n_t)                      # Cov HS, lines
  cov2 <- (MS_t - MS_lt) / (n_r * n_l)                      # Cov HS, testers
  cov3 <- (((n_l - 1) * MS_l + (n_t - 1) * MS_t) / (n_l + n_t - 2) - MS_lt) /
          (n_r * (2 * n_l * n_t - n_l - n_t))               # Cov HS, average
  cov4 <- ((MS_l - Me) + (MS_t - Me) + (MS_lt - Me)) / (3 * n_r)  # Cov FS, average

  Fcoef <- 0
  varA0 <- cov3 * (4 / (1 + Fcoef))
  varD0 <- ((MS_lt - Me) / n_r) * (2 / (1 + Fcoef))

  Fcoef <- 1
  varA1 <- cov3 * (4 / (1 + Fcoef))
  varD1 <- ((MS_lt - Me) / n_r) * (2 / (1 + Fcoef))

  Genetic <- list(
    Cov_HS_line   = cov1,
    Cov_HS_tester = cov2,
    Cov_HS_avg    = cov3,
    Cov_FS_avg    = cov4,
    VarA_F0 = varA0, VarA_F1 = varA1,
    VarD_F0 = varD0, VarD_F1 = varD1
  )

  # ---- proportional contribution to the Crosses sum of squares ----
  Contributions <- list(
    lines   = matriz[6, 2] * 100 / matriz[5, 2],
    testers = matriz[7, 2] * 100 / matriz[5, 2],
    lxt     = matriz[8, 2] * 100 / matriz[5, 2]
  )

  return(list(
    ANOVA_with_parents = matriz1,
    ANOVA_LxT          = matriz4,
    ANOVA_full         = matriz,
    GCA_lines          = GCA.lines,
    GCA_testers        = GCA.testers,
    SCA                = SCA,
    SE                 = SE,
    Genetic            = Genetic,
    Contribution       = Contributions
  ))
}

dataset <- data.frame(
  line = c(rep(c("L1","L2","L3","L4","L5"), each = 12),
           rep(c("L1","L2","L3","L4","L5"), each = 4),
           rep(NA, 12)),
  tester = c(rep(rep(c("T6","T7","T8"), each = 4), 5),
             rep(NA, 20),
             rep(c("T6","T7","T8"), each = 4)),
  Replication = rep(1:4, 23),
  Yield = c(
    74.4, 70.86, 60.94, 68, 91.82, 99.18, 118.88, 120.68, 48.08, 62.1, 58.54, 41.84,
    59.06, 65.62, 81.62, 86.76, 84.16, 109.74, 102.14, 94.52, 96.92, 91.44, 79.86, 74.38,
    109.86, 98.16, 93.26, 102.26, 117.2, 100.28, 116.16, 112.52, 109.68, 116.48, 123.92, 120.86,
    103.14, 109.66, 90.98, 119.4, 53.4, 60.86, 74.46, 69.08, 53.86, 48.3, 40.64, 44.62,
    98.46, 73.1, 89.18, 75.86, 81.36, 72.82, 89.82, 83.74, 86.62, 94.18, 90.32, 108.16,
    104.86, 84.32, 76.92, 76.48,
    88.02, 106.54, 89.82, 108.68,
    77.94, 71.34, 77.52, 69.48,
    80.82, 106.52, 83.28, 95.92,
    59.96, 52.48, 52.98, 50.98,
    96.44, 98.82, 99.14, 107.16,
    91.44, 99.66, 83.28, 89.46,
    91.78, 84.82, 69.92, 81.48
  ),
  stringsAsFactors = FALSE
)

trait_name <- "Yield"

A     <- factor(dataset$line)         # lines
B     <- factor(dataset$tester)       # testers
Block <- factor(dataset$Replication)  # replications

r <- nlevels(Block)   # 4 replications
l <- nlevels(A)       # 5 lines
n_t <- nlevels(B)     # 3 testers

out <- linetester(replications = Block,
                  lines        = A,
                  testers      = B,
                  y            = dataset[[trait_name]])

parent_result <- as.data.frame(out$ANOVA_with_parents)  # Replications, Treatments, Error
full_result   <- as.data.frame(out$ANOVA_full)          # 9 sources + Total

out$GCA_lines       # GCA effects of the five lines
out$GCA_testers     # GCA effects of the three testers
out$SCA             # SCA effects, 5 x 3 matrix of crosses
out$SE              # gca_line, gca_tester, sca, gi_gj_line, gi_gj_tester, sij_skl
out$Contribution    # lines, testers, lxt (% of the Crosses sum of squares)
out$Genetic         # Cov_HS_line/tester/avg, Cov_FS_avg, VarA_F0/F1, VarD_F0/F1
                    # the module reports the F = 1 values

# Error term driving the mean comparisons and the tests of effects
MSE <- as.numeric(full_result[9, 3])
DFE <- as.numeric(full_result[9, 1])

alpha <- 0.05

# --- GCA of lines ---
t_gca_l <- as.numeric(out$GCA_lines) / out$SE$gca_line

# --- GCA of testers ---
t_gca_t <- as.numeric(out$GCA_testers) / out$SE$gca_tester

# --- SCA of crosses ---
t_sca <- as.numeric(out$SCA) / out$SE$sca

datos <- data.frame(Replications = Block,
                    Lines        = A,
                    Testers      = B,
                    Y            = dataset[[trait_name]])

# error mean square and df from the ANOVA WITH PARENTS
Me_par   <- as.numeric(parent_result[3, 3])
df_par   <- as.numeric(parent_result[3, 1])

# cross means, line parent means, tester parent means
cross_mean <- tapply(datos$Y[!is.na(datos$Lines) & !is.na(datos$Testers)],
                     paste(datos$Lines, datos$Testers)[
                       !is.na(datos$Lines) & !is.na(datos$Testers)],
                     mean, na.rm = TRUE)

line_par <- subset(datos, !is.na(Lines) & is.na(Testers))
P1_mean  <- tapply(line_par$Y, droplevels(line_par$Lines), mean, na.rm = TRUE)

test_par <- subset(datos, is.na(Lines) & !is.na(Testers))
P2_mean  <- tapply(test_par$Y, droplevels(test_par$Testers), mean, na.rm = TRUE)

# one row per cross
grid <- expand.grid(line   = levels(A),
                    tester = levels(B),
                    stringsAsFactors = FALSE)
grid$cross <- paste(grid$line, grid$tester)

F1 <- cross_mean[grid$cross]
P1 <- P1_mean[grid$line]
P2 <- P2_mean[grid$tester]

MP <- (P1 + P2) / 2                  # mid-parent
BP <- pmax(P1, P2)                   # better parent

# percentage heterosis
mid_parent   <- ((F1 - MP) / MP) * 100
heterobelt   <- ((F1 - BP) / BP) * 100

# standard errors
SE_H1 <- sqrt((1.5 * Me_par) / r)    # mid-parent
SE_H2 <- sqrt((2   * Me_par) / r)    # better parent
SE_SH <- sqrt((2   * Me_par) / r)    # standard check

# t test, p value and critical differences (shown for heterobeltiosis)
t_val <- (F1 - BP) / SE_H2
p_val <- 2 * pt(-abs(t_val), df = df_par)
CD5   <- SE_H2 * qt(0.975, df = df_par)
CD1   <- SE_H2 * qt(0.995, df = df_par)

MSt <- as.numeric(full_result[2, 3])   # Treatments mean square
MSe <- as.numeric(full_result[9, 3])   # Error mean square

grand_mean <- mean(dataset[[trait_name]], na.rm = TRUE)

sigma2_g <- max((MSt - MSe) / r, 0)    # negative estimates set to zero
sigma2_e <- MSe
sigma2_p <- sigma2_g + sigma2_e

GCV <- (sqrt(sigma2_g) / grand_mean) * 100
PCV <- (sqrt(sigma2_p) / grand_mean) * 100
ECV <- (sqrt(sigma2_e) / grand_mean) * 100

h2_bs <- ifelse(sigma2_p > 0, (sigma2_g / sigma2_p) * 100, 0)  # broad sense (%)

k   <- 2.06                                # 5% selection intensity
GA  <- k * sqrt(sigma2_p) * (h2_bs / 100)  # genetic advance
GAM <- (GA / grand_mean) * 100             # GA as % of mean