The Role of Data Science in Agriculture Sector

Introduction

Data science enables better understanding and management of agricultural price volatility by analyzing trends, predicting fluctuations, and providing actionable insights. This study specifically focuses on leveraging predictive models to assist policymakers and farmers in making informed decisions to optimize production strategies and mitigate risks.

CRISP-DM Framework Analysis

1. Business Understanding

Volatile agricultural prices directly impact farmers’ income and decision-making processes. This study aims to:
1. Predict future agriculture PPI values using time series models.
2. Use other industries’ PPI data to predict agriculture PPI using regression models.

These objectives address critical industry needs, such as stabilizing market conditions and improving income predictability for stakeholders.

2. Data Understanding

Dataset Source: Department of Statistics Malaysia
Fields:
- date: Monthly timestamps
- agriculture: Agriculture PPI
- mining: Mining PPI
- manufacturing: Manufacturing PPI
- electricity: Electricity & Gas Supply PPI
- water: Water Supply PPI

library(dplyr)
library(ggplot2)
library(reshape2)

# Load dataset
data <- read.csv("ppi.csv")

# Summary of the dataset
summary(data)

##      date              overall       agriculture         mining     
##  Length:179         Min.   : 98.4   Min.   : 87.30   Min.   : 44.7  
##  Class :character   1st Qu.:104.5   1st Qu.: 95.05   1st Qu.: 91.4  
##  Mode  :character   Median :108.2   Median :107.90   Median : 99.8  
##                     Mean   :109.2   Mean   :111.77   Mean   :104.7  
##                     3rd Qu.:113.5   3rd Qu.:123.55   3rd Qu.:121.2  
##                     Max.   :124.1   Max.   :178.70   Max.   :157.2  
##  manufacturing    electricity        water      
##  Min.   : 98.8   Min.   : 99.4   Min.   : 99.1  
##  1st Qu.:105.8   1st Qu.:108.4   1st Qu.:106.7  
##  Median :107.6   Median :115.8   Median :112.1  
##  Mean   :109.4   Mean   :113.3   Mean   :111.3  
##  3rd Qu.:110.5   3rd Qu.:117.7   3rd Qu.:114.3  
##  Max.   :122.8   Max.   :123.3   Max.   :128.6

# Check for missing values
cat("Missing values in dataset:", sum(is.na(data)), "\n")

## Missing values in dataset: 0

# Ensure agriculture data is numeric and handle missing values
data$agriculture <- as.numeric(as.character(data$agriculture))
data <- data %>% filter(!is.na(agriculture))
cat("Total rows after cleaning:", nrow(data), "\n")

## Total rows after cleaning: 179

Data Visualizations

# Plot Agriculture PPI over time
ggplot(data, aes(x = as.Date(date), y = agriculture)) +
  geom_line(color = "#2E8B57", size = 1.2) +  
  ggtitle("Agriculture PPI Over Time from 2010 to 2024") +
  xlab("Date") +
  ylab("Agriculture PPI") +
  theme_minimal(base_size = 14) +            
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5), 
    axis.title = element_text(size = 14, face = "bold"),
    axis.text = element_text(size = 12),
    panel.grid.major = element_line(color = "gray80", linetype = "dashed")
  )

# Plot correlation between Agriculture PPI and Mining PPI
correlation_matrix <- cor(
  data[, c("agriculture", "mining", "manufacturing", "electricity", "water")], 
  use = "complete.obs"
)
correlation_data <- melt(correlation_matrix)
ggplot(correlation_data, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile(color = "white") +  # Add white borders for tiles
  scale_fill_gradient2(
    low = "#2b83ba", high = "#d7191c", mid = "#f7f7f7", 
    midpoint = 0, limit = c(-1, 1), name = "Correlation"
  ) +  # Improved color scheme
  ggtitle("Heatmap of Correlations with Agriculture PPI") +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 12, face = "bold"),
    axis.text.y = element_text(size = 12, face = "bold"),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10)
  )

# Plot distribution of Agriculture PPI
ggplot(data, aes(x = agriculture)) +
  geom_histogram(binwidth = 5, fill = "#FF8C00", color = "white", alpha = 0.9) +  # Soft orange with white borders
  ggtitle("Distribution of Agriculture PPI") +
  xlab("Agriculture PPI") +
  ylab("Frequency") +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),  # Centered and bold title
    axis.title = element_text(size = 14, face = "bold"),
    axis.text = element_text(size = 12),
    panel.grid.major = element_line(color = "gray80", linetype = "dashed")
  )

3. Data Preparation

Prepare the dataset for both objectives.

3.1 Time Series Data Preparation

Prepare the agriculture PPI data for time series analysis.

# Create a time series object for agriculture PPI
agriculture_ts <- ts(data$agriculture, start = c(2010, 1), frequency = 12)

3.2 Regression Data Preparation

Prepare the dataset to predict agriculture PPI based on other industries’ PPI data.

# Select relevant columns for regression
data_regression <- data %>%
  select(agriculture, mining, manufacturing, electricity, water) %>%
  na.omit()

# Split data into training and testing sets
set.seed(123)
library(caret)
trainIndex <- createDataPartition(data_regression$agriculture, p = 0.8, list = FALSE)
trainData <- data_regression[trainIndex, ]
testData <- data_regression[-trainIndex, ]

3.3 Feature Selection

Analyze feature importance using Random Forest and select important features.

library(randomForest)

# Train Random Forest model to calculate feature importance
rf_model_importance <- randomForest(
  agriculture ~ ., 
  data = trainData, 
  ntree = 500, 
  mtry = 2, 
  importance = TRUE
)

# Print feature importance
importance <- randomForest::importance(rf_model_importance)
print(importance)

##                %IncMSE IncNodePurity
## mining        21.44922      5756.760
## manufacturing 45.61522     26150.247
## electricity   22.62519      7586.856
## water         28.49505     11854.180

# Visualize feature importance
varImpPlot(rf_model_importance)

# Select important features based on importance scores
selected_features <- c("agriculture", "mining", "manufacturing", "electricity")

# Update training and testing datasets with selected features
trainData_selected <- trainData[, selected_features]
testData_selected <- testData[, selected_features]

4. Modeling

4.1 Time Series Modeling: Multi-Model

Encapsulate ARIMA and ETS models into reusable functions.

library(forecast)

# Function to train ARIMA model
ARIMA_Model <- function(ts_data) {
  model <- auto.arima(ts_data)
  forecast_data <- forecast(model, h = 12)
  rmse <- sqrt(mean((ts_data - fitted(model))^2))
  mae <- mean(abs(ts_data - fitted(model)))
  mape <- mean(abs((ts_data - fitted(model)) / ts_data)) * 100
  r_squared <- 1 - sum((ts_data - fitted(model))^2) / sum((ts_data - mean(ts_data))^2)
  return(list(model = model, forecast = forecast_data, RMSE = rmse,MAE = mae, MAPE = mape, R2 = r_squared))
}

# Function to train ETS model
ETS_Model <- function(ts_data) {
  model <- ets(ts_data)
  forecast_data <- forecast(model, h = 12)
  rmse <- sqrt(mean((ts_data - fitted(model))^2))
  mae <- mean(abs(ts_data - fitted(model)))
  mape <- mean(abs((ts_data - fitted(model)) / ts_data)) * 100
  r_squared <- 1 - sum((ts_data - fitted(model))^2) / sum((ts_data - mean(ts_data))^2)
  return(list(model = model, forecast = forecast_data, RMSE = rmse,MAE = mae, MAPE = mape, R2 = r_squared))
}

# Apply ARIMA and ETS models to the time series
arima_results <- ARIMA_Model(agriculture_ts)
ets_results <- ETS_Model(agriculture_ts)

cat("ARIMA RMSE:", arima_results$RMSE,"\n","ARIMA MAE:",arima_results$MAE,"\n","ARIMA MAPE:", arima_results$MAPE,"%\n","ARIMA R²:",arima_results$R2,"\n")

## ARIMA RMSE: 5.21153 
##  ARIMA MAE: 3.573184 
##  ARIMA MAPE: 3.10699 %
##  ARIMA R²: 0.9269175

cat("ETS RMSE:", ets_results$RMSE,"\n","ETS MAE:",ets_results$MAE,"\n","ETS MAPE:", ets_results$MAPE,"%\n","ETS R²:",ets_results$R2,"\n")

## ETS RMSE: 5.355325 
##  ETS MAE: 3.611136 
##  ETS MAPE: 3.139145 %
##  ETS R²: 0.922829

4.2 Regression Modeling: Multi-Model

Encapsulate Linear Regression, Random Forest, and XGBoost into reusable functions.

library(randomForest)
library(caret)
library(xgboost)
library(class)

# Function to train Linear Regression model
Linear_Regression <- function(train_data, test_data) {
  model <- lm(agriculture ~ ., data = train_data)
  predictions <- predict(model, test_data)
  actual <- test_data$agriculture
  
  rmse <- sqrt(mean((actual - predictions)^2))
  mae <- mean(abs(actual - predictions))
  r2 <- 1 - sum((actual - predictions)^2) / sum((actual - mean(actual))^2)
  mape <- mean(abs((actual - predictions) / actual)) * 100

  return(list(model = model, predictions = predictions, RMSE = rmse, MAE = mae, R2 = r2, MAPE = mape))
}

# Function to train Random Forest model
RF_Model <- function(train_data, test_data, ntree = 1000, mtry = NULL, nodesize = 5) {
  if (is.null(mtry)) {
    mtry <- floor(sqrt(ncol(train_data) - 1))
  }
  model <- randomForest(
    agriculture ~ ., 
    data = train_data, 
    ntree = ntree, 
    mtry = mtry, 
    nodesize = nodesize
  )
  predictions <- predict(model, test_data)
  actual <- test_data$agriculture

  rmse <- sqrt(mean((actual - predictions)^2))
  mae <- mean(abs(actual - predictions))
  r2 <- 1 - sum((actual - predictions)^2) / sum((actual - mean(actual))^2)
  mape <- mean(abs((actual - predictions) / actual)) * 100

  return(list(model = model, predictions = predictions, RMSE = rmse, MAE = mae, R2 = r2, MAPE = mape))
}

# Function to train XGBoost model
XGBoost_Model <- function(train_data, test_data, nrounds = 1000, max_depth = 6, eta = 0.3) {
  # Prepare data for XGBoost
  dtrain <- xgb.DMatrix(data = as.matrix(train_data[, -1]), label = train_data$agriculture)
  dtest <- xgb.DMatrix(data = as.matrix(test_data[, -1]), label = test_data$agriculture)

  # Set XGBoost parameters
  params <- list(
    objective = "reg:squarederror",
    max_depth = max_depth,
    eta = eta
  )

  # Train the XGBoost model
  model <- xgb.train(params = params, data = dtrain, nrounds = nrounds)

  # Predict and calculate metrics
  predictions <- predict(model, dtest)
  actual <- test_data$agriculture

  rmse <- sqrt(mean((actual - predictions)^2))
  mae <- mean(abs(actual - predictions))
  r2 <- 1 - sum((actual - predictions)^2) / sum((actual - mean(actual))^2)
  mape <- mean(abs((actual - predictions) / actual)) * 100

  return(list(model = model, predictions = predictions, RMSE = rmse, MAE = mae, R2 = r2, MAPE = mape))
}

# Function to train KNN Regression model
KNN_Model <- function(data_train, data_test, k = 5) {
  model <- train(agriculture ~ ., data = data_train, method = "knn", tuneGrid = expand.grid(k = k))
  x_test <- data_test[, -1]
  y_test <- data_test$agriculture
  predictions <- predict(model, x_test)

  rmse <- sqrt(mean((y_test - predictions)^2))
  mae <- mean(abs(y_test - predictions))
  r2 <- 1 - sum((y_test - predictions)^2) / sum((y_test - mean(y_test))^2)
  mape <- mean(abs((y_test - predictions) / y_test)) * 100

  return(list(model = model, RMSE = rmse, MAE = mae, R2 = r2, MAPE = mape))
}

# Apply Linear Regression, Random Forest, and XGBoost models
linear_results <- Linear_Regression(trainData_selected, testData_selected)
rf_results <- RF_Model(trainData_selected, testData_selected)
xgboost_results <- XGBoost_Model(trainData_selected, testData_selected, nrounds = 2000, max_depth = 4, eta = 0.01)
knn_results <- KNN_Model(trainData_selected, testData_selected, k = 10)

# Print metrics for all models
cat("Linear Regression Metrics:\n","RMSE:", linear_results$RMSE, "\n","MAE:", linear_results$MAE, "\n","R2:", linear_results$R2, "\n","MAPE:", linear_results$MAPE, "%\n")

## Linear Regression Metrics:
##  RMSE: 12.8489 
##  MAE: 10.03218 
##  R2: 0.5504599 
##  MAPE: 8.716094 %

cat("Random Forest Metrics:\n","RMSE:", rf_results$RMSE, "\n","MAE:", rf_results$MAE, "\n","R2:", rf_results$R2, "\n","MAPE:", rf_results$MAPE, "%\n")

## Random Forest Metrics:
##  RMSE: 10.94292 
##  MAE: 7.2857 
##  R2: 0.6739362 
##  MAPE: 6.254914 %

cat("XGBoost Metrics:\n","RMSE:", xgboost_results$RMSE, "\n","MAE:", xgboost_results$MAE, "\n","R2:", xgboost_results$R2, "\n","MAPE:", xgboost_results$MAPE, "%\n")

## XGBoost Metrics:
##  RMSE: 11.50995 
##  MAE: 6.902211 
##  R2: 0.639269 
##  MAPE: 5.991862 %

cat("KNN Metrics:\n","RMSE:", knn_results$RMSE, "\n","MAE:", knn_results$MAE, "\n","R2:", knn_results$R2, "\n","MAPE:", knn_results$MAPE, "%\n")

## KNN Metrics:
##  RMSE: 11.20911 
##  MAE: 8.794 
##  R2: 0.6578802 
##  MAPE: 7.905431 %

5. Evaluation

5.1 Time Series Model Selection

Evaluate and display results for each time series model, then select the best model.

if (arima_results$RMSE < ets_results$RMSE) {
  best_ts_model <- "ARIMA"
  best_ts_forecast <- arima_results$forecast
} else {
  best_ts_model <- "ETS"
  best_ts_forecast <- ets_results$forecast
}

cat("Best Time Series Model:", best_ts_model, "\n")

## Best Time Series Model: ARIMA

5.2 Regression Model Selection

Evaluate and display results for each regression model, then select the best model.

### Regression Model Selection
if (linear_results$RMSE < rf_results$RMSE & linear_results$RMSE < xgboost_results$RMSE & linear_results$RMSE < knn_results$RMSE) {
  best_reg_model <- "Linear Regression"
  best_reg_predictions <- linear_results$predictions
} else if (rf_results$RMSE < xgboost_results$RMSE & rf_results$RMSE < knn_results$RMSE) {
  best_reg_model <- "Random Forest"
  best_reg_predictions <- rf_results$predictions
} else if (xgboost_results$RMSE < knn_results$RMSE) {
  best_reg_model <- "XGBoost"
  best_reg_predictions <- xgboost_results$predictions
} else {
  best_reg_model <- "KNN"
  best_reg_predictions <- knn_results$predictions
}

cat("Best Regression Model:", best_reg_model, "\n")

## Best Regression Model: Random Forest

6. Deployment

Visualize the insights through trend graphs, regression results, and innovative predictions. Shinny app: Agriculture_Sector_Analysis

Conclusion

By applying the CRISP-DM framework:
- Evaluated ARIMA and ETS models for time series forecasting and selected the best-performing model.
- Evaluated Linear Regression, Random Forest, XGBoost models and K-Nearest Neighbors for predicting agriculture PPI and selected the best-performing model.
- Conducted sensitivity and uncertainty analyses to provide deeper insights into agricultural price trends.
- Developed an alert system to flag significant price volatility, aiding decision-makers.

These models and analyses offer actionable insights for stabilizing agricultural markets, optimizing production, and supporting policy formulation.

References

Dataset: PPI Data CSV
Packages: dplyr, ggplot2, caret, forecast, randomForest, xgboost, class

Agriculture Sector Analysis using CRISP-DM Framework

Wong Yi Ting (S2152880)

2025-01-05