Channel 1: Reduce Phenotyping Error

The Artemis tool improves phenotyping accuracy compared to manual methods, as evidenced by greater agreement among multiple scorers. This higher agreement, measured by the intra-class correlation coefficient (ICC), reflects reduced phenotyping error (residual variance). Lower phenotyping error, in turn, leads to increased heritability estimates, which enhance genetic gain under controlled conditions. Later, we will link this genetic gain to improvements realized in target production environments and the resulting economic outcomes.

The ICC is the portion of the variance in the phenotyping data attributable to real underlying variation in trait expression, given by:

\[ \text{ICC} = \frac{\sigma_{\text{gt}}^2}{\sigma_{\text{gt}}^2 + \sigma_\epsilon^2} \]

Where: \(\sigma_{\text{gt}}^2\) is the variance of the groundtruth data (including both genetic variance (\(\sigma_g^2\)) and environmental variance (\(\sigma_e^2\))), and \(\sigma_\epsilon^2\) is the residual variance (i.e. the phenotyping error, or variance caused by scorer inconsistencies)

This can be approximated using the following estimator: \[ \widehat{\text{ICC}} = \frac{\text{MSB} - \text{MSE}}{\text{MSB} + (k - 1) \cdot \text{MSE}} \]

Consider the following phenotyping datasets where rows represent plots and columns represent scores given by different scorers for a trait, e.g. podcount per plant.

The high-correlation dataset is representative of phenotyping data obtained through the Artemis tool, while the low-correlation dataset represents data obtained manually. The datasets come from the same trial, so the groundtruth is the same. Assume the plots are all in one location, with negligible environmental variation. The plots contain different breeding materials in the same trial.

Dataset from Artemis (High Correlation)
scorer1 scorer2 scorer3
plot 1 20 21 20
plot 2 22 22 23
plot 3 21 22 21
plot 4 23 23 23
plot 5 22 21 22
plot 6 24 25 24
plot 7 25 24 25
plot 8 23 23 22
Dataset from Manual Scoring (Low Correlation)
scorer1 scorer2 scorer3
plot 1 20 15 17
plot 2 22 25 20
plot 3 18 19 22
plot 4 26 21 19
plot 5 21 24 26
plot 6 27 22 24
plot 7 23 28 21
plot 8 19 18 29

Calculating the Phenotyping Error

The first step is to compare the intra-class correlation coefficient (ICC) - or interrater reliability - between the two datasets. The ICC gives the proportion of total variance in the dataset (call this the phenotypic variance) that’s a result of true differences between plots. Assume the plots contain different breeding lines, so the difference in plots is attributable to genetic differences and some random noise in the environment (the ICC doesn’t distinguish between genetic and environmental variation, it just tells us how much variation comes from the scorer). If the ICC is high, it means most of the variation in the data comes from real differences between the plots (that’s what we want!), while a low ICC tells us there’s a lot of inconsistency between scorers. The ICC can range from 0 to 1.

We can calculate \(\widehat{ICC}\) using the above formula using only raw phenotyping datasets like the ones displayed here, no other inputs are needed.

Click to see step-by-step how to estimate ICC from the phenotyping data

First calculate the overall (grand) mean for each dataset:

\[ \bar{Y}_{\cdot\cdot} = \frac{1}{n \cdot k} \sum_{i=1}^n \sum_{j=1}^k Y_{ij} \]

where: \(n\) is the number of plots, \(k\) is the number of scorers, \(Y_{ij}\) is the score for plot \(i\) by scorer \(j\).

# Grand mean for high correlation dataset
grand_mean_high <- mean(as.matrix(data_high))

# Grand mean for low correlation dataset
grand_mean_low <- mean(as.matrix(data_low))

grand_mean_high
## [1] 22.54167
grand_mean_low
## [1] 21.91667

Next calculate the plot means for each dataset:

\[ \bar{Y}_{i\cdot} = \frac{1}{k} \sum_{j=1}^k Y_{ij} \]

# Plot means for high correlation dataset
plot_means_high <- rowMeans(data_high)

# Plot means for low correlation dataset
plot_means_low <- rowMeans(data_low)

plot_means_high
##   plot 1   plot 2   plot 3   plot 4   plot 5   plot 6   plot 7   plot 8 
## 20.33333 22.33333 21.33333 23.00000 21.66667 24.33333 24.66667 22.66667
plot_means_low
##   plot 1   plot 2   plot 3   plot 4   plot 5   plot 6   plot 7   plot 8 
## 17.33333 22.33333 19.66667 22.00000 23.66667 24.33333 24.00000 22.00000

Next for each dataset calculate the variability between plots (difference in plot means). This is the mean square between plots or MSB, a measure of how much of the total variation in scores can be attributed to true differences between plots, rather than noise or measurement error. A higher MSB means the plots are more distinct from each other.

\[ \text{MSB} = \frac{k}{n - 1} \sum_{i=1}^n (\bar{Y}_{i\cdot} - \bar{Y}_{\cdot\cdot})^2 \]

# Number of plots and scorers
n <- nrow(data_high)  # Number of plots
k <- ncol(data_high)  # Number of scorers

# MSB for high correlation dataset
MSB_high <- k * sum((plot_means_high - grand_mean_high)^2) / (n - 1)

# MSB for low correlation dataset
MSB_low <- k * sum((plot_means_low - grand_mean_low)^2) / (n - 1)

MSB_high
## [1] 6.470238
MSB_low
## [1] 16.92857

Next for each dataset calculate the variability of scores within plots, i.e. how much do scores for the same plot vary across scorers. This is the mean square error (MSE). A lower MSE means stronger agreement between scorers.

\[ \text{MSE} = \frac{1}{n \cdot (k - 1)} \sum_{i=1}^n \sum_{j=1}^k (Y_{ij} - \bar{Y}_{i\cdot})^2 \]

# MSE for high correlation dataset
MSE_high <- sum((as.matrix(data_high) - plot_means_high)^2) / (n * (k - 1))

# MSE for low correlation dataset
MSE_low <- sum((as.matrix(data_low) - plot_means_low)^2) / (n * (k - 1))

MSE_high
## [1] 0.2916667
MSE_low
## [1] 11.58333

Final estimation for the ICC of each dataset:

\[ \widehat{\text{ICC}} = \frac{\text{MSB} - \text{MSE}}{\text{MSB} + (k - 1) \cdot \text{MSE}} \]

# ICC for high correlation dataset
ICC_high <- (MSB_high - MSE_high) / (MSB_high + (k - 1) * MSE_high)

# ICC for low correlation dataset
ICC_low <- (MSB_low - MSE_low) / (MSB_low + (k - 1) * MSE_low)

ICC_high
## [1] 0.8759494
ICC_low
## [1] 0.1333135

Note that we can also do this with the psych package in R (I’ve just broken it down manually to show what’s going on). The output displays ICC calculated in a few ways, see https://medium.com/@SalahAssana/a-beginners-guide-to-the-intraclass-correlation-coefficient-icc-288f7fe7bcfc for details, but we want a “single raters” method which means the measure of reliability we care about is that an individual rater is as close as possible to the groundtruth. This matches the manual method worked through above.

## boundary (singular) fit: see help('isSingular')
## Call: ICC(x = data_high)
## 
## Intraclass correlation coefficients 
##                          type  ICC  F df1 df2       p lower bound upper bound
## Single_raters_absolute   ICC1 0.88 22   7  16 4.1e-07        0.66        0.97
## Single_random_raters     ICC2 0.88 22   7  14 1.5e-06        0.66        0.97
## Single_fixed_raters      ICC3 0.88 22   7  14 1.5e-06        0.65        0.97
## Average_raters_absolute ICC1k 0.95 22   7  16 4.1e-07        0.85        0.99
## Average_random_raters   ICC2k 0.95 22   7  14 1.5e-06        0.85        0.99
## Average_fixed_raters    ICC3k 0.95 22   7  14 1.5e-06        0.85        0.99
## 
##  Number of subjects = 8     Number of Judges =  3
## See the help file for a discussion of the other 4 McGraw and Wong estimates,
## boundary (singular) fit: see help('isSingular')
## Call: ICC(x = data_low)
## 
## Intraclass correlation coefficients 
##                          type  ICC   F df1 df2    p lower bound upper bound
## Single_raters_absolute   ICC1 0.13 1.5   7  16 0.25       -0.22        0.65
## Single_random_raters     ICC2 0.13 1.5   7  14 0.26       -0.23        0.65
## Single_fixed_raters      ICC3 0.13 1.5   7  14 0.26       -0.23        0.66
## Average_raters_absolute ICC1k 0.32 1.5   7  16 0.25       -1.20        0.85
## Average_random_raters   ICC2k 0.32 1.5   7  14 0.26       -1.25        0.85
## Average_fixed_raters    ICC3k 0.32 1.5   7  14 0.26       -1.31        0.85
## 
##  Number of subjects = 8     Number of Judges =  3
## See the help file for a discussion of the other 4 McGraw and Wong estimates,


So we have estimated the ICC from the high correlation (Artemis) and low correlation (manual) datasets as follows:

Estimated ICC for High and Low Correlation Datasets
Dataset ICC
High Correlation (Artemis) 0.8759494
Low Correlation (Manual) 0.1333135

From the ICC we can calculate the phenotyping error term, or residual variance, which is a key component of heritability. Reducing the residual variance increases the heritability measure mechanistically as we will see below.

We can think of the ICC as the portion of total variance in the phenotype data that’s explained by real groundtruth variance. Assume for now there is only one replication of each breeding line in the trial. So, if the total variance in the dataset (phenotypic variance, \(\sigma_P^2\)) is given by: \[ \sigma_P^2 = \sigma_{gt}^2 + \sigma_\epsilon^2 \] and the ICC is defined as: \[ \text{ICC} = \frac{\sigma_{gt}^2}{\sigma_P^2} \] Then \[ \sigma_{gt}^2 = \text{ICC} \cdot \sigma_P^2 \]

Substituting this into the equation for \(\sigma_P^2\) (or just intuitively - if ICC is the portion explained by what’s really there, then \(1-ICC\) is the portion explained by measurement error), we get:

\[ \sigma_\epsilon^2 = \sigma_P^2 \cdot (1 - \text{ICC}) \]

We can calculate \(\sigma_P^2\) from the raw data as simply the variance across the whole phenotyping dataset (separately for each dataset). Then from this we can calculate residual variance, \(\sigma_\epsilon^2\)) using the formula above. Note that we have still only used the raw datasets as inputs, no assumptions have been made.

Click to see code for calculating \(\sigma_P^2\), \(\sigma_{gt}^2\), and \(\sigma_\epsilon^2\) 
# Function to calculate the total phenotypic variance in each dataset
calculate_total_variance <- function(data) {
  sigma_P2 <- var(as.vector(as.matrix(data)))   # Flatten the data frame into a single vector of all scores
  return(sigma_P2)
}

sigma_P2_low <- calculate_total_variance(data_low)
sigma_P2_high <- calculate_total_variance(data_high)



# Function to calculate the groundtruth variance from ICC and sigma_P2
calculate_groundtruth_variance <- function(sigma_P2, ICC) {
sigma_gt2 <- sigma_P2 * ICC
  return(sigma_gt2)
}

sigma_gt2_high <- calculate_groundtruth_variance(sigma_P2_high, ICC_high)
sigma_gt2_low <- calculate_groundtruth_variance(sigma_P2_low, ICC_low)


# Function to calculate the phenotyping error, i.e. the residual error variance term
calculate_residual_variance <- function(sigma_gt2, ICC) {
  # sigma_epsilon2 <- sigma_P2 * (1 - ICC)
  sigma_epsilon2 <- sigma_gt2  / ICC * (1 - ICC) 

  return(sigma_epsilon2)
}


sigma_epsilon2_high <- calculate_residual_variance(sigma_gt2_high, ICC_high)
sigma_epsilon2_low <- calculate_residual_variance(sigma_gt2_low, ICC_low)


In the below table we see that the total phenotypic variance is larger for the low correlation dataset, this is due solely to the residual variance caused by phenotyping error. The groundtruth variance in the two datasets is essentially the same (slight difference due to noise). The residual variance term captures the portion of variance caused by scorer error, much larger in the manual low correlation dataset.

Variance components for high correlation (Artemis) and low correlation (manual) phenotyping datasets
Dataset Total Phenotypic Variance Groundtruth Variance Residual Variance
High Correlation (Artemis) 2.172101 1.902651 0.2694506
Low Correlation (Manual) 13.210145 1.761091 11.4490538


Heritability and Phenotyping Error

The heritability of a trait quantifies the proportion of the total phenotypic variance (\(\sigma_P^2\)) that is attributable to genetic variance (\(\sigma_g^2\)). It reflects how much of the observed variability in a trait is due to genetic differences rather than environmental or measurement error. This is what we really care about increasing, because increasing heritability increases the genetic gain we get from a breeding cycle. Reducing the phenotyping error (\(\sigma_\epsilon^2\)) mechanistically increases the heritability estimate because variance is removed from the error term, i.e. some of what appeared to be variance is now “explained” as simply human error in taking the phenotype measurement.

Whereas the ICC combines \(\sigma_g^2\) and \(\sigma_e^2\) to look at the portion of total phenotypic variance explained by groundtruth variance (\(\sigma_{gt}^2\)) (\(ICC = \frac{\sigma_{gt}^2}{\sigma_P^2}\)), heritability measures just the fraction explained by genetic variance (\(h^2 = \frac{\sigma_g^2}{\sigma_P^2}\)). If there is no environmental variance, or infinite replications, they are the same (more on this later.)

Expanding the components of \(\sigma_P^2\), the standard formula for heritability (\(h^2\)) is given by:

\[ h^2 = \frac{\sigma_g^2}{\sigma_g^2 + \frac{\sigma_e^2}{r} + \sigma_\epsilon^2} \]

where: \(\sigma_g^2\) is the true genetic variance, \(\sigma_e^2\) is the environmental variance (variance due to plot-to-plot environmental effects), \(\sigma_\epsilon^2\) is the residual variance (phenotyping error or measurement noise), and \(r\) is the number of replications per genotype. Note that \(\sigma_{gt}^2 = \sigma_g^2 + \frac{\sigma_e^2}{r}\), and \(\sigma_P^2 = \sigma_g^2 + \frac{\sigma_e^2}{r} + \sigma_\epsilon^2\).

So, we have seen that increasing the accuracy of phenotyping (ICC) leads to increased heritability estimates, but by how much exactly?

A key consideration is that the impact of improving scorer agreement depends on the proportion of total variance (or ground truth variance) attributable to environmental variance. Define:

\[\sigma_{gt}^2 = \sigma_g^2 + \sigma_e^2\]

so

\[\sigma_e^2 = \delta(\sigma_{gt}^2)\]

This means that the variance observed in the ground truth data is composed of genetic variance and variance caused by environmental effects, which reflects how plants express traits across environments. Here, \(\delta \in(0, 1)\) represents the proportion of ground truth variance attributable to environmental variance. Estimating \(\delta\) requires an informed assumption, and in other parts of the analysis we can think about the impact channels through which it can be affected by Artemis. Additionally, environmental variance can be reduced by increasing the number of replications (\(r\)) of the same trial.

When environmental variance is minimal (\(\delta\) close to 0) heritability approaches the ICC, and the relationship between them becomes approximately linear (\(ICC \approx h^2\)). However, as environmental variance increases, the relationship becomes asymptotic. This occurs because reducing phenotyping error (i.e., increasing ICC) yields diminishing returns to heritability when environmental variance constitutes a large share of total phenotypic variance. Consequently, under high environmental variance, further improvements in ICC at already high levels have a limited effect on heritability. This relationship is illustrated in the graph below.

Takeaway: the impact of increasing the ICC on heritability depends on the dominance of the term \(\frac{\sigma_e^2}{r}\), which Artemis can reduce through other impact channels.

Returning to the two example datasets, and making an assumption about the value of \(\delta\) and the number of trial replications \(r\), we can now calculate the genetic variance \(\sigma_g^2\), the environmental variance \(\sigma_e^2\), and heritability \(h^2\).

Click to see code for calculating \(\sigma_g^2\), \(\sigma_e^2\), and \(h^2\) 
replications <- 1 # Number of replications of the trial (eg if r = 2, assume the whole trial with n plots is duplicated)
delta <- 0.5

# Function to calculate the environmental variance in each dataset
calculate_environmental_variance <- function(delta, sigma_gt2) {
  sigma_e2 <- delta * sigma_gt2 
  return(sigma_e2)
}

sigma_e2_low <- calculate_environmental_variance(delta, sigma_gt2_low)
sigma_e2_high <- calculate_environmental_variance(delta, sigma_gt2_high)


# Function to calculate the genetic variance in each dataset
calculate_genetic_variance <- function(delta, sigma_gt2) {
  sigma_g2 <- (1 - delta) * sigma_gt2
  return(sigma_g2)
}

sigma_g2_low <- calculate_genetic_variance(delta, sigma_gt2_low)
sigma_g2_high <- calculate_genetic_variance(delta, sigma_gt2_high)


# Function to calculate the heritability in each dataset
calculate_heritability <- function(sigma_g2, sigma_e2, replications, sigma_epsilon2) {
  h2 <- sigma_g2 / (sigma_g2 + 2 * (sigma_e2 / replications) + sigma_epsilon2)
  return(h2)
}

h2_low <- calculate_heritability(sigma_g2_low, sigma_e2_low, replications, sigma_epsilon2_low)
h2_high <- calculate_heritability(sigma_g2_high, sigma_e2_high, replications, sigma_epsilon2_high)


“In the table below, we observe the impact of improved phenotyping (ICC increasing from 0.13 to 0.88) on heritability, assuming \(\delta = 0.5\) and one trial replication. No additional assumptions were made to generate this information (aside from creating the initial datasets, which should eventually be replaced with real data). As shown in the table, the environmental variance and genetic variance remain consistent across the two datasets (with only some random noise), indicating that the entire increase in heritability is due to reduced phenotyping error. It is important to note that the choice of \(\delta\) has significant implications for the effect of increasing scorer reliability. Additionally, this example assumes a very large improvement in scorer reliability, which is likely unrealistic for Artemis.

Variance components for high correlation (Artemis) and low correlation (manual) phenotyping datasets
Dataset Environmental Variance Genetic Variance Heritability
High Correlation (Artemis) 0.9513254 0.9513254 0.3045775
Low Correlation (Manual) 0.8805456 0.8805456 0.0624913


Simulation

To further explore the relationship between phenotyping accuracy and heritability we can simulate additional datasets similar to the high and low correlation examples used earlier, with varying ICCs (scorer agreement). To generate these datasets, we need to assume values for the ground truth variance in a trial (which includes genetic and environmental variance), the mean score (e.g., average pod count or any trait being simulated), the number of plots, and the number of scorers. We can create datasets for ICC values in 0.1 increments from 0 to 1, resulting in a total of 10 datasets. All the datasets should start out the same with the same level of groundtruth variance \(\sigma_{gt}^2\) and mean, and then we add random noise drawn from a normal distribution with a variance equal to the target residual error variance \(\sigma_\epsilon^2\) corresponding to each ICC increment.

Click to see code for simulating datasets 
icc_values <- seq(0.1, 1, by = 0.1) # ICC values for simulation
n_plots <- 100 # Number of plots in the trial, all with the same scoring consistency (ICC)
n_scorers <- 5 # Number of people scoring each plot
mean <- 20 # Average trait value for the trial/dataset
groundtruth_variance <- 1000 # Note that we need this to generate simulated data, but we would NOT need to make an assumption about it if we started with real datasets



# Function to generate simulated datasets
generate_dataset <- function(icc, n_plots, n_scorers, groundtruth_variance) {
  # Fixed genetic variance (true score variance)
  true_scores <- rnorm(n_plots, mean, sd = sqrt(groundtruth_variance))
  
  # Total phenotypic variance based on ICC
  total_phenotypic_variance <- groundtruth_variance / icc
  
  # Residual variance (phenotyping error variance)
  residual_variance <- total_phenotypic_variance * (1 - icc)
  
  # Add phenotyping error for each scorer
  scorer_data <- sapply(1:n_scorers, function(x) {
    true_scores + rnorm(n_plots, sd = sqrt(residual_variance))
  })
  
  return(as.data.frame(scorer_data))
}



# Function to generate and save simulated datasets for each level of ICC 
generate_and_save_datasets <- function(icc_values, n_plots, n_scorers, groundtruth_variance) {
  datasets <- lapply(icc_values, function(icc) {
    generate_dataset(icc, n_plots, n_scorers, groundtruth_variance)
  })
  
  # Assign ICC values as names for easy identification
  names(datasets) <- paste0("ICC_", icc_values)
  
  return(datasets)
}

# Generate and save datasets
set.seed(22) 
simulated_datasets <- generate_and_save_datasets(icc_values, n_plots, n_scorers, groundtruth_variance)


After simulating the datasets, we can calculate heritability for each one to examine how scorer correlation affects heritability. To do this, we need to estimate a few additional parameters: the proportion of total ground truth variance attributable to environmental variance (\(\delta\)) and the number of replications per genotype (\(r\)). Note that the same function used in the simulation can be applied to real datasets to calculate ICC, residual variance, and heritability. For instance, if we had datasets from Artemis and from a manual phenotyping method, we could directly compare heritability between them. While simulations provide an estimate of the theoretical impact of a change in ICC on heritability, real datasets allow for practical validation of these effects.

Click to see code for estimating \(h^2\),\(\sigma_\epsilon^2\), and \(\widehat{ICC}\) from the simulated datasets (this just combines all the same functions used throughout the first example) 
replications <- 1 # Number of replications of the trial (eg if r = 2, assume the whole trial with n plots is duplicated)
delta <- 0.5


calculate_heritability <- function(datasets, delta, replications) {
  results <- lapply(datasets, function(dataset) {
    
    # Calculate ICC using the psych package
    icc_estimate <- ICC(dataset)$results$ICC[1] # Assume one-way ICC model
    
    # Calculate total variance
    sigma_P2 <- var(as.vector(as.matrix(dataset)))
    
    # Calculate genetic and environmental variance
    sigma_g2 <- icc_estimate * sigma_P2 * (1 - delta)
    sigma_e2 <- icc_estimate * sigma_P2 * delta 
    
    # Calculate residual variance
    # sigma_epsilon2 <- sigma_P2 * (1 - icc_estimate)
    sigma_epsilon2 <- sigma_g2  / icc_estimate * (1 - icc_estimate) 

    # Calculate heritability
    h2 <- sigma_g2 / (sigma_g2 + (sigma_e2 / replications) + sigma_epsilon2)
    
    # Return h2, sigma_epsilon2, and the estimated ICC
    return(c(Heritability = h2, ResidualVariance = sigma_epsilon2, ICC = icc_estimate))
  })


  # Convert results to a data frame
  results_df <- as.data.frame(do.call(rbind, results))
  return(results_df)
}


# Calculate heritability from the simulated datasets
heritability_results <- calculate_heritability(simulated_datasets, delta, replications)

The table below shows the heritability, residual variance, and ICC (estimated from the data) under the assumptions of \(\delta = 0.5\) and 1 trial replication. We constructed a dataset for every 0.1 increment of ICC and added some noise, so the calculated ICC values closely follow these increments but do not align perfectly. Reducing \(\delta\) or increasing \(r\) increases the potential impact of improved phenotyping accuracy. The blue line in the first plot would increase linearly from 0 to 1, while the red line in the second plot would decrease asymptotically from 1.

Relationship between ICC, residual variance, and heritability in the simulated datasets
Heritability ResidualVariance ICC
0.135 4439.981 0.156
0.146 2112.331 0.171
0.221 1107.592 0.284
0.295 705.653 0.419
0.325 498.107 0.482
0.405 290.991 0.682
0.425 219.355 0.738
0.433 128.630 0.763
0.472 58.298 0.894
0.500 0.000 1.000

Next steps:

  • Tune the parameters

    • What is a realistic increase in phenotyping accuracy from Artemis?
    • Realistic estimate of \(\delta\)
    • Real or realistic phenotyping datasets from multiple scorers

  • Get a realistic estimate of the heritability gains

  • Set up the rest of the model to show how heritability gains translate into genetic gain on and off station, and eventually into economic gains through increased yields in farmer fields (economic model)

  • Work through other impact channels and combine to see full effects