Getting Started with VeraCrop

Overview

VeraCrop implements Comparative Performance Analysis (CPA) for yield gap estimation: it fits a regression model to field-level agronomic data, identifies which management/environmental variables significantly affect yield, and decomposes the gap between average and attainable yield into the contribution of each variable.

This vignette walks through the full workflow using wheat1, a real 200-field wheat trial dataset bundled with the package. See ?wheat1 for a full description of its columns, and ?wheat2 for a second, smaller and more heterogeneous example dataset used elsewhere in the package documentation to illustrate data-cleaning features such as trim_whitespace().

1. Load the example data

data(wheat1)
str(wheat1)
#> 'data.frame':    200 obs. of  13 variables:
#>  $ Yield              : int  4961 5853 5365 6643 6700 6310 3960 4826 5472 5175 ...
#>  $ Irrigation         : int  315 515 364 553 576 218 411 557 421 383 ...
#>  $ Nitrogen           : int  60 241 150 129 101 220 91 72 43 43 ...
#>  $ Phosphorus         : int  148 21 136 86 59 67 106 12 51 102 ...
#>  $ Potassium          : int  47 137 45 64 35 160 29 165 66 75 ...
#>  $ Soil_pH            : num  6.9 6.6 5.9 5.6 6.3 8.4 7 6.9 7.7 7.5 ...
#>  $ Soil_Organic_Matter: num  1.3 2.3 1 3.1 3 1.9 2.8 1.4 0.7 1.8 ...
#>  $ Plant_density      : int  358 366 247 315 246 397 119 147 336 263 ...
#>  $ Weed_Infestation   : int  14 14 4 30 15 45 33 57 5 57 ...
#>  $ Drought            : int  0 0 0 1 1 0 1 0 0 0 ...
#>  $ Pest               : int  0 1 0 0 0 0 1 0 1 0 ...
#>  $ Cultivar           : chr  "Azar" "Azar" "Azar" "Sardari" ...
#>  $ Sowing_Date        : int  33 29 32 1 7 24 3 28 23 11 ...

2. Preprocess the data

prep_yield_gap() handles variable-type detection, missing-value imputation, encoding of categorical variables, constant/near-zero-variance filtering, and scaling of continuous predictors - all in one call. The response variable (and any binary/dummy columns produced by encoding) are automatically protected from being standardized, since scaling a 0/1 indicator has no meaningful interpretation.

prep <- prep_yield_gap(
  wheat1,
  response_var = "Yield",
  verbose = FALSE
)

str(prep$data)
#> 'data.frame':    200 obs. of  14 variables:
#>  $ Yield                     : int  4961 5853 5365 6643 6700 6310 3960 4826 5472 5175 ...
#>  $ Drought                   : num  0 0 0 1 1 0 1 0 0 0 ...
#>  $ Pest                      : num  0 1 0 0 0 0 1 0 1 0 ...
#>  $ Cultivar_Sardari          : num  0 0 0 1 0 1 0 0 0 0 ...
#>  $ Cultivar_Marvdasht        : num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Irrigation_scaled         : num  -0.799 1.026 -0.352 1.372 1.582 ...
#>  $ Nitrogen_scaled           : num  -0.849 1.6184 0.3779 0.0916 -0.2901 ...
#>  $ Phosphorus_scaled         : num  1.662 -1.241 1.388 0.245 -0.373 ...
#>  $ Potassium_scaled          : num  -0.886 0.657 -0.92 -0.595 -1.092 ...
#>  $ Soil_pH_scaled            : num  -0.1 -0.445 -1.251 -1.596 -0.791 ...
#>  $ Soil_Organic_Matter_scaled: num  -0.808 0.416 -1.175 1.395 1.273 ...
#>  $ Plant_density_scaled      : num  1.3508 1.4439 0.0581 0.85 0.0464 ...
#>  $ Weed_Infestation_scaled   : num  -0.9588 -0.9588 -1.5348 -0.0372 -0.9012 ...
#>  $ Sowing_Date_scaled        : num  1.118 0.782 1.034 -1.573 -1.069 ...

Check that preprocessing produced a usable, well-formed dataset before moving on:

val <- validate_preprocessing(prep, verbose = FALSE)
val$valid
#> [1] TRUE

3. Fit a model and check regression assumptions

model <- lm(Yield ~ ., data = prep$data)
summary(model)$r.squared
#> [1] 0.7605189

diag <- check_assumptions(model, data = prep$data, plot = FALSE, verbose = FALSE)
diag$data_ready
#> [1] TRUE
diag$data_ready_reason
#> [1] "1 warning(s) (<= 2 allowed)"

4. Run the full yield gap analysis

yield_gap_analysis() performs stepwise variable selection, k-fold cross-validation, optimal-value estimation for each significant variable, min-max effect sizes, relative importance, and yield gap decomposition.

Two variable-selection methods are available: "ftest" (the default - a fixed F-to-Enter/F-to-Remove or p-value threshold, similar to software such as SigmaPlot) and "aic" (stepwise selection minimizing AIC via stats::step()). They encode different criteria for what counts as a useful predictor and will often select a different number of variables - neither is universally “correct”.

result <- yield_gap_analysis(prep$data, response = "Yield", verbose = FALSE)

result$metrics
#> $R2
#> [1] 0.5748787
#> 
#> $RMSE
#> [1] 755.5695
#> 
#> $MAE
#> [1] 616.5646
result$yield_mean
#> [1] 5079.56
result$yield_opt
#> [1] 6781.248
result$yield_gap
#> [1] 1701.688

The cpa_table shows, for each significant variable, its coefficient, average value, estimated optimal value, and share of the total yield gap:

result$cpa_table
#>                                              Variable      Beta         Mean
#> (Intercept)                                 Intercept 5042.5149 1.000000e+00
#> Sowing_Date_scaled                 Sowing_Date_scaled -479.1337 1.421048e-16
#> Cultivar_Sardari                     Cultivar_Sardari  985.6168 3.050000e-01
#> Irrigation_scaled                   Irrigation_scaled  378.9352 2.710519e-17
#> Drought                                       Drought -775.2001 3.400000e-01
#> Soil_Organic_Matter_scaled Soil_Organic_Matter_scaled  307.0349 5.740186e-17
#>                                   Opt Contribution_Mean Contribution_Opt
#> (Intercept)                 1.0000000      5.042515e+03        5042.5149
#> Sowing_Date_scaled         -0.6060086     -6.808718e-14         290.3592
#> Cultivar_Sardari            1.0000000      3.006131e+02         985.6168
#> Irrigation_scaled           0.6967774      1.027111e-14         264.0335
#> Drought                     0.0000000     -2.635680e+02           0.0000
#> Soil_Organic_Matter_scaled  0.6472338      1.762438e-14         198.7234
#>                            Gap_Component Share_Percent Share_Percent_Abs
#> (Intercept)                       0.0000            NA                NA
#> Sowing_Date_scaled              290.3592      17.06301          17.06301
#> Cultivar_Sardari                685.0037      40.25437          40.25437
#> Irrigation_scaled               264.0335      15.51598          15.51598
#> Drought                         263.5680      15.48863          15.48863
#> Soil_Organic_Matter_scaled      198.7234      11.67802          11.67802
#>                            Share_Cumulative MinMaxEffect MinMaxPct
#> (Intercept)                              NA           NA        NA
#> Sowing_Date_scaled                 57.31738   -1452.8033  27.05396
#> Cultivar_Sardari                   40.25437     985.6168  18.35406
#> Irrigation_scaled                  72.83336    1179.4826  21.96421
#> Drought                            88.32198    -775.2001  14.43570
#> Soil_Organic_Matter_scaled        100.00000     976.9175  18.19206

Field-level (per-observation) yield gap

In addition to the overall yield gap above, VeraCrop also estimates the gap for each individual field, available in result$predictions:

head(result$predictions)
#>   Observed Predicted Yield_Gap_Obs
#> 1     4961  3956.044      2825.204
#> 2     5853  5184.353      1596.895
#> 3     5365  4053.008      2728.239
#> 4     6643  6955.111         0.000
#> 5     6700  5769.631      1011.617
#> 6     6310  5194.484      1586.764

save_field_level_excel() exports a full per-field breakdown, including each variable’s individual contribution to that field’s gap.

5. Visualize the results

plot_contribution_shares(result)

plot_observed_vs_predicted(result)
#> `geom_smooth()` using formula = 'y ~ x'

All five plots can be generated together with plot_all_graphs(result).

6. Export results

tmp <- tempfile(fileext = ".xlsx")
save_results_excel(result, output_path = tmp, verbose = FALSE)
file.exists(tmp)
#> [1] TRUE

A messier, real-world dataset

wheat2 (40 fields, 29 recorded variables) illustrates several data-quality and small-sample issues that prep_yield_gap() handles automatically:

data(wheat2)

# "Sirvan" and "Sirvan " (trailing space) are automatically merged
length(unique(wheat2$Cultivar))
#> [1] 4
trimmed <- trim_whitespace(wheat2, verbose = FALSE)
length(unique(trimmed$Cultivar))
#> [1] 3

Litracy_Level (farmer education) is a genuinely ordinal variable stored as text. Without an explicit order, it is conservatively treated as nominal; supplying ordinal_level_orders encodes it correctly:

prep_small <- prep_yield_gap(
  wheat2,
  response_var = "Yield",
  ordinal_level_orders = list(
    Litracy_Level = c("Elementary School", "Middle School", "Diploma",
                       "Bachelor's Degree", "Master's Degree")
  ),
  verbose = FALSE
)

Next steps

See ?prep_yield_gap for all preprocessing options (imputation method, encoding thresholds, near-zero-variance / correlation filtering).
See ?check_assumptions for the full set of regression diagnostics.
See ?yield_gap_analysis for variable selection and cross-validation options.
See ?wheat1, ?wheat2, and ?sugarcane1 for details on the bundled example datasets - the latter also illustrates a genuinely ordinal variable (crop_type) and a perfectly collinear (“aliased”) derived column (total_urea_kg_ha), both handled gracefully by the package.