R/augmented_pooled.R

In AugmentedPooledRCBD: Augmented Pooled Randomized Complete Block Design Analysis

#' @importFrom dplyr %>% group_by summarise n
#' @importFrom emmeans emmeans
#' @importFrom rlang .data
#' @importFrom stats anova aov bartlett.test contr.helmert
#'   contrasts<- deviance df.residual lm pf qt residuals sd

augmented_pooled_analysis <- function(
    data,
    trait,
    checks,
    env="env",
    blk="blk",
    trt="trt",
    alpha=0.05
){


  names(data)[names(data)==env] <- "env"
  names(data)[names(data)==blk] <- "blk"
  names(data)[names(data)==trt] <- "trt"

  data$env <- factor(data$env)
  data$blk <- factor(data$blk)
  data$trt <- factor(data$trt)

  data$type <- ifelse(data$trt %in% checks,"control","treatment")

  results <- list()



  results$Descriptive <- data %>%
    group_by(env) %>%
    summarise(
      Mean=mean(.data[[trait]]),
      SD=sd(.data[[trait]]),
      Min=min(.data[[trait]]),
      Max=max(.data[[trait]]),
      N=n(),
      .groups="drop"
    )



  envs <- levels(data$env)

  MSE_vec <- c()
  residual_data <- data.frame()

  for(e in envs){

    sub <- data[data$env == e, ]

    model <- lm(
      sub[[trait]] ~ blk + trt,
      data = sub
    )

    MSE_vec[e] <- anova(model)["Residuals", "Mean Sq"]

    tmp <- data.frame(
      env = e,
      residual = residuals(model)
    )

    residual_data <- rbind(
      residual_data,
      tmp
    )
  }



  bart <- bartlett.test(
    residual ~ env,
    data = residual_data
  )

  y <- residual_data$residual
  grp <- residual_data$env

  means <- tapply(y, grp, mean)

  z <- abs(y - means[grp])

  lev_tab <- summary(aov(z ~ grp))[[1]]

  lev_F <- lev_tab[1, "F value"]
  lev_p <- lev_tab[1, "Pr(>F)"]

  results$Variance <- data.frame(
    Test = c("Bartlett", "Levene"),
    Statistic = c(
      as.numeric(bart$statistic),
      lev_F
    ),
    p_value = c(
      bart$p.value,
      lev_p
    )
  )



  if(bart$p.value > alpha & lev_p > alpha){

    transform <- FALSE
    decision <- "Variances homogeneous"

  }else{

    transform <- TRUE
    decision <- "Variances heterogeneous pooling with transformation"


  }



  if(transform){

    data$y <- data[[trait]]/sqrt(MSE_vec[as.character(data$env)])
    results$Transformed_Data <- data[,c("env","blk","trt","y")]



    model_t <- lm(y ~ env + blk%in%env + trt, data=data)

    adj_t <- suppressMessages(
      as.data.frame(emmeans(model_t,"trt"))
    )

    adj_t$type <- ifelse(adj_t$trt %in% checks,"control","treatment")

    results$Transformed_Means <- adj_t

  }else{

    data$y <- data[[trait]]
    results$Transformed_Data <- NULL
    results$Transformed_Means <- NULL

  }



  for(e in envs){

    sub <- data[data$env==e,]



    m0 <- lm(sub[[trait]] ~ 1, data=sub)
    m1 <- lm(sub[[trait]] ~ trt, data=sub)
    m2 <- lm(sub[[trait]] ~ trt + blk, data=sub)

    ss_trt_I <- deviance(m0) - deviance(m1)
    ss_blk_I <- deviance(m1) - deviance(m2)

    df_trt_I <- df.residual(m0) - df.residual(m1)
    df_blk_I <- df.residual(m1) - df.residual(m2)
    df_err <- df.residual(m2)

    mse <- deviance(m2) / df_err
    ms_trt_I <- ss_trt_I / df_trt_I
    ms_blk_I <- ss_blk_I / df_blk_I

    F_trt_I <- ms_trt_I / mse
    F_blk_I <- ms_blk_I / mse

    p_trt_I <- pf(F_trt_I, df_trt_I, df_err, lower.tail=FALSE)
    p_blk_I <- pf(F_blk_I, df_blk_I, df_err, lower.tail=FALSE)


    old_contrasts <- getOption("contrasts")
    options(contrasts=c("contr.sum","contr.poly"))
    model_type3 <- lm(sub[[trait]] ~ trt + blk, data=sub)
    type3_tab <- car::Anova(model_type3, type=3)
    options(contrasts=old_contrasts)


    model <- model_type3


    type1_df <- data.frame(
      Type="Type I",
      Source=c("trt","blk","Residuals"),
      Df=c(df_trt_I,df_blk_I,df_err),
      SumSq=c(ss_trt_I,ss_blk_I,deviance(m2)),
      MeanSq=c(ms_trt_I,ms_blk_I,mse),
      Fvalue=c(F_trt_I,F_blk_I,NA_real_),
      p_value=c(p_trt_I,p_blk_I,NA_real_)
    )

    type3_df <- data.frame(
      Type="Type III",
      Source=c("trt","blk","Residuals"),
      Df=c(
        type3_tab["trt","Df"],
        type3_tab["blk","Df"],
        type3_tab["Residuals","Df"]
      ),
      SumSq=c(
        type3_tab["trt","Sum Sq"],
        type3_tab["blk","Sum Sq"],
        type3_tab["Residuals","Sum Sq"]
      ),
      MeanSq=c(
        type3_tab["trt","Sum Sq"]/type3_tab["trt","Df"],
        type3_tab["blk","Sum Sq"]/type3_tab["blk","Df"],
        type3_tab["Residuals","Sum Sq"]/type3_tab["Residuals","Df"]
      ),
      Fvalue=c(
        type3_tab["trt","F value"],
        type3_tab["blk","F value"],
        NA_real_
      ),
      p_value=c(
        type3_tab["trt","Pr(>F)"],
        type3_tab["blk","Pr(>F)"],
        NA_real_
      )
    )

    results[[paste0("ANOVA_",e)]] <- rbind(type1_df,type3_df)
    adj <- suppressMessages(
      as.data.frame(emmeans(model,"trt"))
    )

    adj$type <- ifelse(adj$trt %in% checks,"control","treatment")

    results[[paste0("Means_",e)]] <- adj



    n_check <- length(checks)
    n_trt   <- nlevels(sub$trt)

    df_check <- n_check - 1
    df_trt   <- n_trt - 1


    check_levels <- as.character(checks)

    test_levels <- setdiff(
      levels(sub$trt),
      check_levels
    )

    trt_order <- c(
      check_levels,
      test_levels
    )

    trt_aug <- factor(
      sub$trt,
      levels = trt_order
    )



    contr.augmented <- function(n1, n2){

      m1 <- contr.helmert(n1)
      m2 <- contr.helmert(n2)

      m10 <- cbind(
        m1,
        matrix(
          0,
          nrow = nrow(m1),
          ncol = ncol(m2)
        )
      )

      m02 <- cbind(
        matrix(
          0,
          nrow = nrow(m2),
          ncol = ncol(m1)
        ),
        m2
      )

      rbind(m10, m02)
    }



    old_contrasts <- getOption("contrasts")

    contrasts(trt_aug) <- contr.augmented(
      df_check + 1,
      df_trt - df_check
    )



    aug_model <- aov(
      sub[[trait]] ~ trt_aug + blk,
      data = sub
    )



    aug_aov <- summary(
      aug_model,
      split = list(
        trt_aug = list(
          Check = 1:df_check,
          Test = (df_check + 1):(df_trt - 1),
          `Test vs. Check` = df_trt
        )
      )
    )


    options(contrasts = old_contrasts)



    A2 <- aug_aov[[1]]

    results[[paste0("Partition_",e)]] <- data.frame(

      Source = c(
        "Treatment (ignoring Blocks)",
        "Treatment: Check",
        "Treatment: Test",
        "Treatment: Test vs. Check"
      ),

      Df = A2[1:4, "Df"],

      SumSq = A2[1:4, "Sum Sq"],

      MeanSq = A2[1:4, "Mean Sq"],

      Fvalue = A2[1:4, "F value"],

      p_value = A2[1:4, "Pr(>F)"],

      row.names = NULL
    )



    MSE <- mse
    dfE <- df_err


    r <- sum(sub$type == "control") / length(checks)


    grand_mean <- mean(sub[[trait]], na.rm = TRUE)


    CV <- (sqrt(MSE) / grand_mean) * 100


    SEm_check <- sqrt(MSE / r)
    SEm_test  <- sqrt(MSE)


    n_checks <- length(checks)

    SEd_CC <- sqrt(2 * MSE / r)
    SEd_TT_same <- sqrt(2 * MSE)
    SEd_TT_diff <- sqrt(2 * MSE * (1 + 1/n_checks))
    SEd_TC <- sqrt(
      MSE * (1 + 1/r + 1/n_checks + 1/(r * n_checks))
    )


    tval <- qt(
      1 - alpha/2,
      dfE
    )


    CD_CC <- tval * SEd_CC
    CD_TT_same <- tval * SEd_TT_same
    CD_TT_diff <- tval * SEd_TT_diff
    CD_TC <- tval * SEd_TC



    results[[paste0("Precision_",e)]] <- data.frame(
      Statistic = c(
        "MSE",
        "CV (%)",
        "SEm - Check",
        "SEm - Test",
        "SEd - Check vs Check",
        "SEd - Test vs Test (Same Block)",
        "SEd - Test vs Test (Different Blocks)",
        "SEd - Test vs Check",
        "CD - Check vs Check",
        "CD - Test vs Test (Same Block)",
        "CD - Test vs Test (Different Blocks)",
        "CD - Test vs Check"
      ),
      Value = c(
        MSE,
        CV,
        SEm_check,
        SEm_test,
        SEd_CC,
        SEd_TT_same,
        SEd_TT_diff,
        SEd_TC,
        CD_CC,
        CD_TT_same,
        CD_TT_diff,
        CD_TC
      )
    )



    results[[paste0("CD_",e)]] <- data.frame(
      Comparison = c(
        "Control-Control",
        "Treatment-Treatment (Same Block)",
        "Treatment-Treatment (Different Blocks)",
        "Treatment-Control"
      ),
      CD = c(
        CD_CC,
        CD_TT_same,
        CD_TT_diff,
        CD_TC
      )
    )

    chk <- adj[adj$type == "control",]

    best <- max(chk$emmean)

    adj$Status <- "At Par"

    adj$Status[adj$type=="treatment" &
                 adj$emmean > best + results[[paste0("CD_",e)]]$CD[4]] <- "Superior"

    adj$Status[adj$type=="treatment" &
                 adj$emmean < best - results[[paste0("CD_",e)]]$CD[4]] <- "Inferior"


    adj <- adj[order(-adj$emmean), ]


    adj$Rank <- seq_len(nrow(adj))


    adj <- adj[, c(
      "Rank",
      "trt",
      "emmean",
      "SE",
      "df",
      "lower.CL",
      "upper.CL",
      "type",
      "Status"
    )]

    results[[paste0("Ranking_",e)]] <- adj

  }




  m0 <- lm(y ~ 1, data=data)
  m1 <- lm(y ~ env, data=data)
  m2 <- lm(y ~ env + env:blk, data=data)
  m3 <- lm(y ~ env + env:blk + trt, data=data)
  m4 <- lm(y ~ env + env:blk + trt + env:trt, data=data)

  ss_env_I <- deviance(m0) - deviance(m1)
  ss_blk_env_I <- deviance(m1) - deviance(m2)
  ss_trt_I <- deviance(m2) - deviance(m3)
  ss_env_trt_I <- deviance(m3) - deviance(m4)
  ss_error <- deviance(m4)

  df_env_I <- df.residual(m0) - df.residual(m1)
  df_blk_env_I <- df.residual(m1) - df.residual(m2)
  df_trt_I <- df.residual(m2) - df.residual(m3)
  df_env_trt_I <- df.residual(m3) - df.residual(m4)
  df_error <- df.residual(m4)

  ms_env_I <- ss_env_I / df_env_I
  ms_blk_env_I <- ss_blk_env_I / df_blk_env_I
  ms_trt_I <- ss_trt_I / df_trt_I
  ms_env_trt_I <- ss_env_trt_I / df_env_trt_I
  mse <- ss_error / df_error

  F_env_I <- ms_env_I / mse
  F_blk_env_I <- ms_blk_env_I / mse
  F_trt_I <- ms_trt_I / mse
  F_env_trt_I <- ms_env_trt_I / mse

  type1_pooled <- data.frame(
    Type="Type I",
    Source=c("env","env:blk","trt","env:trt","Residuals"),
    Df=c(df_env_I,df_blk_env_I,df_trt_I,df_env_trt_I,df_error),
    SumSq=c(ss_env_I,ss_blk_env_I,ss_trt_I,ss_env_trt_I,ss_error),
    MeanSq=c(ms_env_I,ms_blk_env_I,ms_trt_I,ms_env_trt_I,mse),
    Fvalue=c(F_env_I,F_blk_env_I,F_trt_I,F_env_trt_I,NA_real_),
    p_value=c(
      pf(F_env_I,df_env_I,df_error,lower.tail=FALSE),
      pf(F_blk_env_I,df_blk_env_I,df_error,lower.tail=FALSE),
      pf(F_trt_I,df_trt_I,df_error,lower.tail=FALSE),
      pf(F_env_trt_I,df_env_trt_I,df_error,lower.tail=FALSE),
      NA_real_
    )
  )


  old_contrasts <- getOption("contrasts")
  options(contrasts=c("contr.sum","contr.poly"))
  model_type3 <- lm(y ~ env + env:blk + trt + env:trt, data=data)
  type3_pooled <- car::Anova(model_type3, type=3)
  options(contrasts=old_contrasts)

  type3_pooled_df <- data.frame(
    Type="Type III",
    Source=c("env","env:blk","trt","env:trt","Residuals"),
    Df=c(
      type3_pooled["env","Df"],
      type3_pooled["env:blk","Df"],
      type3_pooled["trt","Df"],
      type3_pooled["env:trt","Df"],
      type3_pooled["Residuals","Df"]
    ),
    SumSq=c(
      type3_pooled["env","Sum Sq"],
      type3_pooled["env:blk","Sum Sq"],
      type3_pooled["trt","Sum Sq"],
      type3_pooled["env:trt","Sum Sq"],
      type3_pooled["Residuals","Sum Sq"]
    ),
    MeanSq=c(
      type3_pooled["env","Sum Sq"]/type3_pooled["env","Df"],
      type3_pooled["env:blk","Sum Sq"]/type3_pooled["env:blk","Df"],
      type3_pooled["trt","Sum Sq"]/type3_pooled["trt","Df"],
      type3_pooled["env:trt","Sum Sq"]/type3_pooled["env:trt","Df"],
      type3_pooled["Residuals","Sum Sq"]/type3_pooled["Residuals","Df"]
    ),
    Fvalue=c(
      type3_pooled["env","F value"],
      type3_pooled["env:blk","F value"],
      type3_pooled["trt","F value"],
      type3_pooled["env:trt","F value"],
      NA_real_
    ),
    p_value=c(
      type3_pooled["env","Pr(>F)"],
      type3_pooled["env:blk","Pr(>F)"],
      type3_pooled["trt","Pr(>F)"],
      type3_pooled["env:trt","Pr(>F)"],
      NA_real_
    )
  )

  results$Pooled_ANOVA <- rbind(type1_pooled,type3_pooled_df)


  model <- m4
  adj <- suppressMessages(
    as.data.frame(emmeans(model,"trt"))
  )

  adj$type <- ifelse(adj$trt %in% checks,"control","treatment")

  results$Pooled_Means <- adj



  n_check <- length(checks)
  n_trt   <- nlevels(data$trt)

  df_check <- n_check - 1
  df_trt   <- n_trt - 1


  check_levels <- as.character(checks)

  test_levels <- setdiff(
    levels(data$trt),
    check_levels
  )

  trt_order <- c(
    check_levels,
    test_levels
  )

  trt_aug <- factor(
    data$trt,
    levels = trt_order
  )



  contr.augmented.pooled <- function(n1, n2) {

    m1 <- contr.helmert(n1)
    m2 <- contr.helmert(n2)

    m10 <- cbind(
      m1,
      matrix(
        0,
        nrow = nrow(m1),
        ncol = ncol(m2)
      )
    )

    m02 <- cbind(
      matrix(
        0,
        nrow = nrow(m2),
        ncol = ncol(m1)
      ),
      m2
    )

    rbind(m10, m02)
  }



  old_contrasts <- getOption("contrasts")

  contrasts(trt_aug) <- contr.augmented.pooled(
    df_check + 1,
    df_trt - df_check
  )



  pooled_aug_model <- aov(
    y ~ env + env:blk + trt_aug + env:trt_aug,
    data = data
  )



  pooled_aug_aov <- summary(
    pooled_aug_model,
    split = list(
      trt_aug = list(
        Check = 1:df_check,
        Test = (df_check + 1):(df_trt - 1),
        `Test vs. Check` = df_trt
      )
    )
  )


  options(contrasts = old_contrasts)



  A_pool <- pooled_aug_aov[[1]]

  rn <- rownames(A_pool)

  i_check <- which(grepl("Check", rn) &
                     grepl("trt_aug", rn))[1]

  i_test <- which(grepl("Test", rn) &
                    !grepl("Check", rn) &
                    !grepl("vs", rn) &
                    grepl("trt_aug", rn))[1]

  i_con <- which(grepl("Test vs. Check", rn))[1]



  trt_total_df <- A_pool[i_check, "Df"] +
    A_pool[i_test, "Df"] +
    A_pool[i_con, "Df"]

  trt_total_ss <- A_pool[i_check, "Sum Sq"] +
    A_pool[i_test, "Sum Sq"] +
    A_pool[i_con, "Sum Sq"]

  trt_total_ms <- trt_total_ss / trt_total_df

  pooled_mse <- type3_pooled["Residuals", "Sum Sq"] /
    type3_pooled["Residuals", "Df"]

  trt_total_F <- trt_total_ms / pooled_mse

  trt_total_p <- pf(
    trt_total_F,
    trt_total_df,
    type3_pooled["Residuals", "Df"],
    lower.tail = FALSE
  )


  rn <- rownames(A_pool)

  i_check <- which(grepl("Check", rn) &
                     grepl("trt_aug", rn))[1]

  i_test <- which(grepl("Test", rn) &
                    !grepl("Check", rn) &
                    !grepl("vs", rn) &
                    grepl("trt_aug", rn))[1]

  i_con <- which(grepl("Test vs. Check", rn))[1]


  results$Pooled_Partition <- data.frame(

    Source = c(
      "Treatment (ignoring Blocks)",
      "Treatment: Check",
      "Treatment: Test",
      "Treatment: Test vs. Check"
    ),

    Df = c(
      trt_total_df,
      A_pool[i_check, "Df"],
      A_pool[i_test, "Df"],
      A_pool[i_con, "Df"]
    ),

    SumSq = c(
      trt_total_ss,
      A_pool[i_check, "Sum Sq"],
      A_pool[i_test, "Sum Sq"],
      A_pool[i_con, "Sum Sq"]
    ),

    MeanSq = c(
      trt_total_ms,
      A_pool[i_check, "Mean Sq"],
      A_pool[i_test, "Mean Sq"],
      A_pool[i_con, "Mean Sq"]
    ),

    Fvalue = c(
      trt_total_F,
      A_pool[i_check, "F value"],
      A_pool[i_test, "F value"],
      A_pool[i_con, "F value"]
    ),

    p_value = c(
      trt_total_p,
      A_pool[i_check, "Pr(>F)"],
      A_pool[i_test, "Pr(>F)"],
      A_pool[i_con, "Pr(>F)"]
    ),

    row.names = NULL
  )

  r <- sum(data$type=="control") / length(checks)

  tval <- qt(
    1 - alpha/2,
    df_error
  )



  MSE_pool <- pooled_mse
  dfE_pool <- df_error


  grand_mean_pool <- mean(data$y, na.rm = TRUE)


  CV_pool <- (sqrt(MSE_pool) / grand_mean_pool) * 100


  SEm_check_pool <- sqrt(MSE_pool / r)
  SEm_test_pool  <- sqrt(MSE_pool)


  n_checks_pool <- length(checks)

  SEd_CC_pool <- sqrt(2 * MSE_pool / r)
  SEd_TT_same_pool <- sqrt(2 * MSE_pool)
  SEd_TT_diff_pool <- sqrt(
    2 * MSE_pool * (1 + 1/n_checks_pool)
  )
  SEd_TC_pool <- sqrt(
    MSE_pool * (
      1 + 1/r + 1/n_checks_pool + 1/(r * n_checks_pool)
    )
  )



  results$Pooled_CD <- data.frame(
    Comparison = c(
      "Control-Control",
      "Treatment-Treatment (Same Block)",
      "Treatment-Treatment (Different Blocks)",
      "Treatment-Control"
    ),
    CD = c(
      tval * SEd_CC_pool,
      tval * SEd_TT_same_pool,
      tval * SEd_TT_diff_pool,
      tval * SEd_TC_pool
    )
  )

  results$Pooled_Precision <- data.frame(
    Statistic = c(
      "MSE",
      "CV (%)",
      "SEm - Check",
      "SEm - Test",
      "SEd - Check vs Check",
      "SEd - Test vs Test (Same Block)",
      "SEd - Test vs Test (Different Blocks)",
      "SEd - Test vs Check",
      "CD - Check vs Check",
      "CD - Test vs Test (Same Block)",
      "CD - Test vs Test (Different Blocks)",
      "CD - Test vs Check"
    ),
    Value = c(
      MSE_pool,
      CV_pool,
      SEm_check_pool,
      SEm_test_pool,
      SEd_CC_pool,
      SEd_TT_same_pool,
      SEd_TT_diff_pool,
      SEd_TC_pool,
      tval * SEd_CC_pool,
      tval * SEd_TT_same_pool,
      tval * SEd_TT_diff_pool,
      tval * SEd_TC_pool
    )
  )
  best <- max(chk$emmean)

  adj$Status <- "At Par"

  adj$Status[adj$type=="treatment" &
               adj$emmean > best + results$Pooled_CD$CD[4]] <- "Superior"

  adj$Status[adj$type=="treatment" &
               adj$emmean < best - results$Pooled_CD$CD[4]] <- "Inferior"


  adj <- adj[order(-adj$emmean), ]


  adj$Rank <- seq_len(nrow(adj))


  adj <- adj[, c(
    "Rank",
    "trt",
    "emmean",
    "SE",
    "df",
    "lower.CL",
    "upper.CL",
    "type",
    "Status"
  )]

  results$Pooled_Ranking <- adj
  results$Decision <- decision

  return(results)

}