{"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"},{"sourceId":193642924,"sourceType":"kernelVersion"}],"dockerImageVersionId":30618,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"options(repr.plot.width = 22, repr.plot.height = 12)\nlibraries <- c(\"here\", \"tidyverse\", \"data.table\", \"zoo\", \"arrow\",\n               # ML\n               \"tidymodels\", \"ranger\")\n\nlapply(libraries, require, character.only = TRUE)","metadata":{"_uuid":"8c2f0f24-fcd4-4e25-9068-538cddad0028","_cell_guid":"36cc19f6-76c6-443d-8cd0-f5fad1099584","collapsed":false,"_kg_hide-output":true,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load, split and preprocess","metadata":{"_uuid":"5dbabdf8-9ec6-412f-8f20-a682a340f69d","_cell_guid":"d1a328da-6c52-48d4-9f38-85942d2dc76e","trusted":true}},{"cell_type":"code","source":"path_input <- \"../input/r-script01-load-save-features-v2/\"\nfiles_l <- list.files(path_input, pattern=\"^(tr|tst).*(.RData)\", full.names=TRUE)\nfor(file_it in files_l){\n    load(file_it)\n    var_name <- str_extract(file_it, \"tr_maize|tr_wheat|tst_maize|tst_wheat\")\n    assign(var_name, tr_file)\n}\n\nrm(tr_file)","metadata":{"_uuid":"4436aa4b-5d9b-4bc9-b520-7cbbe02a93f4","_cell_guid":"98f987e5-9be3-40c4-a40e-88a45ff7fc3d","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train <- rbind(tr_maize, tr_wheat, fill=TRUE)\nhead(train)","metadata":{"_uuid":"ebb2282c-0d7d-4693-86db-e640bf83d49a","_cell_guid":"e4abac0a-2a13-4f10-9d79-99da6f8d428c","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"years_train <- unique(pull(filter(train, real_year <= 2012), real_year))\nyears_test <- unique(pull(filter(train, real_year > 2012), real_year))\n\ntrain_ <- train\ntrain_train <- train_[real_year %in% years_train,]\n#train_test <- train_[real_year %in% years_test,]","metadata":{"_uuid":"148aede2-8da5-43c7-a547-29ae7ac84160","_cell_guid":"f09b838d-575b-41a7-a51c-176234f4fcc0","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_train_avg <- train_train %>%\n    group_by(lat, lon, crop) %>%\n    summarise(min_ = min(yield), \n              max_ = max(yield),\n              sd_ = replace_na(sd(yield), 0),\n              avg_ = mean(yield)\n             ) %>%\n    ungroup()\n\nsummary(train_train_avg)","metadata":{"_uuid":"e7aca1ce-b48a-4834-8a0b-dd4d194fda19","_cell_guid":"b93ecec1-3c5e-413f-9f4b-8f7c2667edfb","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_maize_train <- tr_maize %>%\n    filter(real_year %in% years_train) %>%\n    left_join(train_train_avg) %>%\n    select(-starts_with(\"tt\")) %>%\n    mutate(yield=yield-min_)\n\n#summary(tr_maize_)","metadata":{"_uuid":"867015d1-bb3b-4e2a-aa97-425e76cf6b03","_cell_guid":"68714b2c-dce3-4b30-b759-5dc04463be80","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preproc_m <- recipe(yield ~., data=tr_maize_train) %>%\n    update_role(ID, new_role = \"ID\") %>%\n    update_role(year, new_role = \"fold_sep\") %>%\n    update_role(real_year, new_role = \"plotting\") %>%\n    update_role(crop, new_role = \"plotting\") %>%\n    # Binary\n    #step_integer(crop) %>%\n    # Dummy variables\n    step_mutate(soil = factor(soil, levels = c(1:13))) %>%\n    step_dummy(soil, one_hot=TRUE, \n               naming=function(var, lvl, ordinal) names0(prefix = \"soil\", num = 13)) %>%\n    # Fill missing values in others with -999\n    step_mutate(across(.cols=starts_with(\"flag_period_\"), .fns=~replace_na(.x, 0)),\n                across(.cols=everything(), .fns=~replace_na(.x, -999)),\n                across(.cols=starts_with(\"soil\"), .fns=~as.integer(.x)),\n                across(.cols=starts_with(\"flag_period_\"), .fns=~as.integer(.x)),\n    # Create interaction between soil and nitrogen\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * nitrogen, .names = \"nit_{.col}\"),\n    # Create interaction between soil and precipitation\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * prec_sum_1_v2, .names = \"prec1_{.col}\"),\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * prec_sum_2_v2, .names = \"prec2_{.col}\"),\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * sum_pre, .names = \"prec0_{.col}\"),\n    # Modify real year for scaling\n                real_year = as.numeric(real_year)) %>%\n    # Normalize\n    step_zv(all_predictors()) %>%\n    step_normalize(all_predictors(), -all_integer_predictors())","metadata":{"_uuid":"b5cffa51-b1df-4c76-8a53-4a5f50d9b77f","_cell_guid":"b0a289c9-b673-49d5-8e50-e174421b3392","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Run fixed parameters","metadata":{"_uuid":"edb6f742-e7ec-47d6-b21d-add8a43d2382","_cell_guid":"7e5a58b8-746f-4c5b-9e75-3d448ef0a0d2","trusted":true}},{"cell_type":"markdown","source":"## Maize","metadata":{"_uuid":"6d6de091-da42-43f6-861a-8312a474d4fd","_cell_guid":"a95fcf0e-475f-4bf1-955e-a4f3a8a03f46","trusted":true}},{"cell_type":"code","source":"set.seed(1234)\nrf_spec_m <- rand_forest(trees = 500) %>%\n  set_engine(\"ranger\", importance = \"impurity\", oob.error=FALSE, seed=1234) %>%\n  set_mode(\"regression\") %>%\n  translate()","metadata":{"_uuid":"2c74eee4-6a0b-41b1-927c-2f8efc0c89e9","_cell_guid":"1acada39-71a5-48ac-8995-f089b3bf2467","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_wf <- workflow() %>%\n    add_recipe(preproc_m) %>%\n    add_model(rf_spec_m)\n\nfinal_wf_m <- training_wf\n\nmodel_m <- fit(final_wf_m, tr_maize_train)","metadata":{"_uuid":"0daf06f7-6255-4417-95e5-b65fb7fe3133","_cell_guid":"22eccc36-f4ca-4172-9c2c-42b688ba81ed","collapsed":false,"_kg_hide-input":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fitted_model <- extract_fit_engine(model_m)\n# importance_matrix <- importance(fitted_model)\n\n# importance_matrix <- as.data.frame(importance_matrix) %>%\n#     rownames_to_column() %>%\n#     arrange(desc(importance_matrix))\n\n# importance_matrix[1:50, ]\n# importance_matrix[51:105, ]","metadata":{"_uuid":"1fa777b3-5e2b-4b2f-926f-dd19f15b8496","_cell_guid":"874b381f-5290-4dbd-b3c3-dbf7856a22bd","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Wheat","metadata":{"_uuid":"aab02d79-4837-4744-af29-51b497fe2338","_cell_guid":"6948e224-e768-4689-b40c-a6a491e8b156","trusted":true}},{"cell_type":"code","source":"tr_wheat_train <- tr_wheat %>%\n    filter(real_year %in% years_train) %>%\n    left_join(train_train_avg) %>%\n    select(-starts_with(\"tt\")) %>%\n    mutate(yield=yield-min_)\n\npreproc_w <- recipe(yield ~., data=tr_wheat_train) %>%\n        update_role(ID, new_role = \"ID\") %>%\n    update_role(year, new_role = \"fold_sep\") %>%\n    update_role(real_year, new_role = \"plotting\") %>%\n    update_role(crop, new_role = \"plotting\") %>%\n    # Binary\n    #step_integer(crop) %>%\n    # Dummy variables\n    step_mutate(soil = factor(soil, levels = c(1:13))) %>%\n    step_dummy(soil, one_hot=TRUE, \n               naming=function(var, lvl, ordinal) names0(prefix = \"soil\", num = 13)) %>%\n    # Fill missing values in others with -999\n    step_mutate(across(.cols=starts_with(\"flag_period_\"), .fns=~replace_na(.x, 0)),\n                across(.cols=everything(), .fns=~replace_na(.x, -999)),\n                across(.cols=starts_with(\"soil\"), .fns=~as.integer(.x)),\n                across(.cols=starts_with(\"flag_period_\"), .fns=~as.integer(.x)),\n    # Create interaction between soil and nitrogen\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * nitrogen, .names = \"nit_{.col}\"),\n    # Create interaction between soil and precipitation\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * prec_sum_1_v2, .names = \"prec1_{.col}\"),\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * prec_sum_2_v2, .names = \"prec2_{.col}\"),\n                across(.cols=starts_with(\"soil\"), .fns= ~.x * sum_pre, .names = \"prec0_{.col}\"),\n    # Modify real year for scaling\n                real_year = as.numeric(real_year)) %>%\n    # Normalize\n    step_zv(all_predictors()) %>%\n    step_normalize(all_predictors(), -all_integer_predictors())","metadata":{"_uuid":"22cac786-80f2-4cf6-ae7c-cd3394844fd3","_cell_guid":"82f1cdea-c0eb-420c-8f79-ae938b9e6a52","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"set.seed(1234)\nrf_spec_w <- rand_forest(trees = 500) %>%\n  set_engine(\"ranger\", importance = \"impurity\", oob.error=FALSE, seed=1234) %>%\n  set_mode(\"regression\") %>%\n  translate()","metadata":{"_uuid":"318949b0-8f76-4d72-989e-80362ad011b9","_cell_guid":"483508ef-9488-4168-afc5-41bc088bcc9d","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_wf <- workflow() %>%\n    add_recipe(preproc_w) %>%\n    add_model(rf_spec_w)\n\nfinal_wf_w <- training_wf\n\nmodel_w <- fit(final_wf_w, tr_wheat_train)","metadata":{"_uuid":"3672c1a1-81df-4b97-af03-7af861385702","_cell_guid":"fc770c6d-1b59-40e7-8cc2-9622e07f9bc8","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fitted_model <- extract_fit_engine(model_w)\n# importance_matrix <- importance(fitted_model)\n\n# importance_matrix <- as.data.frame(importance_matrix) %>%\n#     rownames_to_column() %>%\n#     arrange(desc(importance_matrix))\n\n# importance_matrix[1:50, ]\n# importance_matrix[51:105, ]","metadata":{"_uuid":"abf40150-c70b-44c7-9310-d898d6251389","_cell_guid":"cf0babaf-027f-464d-9db0-8dca5d29ffe2","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc()","metadata":{"_uuid":"c52b2fc7-f418-43c9-82bc-48d28b40ab9d","_cell_guid":"8e7707c0-0042-4481-9f48-b0a236007c78","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Eval predictions","metadata":{"_uuid":"f2e9d1d5-0c0d-4f28-8974-d6114e4ea4f6","_cell_guid":"1c686025-e3be-4946-802c-d3691b49753d","trusted":true}},{"cell_type":"markdown","source":"## Training set","metadata":{"_uuid":"e4627678-6fa2-49b9-9b00-49b39892225c","_cell_guid":"fb0d07b8-a13a-4acb-8503-764b421b3d33","trusted":true}},{"cell_type":"markdown","source":"### Ranges","metadata":{"_uuid":"331b6a58-2b29-4856-b46d-8f1bd8c392ff","_cell_guid":"09c8f0ed-b2f4-4450-9298-2fc72882c05f","trusted":true}},{"cell_type":"code","source":"predict_train_m <- model_m %>%\n    predict(tr_maize_train) %>%\n    bind_cols(select(tr_maize_train, real_year, crop, yield, min_, lat, lon)) %>%\n    mutate(yield = yield + min_, .pred = .pred + min_)\n\nsummary(predict_train_m$yield)\nsummary(predict_train_m$.pred)","metadata":{"_uuid":"1b48ef23-61be-49e6-99bf-0bd3cba2151d","_cell_guid":"dd00896f-0960-4f55-b31c-30877fc30d89","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predict_train_w <- model_w %>%\n    predict(tr_wheat_train) %>%\n    bind_cols(select(tr_wheat_train, real_year, crop, yield, min_, lat, lon)) %>%\n    mutate(yield = yield + min_, .pred = .pred + min_)\n\nsummary(predict_train_w$yield)\nsummary(predict_train_w$.pred)","metadata":{"_uuid":"81bbe467-e19e-45b0-94d4-3d73a3d4c569","_cell_guid":"ae66dda1-b572-4e46-97ee-a513827c5f4d","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predict_train <- predict_train_m %>%\n    rbind(predict_train_w)\n\nsummary(predict_train$yield)\nsummary(predict_train$.pred)","metadata":{"_uuid":"39a6f652-4850-4050-b1a5-04de6a82b13f","_cell_guid":"d29b4c37-c416-4168-be12-cea716c87a6d","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Performance","metadata":{"_uuid":"440746c8-c3dc-4276-a2ef-20de616cf7fc","_cell_guid":"5475ddd0-a286-43e1-a7ed-0fbe81794012","trusted":true}},{"cell_type":"code","source":"results_train <- predict_train %>% group_by(lat, lon, crop)\nrsq_ <- results_train %>% rsq_trad(truth=yield, .pred)\nrmse_ <- results_train %>% rmse(truth=yield, .pred)\nsmape_ <- results_train %>% smape(truth=yield, .pred)\nme_ <- results_train %>% msd(truth=yield, .pred)\n\nprint(\"rsq\")\nsummary(rsq_[\".estimate\"])\nprint(\"rmse\")\nsummary(rmse_[\".estimate\"])\nprint(\"smape\")\nsummary(smape_[\".estimate\"])\nprint(\"me\")\nsummary(me_[\".estimate\"])","metadata":{"_uuid":"d316c07b-a7e7-4b8f-b3d5-69b3c303dee5","_cell_guid":"03621ee2-66c7-4ab6-98a7-f48400ab5ad7","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results_train <- predict_train %>% group_by(real_year, crop)\nrmse_ <- results_train %>% rmse(truth=yield, .pred)\nsmape_ <- results_train %>% smape(truth=yield, .pred)\n\nmetrics_train <- rbind(smape_, rmse_)\n\nmetrics_train_sum <- metrics_train %>%\n    group_by(crop, .metric) %>%\n    summarise(median_ = median(.estimate),\n              mean_ = mean(.estimate))\n\nggplot() +\n    facet_grid(.metric ~ crop, scales = \"free\") +\n    geom_point(data=metrics_train, aes(x=real_year, y=.estimate)) +\n    geom_hline(data=metrics_train_sum, aes(yintercept = median_), color = \"blue\") +\n    geom_hline(data=metrics_train_sum, aes(yintercept = mean_), color = \"red\") +\n    theme_bw()","metadata":{"_uuid":"ad2220d9-2515-4a9b-84bb-ad64ed83b8bd","_cell_guid":"b4693191-2ed1-4d3f-8396-7a50f9f68d5a","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test set","metadata":{"_uuid":"17525e9c-a1a5-4e19-894f-b56577f7fd36","_cell_guid":"6b72d589-870f-497f-952e-9eba086e7717","trusted":true}},{"cell_type":"markdown","source":"### Maize","metadata":{"_uuid":"67d42ea0-4afb-4982-b45b-439684431bc3","_cell_guid":"ab222306-b29d-4773-946b-7ad307ac261f","trusted":true}},{"cell_type":"code","source":"tr_maize_tst <- tr_maize %>%\n    filter(real_year %in% years_test) %>%\n    left_join(train_train_avg) %>%\n    select(-starts_with(\"tt\")) %>%\n    mutate(yield=yield-min_)\n\npredict_m_tst <- model_m %>%\n    predict(tr_maize_tst) %>%\n    bind_cols(select(tr_maize_tst, lat, lon, crop, real_year, yield, min_))","metadata":{"_uuid":"c0f03ef9-19e7-4fa1-8200-1d123edb77fd","_cell_guid":"d72ca1e2-df77-4027-adc7-5f41ae64ba19","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predict_m_tst_ <- predict_m_tst %>%\n    mutate(.pred=.pred+min_, yield = yield + min_,\n           rel_error=abs(.pred-yield)*2/(abs(.pred)+abs(yield)),\n           error = .pred-yield,\n           error = if_else(abs(error) < 1, error, error/abs(error)))\n\nsummary(predict_m_tst_$rel_error)\n\nggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=filter(predict_m_tst_, rel_error > 0.2), aes(x=lon, y=lat, color=error)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()\nggsave(\"./error_latlong_m.png\")","metadata":{"_uuid":"e2314301-0747-4e7a-833c-dfcb92b545c0","_cell_guid":"85843ea7-fead-4e8f-9580-60ebf9b624cf","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results_test <- predict_m_tst %>% group_by(lat, lon, crop)\nrsq_ <- results_test %>% rsq_trad(truth=yield, .pred)\nrmse_ <- results_test %>% rmse(truth=yield, .pred)\nsmape_ <- results_test %>% smape(truth=yield, .pred)\nme_ <- results_train %>% msd(truth=yield, .pred)\n\nprint(\"rsq\")\nsummary(rsq_[\".estimate\"])\nprint(\"rmse\")\nsummary(rmse_[\".estimate\"])\nprint(\"smape\")\nsummary(smape_[\".estimate\"])\nprint(\"me\")\nsummary(me_[\".estimate\"])","metadata":{"_uuid":"6ca1e2ab-b157-4b2a-a32f-f7aba7c29c3b","_cell_guid":"d631738f-7fca-4269-80c6-50e610fbf794","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc()","metadata":{"_uuid":"3deebf6c-0727-406f-a49e-2779aa564828","_cell_guid":"a66f4910-2694-4290-a467-02e69be9353e","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Wheat","metadata":{"_uuid":"5f30b3af-59e2-4b75-ae51-7bf489f17a54","_cell_guid":"a4b2f975-3647-4d31-8ed9-692c6d5348d7","trusted":true}},{"cell_type":"code","source":"tr_wheat_tst <- tr_wheat %>%\n    filter(real_year %in% years_test) %>%\n    left_join(train_train_avg) %>%\n    select(-starts_with(\"tt\")) %>%\n    mutate(yield=yield-min_)\n\npredict_w_tst <- model_w %>%\n    predict(tr_wheat_tst) %>%\n    bind_cols(select(tr_wheat_tst, lat, lon, crop, real_year, yield, min_))","metadata":{"_uuid":"c291a842-c610-4e34-8788-fd101555c37b","_cell_guid":"dc92fbee-5564-47d5-8010-e69c4445506b","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predict_w_tst_ <- predict_w_tst %>%\n    mutate(.pred=.pred+min_, yield = yield + min_,\n           rel_error=abs(.pred-yield)*2/(abs(.pred)+abs(yield)),\n           error = .pred-yield,\n           error = if_else(abs(error) < 1, error, error/abs(error)))\n\nsummary(predict_w_tst_$rel_error)\n\nggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=filter(predict_w_tst_, rel_error > 0.2), aes(x=lon, y=lat, color=error)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()\nggsave(\"./error_latlong_w.png\")","metadata":{"_uuid":"3ad242e8-ac37-4ef2-9357-f93cb589aeca","_cell_guid":"ca690040-462f-4c66-a391-3dc913ef4043","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results_test <- predict_w_tst %>% group_by(real_year, crop)\nrsq_ <- results_test %>% rsq_trad(truth=yield, .pred)\nrmse_ <- results_test %>% rmse(truth=yield, .pred)\nsmape_ <- results_test %>% smape(truth=yield, .pred)\nme_ <- results_train %>% msd(truth=yield, .pred)\n\nprint(\"rsq\")\nsummary(rsq_[\".estimate\"])\nprint(\"rmse\")\nsummary(rmse_[\".estimate\"])\nprint(\"smape\")\nsummary(smape_[\".estimate\"])\nprint(\"me\")\nsummary(me_[\".estimate\"])","metadata":{"_uuid":"2fa23b82-30cf-45ba-8ed2-c9893e12e6bb","_cell_guid":"3c020488-d55a-43ae-a6c9-7f38e55c3ec3","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc()","metadata":{"_uuid":"00b234ec-189b-4c5e-aa3e-d8cbce2757f6","_cell_guid":"21cc85e2-786f-42ea-ab54-610857899315","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metrics","metadata":{"_uuid":"c94ff084-21ee-4b95-8556-f5acd2da51cb","_cell_guid":"b0df76cd-bfdd-4d72-9d72-f02cda44ade9","trusted":true}},{"cell_type":"code","source":"predict_train_tst <- rbind(predict_m_tst, predict_w_tst) %>%\n    mutate(yield = yield + min_, .pred = .pred + min_)","metadata":{"_uuid":"3a0ebabc-b0f3-4495-8665-e60d87c475df","_cell_guid":"a37df183-883a-4ca2-9feb-a0875fc6c095","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results_test <- predict_train_tst %>% group_by(lat, lon, crop)\nrsq_ <- results_test %>% rsq_trad(truth=yield, .pred)\nrmse_ <- results_test %>% rmse(truth=yield, .pred)\nsmape_ <- results_test %>% smape(truth=yield, .pred)\n\nprint(\"rsq\")\nsummary(rsq_[\".estimate\"])\nprint(\"rmse\")\nsummary(rmse_[\".estimate\"])\nprint(\"smape\")\nsummary(smape_[\".estimate\"])","metadata":{"_uuid":"1129eb13-1ec2-4b54-bc26-baf9fd63935c","_cell_guid":"b3a65ee7-823c-4211-8101-1c78b72333a9","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results_test <- predict_train_tst %>% group_by(real_year, crop)\nrmse_ <- results_test %>% rmse(truth=yield, .pred)\nsmape_ <- results_test %>% smape(truth=yield, .pred)\n\nmetrics_test <- rbind(smape_, rmse_)\n\nmetrics_test_sum <- metrics_test %>%\n    group_by(crop, .metric) %>%\n    summarise(median_ = median(.estimate),\n              mean_ = mean(.estimate))\n\nggplot() +\n    facet_grid(.metric ~ crop, scales = \"free\") +\n    geom_point(data=metrics_test, aes(x=real_year, y=.estimate)) +\n    geom_hline(data=metrics_test_sum, aes(yintercept = median_), color = \"blue\") +\n    geom_hline(data=metrics_test_sum, aes(yintercept = mean_), color = \"red\") +\n    theme_bw()","metadata":{"_uuid":"b6e90b70-14a0-4f38-abe7-361e23e2b62e","_cell_guid":"80268d84-8058-4d39-9c87-c9c2cda16fc3","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse_ <- predict_train_tst %>% group_by(real_year, crop, lon) %>% rmse(truth=yield, .pred)\nsmape_ <- predict_train_tst %>% group_by(real_year, crop, lon) %>% smape(truth=yield, .pred)\n\nmetrics_test <- rbind(smape_, rmse_) %>%\n    rbind(smape_)\n\nggplot() +\n    facet_grid(.metric ~ crop, scales = \"free\") +\n    geom_point(data=metrics_test, aes(x=lon, y=.estimate, color = as.factor(real_year))) +\n    theme_bw()","metadata":{"_uuid":"de80fe91-d16a-4913-acac-13cc873f5e46","_cell_guid":"ff95d67e-ba41-4733-890d-7e1eac4e147e","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse_ <- predict_train_tst %>% group_by(real_year, crop, lat) %>% rmse(truth=yield, .pred)\nsmape_ <- predict_train_tst %>% group_by(real_year, crop, lat) %>% smape(truth=yield, .pred)\n\nmetrics_test <- rbind(smape_, rmse_) %>%\n    rbind(smape_)\n\nggplot() +\n    facet_grid(.metric ~ crop, scales = \"free\") +\n    geom_point(data=metrics_test, aes(x=lat, y=.estimate, color = as.factor(real_year))) +\n    theme_bw()","metadata":{"_uuid":"998adaf5-4911-4193-aa65-239bd9dcf8f9","_cell_guid":"6467c35a-90fb-4297-ac7e-c716138b4993","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Predicted yields","metadata":{"_uuid":"e563d1e0-93e4-46ee-aa90-52455ce7b8b4","_cell_guid":"20bab9b5-4594-4e9d-8b79-2e98e3214414","trusted":true}},{"cell_type":"code","source":"ggplot() +\n    geom_point(data=filter(predict_train_tst, crop == \"maize\"), aes(x=lon, y=lat, color=.pred)) +\n    scale_color_continuous(type=\"viridis\", limits = c(0, 23)) +\n    theme_bw()\n\nggplot() +\n    geom_point(data=filter(predict_train_tst, crop == \"wheat\"), aes(x=lon, y=lat, color=.pred)) +\n    scale_color_continuous(type=\"viridis\", limits = c(0, 14)) +\n    theme_bw()","metadata":{"_uuid":"65d09a92-9218-42c0-8fec-b9714d088220","_cell_guid":"5b609cec-4cac-487c-b332-04270586c668","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Low r2","metadata":{"_uuid":"0244d97f-023e-4f4b-876c-1614d0f4a50f","_cell_guid":"ab6120e0-426f-435b-b602-e945aaee019f","trusted":true}},{"cell_type":"code","source":"rsq_ <- predict_train_tst %>% group_by(lat, lon, crop) %>% rsq_trad(truth=yield, .pred) %>%\n    pivot_wider(names_from = .metric, values_from = .estimate)\nrmse_ <- predict_train_tst %>% group_by(lat, lon, crop) %>% rmse(truth=yield, .pred) %>%\n    pivot_wider(names_from = .metric, values_from = .estimate)\nsmape_ <- predict_train_tst %>% group_by(lat, lon, crop) %>% smape(truth=yield, .pred) %>%\n    pivot_wider(names_from = .metric, values_from = .estimate)\n\nerror <- left_join(rsq_, rmse_) %>%\n    left_join(smape_) %>%\n    filter(rsq_trad < 0.5)\n\n#head(error)\n\nggplot() +\n    facet_wrap(\"crop\") +\n    geom_point(data=error, aes(x=lon, y=lat, color = rmse)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()\n\nggplot() +\n    facet_wrap(\"crop\") +\n    geom_point(data=error, aes(x=lon, y=lat, color = smape)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"2c888f24-d555-4830-ae8a-a710b3c810d1","_cell_guid":"afa460ef-c8fe-45f8-a0d4-a280f7ff98f2","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### High r2","metadata":{"_uuid":"46991608-45b2-4d2e-935f-7a15d4f57a10","_cell_guid":"febc9085-e90b-411b-bad5-142990509808","trusted":true}},{"cell_type":"code","source":"rsq_ <- predict_train_tst %>% group_by(lat, lon, crop) %>% rsq_trad(truth=yield, .pred) %>%\n    pivot_wider(names_from = .metric, values_from = .estimate)\nrmse_ <- predict_train_tst %>% group_by(lat, lon, crop) %>% rmse(truth=yield, .pred) %>%\n    pivot_wider(names_from = .metric, values_from = .estimate)\nsmape_ <- predict_train_tst %>% group_by(lat, lon, crop) %>% smape(truth=yield, .pred) %>%\n    pivot_wider(names_from = .metric, values_from = .estimate)\n\nerror <- left_join(rsq_, rmse_) %>%\n    left_join(smape_) %>%\n    filter(rsq_trad > 0.5)\n\n#head(error)\n\nggplot() +\n    facet_wrap(\"crop\") +\n    geom_point(data=error, aes(x=lon, y=lat, color = rmse)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()\n\nggplot() +\n    facet_wrap(\"crop\") +\n    geom_point(data=error, aes(x=lon, y=lat, color = smape)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"6c98cc11-e81e-4bdf-9c1b-e674fecec0d3","_cell_guid":"627c6e6b-23e0-4506-b171-995b81dc79b8","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Re-train and submit","metadata":{"_uuid":"e5c216d2-e7c7-48a9-91ce-fcae02cf4036","_cell_guid":"e1060338-c02f-4bda-a265-d46eb60f33e3","trusted":true}},{"cell_type":"code","source":"train_avg <- train %>%\n    group_by(lat, lon, crop) %>%\n    summarise(min_ = min(yield), \n              max_ = max(yield),\n              sd_ = replace_na(sd(yield), 0), \n              avg_ = mean(yield)) %>%\n    ungroup()","metadata":{"_uuid":"66fc2cc3-387a-42c0-a944-400bd92d90cc","_cell_guid":"6f2cb265-501b-494e-9584-461c117952e3","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Maize","metadata":{"_uuid":"90b642e9-064b-4b76-a83e-6a5c3b9cf38a","_cell_guid":"e2ccd79e-9b46-426a-a86c-894ee0ee5448","trusted":true}},{"cell_type":"code","source":"train_m <- tr_maize %>%\n    left_join(train_avg) %>%\n    select(-starts_with(\"tt\")) %>%\n    mutate(yield=yield-min_)\n\nmodel_m <- fit(final_wf_m, train_m)\nsave(model_m, file=\"model_m.Rdata\")","metadata":{"_uuid":"cde931c8-4706-4dd8-bc76-013af4f029ce","_cell_guid":"79e4d047-13af-4291-a4ce-7485db08dfbf","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load(\"../input/archived-r-script02-process-model-predict-xgb-2m/model_m.Rdata\")","metadata":{"_uuid":"41d6e412-53ee-41de-8fe3-00ebb216dfce","_cell_guid":"b03c5522-0dc6-43b3-8e4e-27148f78cb40","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fitted_model <- extract_fit_engine(model_m)\nimportance_matrix <- importance(fitted_model)\n\nimportance_matrix <- as.data.frame(importance_matrix) %>%\n    rownames_to_column() %>%\n    arrange(desc(importance_matrix))\n\nimportance_matrix[1:50, ]\nimportance_matrix[51:105, ]","metadata":{"_uuid":"9ac6053c-c523-468a-89f3-ee508f8d5bf4","_cell_guid":"7d621d78-2500-4afd-9fe3-7c4015cbaa30","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Wheat","metadata":{"_uuid":"69bedf74-de67-4c74-a66a-5d218dcaf299","_cell_guid":"5d9e1ccd-edff-40dc-95ac-bc0fef59f8f4","trusted":true}},{"cell_type":"code","source":"train_w <- tr_wheat %>%\n    left_join(train_avg) %>%\n    select(-starts_with(\"tt\")) %>%\n    mutate(yield=yield-min_)\n\nmodel_w <- fit(final_wf_w, train_w)\nsave(model_w, file=\"model_w.Rdata\")","metadata":{"_uuid":"fdb804ae-0f4d-463c-88df-4dec28b08ec3","_cell_guid":"3f9380d6-789a-440a-b1e9-54c686d37f88","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load(\"../input/archived-r-script02-process-model-predict-xgb-2m/model_w.Rdata\")","metadata":{"_uuid":"4d34198a-2c75-4763-8242-bce95507b07c","_cell_guid":"d787b33f-8837-47c6-a32f-d758dbd4d4cc","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fitted_model <- extract_fit_engine(model_w)\nimportance_matrix <- importance(fitted_model)\n\nimportance_matrix <- as.data.frame(importance_matrix) %>%\n    rownames_to_column() %>%\n    arrange(desc(importance_matrix))\n\nimportance_matrix[1:50, ]\nimportance_matrix[51:105, ]","metadata":{"_uuid":"35d78f51-a533-4628-90e3-0af1cbf02488","_cell_guid":"8aec8875-b4de-4532-9af7-30c6c84e515c","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predict and submit","metadata":{"_uuid":"b100dbe1-dea3-4b8f-81c2-2af8273ff14c","_cell_guid":"d02ffc81-6b06-4e27-a460-8d54e9265201","trusted":true}},{"cell_type":"code","source":"test_m <- tst_maize %>%\n    left_join(train_avg) %>%\n    select(-starts_with(\"tt\"))\n\nyield_ <- pull(predict(model_m, test_m), .pred)\n\nsolution_m <- test_m[, c(\"ID\", \"min_\")]\nsolution_m[, yield := round(as.numeric(yield_ + min_), 6)]\nsolution_m[, min_ := NULL]","metadata":{"_uuid":"51504aef-394d-432b-866a-ccd38ea1fa0f","_cell_guid":"338cd08c-cf29-43f9-a3cc-2d6808fdfd27","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_w <- tst_wheat %>%\n    left_join(train_avg) %>%\n    select(-starts_with(\"tt\"))\n\nyield_ <- pull(predict(model_w, test_w), .pred)\n\nsolution_w <- test_w[, c(\"ID\", \"min_\")]\nsolution_w[, yield := round(as.numeric(yield_ + min_), 6)]\nsolution_w[, min_ := NULL]","metadata":{"_uuid":"b1ed2ce7-f448-4d47-be66-db2e7e25511d","_cell_guid":"1d8fcc07-e56b-4ae2-938a-79610dcaf022","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"solution <- rbind(solution_m, solution_w)\n\nfilename <- \"submission.csv\"\nwrite.csv(solution, file=here(filename), row.names=FALSE, quote=FALSE)","metadata":{"_uuid":"b0bf976f-9942-4368-a26c-411867702c61","_cell_guid":"70774d22-51ac-479b-b4da-cfb5c1c53f3f","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(train$yield)\nsummary(solution$yield)","metadata":{"_uuid":"00cdd7f2-7e66-4b15-8376-0bef8874aaff","_cell_guid":"8c1220b3-2ac7-484b-a006-29b414e12722","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]}]}