{"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":"code","source":"# This R environment comes with many helpful analytics packages installed\n# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats\n# For example, here's a helpful package to load\n\nlibrary(tidyverse) # metapackage of all tidyverse packages\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nlist.files(path = \"../input\")\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2024-06-24T12:48:29.896308Z","iopub.execute_input":"2024-06-24T12:48:29.899241Z","iopub.status.idle":"2024-06-24T12:48:30.967202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(arrow)\nlibrary(matrixStats)\nlibrary(rpart)\nlibrary(rpart.plot)\nlibrary(rattle)\nlibrary(RColorBrewer)\nlibrary(tidymodels)\nlibrary(ranger)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:30.971422Z","iopub.execute_input":"2024-06-24T12:48:31.010644Z","iopub.status.idle":"2024-06-24T12:48:33.942676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Read data and get simple features","metadata":{}},{"cell_type":"markdown","source":"First, we'll open the data. The soil_co2 file contains info on the soil texture, ambient CO2 and fertilization rate, along with the year, latitude and longitude of each datapoint.","metadata":{}},{"cell_type":"code","source":"maize_data <- read_parquet(\"../input/the-future-crop-challenge/soil_co2_maize_train.parquet\") ","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:33.945593Z","iopub.execute_input":"2024-06-24T12:48:33.947151Z","iopub.status.idle":"2024-06-24T12:48:34.088958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(maize_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:34.091978Z","iopub.execute_input":"2024-06-24T12:48:34.093582Z","iopub.status.idle":"2024-06-24T12:48:34.236835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will now calculate some simple predictive features for our model - the total precipitation during six months from sowing, and the mean, minimum and maximum temperature during six months from sowing.","metadata":{}},{"cell_type":"code","source":"# using open_dataset lets us work with the data lazily, using less memory\npr_maize_data <- open_dataset(\"../input/the-future-crop-challenge/pr_maize_train.parquet\")","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:34.254291Z","iopub.execute_input":"2024-06-24T12:48:34.255979Z","iopub.status.idle":"2024-06-24T12:48:34.325058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we have to collect() to read the right columns into memory before aggregating\n# here we are selecting the days from sowing until 6 months after (36:215)\nmaize_data <- pr_maize_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(pr_sum = rowSums(.[, 2:181])) %>% select(ID,pr_sum) %>% merge(maize_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:34.328030Z","iopub.execute_input":"2024-06-24T12:48:34.329632Z","iopub.status.idle":"2024-06-24T12:48:40.174535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(maize_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:40.178955Z","iopub.execute_input":"2024-06-24T12:48:40.180687Z","iopub.status.idle":"2024-06-24T12:48:40.291460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tas_maize_data <- open_dataset(\"../input/the-future-crop-challenge/tas_maize_train.parquet\")","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:40.295584Z","iopub.execute_input":"2024-06-24T12:48:40.297337Z","iopub.status.idle":"2024-06-24T12:48:40.337299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"maize_data <- tas_maize_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(\n        tas_mean = rowMeans(.[, 2:181]),\n        tas_min = rowMins(as.matrix(.[, 2:181])),\n        tas_max = rowMaxs(as.matrix(.[, 2:181])),\n    ) %>% select(ID,tas_mean,tas_min,tas_max) %>% merge(maize_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:40.341670Z","iopub.execute_input":"2024-06-24T12:48:40.343415Z","iopub.status.idle":"2024-06-24T12:48:52.109498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(maize_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:52.113794Z","iopub.execute_input":"2024-06-24T12:48:52.115463Z","iopub.status.idle":"2024-06-24T12:48:52.244968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we'll add the target variable - yield.","metadata":{}},{"cell_type":"code","source":"target_maize <- read_parquet(\"../input/the-future-crop-challenge/train_solutions_maize.parquet\")","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:52.249126Z","iopub.execute_input":"2024-06-24T12:48:52.250672Z","iopub.status.idle":"2024-06-24T12:48:52.316391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"maize_data <- merge(maize_data, target_maize)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:52.320571Z","iopub.execute_input":"2024-06-24T12:48:52.322215Z","iopub.status.idle":"2024-06-24T12:48:52.760573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(maize_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:52.765099Z","iopub.execute_input":"2024-06-24T12:48:52.766956Z","iopub.status.idle":"2024-06-24T12:48:52.903796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Split into train and test\n\nWe have ~40 years of data. We will take the last 10 years as a test set (2010-2020), and use the rest to train our models","metadata":{}},{"cell_type":"code","source":"train <- maize_data %>% filter(real_year <= 2010)\ntest <- maize_data %>% filter(real_year > 2010)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:51:52.186856Z","iopub.execute_input":"2024-06-24T12:51:52.188698Z","iopub.status.idle":"2024-06-24T12:51:52.229883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"head(train)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:52.980974Z","iopub.execute_input":"2024-06-24T12:48:52.982717Z","iopub.status.idle":"2024-06-24T12:48:53.018941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Building a regression tree","metadata":{}},{"cell_type":"code","source":"tree <- rpart(yield ~ tas_mean + tas_max + tas_min + pr_sum + co2 + texture_class + nitrogen + lat + lon,\n              data = train)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:48:53.023254Z","iopub.execute_input":"2024-06-24T12:48:53.024853Z","iopub.status.idle":"2024-06-24T12:49:00.848365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width = 15, repr.plot.height = 10)\nrpart.plot(tree, box.palette =\"GnBu\", type=4, cex=1.2,\n          clip.right.labs = TRUE, branch.lty = 3)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:00.852932Z","iopub.execute_input":"2024-06-24T12:49:00.854618Z","iopub.status.idle":"2024-06-24T12:49:03.351593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"When rpart grows a tree, it uses cross-validation; we can print the results:","metadata":{}},{"cell_type":"code","source":"printcp(tree)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:03.354531Z","iopub.execute_input":"2024-06-24T12:49:03.356068Z","iopub.status.idle":"2024-06-24T12:49:03.394167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plotcp(tree)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:03.396983Z","iopub.execute_input":"2024-06-24T12:49:03.398443Z","iopub.status.idle":"2024-06-24T12:49:03.574301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Generate predictions using our tree on the test set:","metadata":{}},{"cell_type":"code","source":"pred <- predict(tree, newdata = test)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.290087Z","iopub.execute_input":"2024-06-24T12:49:26.291625Z","iopub.status.idle":"2024-06-24T12:49:26.364456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Exercises:","metadata":{}},{"cell_type":"markdown","source":"Calculate how well the tree performed on the test data in terms of MSE.","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What happens if we use less pruning, by setting a very low cp value of 0.001?","metadata":{}},{"cell_type":"code","source":"#tree <- rpart(yield ~ tas_mean + tas_max + tas_min + pr_sum + co2 + texture_class + nitrogen + lat + lon,\n#              data = train, cp=CHOOSE_CP_VALUE_HERE)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.367425Z","iopub.execute_input":"2024-06-24T12:49:26.368998Z","iopub.status.idle":"2024-06-24T12:49:26.380809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot the tree (see code above)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.383707Z","iopub.execute_input":"2024-06-24T12:49:26.385241Z","iopub.status.idle":"2024-06-24T12:49:26.396648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Print the CV results (see code above)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.399515Z","iopub.execute_input":"2024-06-24T12:49:26.401053Z","iopub.status.idle":"2024-06-24T12:49:26.412483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot the curve of X-val error against cp (see code above)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.415400Z","iopub.execute_input":"2024-06-24T12:49:26.416954Z","iopub.status.idle":"2024-06-24T12:49:26.428340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"However, model performance on the training set does not always correspond to skill on the test set. Vary the cp value (amount of pruning) and measure the performance of the resulting trees on the **test** set (in terms of MSE).\n\nWhat is the optimal cp value?","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Optional:** Try altering other hyperparameters when training the tree (e.g., minsplit, minbucket) and try to achieve an even better model performance on the test set.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Random forests\n\nWe'll now train ensembles of trees instead, and see how much better our performance can get.\n\nWe will use 5-fold cross-validation to evaluate the model performance on the training set...","metadata":{}},{"cell_type":"code","source":"# Define a basic random forest with 10 trees\nrf <-   \n  rand_forest(mode = \"regression\", trees=10) %>%\n  set_engine(\"ranger\", seed = 1, importance=\"impurity\")","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.431286Z","iopub.execute_input":"2024-06-24T12:49:26.432802Z","iopub.status.idle":"2024-06-24T12:49:26.453312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rf","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.456548Z","iopub.execute_input":"2024-06-24T12:49:26.458122Z","iopub.status.idle":"2024-06-24T12:49:26.477767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Then create our CV folds:\nfolds <- vfold_cv(train, v = 5)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.481502Z","iopub.execute_input":"2024-06-24T12:49:26.483079Z","iopub.status.idle":"2024-06-24T12:49:26.806792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now we create a workflow, which puts the model and formula together.\n# It could also contain the recipe we made earlier for preprocessing.\nrf_workflow <- workflow() %>% \n  add_model(rf) %>% \n  add_formula(yield ~ tas_mean + tas_max + tas_min + pr_sum + co2 + texture_class + nitrogen + lat + lon)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.810952Z","iopub.execute_input":"2024-06-24T12:49:26.812671Z","iopub.status.idle":"2024-06-24T12:49:26.830513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Then do the fitting and testing for each CV split, and get metrics:\n# Running this will take a couple of minutes!\nrf_cv <- rf_workflow %>%\n  fit_resamples(folds, control = control_resamples(save_pred = TRUE))\ncollect_metrics(rf_cv)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:49:26.834640Z","iopub.execute_input":"2024-06-24T12:49:26.836362Z","iopub.status.idle":"2024-06-24T12:50:42.709671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"collect_predictions(rf_cv) %>% \n      ggplot(aes(yield, .pred)) +\n      geom_point(alpha = 0.3) +\n      labs(\n      x = \"Observed\",\n      y = \"Predicted\",\n      title = 'Our RF predictions on the training set'\n      )\n","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:50:42.714132Z","iopub.execute_input":"2024-06-24T12:50:42.716627Z","iopub.status.idle":"2024-06-24T12:50:59.128918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Adding a categorical variable\n\nAnother way trees are nice is that you can also use categorical variables very easily. We will now include data from wheat, as well as maize, and train our model on data from both crops, using 'crop' as a categorical variable.","metadata":{}},{"cell_type":"code","source":"# First, read and process the wheat data, exactly as we did for maize:\nwheat_data <- read_parquet(\"../input/the-future-crop-challenge/soil_co2_wheat_train.parquet\")\npr_wheat_data <- open_dataset(\"../input/the-future-crop-challenge/pr_wheat_train.parquet\")\nwheat_data <- pr_wheat_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(pr_sum = rowSums(.[, 2:181])) %>% select(ID,pr_sum) %>% merge(wheat_data)\ntas_wheat_data <- open_dataset(\"../input/the-future-crop-challenge/tas_wheat_train.parquet\")\nwheat_data <- tas_wheat_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(\n        tas_mean = rowMeans(.[, 2:181]),\n        tas_min = rowMins(as.matrix(.[, 2:181])),\n        tas_max = rowMaxs(as.matrix(.[, 2:181])),\n    ) %>% select(ID,tas_mean,tas_min,tas_max) %>% merge(wheat_data)\ntarget_wheat <- read_parquet(\"../input/the-future-crop-challenge/train_solutions_wheat.parquet\")\nwheat_data <- merge(wheat_data, target_wheat)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:50:59.131961Z","iopub.execute_input":"2024-06-24T12:50:59.133559Z","iopub.status.idle":"2024-06-24T12:51:12.657310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# add the wheat data to our train and test data\ntrain <- train %>% rbind(wheat_data) %>% filter(real_year <= 2010)\ntest <- test %>% rbind(wheat_data) %>% filter(real_year > 2010)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:52:08.530594Z","iopub.execute_input":"2024-06-24T12:52:08.532236Z","iopub.status.idle":"2024-06-24T12:52:08.714188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define a basic random forest with 10 trees, again:\nrf <-   \n  rand_forest(mode = \"regression\", trees=10) %>%\n  set_engine(\"ranger\", seed = 1, importance=\"impurity\")","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:52:10.623478Z","iopub.execute_input":"2024-06-24T12:52:10.625190Z","iopub.status.idle":"2024-06-24T12:52:10.640198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Then create our CV folds:\nfolds <- vfold_cv(train, v = 5)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:52:11.212387Z","iopub.execute_input":"2024-06-24T12:52:11.214063Z","iopub.status.idle":"2024-06-24T12:52:11.486365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now we will create a recipe for our pre-processing step, where we encode our\n# categorical variables as a dummy variable.\nrf_rec <- \n    recipe(yield ~ tas_mean + tas_max + tas_min + pr_sum + co2 + texture_class + nitrogen + lat + lon,\n         data=train) %>% \n    step_dummy(all_nominal(), -all_outcomes())","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:52:11.761895Z","iopub.execute_input":"2024-06-24T12:52:11.763545Z","iopub.status.idle":"2024-06-24T12:52:11.788826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Our new workflow includes the preprocessing step now,\n# and the formula is already in the recipe, so we don't need to specify that.\nrf_workflow <- workflow() %>% \n  add_recipe(rf_rec) %>%\n  add_model(rf)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:52:12.340250Z","iopub.execute_input":"2024-06-24T12:52:12.342084Z","iopub.status.idle":"2024-06-24T12:52:12.357444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Then do the fitting and testing for each CV split, like before:\n# This takes a couple minutes!\nrf_cv <- rf_workflow %>%\n  fit_resamples(folds, control = control_resamples(save_pred = TRUE))\ncollect_metrics(rf_cv)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:52:12.891283Z","iopub.execute_input":"2024-06-24T12:52:12.893015Z","iopub.status.idle":"2024-06-24T12:54:42.101603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"collect_predictions(rf_cv) %>% \n      ggplot(aes(yield, .pred)) +\n      geom_point(alpha = 0.3) +\n      labs(\n      x = \"Observed\",\n      y = \"Predicted\",\n      title = 'Our RF predictions on the training set, for maize and wheat'\n      )\n","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:54:42.104558Z","iopub.execute_input":"2024-06-24T12:54:42.106069Z","iopub.status.idle":"2024-06-24T12:55:09.949652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not bad!","metadata":{}},{"cell_type":"markdown","source":"### Evaluating on the test set","metadata":{}},{"cell_type":"markdown","source":"Now let's fit the model to the entire training set and look at how it performs on the test set.","metadata":{}},{"cell_type":"code","source":"final_fit <- \n  rf_workflow %>%\n  fit(data = train)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:57:49.601800Z","iopub.execute_input":"2024-06-24T12:57:49.603512Z","iopub.status.idle":"2024-06-24T12:58:27.200806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make predictions:\npred <- predict(final_fit, new_data = test)\npred$truth <- test$yield\nmetrics(pred, truth = truth, estimate = .pred)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:58:27.203825Z","iopub.execute_input":"2024-06-24T12:58:27.205393Z","iopub.status.idle":"2024-06-24T12:58:28.075221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred %>% \n      ggplot(aes(truth, .pred)) +\n      geom_point(alpha = 0.3) +\n      labs(\n      x = \"Observed\",\n      y = \"Predicted\",\n      title = 'Our RF predictions on the test set'\n      )\n","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:58:28.078088Z","iopub.execute_input":"2024-06-24T12:58:28.079679Z","iopub.status.idle":"2024-06-24T12:58:37.641483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Why is our model performance on the test set worse than when we calculated it on the training set using 5-fold cross-validation?**","metadata":{}},{"cell_type":"markdown","source":"### Exercise","metadata":{}},{"cell_type":"markdown","source":"Try to improve the model as much as possible, and then predict yields on the real test set (2020-2100) and submit through Kaggle. There's code cells below to get the data, generate predictions and submit - all you have to do is refine the model and then run the cells.\n\nSome ideas:\n* Tune the hyperparameters (mtry, trees, min_n)\n* Change the features used by the model\n* Use a different cross-validation strategy while tuning and selecting features\n* Generate more sophisticated features (don't forget you will need to do that also for the real test too)","metadata":{}},{"cell_type":"code","source":"## To help you, here's some code that tunes hyperparameters\n## we'll tune mtry and the number of trees - define them as tunable:\n#tuned_rf <- rf %>% \n#  update(mtry = tune(), trees = tune())\n\n## make our new workflow:    \n#rf_workflow <- workflow() %>% \n#  add_model(tuned_rf)\n\n## Now tune:\n#grid_tune <-\n#  rf_workflow %>%\n#  tune_grid(\n#    resamples = folds, \n#    grid = expand.grid( # here we define a grid of potential\n#      mtry = seq(1, 10), # hyperparameter values, with mtry from 1 to 10 and\n#      trees = c(10, 20, 30) # number of trees either 100, 200 or 300.\n#    ),                      # that's 30 different combinations to try.\n#  )\n\n#show_best(grid_tune, n = 3) # print the 3 best combinations.","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:58:37.644412Z","iopub.execute_input":"2024-06-24T12:58:37.645955Z","iopub.status.idle":"2024-06-24T12:58:37.658870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#autoplot(grid_tune, metric = \"rmse\") +\n#    labs(title = \"Results of tuning\")","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:58:37.661745Z","iopub.execute_input":"2024-06-24T12:58:37.663265Z","iopub.status.idle":"2024-06-24T12:58:37.675259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Process data of the real test set","metadata":{}},{"cell_type":"code","source":"maize_test_data <- read_parquet(\"../input/the-future-crop-challenge/soil_co2_maize_test.parquet\") \n\npr_maize_test_data <- open_dataset(\"../input/the-future-crop-challenge/pr_maize_test.parquet\")\n\nmaize_test_data <- pr_maize_test_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(pr_sum = rowSums(.[, 2:181])) %>% select(ID,pr_sum) %>% merge(maize_test_data)\n\ntas_maize_test_data <- open_dataset(\"../input/the-future-crop-challenge/tas_maize_test.parquet\")\n\nmaize_test_data <- tas_maize_test_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(\n        tas_mean = rowMeans(.[, 2:181]),\n        tas_min = rowMins(as.matrix(.[, 2:181])),\n        tas_max = rowMaxs(as.matrix(.[, 2:181])),\n    ) %>% select(ID,tas_mean,tas_min,tas_max) %>% merge(maize_test_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:58:37.678208Z","iopub.execute_input":"2024-06-24T12:58:37.679755Z","iopub.status.idle":"2024-06-24T12:59:09.724146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"wheat_test_data <- read_parquet(\"../input/the-future-crop-challenge/soil_co2_wheat_test.parquet\")\n\npr_wheat_test_data <- open_dataset(\"../input/the-future-crop-challenge/pr_wheat_test.parquet\")\n\nwheat_test_data <- pr_wheat_test_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(pr_sum = rowSums(.[, 2:181])) %>% select(ID,pr_sum) %>% merge(wheat_test_data)\n\ntas_wheat_test_data <- open_dataset(\"../input/the-future-crop-challenge/tas_wheat_test.parquet\")\n\nwheat_test_data <- tas_wheat_test_data %>% select(ID, 36:215) %>% collect() %>%\n    mutate(\n        tas_mean = rowMeans(.[, 2:181]),\n        tas_min = rowMins(as.matrix(.[, 2:181])),\n        tas_max = rowMaxs(as.matrix(.[, 2:181])),\n    ) %>% select(ID,tas_mean,tas_min,tas_max) %>% merge(wheat_test_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T12:59:09.729230Z","iopub.execute_input":"2024-06-24T12:59:09.731439Z","iopub.status.idle":"2024-06-24T12:59:50.984963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"final_test = rbind(maize_test_data, wheat_test_data)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:42:16.738229Z","iopub.execute_input":"2024-06-24T13:42:16.740222Z","iopub.status.idle":"2024-06-24T13:42:16.982290Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"full_train = rbind(train, test)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:42:17.224002Z","iopub.execute_input":"2024-06-24T13:42:17.225878Z","iopub.status.idle":"2024-06-24T13:42:17.383265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# retrain on the entire dataset:\nfinal_model <- \n  rf_workflow %>%\n  fit(data = full_train)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:42:18.648169Z","iopub.execute_input":"2024-06-24T13:42:18.649944Z","iopub.status.idle":"2024-06-24T13:43:15.947385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make predictions:\npred <- predict(final_model, new_data = final_test)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:43:15.952503Z","iopub.execute_input":"2024-06-24T13:43:15.954264Z","iopub.status.idle":"2024-06-24T13:43:21.697196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"my_submission <- data_frame('ID' = final_test$ID, 'yield' = pred$.pred)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:43:21.702414Z","iopub.execute_input":"2024-06-24T13:43:21.704296Z","iopub.status.idle":"2024-06-24T13:43:21.724118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"my_submission","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:43:21.727476Z","iopub.execute_input":"2024-06-24T13:43:21.729176Z","iopub.status.idle":"2024-06-24T13:43:21.792296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# write out the csv file of final predictions\nwrite.csv(my_submission, file=\"submission.csv\", row.names=FALSE, quote=FALSE)","metadata":{"execution":{"iopub.status.busy":"2024-06-24T13:44:09.682528Z","iopub.execute_input":"2024-06-24T13:44:09.684609Z","iopub.status.idle":"2024-06-24T13:44:12.361108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To submit, just press the button on the tab on the right under 'Submit to competition'. Let's see!","metadata":{}}]}