Getting Started with VeraCrop

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.