{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81000,"databundleVersionId":8812083,"sourceType":"competition"}],"dockerImageVersionId":30618,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Goal","metadata":{}},{"cell_type":"markdown","source":"We've already seen that a simple model that predicts the average of each gridcell is not that bad. It is even outperforming the other one I submitted! So the goal of this tutorial/baseline is to exploit another simple model, but to discuss a different concept. In this example, we draw attention to one way to split the dataset. We chose to use splits that reflect the temporal relationship of the problem. We use a sliding window in the cross-validation for hyperparameter tuning and we use a separated validation set, comprised by the last years of the full training data, to estimate model error. We perform splits and tuning manually to go through each step of the processes.\n\nWe will use a regression tree model that is based only on the highest level of information, i.e., crop, location and year to comment on splitting this dataset. We do not expect improved performance, since the tree is going to group several locations and apply the same rule to all of them, losing the granularity of the average per gridcell.\n\n(I am aware my code is a little, err convoluted, to say the least. Bear with me. There may be something useful here anyway. :))","metadata":{}},{"cell_type":"code","source":"# Initialization and libraries\nlibraries <- c(\"here\", \"tidyverse\", \"data.table\",\n               # plots\n               \"zoo\", \"arrow\",\n               # ML\n               \"rpart.plot\")\n\nlapply(libraries, require, character.only = TRUE)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-06-19T20:13:08.138783Z","iopub.execute_input":"2024-06-19T20:13:08.140747Z","iopub.status.idle":"2024-06-19T20:13:09.877394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"read_data <- function(filename, files_, max_rows){\n    \n    if(is.null(max_rows)){\n            \n            x <- read_parquet(files_[grepl(filename, files_)]) %>%\n                as.data.table()\n        } else {\n            \n            x <- open_dataset(files_[grepl(filename, files_)]) %>%\n                head(n=max_rows) %>% \n                compute() %>%\n                as.data.table()\n        }\n    \n    return(x)\n}\n\ngenerate_set <- function(crop, set, max_rows, files_){    \n       \n    # Soil and IDs ------------------------------------------------------------\n    filename <- paste(\"soil_co2\", crop, set, sep = \"_\")\n    soil <- read_data(filename, files_, max_rows)\n    \n    soil <- soil[, texture_class := paste0(\"soil_\", letters[texture_class])] \n    soil <- soil[, crop := ifelse(crop == \"maize\", 0, 1)]\n\n    dataset <- soil[, .(ID, crop, lat, lon, real_year)]\n    \n    if(set==\"train\"){\n        \n        filename <- paste(set, \"solutions\", crop, sep = \"_\")\n        x <- read_data(filename, files_, max_rows)\n        \n        dataset <- cbind(dataset, x[, !c(\"ID\")]) \n        \n    }\n\n    return(dataset)\n    \n}","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2024-06-19T18:55:36.772150Z","iopub.execute_input":"2024-06-19T18:55:36.804133Z","iopub.status.idle":"2024-06-19T18:55:36.819622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Retrieving the training set","metadata":{}},{"cell_type":"code","source":"# Constants ---------------------------------------------------------------\nmax_rows_ <- NULL # Allows for testing a subset of the dataset; NULL reads the whole file\npath_inputs <- \"../input/the-future-crop-challenge/\"\nfiles_ <- list.files(path_inputs, full.names = TRUE, pattern = \".parquet\")","metadata":{"execution":{"iopub.status.busy":"2024-06-19T18:55:36.823581Z","iopub.execute_input":"2024-06-19T18:55:36.824995Z","iopub.status.idle":"2024-06-19T18:55:36.844825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train -------------------------------------------------------------------\ntr_maize <- generate_set(crop=\"maize\", set=\"train\", max_rows=max_rows_, files_)\n\ntr_wheat <- generate_set(crop=\"wheat\", set=\"train\", max_rows=max_rows_, files_)\n\ntrain_full <- rbind(tr_maize, tr_wheat, fill=TRUE)\ntrain_full <- train_full[, !c(\"ID\")]\n\nrm(tr_maize, tr_wheat)\ngc()","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-06-19T18:55:36.848810Z","iopub.execute_input":"2024-06-19T18:55:36.850277Z","iopub.status.idle":"2024-06-19T18:55:38.205488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(train_full)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:12:52.499326Z","iopub.execute_input":"2024-06-19T20:12:52.501272Z","iopub.status.idle":"2024-06-19T20:12:52.619967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dim(train_full)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:12:55.245572Z","iopub.execute_input":"2024-06-19T20:12:55.247830Z","iopub.status.idle":"2024-06-19T20:12:55.270680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Splits and tuning","metadata":{}},{"cell_type":"code","source":"total_years <- unique(train_full$real_year)\nyears_train <- trunc(length(unique(total_years))*0.7)\ntrain <- train_full[real_year <= (1982+years_train-1)]\nvalid <- train_full[real_year > (1982+years_train-1)]","metadata":{"execution":{"iopub.status.busy":"2024-06-19T18:55:38.255454Z","iopub.execute_input":"2024-06-19T18:55:38.256729Z","iopub.status.idle":"2024-06-19T18:55:38.290038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split should preserve the relationship between past and future, based on real_year\n# Temporal sliding years: a fixed number of training years; a fraction of data is not used in all folds\n# See https://www.sciencedirect.com/science/article/pii/S0308521X20308775?via%3Dihub for details\ntotal_years <- unique(train$real_year)\nyears_train <- 15\nyears_test <- 4\nyears_init <- min(total_years) + seq(0, 9)\n\nyears_cv <- data.frame(fold = seq(1, 10), years = years_init) %>%\n    mutate(years = map(years, ~ .x + seq(0, (years_train+years_test-1)))) %>%\n    unnest_wider(years, names_sep = \"_\") %>%\n    pivot_longer(-fold, names_to = \"set\", values_to = \"year\") %>%\n    group_by(fold) %>%\n    mutate(set = if_else(year >= (max(year) - years_test + 1), \"test\", \"train\")) %>%\n    mutate(test_remove = any(!year %in% total_years)) %>%\n    filter(!test_remove) %>%\n    ungroup() %>%\n    select(-test_remove)\n\ntail(years_cv)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:13:36.577822Z","iopub.execute_input":"2024-06-19T20:13:36.580467Z","iopub.status.idle":"2024-06-19T20:13:36.705843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Tuning\n# Organize parameters\nnFolds <- max(years_cv$fold)\nseed_ <- 42\nnSamples <- 100\n\n# Prepare parameters\n# complexity parameter. Any split that does not decrease the \n# overall lack of fit by a factor of cp is not attempted. \n# For instance, with anova splitting, this means that the \n# overall R-squared must increase by cp at each step. \nset.seed(seed_)\ncp <- 10^runif(nSamples, -3, -1)\n\n# the minimum number of observations that must \n# exist in a node in order for a split to be attempted.\nset.seed(seed_)\nminsplit <- runif(nSamples,10,300)\n\n# the minimum number of observations in any terminal <leaf> node. \nset.seed(seed_+1)\nminbucket <- runif(nSamples, 10, 300)\n\nparameters <- data.frame(cp = cp, minsplit = minsplit,\n                         minbucket = minbucket,\n                         id = 1:nSamples) %>%\n    merge(data.frame(fold = 1:nFolds), by = NULL)\n\nparameters$MSE <- NA","metadata":{"execution":{"iopub.status.busy":"2024-06-19T18:55:38.432928Z","iopub.execute_input":"2024-06-19T18:55:38.434410Z","iopub.status.idle":"2024-06-19T18:55:38.477551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tail(parameters)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:13:43.388357Z","iopub.execute_input":"2024-06-19T20:13:43.390011Z","iopub.status.idle":"2024-06-19T20:13:43.417181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (tune_it in 1:nrow(parameters)){\n    \n    parameters_it <- parameters[tune_it, ]\n    \n    years_cv_it <- years_cv[years_cv$fold == parameters_it$fold, ]\n    \n    train_folds <- train[(real_year %in% years_cv_it$year[years_cv_it$set == \"train\"]), ]\n    test_fold <- train[(real_year %in% years_cv_it$year[years_cv_it$set == \"test\"]), ]\n\n    md <- rpart(yield ~ ., \n                data = train_folds,\n                cp = parameters$cp[tune_it],\n                minsplit = parameters$minsplit[tune_it],\n                minbucket = parameters$minbucket[tune_it])\n    y_pred <- predict(md, test_fold)\n    y_real <- pull(test_fold, yield)\n\n    parameters$MSE[tune_it] <- mean((y_pred - y_real)**2)\n        \n }","metadata":{"execution":{"iopub.status.busy":"2024-06-19T18:55:38.521888Z","iopub.execute_input":"2024-06-19T18:55:38.523235Z","iopub.status.idle":"2024-06-19T19:18:35.032933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(arrange(parameters, MSE))","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:13:49.597580Z","iopub.execute_input":"2024-06-19T20:13:49.599243Z","iopub.status.idle":"2024-06-19T20:13:49.630862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"parameters_ <- parameters %>%\n    group_by(across(-all_of(c(\"MSE\", \"fold\")))) %>%\n    summarise(MSE = mean(MSE)) %>%\n    arrange(MSE)\n\nerror_cv <- sqrt(pull(parameters_[1, ], MSE))","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-06-19T19:18:35.072147Z","iopub.execute_input":"2024-06-19T19:18:35.073577Z","iopub.status.idle":"2024-06-19T19:18:35.141757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(parameters_)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:13:56.465450Z","iopub.execute_input":"2024-06-19T20:13:56.467242Z","iopub.status.idle":"2024-06-19T20:13:56.501896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It's noticeable that the lower bound for the cp parameter likely can be reduced, but we didn't so that the tree could still be readable. We are not pursuing performance in this example.","metadata":{}},{"cell_type":"markdown","source":"# Estimate model error","metadata":{}},{"cell_type":"code","source":"model_str <- \"rpart\"\nsuffix <- \"001_tree\"\n\nif(model_str == \"avg\"){\nmodel_ <- round(mean(train$yield), 4)\n} else if (model_str == \"rpart\"){\nmodel_ <- rpart(yield~., data=train, cp = parameters_$cp[1], minsplit = parameters_$minsplit[1], minbucket = parameters_$minbucket[1])\n} else if (model_str == \"rf\"){\nmodel_ <- randomForest(yield~., data=train, ntree=50)\n} else if (model_str == \"xgb\"){\nmodel_ <- xgboost(data=as.matrix(copy(train)[, yield:=NULL]),\n                  label=copy(train)[[\"yield\"]],\n                  nrounds = 20)\n}","metadata":{"execution":{"iopub.status.busy":"2024-06-19T19:18:35.182572Z","iopub.execute_input":"2024-06-19T19:18:35.184040Z","iopub.status.idle":"2024-06-19T19:18:39.681238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Errors\nif(model_str != \"avg\"){\n  if(model_str == \"xgb\"){\n    yield_ <- predict(model_, as.matrix(copy(valid)[, yield:=NULL])) \n  } else {\n    yield_ <- predict(model_, valid) \n  }\n} else {\n  yield_ <- model_\n}\n\n# rm(train)\n\nsolution <- valid[, .(yield)]\nsolution[, pred := round(as.numeric(yield_), 6)]\nsolution[, err := (pred - yield)]\nerror_est <- sqrt(mean(solution$err^2))\n\nfilename <- paste0(\"errorest_\", model_str, suffix, \".csv\")\nwrite.csv(error_est, file=here(filename), row.names=FALSE, quote=FALSE)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T19:18:39.685443Z","iopub.execute_input":"2024-06-19T19:18:39.687030Z","iopub.status.idle":"2024-06-19T19:18:39.849454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(paste0(\"RMSE_CV: \", round(error_cv, 4), \"; RMSE_valid: \", round(error_est, 4)))","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:14:06.497040Z","iopub.execute_input":"2024-06-19T20:14:06.498872Z","iopub.status.idle":"2024-06-19T20:14:06.514579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train model\nObtaining the model, however, will use all data available for training.","metadata":{}},{"cell_type":"code","source":"model_str <- \"rpart\"\nsuffix <- \"001_tree\"\n\nif(model_str == \"avg\"){\nmodel_ <- round(mean(train_full$yield), 4)\n} else if (model_str == \"rpart\"){\nmodel_ <- rpart(yield~., data=train_full, cp = parameters_$cp[1], minsplit = parameters_$minsplit[1], minbucket = parameters_$minbucket[1])\n} else if (model_str == \"rf\"){\nmodel_ <- randomForest(yield~., data=train_full, ntree=50)\n} else if (model_str == \"xgb\"){\nmodel_ <- xgboost(data=as.matrix(copy(train_full)[, yield:=NULL]),\n                  label=copy(train_full)[[\"yield\"]],\n                  nrounds = 20)\n}","metadata":{"execution":{"iopub.status.busy":"2024-06-19T19:18:39.874869Z","iopub.execute_input":"2024-06-19T19:18:39.876327Z","iopub.status.idle":"2024-06-19T19:18:46.949905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rpart.plot(model_)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:14:12.022726Z","iopub.execute_input":"2024-06-19T20:14:12.024496Z","iopub.status.idle":"2024-06-19T20:14:16.399842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Retrieve test set","metadata":{}},{"cell_type":"code","source":"tst_maize <- generate_set(crop=\"maize\", set=\"test\", max_rows=max_rows_, files_)\n\ntst_wheat <- generate_set(crop=\"wheat\", set=\"test\", max_rows=max_rows_, files_)\n\ntest <- rbind(tst_maize, tst_wheat, fill=TRUE) \nids <- test[, c(\"ID\")]\ntest <- test[, !c(\"ID\")]\norder_cols <- colnames(train)[colnames(train) != \"yield\"]\ntest <- test[, ..order_cols]\n\nrm(tst_maize, tst_wheat)\ngc()","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-06-19T19:18:50.114581Z","iopub.execute_input":"2024-06-19T19:18:50.116931Z","iopub.status.idle":"2024-06-19T19:18:51.611937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction","metadata":{}},{"cell_type":"code","source":"if(model_str != \"avg\"){\n  if(model_str == \"xgb\"){\n    yield_ <- predict(model_, as.matrix(copy(test)))\n  } else {\n    yield_ <- predict(model_, test) \n  }\n} else {\n  yield_ <- model_\n}\n\nsolution2 <- ids\nsolution2[, yield := round(as.numeric(yield_), 6)]\nfilename <- paste0(\"submission_\", model_str, suffix, \".csv\")\nwrite.csv(solution2, file=here(filename), row.names=FALSE, quote=FALSE)\n\nfilename <- \"submission.csv\"\nwrite.csv(solution2, file=here(filename), row.names=FALSE, quote=FALSE)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T19:18:51.615616Z","iopub.execute_input":"2024-06-19T19:18:51.617120Z","iopub.status.idle":"2024-06-19T19:18:57.799197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(solution2)","metadata":{"execution":{"iopub.status.busy":"2024-06-19T20:14:17.507349Z","iopub.execute_input":"2024-06-19T20:14:17.509007Z","iopub.status.idle":"2024-06-19T20:14:17.571345Z"},"trusted":true},"execution_count":null,"outputs":[]}]}