Predictive mapping of wholesale grain prices for rural areas in Tanzania
Madaga Lavinia, Jordan Chamberlin
10 November, 2025
Suggested citation: Madaga, Lavinia, Jordan Chamberlin, Bisrat Gebrekidan, Nicholaus Kuboja, Maxwell Mkondiwa, and Silva Joao Vasco, 2025. “Predictive mapping of wholesale crop prices for rural areas in Tanzania.” https://github.com/EiA2030-ex-ante/Tanzania-Price-Data
Introduction
We are interested in estimating the prices of agricultural food commodities across space and time, on the basis of prices as observed at some market locations. Here we explore methods to model these prices, using monthly data from across Tanzania over the period May 2021-July 2025.
This document describes the process of fitting a Random Forest model to predict these prices. The performance of the Random Forest model is evaluated using Root Mean Square Errors (RMSE) and R-Square as test statistics. We also compare the observed prices with the predicted prices using an out of sample dataset.
Load Libraries
library(geodata)## Warning: package 'geodata' was built under R version 4.3.3library(lubridate)
library(terra)
library(data.table)## Warning: package 'data.table' was built under R version 4.3.3library(randomForest)## Warning: package 'randomForest' was built under R version 4.3.3library(httr)
library(caret)## Warning: package 'caret' was built under R version 4.3.3## Warning: package 'ggplot2' was built under R version 4.3.3## Warning: package 'lattice' was built under R version 4.3.3library(Metrics)## Warning: package 'Metrics' was built under R version 4.3.3library(pdp)## Warning: package 'pdp' was built under R version 4.3.3library(gridExtra)
library(stats)
library(dplyr)
library(stringr)
library(fields)## Warning: package 'fields' was built under R version 4.3.3## Warning: package 'spam' was built under R version 4.3.3library(corrplot)## Warning: package 'corrplot' was built under R version 4.3.3library(ggplot2)
library(tidyr)## Warning: package 'tidyr' was built under R version 4.3.3library(gdistance)## Warning: package 'gdistance' was built under R version 4.3.3## Warning: package 'igraph' was built under R version 4.3.3library(sp)
library(raster)
library(foreign)
library(gdistance)
library(foreach)## Warning: package 'foreach' was built under R version 4.3.3library(doParallel)## Warning: package 'doParallel' was built under R version 4.3.3## Warning: package 'iterators' was built under R version 4.3.3library(tidyverse)
library(segmented)## Warning: package 'segmented' was built under R version 4.3.3library(modelbased)## Warning: package 'modelbased' was built under R version 4.3.3library(parameters)## Warning: package 'parameters' was built under R version 4.3.3library(RColorBrewer)
library(effsize)## Warning: package 'effsize' was built under R version 4.3.3library(gghalves)## Warning: package 'gghalves' was built under R version 4.3.3library(viridis)Data
This dataset contains price information for various crops across different regions and markets in Tanzania from 2021 to 2024. The data was acquired from Tanzania’s Ministry of Industry and Trade.
setwd("H:/Tanzania Price data/Datasets")
prices <- fread("Tanzania_Price_Data_AllCrops_with_Coordinates4.csv")
dim(prices)## [1] 11651 21head(prices)## Region Market Maize..min.price. Maize..max.price. Rice..min.price.
## <char> <char> <int> <int> <int>
## 1: Arusha Arusha 43000 45000 140000
## 2: Dar es Salaam Temeke 46000 47000 120000
## 3: Dodoma Majengo 45000 50000 132000
## 4: Dodoma Kibaigwa 30000 42000 NA
## 5: Kagera Bukoba 33000 35000 110000
## 6: Manyara Babati 30000 45000 120000
## Rice..max.price. Sorghum..min.price. Sorghum..max.price.
## <int> <int> <int>
## 1: 170000 60000 65000
## 2: 210000 80000 100000
## 3: 200000 50000 60000
## 4: NA 45000 48000
## 5: 140000 140000 150000
## 6: 170000 60000 80000
## Bulrush.Millet..min.price. Bulrush.Millet..max.price.
## <int> <int>
## 1: 70000 75000
## 2: 80000 100000
## 3: 52000 64000
## 4: NA NA
## 5: 120000 140000
## 6: 80000 100000
## Finger.Millet..min.price. Finger.Millet..max.price. Wheat..min.price.
## <int> <char> <int>
## 1: 120000 125000 70000
## 2: 175000 180000 110000
## 3: 200000 252000 NA
## 4: NA <NA> NA
## 5: 170000 180000 170000
## 6: 120000 130000 100000
## Wheat..max.price. Beans..min.price. Beans..max.price.
## <int> <int> <int>
## 1: 78000 130000 165000
## 2: 120000 220000 260000
## 3: NA 200000 245000
## 4: NA NA NA
## 5: 180000 110000 170000
## 6: 120000 150000 180000
## Irish.Potatoes..min.price. Irish.Potatoes..max.price. Date Latitude
## <int> <int> <char> <num>
## 1: 70000 75000 5/5/2021 -3.36968
## 2: 75000 75000 5/5/2021 -6.85097
## 3: 56000 68000 5/5/2021 -6.17971
## 4: NA NA 5/5/2021 -6.08110
## 5: 60000 75000 5/5/2021 -1.33140
## 6: 60000 60000 5/5/2021 -4.21006
## Longitude
## <num>
## 1: 36.68808
## 2: 39.25672
## 3: 35.74109
## 4: 36.64645
## 5: 31.81293
## 6: 35.74915table(prices$Market)##
## Arusha Babati Bariadi Buguruni Bukoba Igawilo
## 374 365 200 10 511 160
## Ilala Iringa Kibaigwa Kigoma Kilombero Kinondoni
## 325 514 211 369 160 203
## Lindi majengo Majengo Makambako Manyara Mbeya
## 368 61 621 80 116 117
## Mgandini Mnalani Morogoro Moshi Mpanda Mpimbwe
## 160 111 508 374 195 39
## Mtwara Musoma Mwananyamala Mwanza Namfua Njombe
## 522 409 160 471 160 174
## Nyankumbu Pwani Shinyanga Singida Songea Songwe
## 198 89 404 83 483 170
## Sumbawanga Tabora Tandale Tandika Tanga Temeke
## 507 486 177 186 361 291
## Ubungo
## 168sapply(prices, class)## Region Market
## "character" "character"
## Maize..min.price. Maize..max.price.
## "integer" "integer"
## Rice..min.price. Rice..max.price.
## "integer" "integer"
## Sorghum..min.price. Sorghum..max.price.
## "integer" "integer"
## Bulrush.Millet..min.price. Bulrush.Millet..max.price.
## "integer" "integer"
## Finger.Millet..min.price. Finger.Millet..max.price.
## "integer" "character"
## Wheat..min.price. Wheat..max.price.
## "integer" "integer"
## Beans..min.price. Beans..max.price.
## "integer" "integer"
## Irish.Potatoes..min.price. Irish.Potatoes..max.price.
## "integer" "integer"
## Date Latitude
## "character" "numeric"
## Longitude
## "numeric"# Convert to date format
prices$Date <- lubridate::mdy(prices$Date)Basic Data preperation
Check the Region and Market names and coodinates. Make sure the Region and Market names are harmonized and properly geocoded
unique(prices[Region=="Arusha",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Arusha -3.36968 36.68808
## 2: Kilombero -3.37324 36.67870unique(prices[Region=="Dar es Salaam",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Temeke -6.850970 39.25672
## 2: Kinondoni -6.784070 39.27007
## 3: Ilala -6.829410 39.26289
## 4: Tandika -6.867912 39.25480
## 5: Buguruni -6.838380 39.24359
## 6: Tandale -6.795230 39.24085
## 7: Ubungo -6.793620 39.20966
## 8: Mwananyamala -6.788890 39.25840unique(prices[Region=="Dodoma",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Majengo -6.17971 35.74109
## 2: Kibaigwa -6.08110 36.64645unique(prices[Region=="Kagera",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Bukoba -1.3314 31.81293unique(prices[Region=="Manyara",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Babati -4.21006 35.74915
## 2: Manyara -4.46011 37.20217unique(prices[Region=="Rukwa",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Sumbawanga -7.95239 31.62319unique(prices[Region=="Mpanda",.(Market, Latitude, Longitude)])## Empty data.table (0 rows and 3 cols): Market,Latitude,Longitudeunique(prices[Region=="Mtwara",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Mtwara -10.28009 40.18191unique(prices[Region=="Tabora",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Tabora -5.01972 32.80767unique(prices[Region=="Tanga",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Tanga -5.074260 39.09993
## 2: Mgandini -5.086235 39.09561unique(prices[Region=="Iringa",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Iringa -7.78001 35.69671unique(prices[Region=="Kigoma",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Kigoma -4.89697 29.65987unique(prices[Region=="Morogoro",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Morogoro -6.82771 37.66542unique(prices[Region=="Mwanza",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Mwanza -2.51969 32.90144unique(prices[Region=="Mara",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Musoma -1.49982 33.8083unique(prices[Region=="Ruvuma",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Songea -10.67873 35.64836unique(prices[Region=="Shinyanga",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Shinyanga -3.67226 33.43069unique(prices[Region=="Kilimanjaro",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Moshi -3.34865 37.34352unique(prices[Region=="Mbeya",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Mbeya -8.909940 33.45517
## 2: Igawilo -8.924246 33.56788unique(prices[Region=="Katavi",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Mpanda -6.34295 31.07299
## 2: Mpimbwe -7.24425 31.81782
## 3: majengo -6.34809 31.07055
## 4: Majengo -6.34809 31.07055unique(prices[Region=="Njombe",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Njombe -9.33805 34.76949
## 2: Makambako -8.84398 34.82299unique(prices[Region=="Lindi",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Lindi -9.995 39.708unique(prices[Region=="Singida",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Singida -4.812610 34.7428
## 2: Namfua -4.821121 34.7470unique(prices[Region=="Pwani",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Pwani -6.44221 38.90622
## 2: Mnalani -6.44221 38.90622unique(prices[Region=="Simiyu",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Bariadi -2.80355 33.98699unique(prices[Region=="Geita",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Nyankumbu -2.870760 32.23408
## 2: Nyankumbu -2.895955 32.22871unique(prices[Region=="Songwe",.(Market, Latitude, Longitude)])## Market Latitude Longitude
## <char> <num> <num>
## 1: Songwe -8.95179 33.24377setnames(prices, old = "Maize..min.price.", new = "mai.price.min")
setnames(prices, old = "Rice..min.price.", new = "ric.price.min")
setnames(prices, old = "Sorghum..min.price.", new = "sor.price.min")
setnames(prices, old = "Bulrush.Millet..min.price.", new = "bul.price.min")
setnames(prices, old = "Finger.Millet..min.price.", new = "fin.price.min")
setnames(prices, old = "Wheat..min.price.", new = "whe.price.min")
setnames(prices, old = "Beans..min.price.", new = "bea.price.min")
setnames(prices, old = "Irish.Potatoes..min.price.", new = "pot.price.min")
setnames(prices, old = "Maize..max.price.", new = "mai.price.max")
setnames(prices, old = "Rice..max.price.", new = "ric.price.max")
setnames(prices, old = "Sorghum..max.price.", new = "sor.price.max")
setnames(prices, old = "Bulrush.Millet..max.price.", new = "bul.price.max")
setnames(prices, old = "Finger.Millet..max.price.", new = "fin.price.max")
setnames(prices, old = "Wheat..max.price.", new = "whe.price.max")
setnames(prices, old = "Beans..max.price.", new = "bea.price.max")
setnames(prices, old = "Irish.Potatoes..max.price.", new = "pot.price.max")sapply(prices, class)## Region Market mai.price.min mai.price.max ric.price.min
## "character" "character" "integer" "integer" "integer"
## ric.price.max sor.price.min sor.price.max bul.price.min bul.price.max
## "integer" "integer" "integer" "integer" "integer"
## fin.price.min fin.price.max whe.price.min whe.price.max bea.price.min
## "integer" "character" "integer" "integer" "integer"
## bea.price.max pot.price.min pot.price.max Date Latitude
## "integer" "integer" "integer" "Date" "numeric"
## Longitude
## "numeric"#convert prices to numeric
prices$mai.price.min <- as.numeric(prices$mai.price.min)
prices$ric.price.min <- as.numeric(prices$ric.price.min)
prices$sor.price.min <- as.numeric(prices$sor.price.min)
prices$bul.price.min <- as.numeric(prices$bul.price.min)
prices$fin.price.min <- as.numeric(prices$fin.price.min)
prices$whe.price.min <- as.numeric(prices$whe.price.min)
prices$bea.price.min <- as.numeric(prices$bea.price.min)
prices$pot.price.min <- as.numeric(prices$pot.price.min)
prices$mai.price.max <- as.numeric(prices$mai.price.max)
prices$ric.price.max <- as.numeric(prices$ric.price.max)
prices$sor.price.max <- as.numeric(prices$sor.price.max)
prices$bul.price.max <- as.numeric(prices$bul.price.max)
prices$fin.price.max <- as.numeric(prices$fin.price.max)## Warning: NAs introduced by coercionprices$whe.price.max <- as.numeric(prices$whe.price.max)
prices$bea.price.max <- as.numeric(prices$bea.price.max)
prices$pot.price.max <- as.numeric(prices$pot.price.max)
sapply(prices, class)## Region Market mai.price.min mai.price.max ric.price.min
## "character" "character" "numeric" "numeric" "numeric"
## ric.price.max sor.price.min sor.price.max bul.price.min bul.price.max
## "numeric" "numeric" "numeric" "numeric" "numeric"
## fin.price.min fin.price.max whe.price.min whe.price.max bea.price.min
## "numeric" "numeric" "numeric" "numeric" "numeric"
## bea.price.max pot.price.min pot.price.max Date Latitude
## "numeric" "numeric" "numeric" "Date" "numeric"
## Longitude
## "numeric"# convert to price per kg
prices$mai.price.min <- prices$mai.price.min/100
prices$ric.price.min <- prices$ric.price.min/100
prices$sor.price.min <- prices$sor.price.min/100
prices$bul.price.min <- prices$bul.price.min/100
prices$fin.price.min <- prices$fin.price.min/100
prices$whe.price.min <- prices$whe.price.min/100
prices$bea.price.min <- prices$bea.price.min/100
prices$pot.price.min <- prices$pot.price.min/100
prices$mai.price.max <- prices$mai.price.max/100
prices$ric.price.max <- prices$ric.price.max/100
prices$sor.price.max <- prices$sor.price.max/100
prices$bul.price.max <- prices$bul.price.max/100
prices$fin.price.max <- prices$fin.price.max/100
prices$whe.price.max <- prices$whe.price.max/100
prices$bea.price.max <- prices$bea.price.max/100
prices$pot.price.max <- prices$pot.price.max/100# calculate average of min and max
prices$mai.price <- (prices$mai.price.min + prices$mai.price.max) / 2
prices$ric.price <- (prices$ric.price.min + prices$ric.price.max) / 2
prices$sor.price <- (prices$sor.price.min + prices$sor.price.max) / 2
prices$bul.price <- (prices$bul.price.min + prices$bul.price.max) / 2
prices$fin.price <- (prices$fin.price.min + prices$fin.price.max) / 2
prices$whe.price <- (prices$whe.price.min + prices$whe.price.max) / 2
prices$bea.price <- (prices$bea.price.min + prices$bea.price.max) / 2
prices$pot.price <- (prices$pot.price.min + prices$pot.price.max) / 2#We can add dates by using the year and the month names
prices$Day <- day(prices$Date)
prices$Month <- month(prices$Date)
prices$Year <- year(prices$Date)# drop unneccessary columns
prices <- prices[,!c("mai.price.min", "mai.price.max",
"ric.price.min", "ric.price.max",
"sor.price.min", "sor.price.max",
"bul.price.min", "bul.price.max",
"fin.price.min", "fin.price.max",
"whe.price.min", "whe.price.max",
"bea.price.min", "bea.price.max",
"pot.price.min", "pot.price.max")]# calculate monthly mean prices by market
prices.monthly <- prices[, .(mai.price = mean(mai.price, na.rm = TRUE),
ric.price = mean(ric.price, na.rm = TRUE),
sor.price = mean(sor.price, na.rm = TRUE),
bul.price = mean(bul.price, na.rm = TRUE),
fin.price = mean(fin.price, na.rm = TRUE),
whe.price = mean(whe.price, na.rm = TRUE),
bea.price = mean(bea.price, na.rm = TRUE),
pot.price = mean(pot.price, na.rm = TRUE)),
by=.(Region, Market, Month, Year, Latitude, Longitude)]# reshape to long (so that prices for different commodities can be simultaneously estimated)
prices.monthly ## Region Market Month Year Latitude Longitude mai.price
## <char> <char> <int> <int> <num> <num> <num>
## 1: Arusha Arusha 5 2021 -3.36968 36.68808 469.00
## 2: Dar es Salaam Temeke 5 2021 -6.85097 39.25672 485.00
## 3: Dodoma Majengo 5 2021 -6.17971 35.74109 501.50
## 4: Dodoma Kibaigwa 5 2021 -6.08110 36.64645 375.00
## 5: Kagera Bukoba 5 2021 -1.33140 31.81293 426.25
## ---
## 1229: Pwani Mnalani 7 2025 -6.44221 38.90622 1278.00
## 1230: Simiyu Bariadi 7 2025 -2.80355 33.98699 840.00
## 1231: Kigoma Kigoma 7 2025 -4.89697 29.65987 744.75
## 1232: Njombe Makambako 7 2025 -8.84398 34.82299 569.50
## 1233: Rukwa Sumbawanga 7 2025 -7.95239 31.62319 561.00
## ric.price sor.price bul.price fin.price whe.price bea.price pot.price
## <num> <num> <num> <num> <num> <num> <num>
## 1: 1530 695.0 761.0 1333.0 793.0 1545 745
## 2: 1650 900.0 900.0 1775.0 1230.0 2420 730
## 3: 1592 558.0 581.0 1752.0 NaN 2075 618
## 4: NaN 437.5 NaN NaN NaN NaN NaN
## 5: 1260 1375.0 1337.5 1637.5 1637.5 1300 675
## ---
## 1229: 2060 1900.0 1900.0 2150.0 2300.0 2900 1450
## 1230: 2000 1900.0 1900.0 2150.0 2300.0 2900 1450
## 1231: 2175 2325.0 1310.0 2710.0 2850.0 2100 960
## 1232: 2375 850.0 NaN 1625.0 1550.0 2150 1100
## 1233: 2150 NaN NaN 975.0 1080.0 2400 685prices.monthly.long <- melt(prices.monthly, id.vars=c('Region', 'Market', 'Month', 'Year', 'Latitude', 'Longitude'),)# rename columns
setnames(prices.monthly.long, old="variable", new="Crop")
setnames(prices.monthly.long, old="value", new="pkg")# replace crop names
prices.monthly.long[Crop == "mai.price", Crop := "Maize"]
prices.monthly.long[Crop == "ric.price", Crop := "Rice"]
prices.monthly.long[Crop == "sor.price", Crop := "Sorghum"]
prices.monthly.long[Crop == "bul.price", Crop := "B.Millet"]
prices.monthly.long[Crop == "fin.price", Crop := "F.Millet"]
prices.monthly.long[Crop == "whe.price", Crop := "Wheat"]
prices.monthly.long[Crop == "bea.price", Crop := "Beans"]
prices.monthly.long[Crop == "pot.price", Crop := "Potato"]# Reset the factor levels to updated levels
prices.monthly.long[, Crop := factor(Crop)]
# Check the unique values again
unique(prices.monthly.long$Crop)## [1] Maize Rice Sorghum B.Millet F.Millet Wheat Beans Potato
## Levels: Maize Rice Sorghum B.Millet F.Millet Wheat Beans Potato# generate dummies to use in place of factors (for later spatial predictions, which are struggling with factors)
prices.monthly.long[, maize := ifelse(Crop == "Maize",1,0)]
prices.monthly.long[, rice := ifelse(Crop == "Rice",1,0)]
prices.monthly.long[, sorghum := ifelse(Crop == "Sorghum",1,0)]
prices.monthly.long[, bmillet := ifelse(Crop == "B.Millet",1,0)]
prices.monthly.long[, fmillet := ifelse(Crop == "F.Millet",1,0)]
prices.monthly.long[, wheat := ifelse(Crop == "Wheat",1,0)]
prices.monthly.long[, beans := ifelse(Crop == "Beans",1,0)]
prices.monthly.long[, potato := ifelse(Crop == "Potato",1,0)]prices.monthly.long[, jan := ifelse(Month == 1 , 1, 0)]
prices.monthly.long[, feb := ifelse(Month == 2 , 1, 0)]
prices.monthly.long[, mar := ifelse(Month == 3 , 1, 0)]
prices.monthly.long[, apr := ifelse(Month == 4 , 1, 0)]
prices.monthly.long[, may := ifelse(Month == 5 , 1, 0)]
prices.monthly.long[, jun := ifelse(Month == 6 , 1, 0)]
prices.monthly.long[, jul := ifelse(Month == 7 , 1, 0)]
prices.monthly.long[, aug := ifelse(Month == 8 , 1, 0)]
prices.monthly.long[, sep := ifelse(Month == 9 , 1, 0)]
prices.monthly.long[, oct := ifelse(Month == 10 , 1, 0)]
prices.monthly.long[, nov := ifelse(Month == 11 , 1, 0)]
prices.monthly.long[, dec := ifelse(Month == 12 , 1, 0)]# replace NaN with NAs in the price observations
prices.monthly.long[is.nan(pkg), pkg := NA]
# Remove observations with missing observations
prices.monthly.long <- na.omit(prices.monthly.long)# bring in raster stack as predictors
geodata_path("H:/Tanzania Price data/Datasets/geodata")
list.files("H:/Tanzania Price data/Datasets/geodata", recursive=TRUE)## [1] "climate/wc2.1_country/TZA_wc2.1_30s_bio.tif"
## [2] "gadm/gadm41_TZA_0_pk.rds"
## [3] "gadm/gadm41_TZA_1_pk.rds"
## [4] "gadm/gadm41_TZA_2_pk.rds"
## [5] "gadm/gadm41_TZA_3_pk.rds"
## [6] "travel/travel_time_to_cities_u3.tif"
## [7] "travel/travel_time_to_cities_u5.tif"
## [8] "TRUE/gadm/gadm41_TZA_0_pk.rds"
## [9] "TRUE/gadm/gadm41_TZA_1_pk.rds"
## [10] "TRUE/gadm/gadm41_TZA_2_pk.rds"
## [11] "TRUE/gadm/gadm41_TZA_3_pk.rds"
## [12] "TRUE/spam/spam2017v2r1_ssa_area.zip"
## [13] "TRUE/spam/spam2017V2r1_SSA_H_MAIZ_A.tif"
## [14] "TRUE/spam/spam2017V2r1_SSA_H_MAIZ_H.tif"
## [15] "TRUE/spam/spam2017V2r1_SSA_H_MAIZ_I.tif"
## [16] "TRUE/spam/spam2017V2r1_SSA_H_MAIZ_L.tif"
## [17] "TRUE/spam/spam2017V2r1_SSA_H_MAIZ_R.tif"
## [18] "TRUE/spam/spam2017V2r1_SSA_H_MAIZ_S.tif"
## [19] "TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_A.tif"
## [20] "TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_H.tif"
## [21] "TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_I.tif"
## [22] "TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_L.tif"
## [23] "TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_R.tif"
## [24] "TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_S.tif"
## [25] "TRUE/spam/spam2017v2r1_ssa_yield.zip"
## [26] "TRUE/travel/travel_time_to_cities_u5.tif"
## [27] "TRUE/travel/travel_time_to_ports_1.tif"
## [28] "TRUE/wc2.1_country/TZA_wc2.1_30s_bio.tif"# tza0 <- gadm(country="TZA", level=0)
# tza1 <- gadm(country="TZA", level=1)
# tza2 <- gadm(country="TZA", level=2)
# tza3 <- gadm(country="TZA", level=3)
tza0 <- readRDS("H:/Tanzania Price data/Datasets/geodata/TRUE/gadm/gadm41_TZA_0_pk.rds")
tza1 <- readRDS("H:/Tanzania Price data/Datasets/geodata/TRUE/gadm/gadm41_TZA_1_pk.rds")
tza2 <- readRDS("H:/Tanzania Price data/Datasets/geodata/TRUE/gadm/gadm41_TZA_2_pk.rds")
tza3 <- readRDS("H:/Tanzania Price data/Datasets/geodata/TRUE/gadm/gadm41_TZA_3_pk.rds")# convert prices observations to vector for mapping
mypts <- vect(prices.monthly.long, geom=c("Longitude", "Latitude"), crs=crs(tza0), keepgeom=TRUE)
# see if these show up correctly
plot(tza1)
plot(mypts, col="Red", add=TRUE)
# text(mypts, label="Market")# Extract Region information from GADM
region_info <- terra::extract(tza1, mypts)
head(region_info)## id.y GID_1 GID_0 COUNTRY NAME_1 VARNAME_1 NL_NAME_1
## 1 1 TZA.1_1 TZA Tanzania Arusha <NA> <NA>
## 2 2 TZA.2_1 TZA Tanzania Dar es Salaam Dar-es-salaam|Pwani <NA>
## 3 3 TZA.3_1 TZA Tanzania Dodoma <NA> <NA>
## 4 4 TZA.3_1 TZA Tanzania Dodoma <NA> <NA>
## 5 5 TZA.6_1 TZA Tanzania Kagera <NA> <NA>
## 6 6 TZA.11_1 TZA Tanzania Manyara <NA> <NA>
## TYPE_1 ENGTYPE_1 CC_1 HASC_1 ISO_1
## 1 Mkoa Region 02 TZ.AS TZ-01
## 2 Mkoa Region 07 TZ.DS TZ-02
## 3 Mkoa Region 01 TZ.DO TZ-03
## 4 Mkoa Region 01 TZ.DO TZ-03
## 5 Mkoa Region 18 TZ.KG TZ-05
## 6 Mkoa Region 21 TZ.MY TZ-26# merge district names to the dataset
prices.monthly.long$Region_GADM <- region_info$NAME_1
#prices.monthly.long$District_GADM <- region_info$NAME_2
head(prices.monthly.long)## Region Market Month Year Latitude Longitude Crop pkg maize
## <char> <char> <int> <int> <num> <num> <fctr> <num> <num>
## 1: Arusha Arusha 5 2021 -3.36968 36.68808 Maize 469.00 1
## 2: Dar es Salaam Temeke 5 2021 -6.85097 39.25672 Maize 485.00 1
## 3: Dodoma Majengo 5 2021 -6.17971 35.74109 Maize 501.50 1
## 4: Dodoma Kibaigwa 5 2021 -6.08110 36.64645 Maize 375.00 1
## 5: Kagera Bukoba 5 2021 -1.33140 31.81293 Maize 426.25 1
## 6: Manyara Babati 5 2021 -4.21006 35.74915 Maize 375.00 1
## rice sorghum bmillet fmillet wheat beans potato jan feb mar apr
## <num> <num> <num> <num> <num> <num> <num> <num> <num> <num> <num>
## 1: 0 0 0 0 0 0 0 0 0 0 0
## 2: 0 0 0 0 0 0 0 0 0 0 0
## 3: 0 0 0 0 0 0 0 0 0 0 0
## 4: 0 0 0 0 0 0 0 0 0 0 0
## 5: 0 0 0 0 0 0 0 0 0 0 0
## 6: 0 0 0 0 0 0 0 0 0 0 0
## may jun jul aug sep oct nov dec Region_GADM
## <num> <num> <num> <num> <num> <num> <num> <num> <char>
## 1: 1 0 0 0 0 0 0 0 Arusha
## 2: 1 0 0 0 0 0 0 0 Dar es Salaam
## 3: 1 0 0 0 0 0 0 0 Dodoma
## 4: 1 0 0 0 0 0 0 0 Dodoma
## 5: 1 0 0 0 0 0 0 0 Kagera
## 6: 1 0 0 0 0 0 0 0 Manyaraunique(prices.monthly.long[, c("Region", "Market", "Region_GADM")])## Region Market Region_GADM
## <char> <char> <char>
## 1: Arusha Arusha Arusha
## 2: Dar es Salaam Temeke Dar es Salaam
## 3: Dodoma Majengo Dodoma
## 4: Dodoma Kibaigwa Dodoma
## 5: Kagera Bukoba Kagera
## 6: Manyara Babati Manyara
## 7: Rukwa Sumbawanga Rukwa
## 8: Katavi Mpanda Katavi
## 9: Mtwara Mtwara Mtwara
## 10: Tabora Tabora Tabora
## 11: Tanga Tanga Tanga
## 12: Dar es Salaam Kinondoni Dar es Salaam
## 13: Dar es Salaam Ilala Dar es Salaam
## 14: Iringa Iringa Iringa
## 15: Kigoma Kigoma Kigoma
## 16: Morogoro Morogoro Morogoro
## 17: Mwanza Mwanza Mwanza
## 18: Mara Musoma Mara
## 19: Ruvuma Songea Ruvuma
## 20: Shinyanga Shinyanga Shinyanga
## 21: Kilimanjaro Moshi Kilimanjaro
## 22: Mbeya Mbeya Mbeya
## 23: Njombe Njombe Njombe
## 24: Lindi Lindi Lindi
## 25: Manyara Manyara Manyara
## 26: Dar es Salaam Tandika Dar es Salaam
## 27: Dar es Salaam Buguruni Dar es Salaam
## 28: Dar es Salaam Tandale Dar es Salaam
## 29: Singida Singida Singida
## 30: Pwani Pwani Pwani
## 31: Simiyu Bariadi Simiyu
## 32: Katavi Mpimbwe Katavi
## 33: Geita Nyankumbu Geita
## 34: Songwe Songwe Mbeya
## 35: Dar es Salaam Mwananyamala Dar es Salaam
## 36: Singida Namfua Singida
## 37: Arusha Kilombero Arusha
## 38: Tanga Mgandini Tanga
## 39: Katavi majengo Katavi
## 40: Mbeya Igawilo Mbeya
## 41: Pwani Mnalani Pwani
## 42: Katavi Majengo Katavi
## 43: Njombe Makambako Njombe
## 44: Dar es Salaam Ubungo Dar es Salaam
## Region Market Region_GADM# Check if Region matches Region_GADM
prices.monthly.long$region_match <- prices.monthly.long$Region == prices.monthly.long$Region_GADM
table(prices.monthly.long$region_match)##
## FALSE TRUE
## 129 8573# Show mismatched rows
mis <- subset(prices.monthly.long, region_match == FALSE)
mis## Region Market Month Year Latitude Longitude Crop pkg maize rice
## <char> <char> <int> <int> <num> <num> <fctr> <num> <num> <num>
## 1: Songwe Songwe 3 2024 -8.95179 33.24377 Maize 366.5000 1 0
## 2: Songwe Songwe 4 2024 -8.95179 33.24377 Maize 435.3889 1 0
## 3: Songwe Songwe 5 2024 -8.95179 33.24377 Maize 438.4615 1 0
## 4: Songwe Songwe 6 2024 -8.95179 33.24377 Maize 433.0000 1 0
## 5: Songwe Songwe 7 2024 -8.95179 33.24377 Maize 433.0000 1 0
## ---
## 125: Songwe Songwe 3 2025 -8.95179 33.24377 Potato 1500.0000 0 0
## 126: Songwe Songwe 4 2025 -8.95179 33.24377 Potato 1812.5000 0 0
## 127: Songwe Songwe 5 2025 -8.95179 33.24377 Potato 1769.2308 0 0
## 128: Songwe Songwe 6 2025 -8.95179 33.24377 Potato 2000.0000 0 0
## 129: Songwe Songwe 7 2025 -8.95179 33.24377 Potato 1000.0000 0 0
## sorghum bmillet fmillet wheat beans potato jan feb mar apr may
## <num> <num> <num> <num> <num> <num> <num> <num> <num> <num> <num>
## 1: 0 0 0 0 0 0 0 0 1 0 0
## 2: 0 0 0 0 0 0 0 0 0 1 0
## 3: 0 0 0 0 0 0 0 0 0 0 1
## 4: 0 0 0 0 0 0 0 0 0 0 0
## 5: 0 0 0 0 0 0 0 0 0 0 0
## ---
## 125: 0 0 0 0 0 1 0 0 1 0 0
## 126: 0 0 0 0 0 1 0 0 0 1 0
## 127: 0 0 0 0 0 1 0 0 0 0 1
## 128: 0 0 0 0 0 1 0 0 0 0 0
## 129: 0 0 0 0 0 1 0 0 0 0 0
## jun jul aug sep oct nov dec Region_GADM region_match
## <num> <num> <num> <num> <num> <num> <num> <char> <lgcl>
## 1: 0 0 0 0 0 0 0 Mbeya FALSE
## 2: 0 0 0 0 0 0 0 Mbeya FALSE
## 3: 0 0 0 0 0 0 0 Mbeya FALSE
## 4: 1 0 0 0 0 0 0 Mbeya FALSE
## 5: 0 1 0 0 0 0 0 Mbeya FALSE
## ---
## 125: 0 0 0 0 0 0 0 Mbeya FALSE
## 126: 0 0 0 0 0 0 0 Mbeya FALSE
## 127: 0 0 0 0 0 0 0 Mbeya FALSE
## 128: 1 0 0 0 0 0 0 Mbeya FALSE
## 129: 0 1 0 0 0 0 0 Mbeya FALSEunique(mis[, c("Region", "Market", "Region_GADM")])## Region Market Region_GADM
## <char> <char> <char>
## 1: Songwe Songwe Mbeya# Songwe is mismatched, rename correctly
prices.monthly.long[Region_GADM == "Mbeya" & Region == "Songwe", Region_GADM := "Songwe"]
unique(prices.monthly.long[, c("Region", "Market", "Region_GADM")])## Region Market Region_GADM
## <char> <char> <char>
## 1: Arusha Arusha Arusha
## 2: Dar es Salaam Temeke Dar es Salaam
## 3: Dodoma Majengo Dodoma
## 4: Dodoma Kibaigwa Dodoma
## 5: Kagera Bukoba Kagera
## 6: Manyara Babati Manyara
## 7: Rukwa Sumbawanga Rukwa
## 8: Katavi Mpanda Katavi
## 9: Mtwara Mtwara Mtwara
## 10: Tabora Tabora Tabora
## 11: Tanga Tanga Tanga
## 12: Dar es Salaam Kinondoni Dar es Salaam
## 13: Dar es Salaam Ilala Dar es Salaam
## 14: Iringa Iringa Iringa
## 15: Kigoma Kigoma Kigoma
## 16: Morogoro Morogoro Morogoro
## 17: Mwanza Mwanza Mwanza
## 18: Mara Musoma Mara
## 19: Ruvuma Songea Ruvuma
## 20: Shinyanga Shinyanga Shinyanga
## 21: Kilimanjaro Moshi Kilimanjaro
## 22: Mbeya Mbeya Mbeya
## 23: Njombe Njombe Njombe
## 24: Lindi Lindi Lindi
## 25: Manyara Manyara Manyara
## 26: Dar es Salaam Tandika Dar es Salaam
## 27: Dar es Salaam Buguruni Dar es Salaam
## 28: Dar es Salaam Tandale Dar es Salaam
## 29: Singida Singida Singida
## 30: Pwani Pwani Pwani
## 31: Simiyu Bariadi Simiyu
## 32: Katavi Mpimbwe Katavi
## 33: Geita Nyankumbu Geita
## 34: Songwe Songwe Songwe
## 35: Dar es Salaam Mwananyamala Dar es Salaam
## 36: Singida Namfua Singida
## 37: Arusha Kilombero Arusha
## 38: Tanga Mgandini Tanga
## 39: Katavi majengo Katavi
## 40: Mbeya Igawilo Mbeya
## 41: Pwani Mnalani Pwani
## 42: Katavi Majengo Katavi
## 43: Njombe Makambako Njombe
## 44: Dar es Salaam Ubungo Dar es Salaam
## Region Market Region_GADM# Check if Region matches Region_GADM
prices.monthly.long$region_match <- prices.monthly.long$Region == prices.monthly.long$Region_GADM
table(prices.monthly.long$region_match)##
## TRUE
## 8702# convert prices observations to vector for mapping
mypts <- vect(prices.monthly.long, geom=c("Longitude", "Latitude"), crs=crs(tza0), keepgeom=TRUE)
mkt_pts <- unique(mypts[, c("Market", "Latitude", "Longitude")])
head(mkt_pts)## Market Latitude Longitude
## 1 Arusha -3.36968 36.68808
## 2 Temeke -6.85097 39.25672
## 3 Majengo -6.17971 35.74109
## 4 Kibaigwa -6.08110 36.64645
## 5 Bukoba -1.33140 31.81293
## 6 Babati -4.21006 35.74915# create reference raster
tza_extent <- ext(tza1) |> floor()
r <- crop(rast(res=1/12), tza_extent)Interpolate
## Interpolate
#xy <- as.matrix(mypts[,c("Longitude", "Latitude")])
xy <- geom(mypts)[,c("y","x")]
#tps <- Tps(xy, p$spatial)
tps <- Tps(xy, mypts$pkg)
sp <- interpolate(r, tps)
sp <- mask(sp, tza1)
plot(sp)
lines(tza1)
Predict prices with coodinates only
## Predict Maize prices with coodinates only
maize_mypts <- mypts[mypts$Crop == "Maize", ]
rf <- randomForest(pkg ~ Longitude + Latitude , data=maize_mypts)
sp3 <- interpolate(r, rf, xyNames=c("Longitude", "Latitude"))
sp3 <- mask(sp3, tza1)
plot(sp3)
lines(tza1)
Covariates
The covariates used include a mix of crop-specific indicators, temporal variables to capture monthly and yearly effects, geographical coordinates, accessibility measures, climatic conditions, and lagged rainfall to account for delayed effects of weather on crop prices.
ttcity <- rast("H:/Tanzania Price data/Datasets/geodata/TRUE/travel/travel_time_to_cities_u5.tif")
ttport <- rast("H:/Tanzania Price data/Datasets/geodata/TRUE/travel/travel_time_to_ports_1.tif")
clm <- rast("H:/Tanzania Price data/Datasets/geodata/TRUE/wc2.1_country/TZA_wc2.1_30s_bio.tif")
area <- rast("H:/Tanzania Price data/Datasets/geodata/TRUE/spam/spam2017V2r1_SSA_H_MAIZ_A.tif")
yield <- rast("H:/Tanzania Price data/Datasets/geodata/TRUE/spam/spam2017V2r1_SSA_Y_MAIZ_R.tif")
popd <- rast("gpw_v4_population_density_rev11_2020_10m.tif")names(ttcity) <- c("ttcity_u5") ## travel time cities of 100k or more
names(ttport) <- c("ttport_1") ## travel time to major ports
names(clm) <- gsub("wc2.1_30s_", "", names(clm))
names(area) <- c("MAI_ARE") # SPAM maize area 2010
names(yield) <- c("MAI_YLD") # SPAM maize yield 2010
names(popd) <- c("popdens") # GPW4
comment(ttcity) <- "travel time to cities 100k or more"
comment(ttport) <- "travel time to major ports"
comment(popd) <- "population density 2020 (GPW4 @ 10dm)"
comment(clm)[1] <-"BIO1 = Annual Mean Temperature"
comment(clm)[2] <-"BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))"
comment(clm)[3] <-"BIO3 = Isothermality (BIO2/BIO7) (×100)"
comment(clm)[4] <-"BIO4 = Temperature Seasonality (standard deviation ×100)"
comment(clm)[5] <-"BIO5 = Max Temperature of Warmest Month"
comment(clm)[6] <-"BIO6 = Min Temperature of Coldest Month"
comment(clm)[7] <-"BIO7 = Temperature Annual Range (BIO5-BIO6)"
comment(clm)[8] <-"BIO8 = Mean Temperature of Wettest Quarter"
comment(clm)[9] <-"BIO9 = Mean Temperature of Driest Quarter"
comment(clm)[10] <-"BIO10 = Mean Temperature of Warmest Quarter"
comment(clm)[11] <-"BIO11 = Mean Temperature of Coldest Quarter"
comment(clm)[12] <-"BIO12 = Annual Precipitation"
comment(clm)[13] <-"BIO13 = Precipitation of Wettest Month"
comment(clm)[14] <-"BIO14 = Precipitation of Driest Month"
comment(clm)[15] <-"BIO15 = Precipitation Seasonality (Coefficient of Variation)"
comment(clm)[16] <-"BIO16 = Precipitation of Wettest Quarter"
comment(clm)[17] <-"BIO17 = Precipitation of Driest Quarter"
comment(clm)[18] <-"BIO18 = Precipitation of Warmest Quarter"
comment(clm)[19] <-"BIO19 = Precipitation of Coldest Quarter"Harmonize rasters to national boundaries and common resolution
ttcity <- resample(ttcity, r)
ttport <- resample(ttport, r)
clm <- resample(clm, r)
area <- resample(area, r)
popd <- resample(popd, r)
freq(is.na(area))## layer value count
## 1 1 0 14393
## 2 1 1 6343area <- classify(area, cbind(NA,0))
yield <- resample(yield, r)
freq(is.na(yield))## layer value count
## 1 1 0 14393
## 2 1 1 6343yield <- classify(yield, cbind(NA,0))
# check again
compareGeom(ttcity, ttport, clm, area, yield, popd)## [1] TRUEGenerate Latitude and Longitude grid
latgrd <- longrd <- r
latgrd[] <- yFromCell(latgrd, 1:ncell(latgrd))
longrd[] <- xFromCell(longrd, 1:ncell(longrd))
names(latgrd) <- c("latitude")
names(longrd) <- c("longitude")Prepare Predictor Stack
rstack <- c(ttcity, ttport, clm, area, yield, popd, latgrd, longrd)
names(rstack)## [1] "ttcity_u5" "ttport_1" "bio_1" "bio_2" "bio_3" "bio_4"
## [7] "bio_5" "bio_6" "bio_7" "bio_8" "bio_9" "bio_10"
## [13] "bio_11" "bio_12" "bio_13" "bio_14" "bio_15" "bio_16"
## [19] "bio_17" "bio_18" "bio_19" "MAI_ARE" "MAI_YLD" "popdens"
## [25] "latitude" "longitude"# create focal mean to extract from (as alternative to using buffers for extraction, which are not supported in terra)
fm <- focalMat(r, d=0.18, type='circle', fillNA=FALSE)
rstack2 <- focal(rstack, w=fm, fun="mean", na.policy="all", fillvalue=NA, # na.rm=TRUE,
expand=TRUE, silent=FALSE) #, filename="", overwrite=FALSE) # extract values to dataset -- use a 20km buffer
# do a focal sum of 20km radius - this is about 0.18 of a decimal degree... 0.18*112=20.16
fm <- focalMat(r, d=0.18, type='circle', fillNA=FALSE)
rstack2 <- focal(rstack2, w=fm, fun="sum", na.policy="all", fillvalue=NA, na.rm=TRUE,
expand=TRUE, silent=FALSE) #, filename="", overwrite=FALSE) Lag Rainfall
Bring in Rainfall Data from Chirps
chirps_path <- "H:/Tanzania Price data/chirps_data"
chirps_files <- list.files(chirps_path, pattern = ".tif$", full.names = TRUE)
# Read all CHIRPS data files into a SpatRaster collection
chirps_rasters <- rast(chirps_files)
#crop to Tanzania boundary
Chirps_Tz <- crop(chirps_rasters, tza1)
writeRaster(Chirps_Tz, "Tz_chirps_monthly_croped.tif", overwrite=TRUE)
Tz_chirps_monthly <- terra::rast("Tz_chirps_monthly_croped.tif")
Tz_chirps_monthly## class : SpatRaster
## dimensions : 215, 222, 57 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : 29.35, 40.45, -11.75, -1.000001 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : Tz_chirps_monthly_croped.tif
## names : chirp~20.11, chirp~20.12, chirp~21.01, chirp~21.02, chirp~21.03, chirp~21.04, ...
## min values : -9999.0000, -9999.0000, -9999.000, -9999.0000, -9999.0000, -9999.0000, ...
## max values : 459.7112, 504.3329, 508.448, 585.3241, 750.1949, 960.3376, ...#Replace -9999 with NA
Tz_chirps_monthly <- classify(Tz_chirps_monthly, cbind(-9999,NA))
#extract layer names
layer_names <- names(Tz_chirps_monthly)
layer_names## [1] "chirps-v2.0.2020.11" "chirps-v2.0.2020.12" "chirps-v2.0.2021.01"
## [4] "chirps-v2.0.2021.02" "chirps-v2.0.2021.03" "chirps-v2.0.2021.04"
## [7] "chirps-v2.0.2021.05" "chirps-v2.0.2021.06" "chirps-v2.0.2021.07"
## [10] "chirps-v2.0.2021.08" "chirps-v2.0.2021.09" "chirps-v2.0.2021.10"
## [13] "chirps-v2.0.2021.11" "chirps-v2.0.2021.12" "chirps-v2.0.2022.01"
## [16] "chirps-v2.0.2022.02" "chirps-v2.0.2022.03" "chirps-v2.0.2022.04"
## [19] "chirps-v2.0.2022.05" "chirps-v2.0.2022.06" "chirps-v2.0.2022.07"
## [22] "chirps-v2.0.2022.08" "chirps-v2.0.2022.09" "chirps-v2.0.2022.10"
## [25] "chirps-v2.0.2022.11" "chirps-v2.0.2022.12" "chirps-v2.0.2023.01"
## [28] "chirps-v2.0.2023.02" "chirps-v2.0.2023.03" "chirps-v2.0.2023.04"
## [31] "chirps-v2.0.2023.05" "chirps-v2.0.2023.06" "chirps-v2.0.2023.07"
## [34] "chirps-v2.0.2023.08" "chirps-v2.0.2023.09" "chirps-v2.0.2023.10"
## [37] "chirps-v2.0.2023.11" "chirps-v2.0.2023.12" "chirps-v2.0.2024.01"
## [40] "chirps-v2.0.2024.02" "chirps-v2.0.2024.03" "chirps-v2.0.2024.04"
## [43] "chirps-v2.0.2024.05" "chirps-v2.0.2024.06" "chirps-v2.0.2024.07"
## [46] "chirps-v2.0.2024.08" "chirps-v2.0.2024.09" "chirps-v2.0.2024.10"
## [49] "chirps-v2.0.2024.11" "chirps-v2.0.2024.12" "chirps-v2.0.2025.01"
## [52] "chirps-v2.0.2025.02" "chirps-v2.0.2025.03" "chirps-v2.0.2025.04"
## [55] "chirps-v2.0.2025.05" "chirps-v2.0.2025.06" "chirps-v2.0.2025.07"# We need to create a sequence of dates from the layer names
# Extract year and month from layer names and convert to Date
dates <- as.Date(paste0(sub("chirps-v2.0\\.", "", layer_names), "-01"), format = "%Y.%m-%d")
dates## [1] "2020-11-01" "2020-12-01" "2021-01-01" "2021-02-01" "2021-03-01"
## [6] "2021-04-01" "2021-05-01" "2021-06-01" "2021-07-01" "2021-08-01"
## [11] "2021-09-01" "2021-10-01" "2021-11-01" "2021-12-01" "2022-01-01"
## [16] "2022-02-01" "2022-03-01" "2022-04-01" "2022-05-01" "2022-06-01"
## [21] "2022-07-01" "2022-08-01" "2022-09-01" "2022-10-01" "2022-11-01"
## [26] "2022-12-01" "2023-01-01" "2023-02-01" "2023-03-01" "2023-04-01"
## [31] "2023-05-01" "2023-06-01" "2023-07-01" "2023-08-01" "2023-09-01"
## [36] "2023-10-01" "2023-11-01" "2023-12-01" "2024-01-01" "2024-02-01"
## [41] "2024-03-01" "2024-04-01" "2024-05-01" "2024-06-01" "2024-07-01"
## [46] "2024-08-01" "2024-09-01" "2024-10-01" "2024-11-01" "2024-12-01"
## [51] "2025-01-01" "2025-02-01" "2025-03-01" "2025-04-01" "2025-05-01"
## [56] "2025-06-01" "2025-07-01"# Assign these dates to the SpatRaster object
time(Tz_chirps_monthly) <- dates
#rename the layers to the formatted dates
names(Tz_chirps_monthly) <- dates
# Check the SpatRaster object
print(Tz_chirps_monthly)## class : SpatRaster
## dimensions : 215, 222, 57 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : 29.35, 40.45, -11.75, -1.000001 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : Tz_chirps_monthly_croped
## names : 2020-11-01, 2020-12-01, 2021-01-01, 2021-02-01, 2021-03-01, 2021-04-01, ...
## min values : 5.859746, 1.903993, 0.8161677, 0.3187093, 0.9539621, 1.17571, ...
## max values : 459.711212, 504.332947, 508.4479980, 585.3240967, 750.1949463, 960.33765, ...
## time (days) : 2020-11-01 to 2025-07-01# do a focal mean of 100km radius - this is about 0.9 of a decimal degree... 0.9009*112=100.9008
# Calculate the focal mean for each layer (month)
fm_r <- focalMat(Tz_chirps_monthly, d=0.9, type='circle', fillNA=FALSE)
Rainfall_focal_sum_100km <- focal(Tz_chirps_monthly, w=fm_r, fun="mean", na.policy="all", fillvalue=NA, na.rm=TRUE,
expand=TRUE, silent=FALSE)
# Check the result
Rainfall_focal_sum_100km## class : SpatRaster
## dimensions : 215, 222, 57 (nrow, ncol, nlyr)
## resolution : 0.05, 0.05 (x, y)
## extent : 29.35, 40.45, -11.75, -1.000001 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## varname : Tz_chirps_monthly_croped
## names : 2020-11-01, 2020-12-01, 2021-01-01, 2021-02-01, 2021-03-01, 2021-04-01, ...
## min values : 20.79649, 27.27837, 8.092316, 4.489223, 11.30184, 15.74836, ...
## max values : 313.14248, 326.88227, 416.590468, 395.422070, 464.67840, 569.81990, ...
## time (days) : 2020-11-01 to 2025-07-01Rainfall <- Rainfall_focal_sum_100km#Resample
Rainfall_res <- resample(Rainfall, r)
Rainfall_res## class : SpatRaster
## dimensions : 144, 144, 57 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : 29, 41, -12, 0 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## names : 2020-11-01, 2020-12-01, 2021-01-01, 2021-02-01, 2021-03-01, 2021-04-01, ...
## min values : 20.89972, 32.9472, 8.116516, 4.527868, 11.83634, 15.90436, ...
## max values : 310.41034, 326.4339, 416.389160, 395.367279, 462.19266, 550.26270, ...
## time (days) : 2020-11-01 to 2025-07-01Define a function to calculate the 6-month lagged sum of rainfall values.
Here we define function that calculates 6 month lag sum of rainfall for each month in our data. The output is raster datsets.
calculate_lagged_sum <- function(raster_stack, num_months = 6) {
# Get the time vector from the raster stack
time_vector <- time(raster_stack)
# Initialize list to store lagged sum rasters
lagged_sum_rasters <- vector("list", length(time_vector))
# Loop through each layer in the raster stack
for (i in seq_along(time_vector)) {
if (i > num_months) { # We need at least 'num_months' previous layers to calculate the lagged sum
# Determine the start and end dates for the lag period
end_date <- time_vector[i] # Date of the current layer being processed
start_date <- end_date %m-% months(num_months) #The date num_months before the end_date
# Select the layers that fall within the lag period
lag_period_layers <- raster_stack[[which(time_vector > start_date & time_vector <= end_date)]]
# Calculate the sum of the selected layers
if (nlyr(lag_period_layers) == num_months) {
lagged_sum_rasters[[i]] <- sum(lag_period_layers, na.rm = TRUE)
} else {
lagged_sum_rasters[[i]] <- rast(nrow = nrow(raster_stack), ncol = ncol(raster_stack),
crs = crs(raster_stack), ext = ext(raster_stack),
vals = NA) # Use an empty raster with NA values
}
} else {
lagged_sum_rasters[[i]] <- rast(nrow = nrow(raster_stack), ncol = ncol(raster_stack),
crs = crs(raster_stack), ext = ext(raster_stack),
vals = NA) # Use an empty raster with NA values
}
}
# Combine the lagged sum rasters into a single raster stack, excluding empty rasters
lagged_sum_stack <- rast(lagged_sum_rasters)
# Set the names for the layers in the lagged sum stack
names(lagged_sum_stack) <- names(raster_stack)[!is.na(lagged_sum_rasters)]
return(lagged_sum_stack)
}
# Calculate the 6-month lagged sum for each period in the raster stack
lagged_rainfall_sum <- calculate_lagged_sum(Rainfall_res, num_months = 6)
# Remove the first 6 layers from the raster stack since they are empty
lagged_rainfall_sum_filtered <- lagged_rainfall_sum[[7:nlyr(lagged_rainfall_sum)]]
# check the result
print(lagged_rainfall_sum_filtered)## class : SpatRaster
## dimensions : 144, 144, 51 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : 29, 41, -12, 0 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## names : 2021-05-01, 2021-06-01, 2021-07-01, 2021-08-01, 2021-09-01, 2021-10-01, ...
## min values : 172.7856, 105.5447, 105.4121, 105.1864, 26.20571, 7.070483, ...
## max values : 1512.4689, 1355.9669, 1052.9220, 1051.9214, 1038.13874, 643.716021, ...names(lagged_rainfall_sum_filtered) <- paste0(names(lagged_rainfall_sum_filtered), "_rain.sum.lag")
plot(lagged_rainfall_sum_filtered)
## We'll have to include lagged_rainfall_sum_filtered in the predictor stack.
rstack## class : SpatRaster
## dimensions : 144, 144, 26 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : 29, 41, -12, 0 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## names : ttcity_u5, ttport_1, bio_1, bio_2, bio_3, bio_4, ...
## min values : 3.350081e-02, 997.9427, 4.333545, 6.546645, 53.85881, 19.75255, ...
## max values : 1.166560e+03, 3043.3459, 28.472235, 15.285786, 93.14099, 266.33594, ...#names(rstack)rstack3 <- c(rstack, lagged_rainfall_sum_filtered)
names(rstack3)## [1] "ttcity_u5" "ttport_1"
## [3] "bio_1" "bio_2"
## [5] "bio_3" "bio_4"
## [7] "bio_5" "bio_6"
## [9] "bio_7" "bio_8"
## [11] "bio_9" "bio_10"
## [13] "bio_11" "bio_12"
## [15] "bio_13" "bio_14"
## [17] "bio_15" "bio_16"
## [19] "bio_17" "bio_18"
## [21] "bio_19" "MAI_ARE"
## [23] "MAI_YLD" "popdens"
## [25] "latitude" "longitude"
## [27] "2021-05-01_rain.sum.lag" "2021-06-01_rain.sum.lag"
## [29] "2021-07-01_rain.sum.lag" "2021-08-01_rain.sum.lag"
## [31] "2021-09-01_rain.sum.lag" "2021-10-01_rain.sum.lag"
## [33] "2021-11-01_rain.sum.lag" "2021-12-01_rain.sum.lag"
## [35] "2022-01-01_rain.sum.lag" "2022-02-01_rain.sum.lag"
## [37] "2022-03-01_rain.sum.lag" "2022-04-01_rain.sum.lag"
## [39] "2022-05-01_rain.sum.lag" "2022-06-01_rain.sum.lag"
## [41] "2022-07-01_rain.sum.lag" "2022-08-01_rain.sum.lag"
## [43] "2022-09-01_rain.sum.lag" "2022-10-01_rain.sum.lag"
## [45] "2022-11-01_rain.sum.lag" "2022-12-01_rain.sum.lag"
## [47] "2023-01-01_rain.sum.lag" "2023-02-01_rain.sum.lag"
## [49] "2023-03-01_rain.sum.lag" "2023-04-01_rain.sum.lag"
## [51] "2023-05-01_rain.sum.lag" "2023-06-01_rain.sum.lag"
## [53] "2023-07-01_rain.sum.lag" "2023-08-01_rain.sum.lag"
## [55] "2023-09-01_rain.sum.lag" "2023-10-01_rain.sum.lag"
## [57] "2023-11-01_rain.sum.lag" "2023-12-01_rain.sum.lag"
## [59] "2024-01-01_rain.sum.lag" "2024-02-01_rain.sum.lag"
## [61] "2024-03-01_rain.sum.lag" "2024-04-01_rain.sum.lag"
## [63] "2024-05-01_rain.sum.lag" "2024-06-01_rain.sum.lag"
## [65] "2024-07-01_rain.sum.lag" "2024-08-01_rain.sum.lag"
## [67] "2024-09-01_rain.sum.lag" "2024-10-01_rain.sum.lag"
## [69] "2024-11-01_rain.sum.lag" "2024-12-01_rain.sum.lag"
## [71] "2025-01-01_rain.sum.lag" "2025-02-01_rain.sum.lag"
## [73] "2025-03-01_rain.sum.lag" "2025-04-01_rain.sum.lag"
## [75] "2025-05-01_rain.sum.lag" "2025-06-01_rain.sum.lag"
## [77] "2025-07-01_rain.sum.lag"#Extract to the point dataset
extr1 <- terra::extract(rstack3, mypts, method = "bilinear")
mypts <- cbind(mypts, extr1)
# Remove the ID column from the dataset
mypts <- mypts[, !names(mypts) %in% "ID"]Here we extract sum of lag rainfall for each row under the column rain.sum.lag
mypts_df <- as.data.frame(mypts)
# Define the function to obtain sum of lag rainfall from corresponding rasters to mypts under rain.sum.lag column (for each row)
# Each extraction has to match the month and year
get_rain_sum_row <- function(current_date, mypts_row) {
# Extract the rainfall value for the current date
rain_sum <- mypts_row[[paste0(current_date, "_rain.sum.lag")]]
return(rain_sum)
}
# Loop through each row and obtain the rainfall sum for each month and year
for (i in 1:nrow(mypts_df)) {
# Extract relevant data for the current row
month <- mypts_df$Month[i]
year <- mypts_df$Year[i]
current_date <- paste0(year, "-", sprintf("%02d", month), "-01") # Format date correctly
# Pass necessary data to the function
rain_sum <- get_rain_sum_row(current_date, mypts_df[i, ])
# Update the rain.sum.lag column
mypts_df$rain.sum.lag[i] <- rain_sum
}
# Update the SpatVector with the new rain.avg column
mypts$rain.sum.lag <- mypts_df$rain.sum.lag# I'll drop the dates with rain.sum.lag from mypts, seems redundant
column_indices <- grep("^202[0-4]-", names(mypts))
mypts <- mypts[, -column_indices]
names(mypts)## [1] "Region" "Market"
## [3] "Month" "Year"
## [5] "Latitude" "Longitude"
## [7] "Crop" "pkg"
## [9] "maize" "rice"
## [11] "sorghum" "bmillet"
## [13] "fmillet" "wheat"
## [15] "beans" "potato"
## [17] "jan" "feb"
## [19] "mar" "apr"
## [21] "may" "jun"
## [23] "jul" "aug"
## [25] "sep" "oct"
## [27] "nov" "dec"
## [29] "Region_GADM" "region_match"
## [31] "ttcity_u5" "ttport_1"
## [33] "bio_1" "bio_2"
## [35] "bio_3" "bio_4"
## [37] "bio_5" "bio_6"
## [39] "bio_7" "bio_8"
## [41] "bio_9" "bio_10"
## [43] "bio_11" "bio_12"
## [45] "bio_13" "bio_14"
## [47] "bio_15" "bio_16"
## [49] "bio_17" "bio_18"
## [51] "bio_19" "MAI_ARE"
## [53] "MAI_YLD" "popdens"
## [55] "latitude" "longitude"
## [57] "2025-01-01_rain.sum.lag" "2025-02-01_rain.sum.lag"
## [59] "2025-03-01_rain.sum.lag" "2025-04-01_rain.sum.lag"
## [61] "2025-05-01_rain.sum.lag" "2025-06-01_rain.sum.lag"
## [63] "2025-07-01_rain.sum.lag" "rain.sum.lag"#define Month as a factor
mypts$Month <- as.factor(mypts$Month)
levels(mypts$Month)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12"#We'll define month as an interger instead.
# Check to make sure Month is interger
sapply(mypts, class)## Region Market Month
## "character" "character" "factor"
## Year Latitude Longitude
## "integer" "numeric" "numeric"
## Crop pkg maize
## "factor" "numeric" "numeric"
## rice sorghum bmillet
## "numeric" "numeric" "numeric"
## fmillet wheat beans
## "numeric" "numeric" "numeric"
## potato jan feb
## "numeric" "numeric" "numeric"
## mar apr may
## "numeric" "numeric" "numeric"
## jun jul aug
## "numeric" "numeric" "numeric"
## sep oct nov
## "numeric" "numeric" "numeric"
## dec Region_GADM region_match
## "numeric" "character" "logical"
## ttcity_u5 ttport_1 bio_1
## "numeric" "numeric" "numeric"
## bio_2 bio_3 bio_4
## "numeric" "numeric" "numeric"
## bio_5 bio_6 bio_7
## "numeric" "numeric" "numeric"
## bio_8 bio_9 bio_10
## "numeric" "numeric" "numeric"
## bio_11 bio_12 bio_13
## "numeric" "numeric" "numeric"
## bio_14 bio_15 bio_16
## "numeric" "numeric" "numeric"
## bio_17 bio_18 bio_19
## "numeric" "numeric" "numeric"
## MAI_ARE MAI_YLD popdens
## "numeric" "numeric" "numeric"
## latitude longitude 2025-01-01_rain.sum.lag
## "numeric" "numeric" "numeric"
## 2025-02-01_rain.sum.lag 2025-03-01_rain.sum.lag 2025-04-01_rain.sum.lag
## "numeric" "numeric" "numeric"
## 2025-05-01_rain.sum.lag 2025-06-01_rain.sum.lag 2025-07-01_rain.sum.lag
## "numeric" "numeric" "numeric"
## rain.sum.lag
## "numeric"# drop levels that don't exist in Crop field
mypts$Crop <- mypts$Crop[,drop=TRUE]
levels(mypts$Crop)## [1] "Maize" "Rice" "Sorghum" "B.Millet" "F.Millet" "Wheat" "Beans"
## [8] "Potato"Linear model price Prediction
Coefficient estimates from the linear model provide a detailed insight into the relationship between each predictor and the response variable.
# Fit the linear model
lm_model <- lm(pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato +
Month +
Year +
ttcity_u5 + ttport_1 +
bio_1 + bio_2 + bio_3 + bio_4 + bio_5 + bio_6 +
bio_7 + bio_8 + bio_9 + bio_10 + bio_11 + bio_12 + bio_13
+ bio_14 + bio_15 + bio_16 + bio_17 + bio_18 + bio_19 +
MAI_ARE + MAI_YLD +
Longitude + Latitude +
rain.sum.lag,
data = mypts)
summary(lm_model)##
## Call:
## lm(formula = pkg ~ maize + rice + sorghum + bmillet + fmillet +
## wheat + beans + potato + Month + Year + ttcity_u5 + ttport_1 +
## bio_1 + bio_2 + bio_3 + bio_4 + bio_5 + bio_6 + bio_7 + bio_8 +
## bio_9 + bio_10 + bio_11 + bio_12 + bio_13 + bio_14 + bio_15 +
## bio_16 + bio_17 + bio_18 + bio_19 + MAI_ARE + MAI_YLD + Longitude +
## Latitude + rain.sum.lag, data = mypts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1362.94 -267.39 -25.44 247.62 2292.08
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.812e+05 8.096e+03 -34.735 < 2e-16 ***
## maize -1.090e+02 1.683e+01 -6.473 1.01e-10 ***
## rice 1.334e+03 1.679e+01 79.452 < 2e-16 ***
## sorghum 3.952e+02 1.758e+01 22.485 < 2e-16 ***
## bmillet 5.070e+02 1.857e+01 27.304 < 2e-16 ***
## fmillet 8.862e+02 1.733e+01 51.126 < 2e-16 ***
## wheat 9.451e+02 1.868e+01 50.586 < 2e-16 ***
## beans 1.593e+03 1.681e+01 94.814 < 2e-16 ***
## potato NA NA NA NA
## Month2 2.953e+01 2.269e+01 1.302 0.192997
## Month3 2.289e+01 2.341e+01 0.978 0.328306
## Month4 4.101e+01 2.465e+01 1.664 0.096236 .
## Month5 2.217e+01 2.356e+01 0.941 0.346656
## Month6 7.826e+00 2.221e+01 0.352 0.724572
## Month7 -2.822e+01 2.157e+01 -1.308 0.190949
## Month8 -5.102e+01 2.326e+01 -2.193 0.028304 *
## Month9 -3.504e+01 2.370e+01 -1.479 0.139232
## Month10 -1.211e+01 2.490e+01 -0.487 0.626605
## Month11 3.534e+01 2.436e+01 1.451 0.146925
## Month12 7.585e+01 2.305e+01 3.291 0.001001 **
## Year 1.328e+02 3.819e+00 34.775 < 2e-16 ***
## ttcity_u5 1.461e+00 1.533e-01 9.528 < 2e-16 ***
## ttport_1 -1.547e+00 1.697e-01 -9.116 < 2e-16 ***
## bio_1 -2.407e+03 1.825e+02 -13.193 < 2e-16 ***
## bio_2 -1.672e+03 1.545e+02 -10.819 < 2e-16 ***
## bio_3 1.743e+02 1.435e+01 12.145 < 2e-16 ***
## bio_4 9.592e+00 3.429e+00 2.797 0.005164 **
## bio_5 1.347e+03 1.239e+02 10.872 < 2e-16 ***
## bio_6 -1.142e+03 1.163e+02 -9.818 < 2e-16 ***
## bio_7 NA NA NA NA
## bio_8 8.984e+02 7.246e+01 12.398 < 2e-16 ***
## bio_9 -2.292e+01 5.912e+01 -0.388 0.698180
## bio_10 -4.150e+02 1.510e+02 -2.748 0.006013 **
## bio_11 1.775e+03 2.288e+02 7.757 9.65e-15 ***
## bio_12 1.370e+00 2.436e-01 5.624 1.93e-08 ***
## bio_13 1.769e+00 7.849e-01 2.254 0.024223 *
## bio_14 3.591e+01 9.888e+00 3.631 0.000284 ***
## bio_15 1.387e+01 1.862e+00 7.449 1.03e-13 ***
## bio_16 -1.646e+00 5.187e-01 -3.174 0.001510 **
## bio_17 -8.649e+00 2.677e+00 -3.231 0.001237 **
## bio_18 -7.503e-01 1.888e-01 -3.974 7.14e-05 ***
## bio_19 -4.436e-01 2.098e+00 -0.211 0.832550
## MAI_ARE 4.414e-02 1.775e-02 2.486 0.012921 *
## MAI_YLD 7.496e-02 1.006e-02 7.450 1.02e-13 ***
## Longitude 3.727e+01 2.039e+01 1.828 0.067578 .
## Latitude -7.499e+01 2.176e+01 -3.447 0.000570 ***
## rain.sum.lag -9.795e-02 2.445e-02 -4.005 6.24e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 413.6 on 8657 degrees of freedom
## Multiple R-squared: 0.705, Adjusted R-squared: 0.7035
## F-statistic: 470.2 on 44 and 8657 DF, p-value: < 2.2e-16RandomForest and TPS
Random Forest price prediction
First check if there are any predictors with NA values
for(column in seq_along(mypts)){
if(any(is.na(mypts[column]))){
print(paste0("Column: ", colnames(mypts)[column], " has at least one NA value"))
}
}
#There are no columns with missing valuesSplit data the to be used for Training and validation
Data from May 2021 - Dec 2023 will be used for model training while more recent data from Jan - July 2025 will be used for Validation.
mypts <- as.data.frame(mypts)
# Filter the data for training (May 2021 - Dec 2023)
training_data <- mypts[mypts$Year %in% c(2021, 2022, 2023), ]
# Check training data
#head(training_data)# Filter the data for validation (Jan 2024 - July 2025)
validation_data <- mypts[mypts$Year %in% c(2024, 2025), ]
# Check validation data
#head(validation_data)Random Forest for generating variable of importance
Tune The Forest
The tuneRF function in the randomForest package is used to tune the mtry parameter, which is the number of variables randomly sampled as candidates at each split in the random forest. The function requires a data frame of predictor variables and a response variable.
# Convert training_data data to data frame
mypts_df <- as.data.frame(training_data)
trf <- tuneRF(x=mypts_df[,1:ncol(mypts_df)], # Prediction variables
y=mypts_df$pkg) # Response variable## mtry = 21 OOB error = 1863.607
## Searching left ...
## mtry = 11 OOB error = 8910.151
## -3.781133 0.05
## Searching right ...
## mtry = 42 OOB error = 188.1819
## 0.8990227 0.05
## mtry = 64 OOB error = 98.48041
## 0.4766745 0.05
Based on the output from tuneRF, you can observe that the mtry value that gives the lowest Out-of-Bag (OOB) error. To build the first random forest model, we will use this mtry value.
(mintree <- trf[which.min(trf[,2]),1])## [1] 64Fit The Random Forest Model (1)
Random Forest for generating variable of importance
Here, the model is fitted using the randomForest function to build a predictive model for food commodity prices. The model takes the response variable, the prediction variables and the optimal number of variables to consider at each split. The goal is to generate Variable importance scores which will help us understand which variables have the most significant impact on the response variable, enabling us to interpret the model and possibly simplify it by focusing on the most important predictors.
# Create the random forest model
rf1 <- randomForest(pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato +
Month +
Year +
ttcity_u5 + ttport_1 +
bio_1 + bio_2 + bio_3 + bio_4 + bio_5 + bio_6 +
bio_7 + bio_8 + bio_9 + bio_10 + bio_11 + bio_12 +
bio_13 + bio_14 + bio_15 + bio_16 + bio_17 +
bio_18 + bio_19 + MAI_ARE + MAI_YLD +
Longitude + Latitude +
rain.sum.lag,
data = training_data,mtry=mintree,
importance=TRUE,na.rm=TRUE)## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## rangerf1##
## Call:
## randomForest(formula = pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato + Month + Year + ttcity_u5 + ttport_1 + bio_1 + bio_2 + bio_3 + bio_4 + bio_5 + bio_6 + bio_7 + bio_8 + bio_9 + bio_10 + bio_11 + bio_12 + bio_13 + bio_14 + bio_15 + bio_16 + bio_17 + bio_18 + bio_19 + MAI_ARE + MAI_YLD + Longitude + Latitude + rain.sum.lag, data = training_data, mtry = mintree, importance = TRUE, na.rm = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 36
##
## Mean of squared residuals: 32910.31
## % Var explained: 93.69varImpPlot(rf1)
## evaluate
(oob <- sqrt(rf1$mse[which.min(rf1$mse)]))## [1] 181.0919This calculates the RMSE of the tree in the Random Forest model that has the lowest OOB error.
importance_metrics <- importance(rf1, type=1) # %IncMSE
impvar <- rownames(importance_metrics)[order(importance_metrics[, 1], decreasing=TRUE)]
# Get the top 20 variables
top_20_vars <- impvar[1:20]
top_20_vars## [1] "Year" "Month" "beans" "rice" "Longitude"
## [6] "sorghum" "fmillet" "wheat" "bmillet" "rain.sum.lag"
## [11] "bio_12" "bio_18" "MAI_YLD" "bio_15" "MAI_ARE"
## [16] "Latitude" "bio_3" "ttport_1" "bio_4" "ttcity_u5"node_purity <- importance(rf1, type=2) # IncNodePurity
# Sort variables by importance (IncNodePurity)
node_purity_sorted <- sort(node_purity[,1], decreasing = TRUE)
# Select the top 20 important variables
top_vars <- names(node_purity_sorted)[1:20]
print(top_vars)## [1] "rice" "beans" "Year" "Month" "wheat"
## [6] "fmillet" "potato" "bio_12" "maize" "Longitude"
## [11] "sorghum" "bmillet" "rain.sum.lag" "bio_3" "MAI_YLD"
## [16] "bio_7" "bio_18" "ttcity_u5" "Latitude" "MAI_ARE"rf1$importanceSD## maize rice sorghum bmillet fmillet wheat
## 1790.6975 2755.7963 457.2943 587.0560 1633.1918 1571.4940
## beans potato Month Year ttcity_u5 ttport_1
## 2970.6954 1785.4471 294.1134 527.9516 231.8652 189.1793
## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6
## 119.7947 230.9594 564.9766 150.5680 173.8270 281.5453
## bio_7 bio_8 bio_9 bio_10 bio_11 bio_12
## 696.0512 108.5006 200.2282 142.8832 127.6657 762.8554
## bio_13 bio_14 bio_15 bio_16 bio_17 bio_18
## 102.8958 126.9305 241.4402 121.1899 138.4660 234.0913
## bio_19 MAI_ARE MAI_YLD Longitude Latitude rain.sum.lag
## 176.7645 221.9445 429.8202 314.9773 270.0510 155.9511Fit The Random Forest Model (2)
Estimate more parsimonious specification
In this section, we aim to refine our model by selecting the most important variables. We will review the importance metrics (%IncMSE and IncNodePurity) to identify the variables that contribute the most to the model’s predictive power. Our focus will be on variables with higher importance values to ensure a more efficient and interpretable model.
### Fit The Random Forest Model (2)
rf <- randomForest(pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato +
Month + Year +
#latitude + longitude +
ttport_1 + ttcity_u5 + popdens +
bio_3 + bio_6 + bio_9 + bio_12 + bio_18 +
rain.sum.lag,
data = training_data,
ntree = 500,
importance = TRUE,
na.action = na.omit)
rf##
## Call:
## randomForest(formula = pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato + Month + Year + ttport_1 + ttcity_u5 + popdens + bio_3 + bio_6 + bio_9 + bio_12 + bio_18 + rain.sum.lag, data = training_data, ntree = 500, importance = TRUE, na.action = na.omit)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 6
##
## Mean of squared residuals: 36438.71
## % Var explained: 93.01# evaluate
varImpPlot(rf)
(oob <- sqrt(rf$mse[which.min(rf$mse)]))## [1] 190.8051partialPlot(rf, as.data.frame(training_data), "rain.sum.lag")
spatial prediction
These are prediction plots for each food commodities (maize, beans, rice, sorghum, bmillet, fmillet, wheat, potato) with their respective month of the year 2023.
year <- 2024Note: we must set the rain.sum.lag variables for each month we are predicting
Predicted Maize Prices
#Maize
# Create an empty list to store predictions for maize
predict_for_month <- function(month){
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")] # Remember to change depending on year
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_maize <- data.frame(
maize = 1, rice = 0, sorghum = 0, bmillet = 0, fmillet = 0, wheat = 0, beans = 0, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_maize, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_maize <- lapply(1:12, predict_for_month)
# Extract pixel values from predictions_maize
maize_values <- unlist(lapply(predictions_maize, values))
# Get min and max values
min_maize <- min(maize_values, na.rm = TRUE)
max_maize <- max(maize_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
break_interval <- 100
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
break_interval <- 100
# Create a 3x4 matrix of plots
par(mar = c(0, 0, 0, 0)) # Set margins to 0 for inner plots
for (i in 1:12) {
plot(predictions_maize[[i]], main = paste("Maize prices", toupper(i), year),
zlim = c(min_maize, max_maize), col = color_palette, breaks = seq(min_maize, max_maize, by = break_interval), legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_maize, max_maize), n = 4)
# Reset plot layout for the legend
layout(matrix(1))
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_maize, max_maize), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Maize Price (TZS Per Kg)", side = 1, line = 2, cex = 0.9))
Predicted Beans Prices
# Beans
# Function to predict beans for a given month
predict_for_beans <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_beans <- data.frame(
maize = 0, rice = 0, sorghum = 0, bmillet = 0, fmillet = 0, wheat = 0, beans = 1, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_beans, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_beans <- lapply(1:12, predict_for_beans)
# Extract pixel values from predictions_beans
bean_values <- unlist(lapply(predictions_beans, values))
min_bean <- min(bean_values, na.rm = TRUE)
max_bean <- max(bean_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
break_interval <- 150
# Loop through each month to plot beans prices
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_beans[[i]], main = paste("Beans prices", toupper(i), year),
zlim = c(min_bean, max_bean), col = color_palette, breaks = seq(min_bean, max_bean, by = break_interval), legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_bean, max_bean), n = 4)
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_bean, max_bean), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Beans Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Predicted Rice Prices
# Rice
# Function to predict rice prices for a given month
predict_for_rice <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_rice <- data.frame(
maize = 0, rice = 1, sorghum = 0, bmillet = 0, fmillet = 0, wheat = 0, beans = 0, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_rice, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_rice <- lapply(1:12, predict_for_rice)
# Extract pixel values from predictions_rice
rice_values <- unlist(lapply(predictions_rice, values))
min_rice <- min(rice_values, na.rm = TRUE)
max_rice <- max(rice_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
break_interval <- 150
# Loop through each month to plot rice prices
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_rice[[i]], main = paste("Rice prices", toupper(i), year),
zlim = c(min_rice, max_rice), col = color_palette, breaks = seq(min_rice, max_rice, by = break_interval), legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_rice, max_rice), n = 5)
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_rice, max_rice), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Rice Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Predicted Sorghum Prices
# Function to predict sorghum prices for a given month
predict_for_sorghum <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_sorghum <- data.frame(
maize = 0, rice = 0, sorghum = 1, bmillet = 0, fmillet = 0, wheat = 0, beans = 0, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_sorghum, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_sorghum <- lapply(1:12, predict_for_sorghum)
# Extract pixel values from predictions_sorghum
sorghum_values <- unlist(lapply(predictions_sorghum, values))
min_sorghum <- min(sorghum_values, na.rm = TRUE)
max_sorghum <- max(sorghum_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
break_interval <- 300
# Loop through each month to plot sorghum prices
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_sorghum[[i]], main = paste("Sorghum prices", toupper(i), year),
zlim = c(min_sorghum, max_sorghum), col = color_palette, breaks = seq(min_sorghum, max_sorghum, by = break_interval), legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_sorghum, max_sorghum), n = 5)
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_sorghum, max_sorghum), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Sorghum Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Predicted Bulrush Millet Prices
# bmillet
# Function to predict bmillet prices for a given month
predict_for_bmillet <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_bmillet <- data.frame(
maize = 0, rice = 0, sorghum = 0, bmillet = 1, fmillet = 0, wheat = 0, beans = 0, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_bmillet, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_bmillet <- lapply(1:12, predict_for_bmillet)
# Extract pixel values from predictions_bmillet
bmillet_values <- unlist(lapply(predictions_bmillet, values))
min_bmillet <- min(bmillet_values, na.rm = TRUE)
max_bmillet <- max(bmillet_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout with an extra row for the legend
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
# Loop through each month to plot bmillet prices
break_interval <- 150
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_bmillet[[i]], main = paste("Bmillet prices", toupper(i), year),
zlim = c(min_bmillet, max_bmillet), col = color_palette,
breaks = seq(min_bmillet, max_bmillet, by = break_interval),
legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_bmillet, max_bmillet), n = 5)
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_bmillet, max_bmillet), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Bulrush Millet Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Predicted Finger Millet Prices
#fmillet
# Function to predict fmillet prices for a given month
predict_for_fmillet <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_fmillet <- data.frame(
maize = 0, rice = 0, sorghum = 0, bmillet = 0, fmillet = 1, wheat = 0, beans = 0, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_fmillet, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_fmillet <- lapply(1:12, predict_for_fmillet)
# Extract pixel values from predictions_fmillet
fmillet_values <- unlist(lapply(predictions_fmillet, values))
min_fmillet <- min(fmillet_values, na.rm = TRUE)
max_fmillet <- max(fmillet_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout with an extra row for the legend
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
# Loop through each month to plot fmillet prices
break_interval <- 300
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_fmillet[[i]], main = paste("Fmillet prices", toupper(i), year),
zlim = c(min_fmillet, max_fmillet), col = color_palette,
breaks = seq(min_fmillet, max_fmillet, by = break_interval),
legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_fmillet, max_fmillet), n = 3)
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_fmillet, max_fmillet), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Finger Millet Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Predicted Wheat Prices
#Wheat
# Function to predict wheat prices for a given month
predict_for_wheat <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_wheat <- data.frame(
maize = 0, rice = 0, sorghum = 0, bmillet = 0, fmillet = 0, wheat = 1, beans = 0, potato = 0,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_wheat, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_wheat <- lapply(1:12, predict_for_wheat)
# Extract pixel values from predictions_wheat
wheat_values <- unlist(lapply(predictions_wheat, values))
min_wheat <- min(wheat_values, na.rm = TRUE)
max_wheat <- max(wheat_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
# Define the break interval for both plot and legend
break_interval <- 100
# Loop through each month to plot wheat prices
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_wheat[[i]], main = paste("Wheat prices", toupper(i), year),
zlim = c(min_wheat, max_wheat), col = color_palette,
breaks = seq(min_wheat, max_wheat, by = break_interval),
legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_wheat, max_wheat), n = 5)
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_wheat, max_wheat), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Wheat Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Predicted potato prices
#potatoes
# Function to predict potato prices for a given month
predict_for_potato <- function(month) {
rain_sum_lag <- rstack3[paste0("2024-", sprintf("%02d", month), "-01_rain.sum.lag")]
names(rain_sum_lag) <- "rain.sum.lag"
newstack <- c(rstack, rain_sum_lag)
const_potato <- data.frame(
maize = 0, rice = 0, sorghum = 0, bmillet = 0, fmillet = 0, wheat = 0, beans = 0, potato = 1,
Month = factor(month, levels = levels(training_data$Month)), # match factor levels
Year = year
)
pred <- predict(newstack, rf, const = const_potato, na.rm = TRUE)
pred <- mask(pred, tza1)
return(pred)
}
# Create predictions for all months
predictions_potato <- lapply(1:12, predict_for_potato)
# Extract pixel values from predictions_potato
potato_values <- unlist(lapply(predictions_potato, values))
min_potato <- min(potato_values, na.rm = TRUE)
max_potato <- max(potato_values, na.rm = TRUE)
# Define the continuous color palette and reverse it
color_palette <- colorRampPalette(c("blue", "wheat", "red"))(100)
# Set up plot layout with an extra row for the legend
layout_matrix <- matrix(c(1, 2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 13, 13, 13), nrow = 4, byrow = TRUE)
# Set up layout
layout(layout_matrix, heights = c(1, 1, 1, 0.5))
# Loop through each month to plot potato prices
break_interval <- 150
par(mar = c(0, 0, 0, 0))
for (i in 1:12) {
plot(predictions_potato[[i]], main = paste("Potato prices", toupper(i), year),
zlim = c(min_potato, max_potato), col = color_palette,
breaks = seq(min_potato, max_potato, by = break_interval),
legend = FALSE, axes = FALSE)
# Add region boundaries
plot(tza1, add = TRUE, border = "black", lwd = 0.1)
points(training_data, pch = 20, col = "red", cex = 0.5)
}
# Generate pretty breaks for the legend
legend_breaks <- pretty(c(min_potato, max_potato), n = 5) # Adjust n as needed
# Reset plot layout to 1x1 for the legend
layout(matrix(1))
# Set margins for the legend
par(mar = c(5, 4, 2, 1))
# Plot the legend in the bottom row, centered horizontally
image.plot(zlim = c(min_potato, max_potato), legend.only = TRUE,
col = color_palette, horizontal = TRUE,
legend.width = 0.7, legend.shrink = 0.9,
axis.args = list(at = legend_breaks, labels = legend_breaks, cex.axis = 0.9),
legend.args = list(text = "Predicted Potato Price (Tsh) Per Kg", side = 1, line = 2, cex = 0.9))
Cumulative probability plots
# List of crops excluding maize
crops <- c("beans", "rice", "sorghum", "bmillet", "fmillet", "wheat", "potato")
# Store CDF data for plotting
cdf_list <- list()
for (crop in crops) {
pred_list <- get(paste0("predictions_", crop)) # e.g., predictions_beans
df_crop <- data.frame()
for (month in 1:12) {
vals <- values(pred_list[[month]])
vals <- vals[!is.na(vals)]
vals_sorted <- sort(vals)
probs <- seq_along(vals_sorted) / length(vals_sorted)
df_month <- data.frame(
price = vals_sorted,
cum_prob = probs,
month = month,
crop = crop
)
df_crop <- rbind(df_crop, df_month)
}
cdf_list[[crop]] <- df_crop
}
cdf_data <- bind_rows(cdf_list)
cdf_data$month <- factor(cdf_data$month, labels = month.name)
# Define nicer labels
crop_labels <- c(
beans = "Bean",
bmillet = "Bulrush Millet",
fmillet = "Finger Millet",
potato = "Potato",
rice = "Rice",
sorghum = "Sorghum",
wheat = "Wheat"
)
# 12 colors
month_colors <- c("#E31A1C", "#6A3D9A", "#FF7F00", "#1F78B4", "#33A02C", "#FB9A99","#A6CEE3", "#B2DF8A", "#FDBF6F", "#CAB2D6", "#FFFF99", "#B15928")
ggplot(cdf_data, aes(x = price, y = cum_prob, color = month)) +
geom_line(size = 1) +
facet_wrap(~crop, scales = "free", ncol = 2,
labeller = labeller(crop = crop_labels)) +
scale_color_manual(values = month_colors) +
labs(x = "Predicted price (TZs/kg)", y = "Cumulative probability", color = "") +
theme_minimal(base_size = 12) +
theme(
legend.justification = c("right", "top"),
legend.title = element_text(size = 10),
legend.text = element_text(size = 10),
axis.text.x = element_text(angle = 0, vjust = 1, size = 10),
axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 10),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
legend.key.width = unit(1.5, "lines")
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave("cumulative_probability_plots1.png", width = 12, height = 10, dpi = 300)Variable importance plots
# Extract variable importance
imp <- importance(rf) %>% as.data.frame()
imp$Variable <- rownames(imp)
# Choose importance column
importance_col <- if ("IncMSE" %in% names(imp)) {
"IncMSE"
} else if ("%IncMSE" %in% names(imp)) {
"%IncMSE"
} else {
"IncNodePurity"
}
# Create nicer labels for publication
pretty_labels <- c(
maize = "Maize",
rice = "Rice",
sorghum = "Sorghum",
bmillet = "Bulrush millet",
fmillet = "Finger millet",
wheat = "Wheat",
beans = "Common bean",
potato = "Irish potato",
Month = "Month",
Year = "Year",
ttport_1 = "Travel time to port",
ttcity_u5 = "Travel time to secondary city",
popdens = "Population density",
bio_3 = "Isothermality (bio3)",
bio_6 = "Min temp coldest month (bio6)",
bio_9 = "Mean temp driest quarter (bio9)",
bio_12 = "Annual precipitation (bio12)",
bio_18 = "Precip warmest quarter (bio18)",
rain.sum.lag = "Lagged rainfall"
)
# Arrange and keep top 10
imp_plot <- imp %>%
arrange(desc(.data[[importance_col]])) %>%
mutate(
PrettyVar = recode(Variable, !!!pretty_labels),
PrettyVar = factor(PrettyVar, levels = PrettyVar)
)
# Plot
ggplot(imp_plot, aes(x = PrettyVar, y = .data[[importance_col]])) +
geom_col(fill = "steelblue") +
coord_flip() +
labs(
x = "Predictor",
y = "Importance",
title = ""
) +
theme_minimal(base_size = 14) +
theme(
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 10),
axis.title.x = element_text(size = 10),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5)
)
Prediction Evaluation
1. Using 2024/2025 Validation data
pred<-predict(object=rf, newdata=validation_data)
actual<-validation_data$pkg
result<-data.frame(actual=actual, predicted=pred)mse <- mean((actual - pred)^2, na.rm=TRUE)
paste('Mean Squared Error:', mse)## [1] "Mean Squared Error: 216507.550059363"rmse <- sqrt(mse)
paste('Root Mean Squared error: ',mean(sqrt(rf$mse)))## [1] "Root Mean Squared error: 194.21431712721"#Save predicted & observed price
write.csv(result, "result.csv")#reading result.csv file (predicted vs observed)
rslt <- read.csv("result.csv", header=T)
print(names(rslt))## [1] "X" "actual" "predicted"#RMSE predicting from rf - predicited vs observed
rf.rmse<-round(sqrt(mean( (rslt$actual-rslt$predicted)^2 , na.rm = TRUE )),2)
print(rf.rmse)## [1] 465.3#R-square
rf.r2<-round(summary(lm(actual~predicted, rslt))$r.squared,2)
print(rf.r2)## [1] 0.66range(actual)## [1] 327.5 4050.0range(pred)## [1] 637.0127 3378.3272#plotting predicted Vs observed
ggplot(result, aes(x=actual, y=predicted), alpha=0.6) +
geom_point(colour = "blue", size = 1.4, alpha=0.6) +
ggtitle('Random Forest "Wholesale Grain Prices in Tanzania"') +
scale_x_continuous("Observed Price (Tsh) Per Kg",
limits = c(0, 5000),
breaks = seq(0, 5000, 1000)) +
scale_y_continuous("Predicted Price (Tsh) Per Kg",
limits = c(0, 5000),
breaks = seq(0, 5000, 1000)) +
theme(axis.line = element_line(colour = "black"),
axis.text.y = element_text(size = 8, angle = 90, hjust = 0.5, vjust = 1),
axis.text.x = element_text(size = 8)) +
geom_abline(intercept = 0, slope = 1, linewidth = 0.5) +
geom_smooth(aes(x = actual, y = predicted), formula = y ~ x, method = "lm", se = FALSE, colour = "red", linetype = 2, size = 0.9) +
annotate("text", x = 300, y = 4500, label = paste("RMSE:", rf.rmse)) +
annotate("text", x = 300, y = 4200, label = paste("R^2: ", rf.r2), parse = TRUE)
2. Compare the observed Prices (the training data) with the predicted Prices (predicted using the training data) using stats package
library(stats)mypts_df$pred <- stats::predict(rf)rsq <- function (obs, pred) cor(obs, pred, use = 'complete.obs') ^ 2
RMSE <- function(obs, pred){sqrt(mean((pred - obs)^2, na.rm = TRUE))}fr2_rsq <- rsq(mypts_df$pkg, mypts_df$pred) %>% round(digits = 2)
fr2_rmse <- RMSE(mypts_df$pkg, mypts_df$pred) %>% round(digits = 0)Price_fit_plot <- ggplot(data = mypts_df, aes(x = pkg, y = pred)) +
geom_point(colour = "blue", size = 1.4 ,alpha=0.6) +
ggtitle('Observed vs Predicted "Wholesale Grain Prices in Tanzania"') +
geom_abline(slope = 1, alpha=0.3) +
annotate('text', x = 150, y = 4500, label = paste0("R^{2}==", fr2_rsq), parse = TRUE, size=3) +
annotate('text', x = 150, y = 4200, label = paste0("RMSE==", fr2_rmse), parse = TRUE, size=3) +
labs(x = "Observed Price (Tsh) Per Kg", y = "Predicted Price (Tsh) Per Kg") +
xlim(0, 5000) + ylim(0, 5000)
Price_fit_plot
Partial dependence plots
Plot Partial dependence of all the variables used except food commodities and months.
library(caret)
var_importance <- varImp(rf)
impvar <- rownames(var_importance)[order(var_importance[, 1], decreasing=TRUE)]op <- par(mfrow=c(2, 4))
# exclude food commodities and months
predictors_to_plot <- setdiff(impvar, c("maize", "rice", "sorghum", "bmillet", "fmillet", "wheat", "beans", "potato", "jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"))
for (i in seq_along(predictors_to_plot)) {
partialPlot(rf, as.data.frame(training_data), predictors_to_plot[i], xlab=predictors_to_plot[i],
main="Partial Dependence")
}
Descriptive statistics for each crop
# Compute descriptive statistics for each crop
stats <- mypts %>%
filter(pkg > 0) %>%
group_by(Crop) %>%
summarise(
Mean = mean(pkg, na.rm = TRUE),
Median = median(pkg, na.rm = TRUE),
Minimum = min(pkg, na.rm = TRUE),
Maximum = max(pkg, na.rm = TRUE),
Std_Dev = sd(pkg, na.rm = TRUE),
IQR = IQR(pkg, na.rm = TRUE),
Observations = n()
) %>%
arrange(Crop)
print(stats)## # A tibble: 8 × 8
## Crop Mean Median Minimum Maximum Std_Dev IQR Observations
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Maize 784. 750 268. 2388. 263. 389. 1216
## 2 Rice 2230. 2250 800 3518. 529. 814. 1227
## 3 Sorghum 1311. 1300 328. 3250 463. 650 1041
## 4 B.Millet 1421. 1300 469. 4000 587. 800 865
## 5 F.Millet 1784. 1756. 700 4000 501. 500 1087
## 6 Wheat 1839. 1811. 368. 4050 584. 646. 842
## 7 Beans 2489. 2550 985 3850 554. 876. 1223
## 8 Potato 899. 871. 300 2045. 266. 275 1201correlation Plots
These are correlation plots for pooled sample in two periods: post-harvest (May-Oct) and lean season (Nov-April).
#We'll use mean Monthly data in wide format
prices.monthly## Region Market Month Year Latitude Longitude mai.price
## <char> <char> <int> <int> <num> <num> <num>
## 1: Arusha Arusha 5 2021 -3.36968 36.68808 469.00
## 2: Dar es Salaam Temeke 5 2021 -6.85097 39.25672 485.00
## 3: Dodoma Majengo 5 2021 -6.17971 35.74109 501.50
## 4: Dodoma Kibaigwa 5 2021 -6.08110 36.64645 375.00
## 5: Kagera Bukoba 5 2021 -1.33140 31.81293 426.25
## ---
## 1229: Pwani Mnalani 7 2025 -6.44221 38.90622 1278.00
## 1230: Simiyu Bariadi 7 2025 -2.80355 33.98699 840.00
## 1231: Kigoma Kigoma 7 2025 -4.89697 29.65987 744.75
## 1232: Njombe Makambako 7 2025 -8.84398 34.82299 569.50
## 1233: Rukwa Sumbawanga 7 2025 -7.95239 31.62319 561.00
## ric.price sor.price bul.price fin.price whe.price bea.price pot.price
## <num> <num> <num> <num> <num> <num> <num>
## 1: 1530 695.0 761.0 1333.0 793.0 1545 745
## 2: 1650 900.0 900.0 1775.0 1230.0 2420 730
## 3: 1592 558.0 581.0 1752.0 NaN 2075 618
## 4: NaN 437.5 NaN NaN NaN NaN NaN
## 5: 1260 1375.0 1337.5 1637.5 1637.5 1300 675
## ---
## 1229: 2060 1900.0 1900.0 2150.0 2300.0 2900 1450
## 1230: 2000 1900.0 1900.0 2150.0 2300.0 2900 1450
## 1231: 2175 2325.0 1310.0 2710.0 2850.0 2100 960
## 1232: 2375 850.0 NaN 1625.0 1550.0 2150 1100
## 1233: 2150 NaN NaN 975.0 1080.0 2400 685# Create a function to determine the season
get_season <- function(month) {
if (month %in% 5:10) {
return("Post-Harvest")
} else {
return("Lean Season")
}
}# Add a 'Season' column to the data
prices.monthly$Season <- sapply(prices.monthly$Month, get_season)# Post Harvest Data
post_harvest_data <- prices.monthly[prices.monthly$Season == "Post-Harvest", ]
# Remove rows with NaNs from the Post-Harvest data
post_harvest_data <- post_harvest_data[complete.cases(post_harvest_data[, 7:14]), ]
# Lean Season data
lean_season_data <- prices.monthly[prices.monthly$Season == "Lean Season", ]
# Remove rows with NaNs from the Lean Season data
lean_season_data <- lean_season_data[complete.cases(lean_season_data[, 7:14]), ]# Calculate correlation matrix for Post-Harvest season
post_harvest_corr <- cor(post_harvest_data[, 7:14])
# Calculate correlation matrix for Lean Season
lean_season_corr <- cor(lean_season_data[, 7:14])# Plot correlation matrix for Post-Harvest season
corrplot(post_harvest_corr,
method = "color",
title = "",
tl.col = "black",
tl.cex = 0.5,
addCoef.col = "black",
number.cex = 0.5,
number.digits = 2)
# Add title to the plot
title(main = "Post-Harvest Correlation Matrix",
line = 3,
cex.main = 0.9)
# Plot correlation matrix for Lean Season
corrplot(lean_season_corr,
method = "color",
title = "",
tl.col = "black",
tl.cex = 0.5,
addCoef.col = "black",
number.cex = 0.5,
number.digits = 2)
# Add title to the plot
title(main = "Lean Season Correlation Matrix",
line = 3,
cex.main = 0.9) 
Comparing Pooled Vs Crop Specific Price Predictions Using different validation methods
We are comparing pooled vs Crop Specific Price Prediction models to determine which model performs better at prediction. This is done by comparing their prediction’s r2 or rmse respectively.
“train-test split” validation
Here the dataset is split once (70:30). The model is trained on one part of the data (70%) and evaluated on the remaining set of the data (30%).
# Pooled RF model
##Split the data into train and test datasets
set.seed(1983)
mypts <- as.data.frame(mypts)
rows <- sample(x=1:nrow(mypts), size = 0.70* nrow(mypts))
train <- mypts[rows, ]
test <- mypts[! rownames(mypts) %in% rows, ]
# Hyperparameter tuning
x_vars1 <- train %>%
dplyr::select(-pkg)
y_var1 <- train$pkg
trf1 <- tuneRF(x = x_vars1,
y = y_var1,
stepFactor = 1.5,
improve = 0.01,
ntreeTry = 500,
trace = TRUE,
plot = TRUE)## mtry = 21 OOB error = 58568.45
## Searching left ...
## mtry = 14 OOB error = 66760.4
## -0.1398697 0.01
## Searching right ...
## mtry = 31 OOB error = 53474.18
## 0.08697979 0.01
## mtry = 46 OOB error = 51568.84
## 0.03563095 0.01
## mtry = 63 OOB error = 51022.28
## 0.0105987 0.01
# number of valid predictors (excluding the target variable)
n_predictors1 <- length(setdiff(names(train), c("pkg")))
n_predictors1## [1] 63# Extract the mtry value that gave the lowest OOB error
best_mtry1 <- trf1[which.min(trf1[, 2]), 1]
best_mtry1## [1] 63# Ensure best_mtry is within valid range
best_mtry1 <- min(best_mtry1, n_predictors1)
best_mtry1## [1] 63RF <- randomForest(pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato +
Month + Year +
#latitude + longitude +
ttport_1 + ttcity_u5 + popdens +
bio_3 + bio_6 + bio_9 + bio_12 + bio_18 +
rain.sum.lag,
data = train,
mtry = best_mtry1,
ntree = 500,
importance = TRUE,
na.action = na.omit)## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## rangeRF##
## Call:
## randomForest(formula = pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato + Month + Year + ttport_1 + ttcity_u5 + popdens + bio_3 + bio_6 + bio_9 + bio_12 + bio_18 + rain.sum.lag, data = train, mtry = best_mtry1, ntree = 500, importance = TRUE, na.action = na.omit)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 19
##
## Mean of squared residuals: 50532.19
## % Var explained: 91.15evaluate_models <- function(crop) {
# Filter data for the specific crop
crop_data <- mypts[mypts$Crop == crop, ]
set.seed(1983)
rownames(crop_data)<-1:nrow(crop_data)
rows <- sample(x=1:nrow(crop_data), size = 0.70* nrow(crop_data))
training_data_crop <- crop_data[rows, ]
validation_data_crop <- crop_data[! rownames(crop_data) %in% rows, ]
# Predictions for pooled model on specific crop data
pooled_pred_crop <- predict(RF, newdata = validation_data_crop)
actual_pooled <- validation_data_crop$pkg
result_pooled <-data.frame(actual=actual_pooled, predicted=pooled_pred_crop)
# Crop-Specific Model
rf_crop <- randomForest(pkg ~ Month + Year +
#latitude + longitude +
ttport_1 + ttcity_u5 + popdens +
bio_3 + bio_6 + bio_9 + bio_12 + bio_18 +
rain.sum.lag,
data = training_data_crop,
mtry = best_mtry1,
ntree = 500,
importance = TRUE,
na.action = na.omit)
# Predictions for crop-specific model
predictions_crop <- predict(rf_crop, newdata = validation_data_crop)
actual_crop <- validation_data_crop$pkg
result_crop_specific <- data.frame(actual=actual_crop, predicted=predictions_crop)
# Calculate performance metrics
rmse_pooled <- round(sqrt(mean((result_pooled$actual-result_pooled$predicted)^2 , na.rm = TRUE )),2)
r2_pooled <- round(summary(lm(actual~predicted, result_pooled))$r.squared,2)
rmse_crop <- round(sqrt(mean((result_crop_specific$actual-result_crop_specific$predicted)^2 , na.rm = TRUE )),2)
r2_crop <- round(summary(lm(actual~predicted, result_crop_specific))$r.squared,2)
return(data.frame(Crop = crop,
Model = c("Pooled", "Crop-Specific"),
RMSE = c(rmse_pooled, rmse_crop),
R_squared = c(r2_pooled, r2_crop)))
}# Apply function to all crops
crop <- unique(mypts$Crop)
comparison_df3 <- do.call(rbind, lapply(crop, evaluate_models))## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## range
## Warning in randomForest.default(m, y, ...): invalid mtry: reset to within valid
## rangecomparison_df3## Crop Model RMSE R_squared
## 1 Maize Pooled 67.75 0.93
## 2 Maize Crop-Specific 120.02 0.77
## 3 Rice Pooled 140.25 0.93
## 4 Rice Crop-Specific 207.79 0.84
## 5 Sorghum Pooled 146.59 0.91
## 6 Sorghum Crop-Specific 228.74 0.77
## 7 B.Millet Pooled 172.88 0.91
## 8 B.Millet Crop-Specific 312.89 0.66
## 9 F.Millet Pooled 181.18 0.89
## 10 F.Millet Crop-Specific 282.75 0.71
## 11 Wheat Pooled 182.76 0.92
## 12 Wheat Crop-Specific 317.81 0.68
## 13 Beans Pooled 150.27 0.92
## 14 Beans Crop-Specific 204.52 0.86
## 15 Potato Pooled 107.88 0.84
## 16 Potato Crop-Specific 164.05 0.62comparison_df_wide3 <- comparison_df3 %>%
pivot_wider(names_from = Model, values_from = c(RMSE, R_squared)) %>%
rename(RMSE_Crop_Specific = `RMSE_Crop-Specific`,
RMSE_Pooled = `RMSE_Pooled`,
R_squared_Crop_Specific = `R_squared_Crop-Specific`,
R_squared_Pooled = `R_squared_Pooled`)knitr::kable(comparison_df_wide3, caption = "Model Comparison: R² and RMSE Based on Train-Test Split Validation")| Crop | RMSE_Pooled | RMSE_Crop_Specific | R_squared_Pooled | R_squared_Crop_Specific |
|---|---|---|---|---|
| Maize | 67.75 | 120.02 | 0.93 | 0.77 |
| Rice | 140.25 | 207.79 | 0.93 | 0.84 |
| Sorghum | 146.59 | 228.74 | 0.91 | 0.77 |
| B.Millet | 172.88 | 312.89 | 0.91 | 0.66 |
| F.Millet | 181.18 | 282.75 | 0.89 | 0.71 |
| Wheat | 182.76 | 317.81 | 0.92 | 0.68 |
| Beans | 150.27 | 204.52 | 0.92 | 0.86 |
| Potato | 107.88 | 164.05 | 0.84 | 0.62 |
write.csv(comparison_df_wide3, "Train-Test-Validation-results.csv", row.names = FALSE)Temporal Split validation
Function to evaluate models for each crop using 2024 Data as the validation data
Here we are using 2024 and 2025 data for validation.
# Comparison between Pooled model and Crop Specific Model
set.seed(1983)
evaluate_models <- function(crop) {
# Filter data for the specific crop
crop_data <- mypts[mypts$Crop == crop, ]
training_data_crop <- crop_data[crop_data$Year %in% c(2021, 2022, 2023), ]
training_data_crop <- as.data.frame(training_data_crop)
validation_data_crop <- crop_data[crop_data$Year %in% c(2024, 2025), ]
validation_data_crop <- as.data.frame(validation_data_crop)
# Predictions for pooled model on specific crop data
pooled_pred_crop <- predict(rf, newdata = validation_data_crop)
actual_pooled <- validation_data_crop$pkg
result_pooled <-data.frame(actual=actual_pooled, predicted=pooled_pred_crop)
# Crop-Specific RF Model
rf_crop <- randomForest(pkg ~ Month + Year + ttport_1 + ttcity_u5 + popdens + bio_3 + bio_6 + bio_9 + bio_12 + bio_18 + rain.sum.lag,
data = training_data_crop, na.rm = TRUE)
# Predictions for crop-specific model
predictions_crop <- predict(rf_crop, newdata = validation_data_crop)
actual_crop <- validation_data_crop$pkg
result_crop_specific <- data.frame(actual=actual_crop, predicted=predictions_crop)
# Calculate performance metrics
rmse_pooled <- round(sqrt(mean((result_pooled$actual-result_pooled$predicted)^2 , na.rm = TRUE )),2)
r2_pooled <- round(summary(lm(actual~predicted, result_pooled))$r.squared,2)
rmse_crop <- round(sqrt(mean((result_crop_specific$actual-result_crop_specific$predicted)^2 , na.rm = TRUE )),2)
r2_crop <- round(summary(lm(actual~predicted, result_crop_specific))$r.squared,2)
return(data.frame(Crop = crop,
Model = c("Pooled", "Crop-Specific"),
RMSE = c(rmse_pooled, rmse_crop),
R_squared = c(r2_pooled, r2_crop)))
}# Apply the function to all crops
crop <- unique(mypts$Crop)
set.seed(1983)
comparison_df <- do.call(rbind, lapply(crop, evaluate_models))
comparison_df## Crop Model RMSE R_squared
## 1 Maize Pooled 325.34 0.21
## 2 Maize Crop-Specific 308.30 0.20
## 3 Rice Pooled 554.02 0.20
## 4 Rice Crop-Specific 534.72 0.20
## 5 Sorghum Pooled 370.53 0.40
## 6 Sorghum Crop-Specific 373.60 0.40
## 7 B.Millet Pooled 562.32 0.25
## 8 B.Millet Crop-Specific 562.76 0.26
## 9 F.Millet Pooled 620.49 0.10
## 10 F.Millet Crop-Specific 622.74 0.09
## 11 Wheat Pooled 538.15 0.31
## 12 Wheat Crop-Specific 539.04 0.30
## 13 Beans Pooled 345.19 0.35
## 14 Beans Crop-Specific 343.04 0.35
## 15 Potato Pooled 318.71 0.00
## 16 Potato Crop-Specific 316.85 0.01#write.csv(comparison_df, "model_comparison_df.csv")#comparison_df <- read.csv("model_comparison_df.csv")
comparison_df_wide <- comparison_df %>%
pivot_wider(names_from = Model, values_from = c(RMSE, R_squared)) %>%
rename(
RMSE_Crop_Specific = `RMSE_Crop-Specific`,
RMSE_Pooled = `RMSE_Pooled`,
R_squared_Crop_Specific = `R_squared_Crop-Specific`,
R_squared_Pooled = `R_squared_Pooled`
) %>%
group_by(Crop) %>%
summarize(
RMSE_Pooled = max(RMSE_Pooled, na.rm = TRUE),
RMSE_Crop_Specific = max(RMSE_Crop_Specific, na.rm = TRUE),
R_squared_Pooled = max(R_squared_Pooled, na.rm = TRUE),
R_squared_Crop_Specific = max(R_squared_Crop_Specific, na.rm = TRUE)
)knitr::kable(comparison_df_wide, caption = "Model Comparison: R² and RMSE")| Crop | RMSE_Pooled | RMSE_Crop_Specific | R_squared_Pooled | R_squared_Crop_Specific |
|---|---|---|---|---|
| Maize | 325.34 | 308.30 | 0.21 | 0.20 |
| Rice | 554.02 | 534.72 | 0.20 | 0.20 |
| Sorghum | 370.53 | 373.60 | 0.40 | 0.40 |
| B.Millet | 562.32 | 562.76 | 0.25 | 0.26 |
| F.Millet | 620.49 | 622.74 | 0.10 | 0.09 |
| Wheat | 538.15 | 539.04 | 0.31 | 0.30 |
| Beans | 345.19 | 343.04 | 0.35 | 0.35 |
| Potato | 318.71 | 316.85 | 0.00 | 0.01 |
write.csv(comparison_df_wide, "Temporal-split-results.csv", row.names = FALSE)Autoregressive Random Forest
This approach involves incorporating lagged Prices as predictors to capture temporal dependencies in the data.
mypts <- mypts %>%
arrange(Crop, Market, Year, Month) %>%
group_by(Crop, Market) %>%
mutate(pkg_lag1 = lag(pkg, 1), # 1-month lag
pkg_lag2 = lag(pkg, 2), # 2-month lag
pkg_lag3 = lag(pkg, 3)) %>% # 3-month lag
ungroup()
mypts <- na.omit(mypts)# Filter the data for training (May 2021 - Dec 2023)
training_data <- mypts[mypts$Year %in% c(2021, 2022, 2023), ]
# Check training data
#head(training_data)# Filter the data for validation (Jan 2024 - July 2025)
validation_data <- mypts[mypts$Year %in% c(2024, 2025), ]
# Check validation data
#head(validation_data)rf_ar <- randomForest(pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans + potato +
Month + Year +
ttport_1 + ttcity_u5 + popdens +
bio_3 + bio_6 + bio_9 + bio_12 + bio_18 +
rain.sum.lag +
pkg_lag1 + pkg_lag2 + pkg_lag3, # Adding lagged prices
data=training_data, na.rm=TRUE)
varImpPlot(rf_ar)
evaluate_models_ar <- function(crop) {
crop_data <- mypts[mypts$Crop == crop, ]
training_data_crop <- crop_data[crop_data$Year %in% c(2021, 2022, 2023), ]
validation_data_crop <- crop_data[crop_data$Year %in% c(2024, 2025), ]
# Predictions for pooled model with lags
pooled_pred_crop <- predict(rf_ar, newdata = validation_data_crop)
actual_pooled <- validation_data_crop$pkg
result_pooled <- data.frame(actual=actual_pooled, predicted=pooled_pred_crop)
# Train crop-specific model with lags
rf_crop_ar <- randomForest(pkg ~ Month + Year + ttport_1 + ttcity_u5 + popdens +
bio_3 + bio_6 + bio_9 + bio_12 + bio_18 + rain.sum.lag +
pkg_lag1 + pkg_lag2 + pkg_lag3, # Adding lagged prices
data = training_data_crop, na.rm = TRUE)
# Predictions for crop-specific model with lags
predictions_crop <- predict(rf_crop_ar, newdata = validation_data_crop)
actual_crop <- validation_data_crop$pkg
result_crop_specific <- data.frame(actual=actual_crop, predicted=predictions_crop)
# Calculate performance metrics
rmse_pooled <- round(sqrt(mean((result_pooled$actual - result_pooled$predicted)^2, na.rm = TRUE)), 2)
r2_pooled <- round(summary(lm(actual ~ predicted, result_pooled))$r.squared, 2)
rmse_crop <- round(sqrt(mean((result_crop_specific$actual - result_crop_specific$predicted)^2, na.rm = TRUE)), 2)
r2_crop <- round(summary(lm(actual ~ predicted, result_crop_specific))$r.squared, 2)
return(data.frame(Crop = crop,
Model = c("Pooled_AR", "Crop-Specific_AR"),
RMSE = c(rmse_pooled, rmse_crop),
R_squared = c(r2_pooled, r2_crop)))
}# Run model evaluation for all crops
set.seed(1983)
comparison_df_ar <- do.call(rbind, lapply(unique(mypts$Crop), evaluate_models_ar))
print(comparison_df_ar)## Crop Model RMSE R_squared
## 1 Maize Pooled_AR 190.25 0.44
## 2 Maize Crop-Specific_AR 175.28 0.47
## 3 Rice Pooled_AR 231.90 0.71
## 4 Rice Crop-Specific_AR 268.08 0.69
## 5 Sorghum Pooled_AR 254.32 0.72
## 6 Sorghum Crop-Specific_AR 283.02 0.69
## 7 B.Millet Pooled_AR 269.39 0.84
## 8 B.Millet Crop-Specific_AR 393.85 0.71
## 9 F.Millet Pooled_AR 318.01 0.76
## 10 F.Millet Crop-Specific_AR 452.56 0.60
## 11 Wheat Pooled_AR 331.52 0.75
## 12 Wheat Crop-Specific_AR 371.37 0.71
## 13 Beans Pooled_AR 220.59 0.71
## 14 Beans Crop-Specific_AR 224.35 0.71
## 15 Potato Pooled_AR 158.59 0.66
## 16 Potato Crop-Specific_AR 188.10 0.58#comparison_df <- read.csv("model_comparison_df.csv")
comparison_df_wide_ar <- comparison_df_ar %>%
pivot_wider(names_from = Model, values_from = c(RMSE, R_squared)) %>%
rename(
RMSE_Crop_Specific = `RMSE_Crop-Specific_AR`,
RMSE_Pooled = `RMSE_Pooled_AR`,
R_squared_Crop_Specific = `R_squared_Crop-Specific_AR`,
R_squared_Pooled = `R_squared_Pooled_AR`
) %>%
group_by(Crop) %>%
summarize(
RMSE_Pooled = max(RMSE_Pooled, na.rm = TRUE),
RMSE_Crop_Specific = max(RMSE_Crop_Specific, na.rm = TRUE),
R_squared_Pooled = max(R_squared_Pooled, na.rm = TRUE),
R_squared_Crop_Specific = max(R_squared_Crop_Specific, na.rm = TRUE)
)knitr::kable(comparison_df_wide_ar, caption = "Model Comparison: R² and RMSE for Autoregressive Random Forest")| Crop | RMSE_Pooled | RMSE_Crop_Specific | R_squared_Pooled | R_squared_Crop_Specific |
|---|---|---|---|---|
| Maize | 190.25 | 175.28 | 0.44 | 0.47 |
| Rice | 231.90 | 268.08 | 0.71 | 0.69 |
| Sorghum | 254.32 | 283.02 | 0.72 | 0.69 |
| B.Millet | 269.39 | 393.85 | 0.84 | 0.71 |
| F.Millet | 318.01 | 452.56 | 0.76 | 0.60 |
| Wheat | 331.52 | 371.37 | 0.75 | 0.71 |
| Beans | 220.59 | 224.35 | 0.71 | 0.71 |
| Potato | 158.59 | 188.10 | 0.66 | 0.58 |
write.csv(comparison_df_wide_ar, "Temporal-split-results_AR.csv", row.names = FALSE)Validation and robustness checks plots
# Relabel Validation more explicitly
df1 <- comparison_df_wide3 %>% mutate(Validation = "Train–Test Split")
df2 <- comparison_df_wide %>% mutate(Validation = "Temporal Split (standard RF)")
df3 <- comparison_df_wide_ar %>% mutate(Validation = "Temporal Split (autoregressive RF)")
# Combine into one dataset
all_df <- bind_rows(df1, df2, df3)
# Reshape to long for RMSE plotting
rmse_long <- all_df %>%
dplyr::select(Crop, Validation, RMSE_Pooled, RMSE_Crop_Specific) %>%
pivot_longer(cols = starts_with("RMSE"),
names_to = "Model", values_to = "RMSE") %>%
mutate(Model = recode(Model,
"RMSE_Pooled" = "Pooled",
"RMSE_Crop_Specific" = "Crop-Specific"))
# Plot RMSE
p1 <- ggplot(rmse_long, aes(x = Crop, y = RMSE, fill = Model)) +
geom_col(position = "dodge") +
facet_wrap(~Validation, ncol = 1, scales = "free_y") +
theme_minimal(base_size = 14) +
labs(y = "RMSE (TZS/kg)", x = "Crop") +
theme(
legend.position = "bottom",
legend.title = element_blank(),
axis.text.x = element_text(angle = 0, hjust = 1, size = 10),
axis.text.y = element_text(size = 10),
axis.title.x = element_text(size = 10),
axis.title.y = element_text(size = 10),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5)
)
# For R²: reshape
r2_long <- all_df %>%
dplyr::select(Crop, Validation, R_squared_Pooled, R_squared_Crop_Specific) %>%
pivot_longer(cols = starts_with("R_squared"),
names_to = "Model", values_to = "R2") %>%
mutate(Model = recode(Model,
"R_squared_Pooled" = "Pooled",
"R_squared_Crop_Specific" = "Crop-Specific"))
# Plot R²
p2 <- ggplot(r2_long, aes(x = Crop, y = R2, color = Model, group = Model)) +
geom_point(size = 3) +
geom_line() +
facet_wrap(~Validation, ncol = 1, scales = "free_y") +
theme_minimal(base_size = 14) +
labs(y = expression(R^2), x = "Crop") +
theme(
legend.position = "bottom",
legend.title = element_blank(),
axis.text.x = element_text(angle = 0, hjust = 1, size = 10),
axis.text.y = element_text(size = 10),
axis.title.x = element_text(size = 10),
axis.title.y = element_text(size = 10),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5)
)
# Combine panels (RMSE + R² side by side)
library(patchwork)## Warning: package 'patchwork' was built under R version 4.3.3##
## Attaching package: 'patchwork'## The following object is masked from 'package:MASS':
##
## area## The following object is masked from 'package:raster':
##
## area## The following object is masked from 'package:terra':
##
## areacombined_plot <- p1 | p2
combined_plot
Market Holdout Validation Analysis
This analysis evaluates the spatial generalizability of a random forest model for predicting crop prices (pkg) across 44 markets. Using a leave-N-markets-out cross-validation approach, we assess how model performance (measured by R²) changes as we incrementally increase the number of markets held out from model training.
For each number of held-out markets (num_holdout from 1 to 43), the process is repeated 150 times with different random market combinations. A random forest is trained on the remaining markets and evaluated on the held-out ones. Average and standard deviation of R² values across repetitions are computed for each holdout size.
A segmented (piecewise) regression is then fitted to the resulting data to visualize trends in prediction accuracy as market coverage declines.
#registerDoParallel(cores = parallel::detectCores() - 1)
#set.seed(1983)
#unique_markets <- unique(mypts$Market)
#total_markets <- length(unique_markets)
#market_holdout_results <- list()
#for (num_holdout in 1:(total_markets - 1)) {
# r2_vals <- foreach(rep = 1:150, .combine = c, .packages = c("randomForest", #"tidyverse")) %dopar% {
# heldout_markets <- sample(unique_markets, num_holdout)
# train_subset <- mypts %>% filter(!Market %in% heldout_markets)
# test_subset <- mypts %>% filter(Market %in% heldout_markets)
# model_rf <- randomForest(pkg ~ maize + rice + sorghum + bmillet + fmillet + wheat + beans +
# Month + Year +
# latitude + longitude +
# ttport_1 + ttcity_u5 + popdens +
# bio_3 + bio_6 + bio_9 + bio_12 + bio_18 +
# rain.sum.lag,
# data = train_subset,
# ntree = 500,
# importance = TRUE,
# na.action = na.omit)
#
# preds <- predict(model_rf, newdata = test_subset)
# actual_vals <- test_subset$pkg
# if (length(unique(actual_vals)) > 1) {
# return(summary(lm(actual_vals ~ preds))$r.squared)
# } else {
# return(NA_real_)
# }
# }
# market_holdout_results[[num_holdout]] <- data.frame(
# Markets_Heldout = num_holdout,
# Avg_R2 = mean(r2_vals, na.rm = TRUE),
# SD_R2 = sd(r2_vals, na.rm = TRUE)
# )
# cat("Done with num_holdout =", num_holdout, "\n")
#}#market_holdout_df <- bind_rows(market_holdout_results)
#market_holdout_df#write.csv(market_holdout_df, "market_holdout_df2.csv", row.names = FALSE)# To better show the downward trend as you hold out more markets
# Fit linear model
#lm_fit <- lm(Avg_R2 ~ Markets_Heldout, data = market_holdout_df)
# Fit segmented (piecewise) model with breakpoint estimation
#seg_fit <- segmented(lm_fit, seg.Z = ~Markets_Heldout)
# Add predicted values to data
#market_holdout_df$seg_fit <- predict(seg_fit)# Plot
#ggplot(market_holdout_df, aes(x = Markets_Heldout, y = Avg_R2)) +
# geom_ribbon(aes(ymin = Avg_R2 - SD_R2, ymax = Avg_R2 + SD_R2), fill = "grey", alpha = 0.4) +
# geom_line(aes(y = seg_fit), color = "black", size = 1) +
# geom_point(color = "black") +
# labs(title = "",
# x = "Number of Markets Held Out",
# y = "Mean R²") +
# theme_minimal() +
# theme(
# axis.line = element_line(color = "black", size = 0.6),
# axis.line.x = element_line(color = "black", size = 0.6),
# axis.line.y = element_line(color = "black", size = 0.6),
# panel.border = element_blank()
# )4.3. Comparison of national, regional, and predicted rural maize prices
dta_reg <- mypts %>%
filter(Crop == "Maize", Month == 8, Year == 2024) %>%
group_by(Region_GADM) %>%
summarise(Reg_max_july_pkg = max(pkg, na.rm = TRUE))
head(dta_reg)## # A tibble: 6 × 2
## Region_GADM Reg_max_july_pkg
## <chr> <dbl>
## 1 Arusha 639.
## 2 Dar es Salaam 1000
## 3 Dodoma 532.
## 4 Geita 650
## 5 Iringa 471.
## 6 Kagera 662.We use Dar es Salaam price as the national price bench mark
# Assign national Average (Dar es Salaam price) value to Raster
# National Average value for maize in August is 1000
# create reference raster
tza_extent <- ext(tza1) |> floor()
rnat <- crop(rast(res=1/12), tza_extent)
rnat <- project(rnat, crs(tza1))
values(rnat) <- 1000
rnatavg <- terra::mask(rnat, tza1)
plot(rnatavg) 
names(rnatavg) <- "national_avg"# Njombes price is missing, use Iringa's price
# Extract Iringa's price
njombe_price <- dta_reg %>%
filter(Region_GADM == "Iringa") %>%
mutate(Region_GADM = "Njombe")
njombe_price## # A tibble: 1 × 2
## Region_GADM Reg_max_july_pkg
## <chr> <dbl>
## 1 Njombe 471.# Add to regional summary
dta_reg <- bind_rows(dta_reg, njombe_price)
filter(dta_reg, Region_GADM == "Njombe")## # A tibble: 1 × 2
## Region_GADM Reg_max_july_pkg
## <chr> <dbl>
## 1 Njombe 471.#spatial join with the tza1
tza1_with_reg_prices <- merge(tza1, dta_reg, by.x = "NAME_1", by.y = "Region_GADM", all.x = TRUE)
head(tza1_with_reg_prices)## NAME_1 GID_1 GID_0 COUNTRY VARNAME_1 NL_NAME_1 TYPE_1
## 1 Arusha TZA.1_1 TZA Tanzania <NA> <NA> Mkoa
## 2 Dar es Salaam TZA.2_1 TZA Tanzania Dar-es-salaam|Pwani <NA> Mkoa
## 3 Dodoma TZA.3_1 TZA Tanzania <NA> <NA> Mkoa
## 4 Geita TZA.4_1 TZA Tanzania <NA> <NA> Mkoa
## 5 Iringa TZA.5_1 TZA Tanzania <NA> <NA> Mkoa
## 6 Kagera TZA.6_1 TZA Tanzania <NA> <NA> Mkoa
## ENGTYPE_1 CC_1 HASC_1 ISO_1 Reg_max_july_pkg
## 1 Region 02 TZ.AS TZ-01 639.0909
## 2 Region 07 TZ.DS TZ-02 1000.0000
## 3 Region 01 TZ.DO TZ-03 531.8182
## 4 Region 25 TZ.GE TZ-27 650.0000
## 5 Region 11 TZ.IG TZ-04 471.3636
## 6 Region 18 TZ.KG TZ-05 661.8182unique(tza1_with_reg_prices$NAME_1)## [1] "Arusha" "Dar es Salaam" "Dodoma" "Geita"
## [5] "Iringa" "Kagera" "Kaskazini Pemba" "Kaskazini Unguja"
## [9] "Katavi" "Kigoma" "Kilimanjaro" "Kusini Pemba"
## [13] "Kusini Unguja" "Lindi" "Manyara" "Mara"
## [17] "Mbeya" "Mjini Magharibi" "Morogoro" "Mtwara"
## [21] "Mwanza" "Njombe" "Pwani" "Rukwa"
## [25] "Ruvuma" "Shinyanga" "Simiyu" "Singida"
## [29] "Songwe" "Tabora" "Tanga"#plot(tza1_with_reg_prices)
#text(tza1_with_reg_prices, label="NAME_1")# convert the shapefile with the price data to a raster
regional_raster <- rasterize(tza1_with_reg_prices, rast(res=1/12, extent=ext(tza1), crs=crs(tza1)), field = "Reg_max_july_pkg")
regional_raster <- resample(regional_raster, r)
plot(regional_raster)
# Bring in predicted june prices for maize 2024
pred_maize <- rast("H:/Tanzania Price data/Datasets/Pred-plots/Maize/maize_price_rf_pred_08.tif")
pred_maize## class : SpatRaster
## dimensions : 144, 144, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : 29, 41, -12, 0 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : maize_price_rf_pred_08.tif
## name : lyr1
## min value : 763.1077
## max value : 1201.0100pred_maize <- project(pred_maize, crs(tza1))
names(pred_maize) <- "pred_maize"# 3 plots in one row
par(mfrow = c(1, 3), mar = c(4, 4, 4, 1))
color_palette <- brewer.pal(9, "YlGnBu")
# Dar es Salaam Price Map
plot(rnatavg,
col = color_palette,
main = "Dar es Salaam Price",
axes = TRUE, box = TRUE)
plot(tza1, add = TRUE, border = "gray30", lwd = 0.5)
# -------------------------------------------
# Regional Prices
plot(regional_raster,
col = color_palette,
main = "Regional Prices",
axes = TRUE, box = TRUE)
plot(tza1, add = TRUE, border = "gray30", lwd = 0.5)
# -------------------------------------------
# Predicted Prices Map
plot(pred_maize,
col = color_palette,
main = "Predicted Prices",
axes = TRUE, box = TRUE)
plot(tza1, add = TRUE, border = "gray30", lwd = 0.5)
# Pairwise comparison
# Extract values from rasters
pred_values <- values(pred_maize)
nat_values <- values(rnatavg)
reg_values <- values(regional_raster)
# Identify rows with all finite values
valid_idx <- which(is.finite(pred_values) & is.finite(nat_values) & is.finite(reg_values))
# Filter vectors
pred_values_f <- pred_values[valid_idx]
nat_values_f <- nat_values[valid_idx]
reg_values_f <- reg_values[valid_idx]
#dist_values_f <- dist_values[valid_idx]# Create a data frame from the filtered values
comparison_df <- data.frame(
Predicted = pred_values_f,
National_Avg = nat_values_f,
Regional_Avg = reg_values_f
)
# Define pairwise differences
comparison_df <- comparison_df %>%
mutate(
Pred_vs_Nat = National_Avg - Predicted,
Pred_vs_Reg = Regional_Avg - Predicted,
#pred_vs_dist = Predicted - District_Avg
)
head(comparison_df)## Predicted National_Avg Regional_Avg Pred_vs_Nat Pred_vs_Reg
## 1 996.7144 1000 661.8182 3.2855835 -334.8962
## 2 979.4234 1000 661.8182 20.5765991 -317.6052
## 3 992.1823 1000 661.8182 7.8177490 -330.3641
## 4 1005.4178 1000 661.8182 -5.4178467 -343.5997
## 5 1000.8348 1000 661.8182 -0.8347778 -339.0166
## 6 1001.3495 1000 661.8182 -1.3495483 -339.5314# % of predicted prices compared to national and regional
mean(comparison_df$Pred_vs_Nat > 0) ## [1] 0.7849723mean(comparison_df$Pred_vs_Reg > 0) ## [1] 0.02608495# Difference: National Average - Predicted
diff_nat <- rnatavg - pred_maize
# Difference: Regional Average - Predicted
diff_reg <- regional_raster - pred_maize# Plot National vs Predicted difference
color_palette <- colorRampPalette(c("darkred", "wheat", "darkgreen"))(100)
plot(diff_nat, main = "National Avg - Predicted", col = color_palette)
plot(diff_reg, main = "Regional Avg - Predicted", col = color_palette)