## ----------------------------------------------------------------------------- # install packages install.packages('randomForest', repos="http://cran.us.r-project.org", dependencies=T, quiet=T) install.packages('caret', repos="http://cran.us.r-project.org", dependencies=T, quiet=T) install.packages('caretEnsemble', repos="http://cran.us.r-project.org", dependencies=T, quiet=T) install.packages('iml', repos="http://cran.us.r-project.org", dependencies=T, quiet=T) install.packages('tidyverse', repos="http://cran.us.r-project.org", dependencies=T, quiet=T) install.packages('dplyr', repos="http://cran.us.r-project.org", dependencies=T, quiet=T) # # install packages packages <- c("elasticnet", "partykit", "rpart", "rpart.plot") installed_packages <- packages %in% rownames(installed.packages()) if(any(installed_packages == FALSE)){ install.packages(packages[!installed_packages], repos="http://cran.us.r-project.org", quiet=T)} # # load packages packages <- c("elasticnet", "partykit", "rpart", "rpart.plot", "randomForest", "caret", "caretEnsemble", "iml", "tidyverse", "dplyr") invisible(lapply(packages, function(x) suppressMessages(require(x, character.only=T, quietly=T, warn.conflicts=F)))) ## ----------------------------------------------------------------------------- # read .csv file with data file <- 'https://raw.githubusercontent.com/jvasco323/eia-yg-training-ppt/master/pooled_data_wheat_phd_fcr.csv' pooled_data <- read.csv(url(file), header=TRUE) # # list variables of interest for visualization str(pooled_data) ## ----------------------------------------------------------------------------- set.seed(1234) pooled_data1 <- pooled_data %>% sample_frac(size = 0.3) ## ----------------------------------------------------------------------------- # leveling categorical variables pooled_data1$lodging_cat <- as.factor(pooled_data1$lodging_cat) pooled_data1$lodging_cat <- relevel(pooled_data1$lodging_cat, ref="No") pooled_data1$irri_num <- as.factor(pooled_data1$irri_num) pooled_data1$irri_num <- relevel(pooled_data1$irri_num, ref="<=2") pooled_data1$rice_till <- as.factor(pooled_data1$rice_till) pooled_data1$rice_till <- relevel(pooled_data1$rice_till, ref="<=4") pooled_data1$residue <- as.factor(pooled_data1$residue) pooled_data1$residue <- relevel(pooled_data1$residue, ref="No") pooled_data1$till_int <- as.factor(pooled_data1$till_int) pooled_data1$till_int <- relevel(pooled_data1$till_int, ref="Intensive") pooled_data1$Texture <- as.factor(pooled_data1$Texture) pooled_data1$Texture <- relevel(pooled_data1$Texture, ref="clay") pooled_data1$insect_severity <- as.factor(pooled_data1$insect_severity) pooled_data1$insect_severity <- relevel(pooled_data1$insect_severity, ref="None") pooled_data1$disease_severity <- as.factor(pooled_data1$disease_severity) pooled_data1$disease_severity <- relevel(pooled_data1$disease_severity, ref="None") pooled_data1$weed_severity <- as.factor(pooled_data1$weed_severity) pooled_data1$weed_severity <- relevel(pooled_data1$weed_severity, ref="No") pooled_data1$X = NULL ## ----------------------------------------------------------------------------- # partitioning data into train and test set set.seed(458) train_seq <- createDataPartition(pooled_data1$cc_yield, p=0.7, list=F) train_data <- pooled_data1[train_seq,] test_data <- pooled_data1[-train_seq,] # # check near zero variance predictor (nzv <- nearZeroVar(train_data, saveMetrics=T)) ## ----------------------------------------------------------------------------- # five-fold cross-validation repeated once with random sampling control <- trainControl(method="repeatedcv", number=5, repeats=1, search="random") algorithmList <- c("ranger", "rpart", "lasso", "lm") # # fit the models models <- caretList(cc_yield ~ ., data=train_data, trControl=control, methodList=algorithmList) # # check model performance and select the final model model_selection <- summary(resamples(models)) model_selection ## ---- fig.height=4.5, fig.width=4.5------------------------------------------- # define five-fold cross-validation scheme control_feature <- rfeControl(functions=rfFuncs, method="cv", number=5) # # recursive feature elimination; 'sizes' stands for the number of variables feature_selection <- rfe(train_data[,2:23], train_data[,1], sizes=c(1:23), rfeControl=control_feature) # # rmse vs. number of variables (ordered from very to not important) plot(feature_selection) # # r2 vs. number of variables (ordered from very to not important) # plot(feature_selection$results$Rsquared) ## ----------------------------------------------------------------------------- # list all variables list_of_var <- feature_selection$optVariables list_of_var <- append("cc_yield", list_of_var) # # select 10 most important variables only train_data <- train_data[,list_of_var[1:11]] test_data <- test_data[,list_of_var[1:11]] ## ----------------------------------------------------------------------------- # parameter combinations for grid search approach ntrees <- c(100, 500) nodesize <- seq(10,130, 15) tuneGrid <- expand.grid(.mtry = c(4,6,8,10,12)) params <- expand.grid(ntrees = ntrees, nodesize = nodesize) # # cross-validation scheme for the hyper-parameter tuning control <- trainControl(method="repeatedcv", number=5, repeats=1, search='grid', verboseIter=F) # # create vector to store the 18 combinations of nodesize x ntree, at the best mtry value all_model <- vector("list", nrow(params)) # # fit the model for different combinations of hyper-parameters for(i in 1:nrow(params)){ nodesize <- params[i,2] ntree <- params[i,1] set.seed(65) rf_model <- train(cc_yield ~ ., # y = f(x) data=train_data, # select training data method="rf", # choose random forest tuneGrid=tuneGrid, # evaluate mtry for each combination trControl=control, # set cross-validation scheme ntree=ntree, # set number of trees nodesize=nodesize) # set node size all_model[[i]] <- rf_model # store model outputs } # # set dataframe names names(all_model) <- paste("ntrees:", params$ntrees, "nodesize:", params$nodesize) ## ---- fig.height=5, fig.width=5----------------------------------------------- # extract cross-validated r2 (proxy of test set r2) rsquared <- sapply(all_model, function(object) object$results["Rsquared"]) names(rsquared) <- paste("ntrees:", params$ntrees, "nodesize:", params$nodesize) rsq_df = data.frame(matrix(unlist(rsquared), nrow=length(rsquared), byrow=T)) rsq_df$nodesize = paste(params$nodesize) rsq_df$ntree = paste(params$ntrees) rsq_df$mean_r2 = rowMeans(rsq_df[1:5]) # # extract the training r2 (optional) store_validation <- vector("list", nrow(params)) for(i in 1:nrow(params)){ store_validation[[i]] <- list(predict(all_model[[i]], newdata=all_model[[i]][["trainingData"]][2:11]), all_model[[i]][["trainingData"]]$.outcome) } rsq_train = vector("list", nrow(params)) for(i in 1:nrow(params)){ rsq_train[[i]] <- cor(store_validation[[i]][[1]], store_validation[[i]][[2]])^2 } rsq_train_df = data.frame(matrix(unlist(rsq_train))) rsq_train_df$nodesize = paste(params$nodesize) rsq_train_df$ntree = paste(params$ntrees) names(rsq_train_df)[1] <- 'r2_train' # # merge training and cross-validated r2 into single data frame r2 <- merge(rsq_train_df, rsq_df, by=c('nodesize', 'ntree')) r2 <- r2[order(r2$r2_train, decreasing=T),] r2$row <- seq(1,18,1) plot(r2$row, r2$r2_train, col='red', ylim=c(0,1)) points(r2$row, r2$mean_r2, col='blue') ## ----------------------------------------------------------------------------- # set final hyper-parameter values ntrees_model <- c(500) nodesize_model <- c(55) tuneGrid_model <- expand.grid(.mtry=c(4)) # # define cross-validation scheme control_model <- trainControl(method="repeatedcv", number=10, repeats=5, verboseIter=F) # # fit the final random forest model rf_model_final <- train(cc_yield ~ ., data=train_data, method="rf", tuneGrid=tuneGrid_model, trControl=control_model, ntree=ntrees_model, nodesize=nodesize_model) ## ---- fig.height=5, fig.width=5----------------------------------------------- # create data frame for the tree train_data_tree <- train_data # # change variable names names(train_data_tree)[1] = c("Yield") # # fit the regression tree yield.model <- rpart(Yield ~ ., data=train_data_tree) # # get a good cp value by ploting the relative error vs. tree size # the figure shows cp of 0.015 or 0.019 being optimum plotcp(yield.model) # # fit the regression tree with the best cp value yield.model.opt = rpart(Yield ~ ., data=train_data_tree, cp=0.015) ## ---- fig.height=6, fig.width=7----------------------------------------------- # plot the regression tree rpart.plot(yield.model.opt, type=2, extra="auto", round=0, under=T, box.palette="BlGnYl") ## ---- fig.height=4.5, fig.width=5--------------------------------------------- # select xy variables X <- train_data[which(names(train_data) != "cc_yield")] y <- train_data$cc_yield # # check variable importance based on model rmse model_imp <- Predictor$new(rf_model_final, data=X, y=train_data$cc_yield) imp_rf <- FeatureImp$new(model_imp, loss="rmse") plot(imp_rf) ## ---- fig.height=4.5, fig.width=4.5------------------------------------------- # example partial dependency plot for nitrogen applied pdp_totN <- FeatureEffect$new(model_imp, feature=c("total_nitrogen_ha"), method="pdp") %>% plot() + ggtitle("PDP of Total N applied") pdp_totN ## ---- fig.height=4.5, fig.width=6--------------------------------------------- # overall two-way interaction strength ia <- Interaction$new(model_imp) # # interaction between production practices ia_disease_severity <- Interaction$new(model_imp,feature="insect_severity") # changed from disease severity... # # continuous x categorical variable pdp_dis_N <- FeatureEffect$new(model_imp, c("insect_severity", "total_nitrogen_ha"), method="pdp") %>% plot() + ggtitle("Disease incidence x N applied") pdp_dis_N # # continuous x continuous variable pdp_N_temp <- FeatureEffect$new(model_imp, c("total_nitrogen_ha", "Min_temperature"), method="pdp") %>% plot() + ggtitle("N applied x Minimum temperature") pdp_N_temp ## ---- fig.height=8, fig.width=8----------------------------------------------- surrigate_tree <- TreeSurrogate$new(model_imp) plot(surrigate_tree) ## ----------------------------------------------------------------------------- # fit model with local importance yield_model_localImp <- ranger::ranger(cc_yield ~ ., data=train_data, importance='permutation', local.importance=T, scale.permutation.importance=T, mtry=3) # # transform to data frame and print summary local_imp <- data.frame(yield_model_localImp$variable.importance.local) # # global importance data of each variable global_imp <- data.frame(ranger::importance(yield_model_localImp)) global_imp ## ---- fig.height=5, fig.width=5----------------------------------------------- plot(train_data$total_nitrogen_ha, yield_model_localImp$variable.importance.local[,"total_nitrogen_ha"], xlab="Nitrogen applied", ylab="Importance of Nitrogen", pch = 20) abline(h=ranger::importance(yield_model_localImp)["total_nitrogen_ha"], col='red') ## ---- fig.height=4.5, fig.width=6--------------------------------------------- # used model prediction X <- train_data[which(names(train_data) != "cc_yield")] y <- train_data$cc_yield model_imp <- Predictor$new(rf_model_final, data=X, y=train_data$cc_yield) # # get shapely values shap.explain <- iml::Shapley$new(model_imp, x.interest=X[1,]) plot(shap.explain) ## ---- fig.height=4.5, fig.width=6--------------------------------------------- ix.minyld <- which.min(y) shap.explain.minyld <- iml::Shapley$new(model_imp, x.interest=X[ix.minyld,]) plot(shap.explain.minyld)