R/Phenotypic_Path.R

In TBA: Collection of 'shiny' Apps for Tree Breeding Analysis

Phenotypic_Path<- function(data) {
  old_options <- options(scipen = 999)  # Save current options
  on.exit(options(old_options))         # Restore options when function exits


  # Convert the first two columns to factor type
  data[, 1:2] <- lapply(data[, 1:2], as.factor)

  # Convert the remaining columns to numeric
  data[, -c(1, 2)] <- lapply(data[, -c(1, 2)], as.numeric)

  # Extract trait names (excluding the first two columns)
  traits <- names(data)[-c(1, 2)][sapply(data[-c(1, 2)], is.numeric)]

  # Prepare a matrix to store correlations
  correlation_matrix <- matrix(NA, nrow = length(traits), ncol = length(traits))
  formatted_correlation_matrix <- matrix(NA, nrow = length(traits), ncol = length(traits))

  # Calculate correlations for each pair of traits
  for (i in 1:length(traits)) {
    for (j in 1:length(traits)) {
      trait1 <- traits[i]
      trait2 <- traits[j]

      if (i == j) {
        correlation_matrix[i, j] <- 1 # Set correlation to 1 if it's the same trait
        formatted_correlation_matrix[i, j] <- 1  # Set formatted correlation value to 1
      } else {
        # Perform linear regression for trait1
        formula1 <- as.formula(paste0("`", trait1, "` ~ `", names(data)[1], "` + `", names(data)[2], "`"))
        model1 <- lm(formula1, data = data)
        anova_result1 <- anova(model1)

        # Perform linear regression for trait2
        formula2 <- as.formula(paste0("`", trait2, "` ~ `", names(data)[1], "` + `", names(data)[2], "`"))
        model2 <- lm(formula2, data = data)
        anova_result2 <- anova(model2)

        # Calculate phenotypic variance for trait1 and trait2
        replication_levels <- nlevels(data[[1]])
        genotypic_variance1 <- round((anova_result1$`Mean Sq`[2] - anova_result1$`Mean Sq`[3]) / replication_levels,4)

        phenotypic_variance1<-round(genotypic_variance1+anova_result1$`Mean Sq`[3],4)
        genotypic_variance2 <- round((anova_result2$`Mean Sq`[2] - anova_result2$`Mean Sq`[3]) / replication_levels,4)
        phenotypic_variance2<-round(genotypic_variance2+anova_result2$`Mean Sq`[3],4)

        # Calculate covariance sums
        total_of_genotypes_trait1 <- tapply(data[[trait1]], data[[2]], sum)
        total_of_genotypes_trait2 <- tapply(data[[trait2]], data[[2]], sum)
        total_of_replication_trait1 <- tapply(data[[trait1]], data[[1]], sum)
        total_of_replication_trait2 <- tapply(data[[trait2]], data[[1]], sum)

        number_of_replication <- nlevels(data[[1]])
        number_of_genotype <- nlevels(data[[2]])
        Grand_total_trait1 <- sum(data[[trait1]])
        Grand_total_trait2 <- sum(data[[trait2]])
        CF <- (Grand_total_trait1 * Grand_total_trait2) / (number_of_replication * number_of_genotype)

        Total_SP <- round(sum(data[[trait1]] * data[[trait2]]) - CF,4)
        Genotypic_SP <- round((sum(total_of_genotypes_trait1 * total_of_genotypes_trait2) / number_of_replication) - CF,4)
        Replication_SP <- round((sum(total_of_replication_trait1 * total_of_replication_trait2) / number_of_genotype) - CF,4)
        Error_SP <- Total_SP - Genotypic_SP - Replication_SP
        DF_Replication <- number_of_replication - 1
        DF_Genotypes <- number_of_genotype - 1
        DF_Error <- DF_Replication * DF_Genotypes
        Replication_MP <- round(Replication_SP / DF_Replication,4)
        Genotypic_MP <- round(Genotypic_SP / DF_Genotypes,4)
        Error_MP <- Error_SP / DF_Error

        Genotypic_Covariance <- round((Genotypic_MP - Error_MP) / number_of_replication,4)

        Phenotypic_Covariance <- round(Error_MP + Genotypic_Covariance,4)

        # Calculate correlation
        correlation <- round(Phenotypic_Covariance / sqrt(phenotypic_variance1 * phenotypic_variance2), 4)

        # Perform significance test
        n <- number_of_replication * number_of_genotype # Number of observations
        df <- n - 2  # Degrees of freedom for Pearson correlation
        if (!is.nan(correlation)&& !is.na(correlation)) {
          t_stat <- (correlation) * (sqrt(df / (1 - (correlation)^2)))  # Calculate t-statistic
          p_value <- 2 * pt(abs(t_stat), df = df, lower.tail = FALSE)  # Calculate two-tailed p-value
        } else {
          t_stat <- NA
          p_value <- NA
        }

        # Determine significance level symbol
        if (!is.nan(t_stat) && !is.na(t_stat)) {
          if (p_value < 0.05) {
            significance_symbol <- "*"  # Significant at 5%
          } else {
            significance_symbol <- "NS"  # Non-significant
          }
        } else {
          significance_symbol <- ""  # No significance symbol if t_stat is NA and NaN
        }

        # Store correlation value in the matrices

        formatted_correlation_matrix[i, j] <- correlation
        correlation_matrix[i, j] <- paste0(format(correlation, scientific = FALSE), significance_symbol)
      }
    }
  }
  phenotypic_correlation_matrix <- noquote(correlation_matrix)
  correlation_only <- noquote(formatted_correlation_matrix)

  # Path Analysis
  dependent_variable <- correlation_only[1:(length(traits) - 1), length(traits)]
  dependent_variable_matrix <- matrix(dependent_variable, ncol = 1)
  independent_variable <- correlation_only[1:(length(traits) - 1), 1:(length(traits) - 1)]
  direct_effect <- solve(independent_variable,dependent_variable_matrix)
  Direct_and_indirect_effect <- matrix(nrow = (length(traits) - 1), ncol = (length(traits) - 1))

  for (i in 1:(length(traits) - 1)) {
    for (j in 1:(length(traits) - 1)) {
      Direct_and_indirect_effect[i, j] <- round(direct_effect[j] * independent_variable[i, j], 4)
    }
  }

  Path_effects <- cbind(Direct_and_indirect_effect, phenotypic_correlation_matrix[1:(length(traits) - 1), length(traits)])
  rownames(Path_effects) <- traits[1:(length(traits) - 1)]
  colnames(Path_effects) <- traits[1:(length(traits))]

  residual <- 1 - t(direct_effect) %*% dependent_variable_matrix
  Residual_effect <- round(sqrt(residual), 4)

  rownames(Direct_and_indirect_effect) <- traits[1:(length(traits) - 1)]
  colnames(Direct_and_indirect_effect) <- traits[1:(length(traits) - 1)]
  rownames(Residual_effect)<-"Residual Effect"

  # Convert Path_effects to data frame
  Path_effects <- as.data.frame(Path_effects)

  # Return the data frame with row names and residual effect
  return(list(Path_effects = Path_effects, Residual_effect = Residual_effect))
}