# Script for 'Impact of rising temperature and management practices on global crop yields and food security', Agnolucci et al. (2020) # ----------------------------------------------------------------------------------------------------------------------------------- rm(list=ls()) library(xlsx) library(data.table) library("Formula") library("plm") library("lmtest") library(ggplot2) library(tools) setwd("ADD HERE") # # # Barley --------------------------------------------------------------------------------- myCrop <- "Barley" load(paste(myCrop, "YieldDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + as.factor(iso3_o):(as.numeric(year)), data = YieldWeatherUse, index=c("iso3_o", "year"), effect = "individual", model = "within" ) sel_mod_RobCoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2) + p ) ) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- sel_mod$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Cassava --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Cassava" load(paste(myCrop, "YieldIrrDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p +t:AAI + as.factor(iso3_o):(as.numeric(year)) + as.factor(iso3_o):I((as.numeric(year))^2), data = YieldWeatherUseIrr, index=c("iso3_o", "year"), model = "random", random.method = "amemiya", effect = "individual" ) sel_mod_robcoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_t_irr <- - ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]] ) / (2*sel_mod$coefficients[["I(t^2)"]]) sel_mod_sum <- data.table( Vars = c("T_Vertex", "T_Vertex_Irr", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_t_irr, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2) + p +t:AAI)), vcov=vcovBK) # sdev sdev_t <- mean(aggregate(YieldWeatherUseIrr$t, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUseIrr$p, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_t_global_irr <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]]) me_p_global <- sel_mod$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_t_global_irr_sd <- me_t_global_irr * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Cotton --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Cotton" load(paste(myCrop, "YieldFertDataset.RData",sep="")) sel_mod_ <- plm(log(yield) ~ t + I(t^2) + p + I(p^2) + as.numeric(year), data = YieldWeatherUseFert, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod_)$coeff v_t <- - sel_mod_coef["t",1] / (2*sel_mod_coef["I(t^2)",1] ) v_p <- - sel_mod_coef["p",1] / (2*sel_mod_coef["I(p^2)",1] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "P_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod_)$r.squared, summary(sel_mod_)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod_, update(sel_mod_, .~.-log(yield)- (t + I(t^2) + p + I(p^2)))) # sdev sdev_t <- mean(aggregate(YieldWeatherUseFert$t, list(YieldWeatherUseFert$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUseFert$p, list(YieldWeatherUseFert$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod_$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod_$coefficients[["t"]] me_p_global <- ( (2*sel_mod_$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod_$coefficients[["p"]]) *10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Groundnuts --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Groundnuts" load(paste(myCrop, "YieldIrrDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ + t + I(t^2) + p +t:AAI + p:AAI + as.factor(iso3_o):(as.numeric(year)) + as.factor(iso3_o):I((as.numeric(year))^2), data = YieldWeatherUseIrr, index=c("iso3_o", "year"), model = "random", random.method = "amemiya", effect = "individual" ) sel_mod_robcoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_t_irr <- - ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]] ) / (2*sel_mod$coefficients[["I(t^2)"]]) sel_mod_sum <- data.table( Vars = c("T_Vertex", "T_Vertex_Irr", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_t_irr, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p +t:AAI + p:AAI)), vcov=vcovBK) # sdev sdev_t <- mean(aggregate(YieldWeatherUseIrr$t, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUseIrr$p, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_t_global_irr <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]]) me_p_global <- sel_mod$coefficients[["p"]]*10 me_p_global_irr <- sel_mod$coefficients[["p"]]*10 + sel_mod$coefficients[["p:AAI"]]*10 me_t_global_sd <- me_t_global * sdev_t me_t_global_irr_sd <- me_t_global_irr * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 me_p_global_irr_sd <- me_p_global_irr * sdev_p / 10 # # # Maize --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Maize" load(paste(myCrop, "YieldIrrDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p +t:AAI + as.factor(iso3_o):(as.numeric(year)) + as.factor(iso3_o):I((as.numeric(year))^2), data = YieldWeatherUseIrr, index=c("iso3_o", "year"), effect = "individual", model = "random", random.method = "amemiya") sel_mod_robcoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_t_irr <- - ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]] ) / (2*sel_mod$coefficients[["I(t^2)"]]) sel_mod_sum <- data.table( Vars = c("T_Vertex", "T_Vertex_Irr", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_t_irr, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p +t:AAI)), vcov=vcovBK) # sdev sdev_t <- mean(aggregate(YieldWeatherUseIrr$t, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUseIrr$p, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_t_global_irr <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]]) me_p_global <- sel_mod$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_t_global_irr_sd <- me_t_global_irr * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # Millet --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Millet" load(paste(myCrop, "YieldDataset.RData",sep="")) sel_mod_ <- plm(log(yield) ~ t + I(t^2) + p + (as.numeric(year)) + I((as.numeric(year))^2), data = YieldWeatherUse, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod_)$coeff v_t <- - sel_mod_$coefficients[["t"]]/(2*sel_mod_$coefficients[["I(t^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod_)$r.squared, summary(sel_mod_)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod_, update(sel_mod_, .~.-log(yield)- ( t + I(t^2) + p ) ) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod_$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod_$coefficients[["t"]] me_p_global <- sel_mod_$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Oats --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Oats" load(paste(myCrop, "YieldDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ + t + I(t^2) + p + I(p^2) + as.factor(iso3_o):(as.numeric(year)) + as.factor(iso3_o):I((as.numeric(year))^2), data = YieldWeatherUse, index=c("iso3_o", "year"), model = "random", random.method = "walhus", effect = "individual" ) sel_mod_robcoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]]/(2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "P_Vertex" , "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2) + p + I(p^2))), vcov=vcovBK) # sdev sdev_t <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Potatoes --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Potatoes" load(paste(myCrop, "YieldPestDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + pest , data = YieldWeatherUsePest, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod)$coeff v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2) + p ))) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePest$t, list(YieldWeatherUsePest$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePest$p, list(YieldWeatherUsePest$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- sel_mod$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Pulses --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Pulses" load(paste(myCrop, "YieldPestDataset.RData",sep="")) sel_mod_ <- plm(log(yield) ~ t + I(t^2) + p + I(p^2) + pest + (as.numeric(year)) + I((as.numeric(year))^2), data = YieldWeatherUsePest, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod_)$coeff v_t <- - sel_mod_$coefficients[["t"]]/(2*sel_mod_$coefficients[["I(t^2)"]] ) v_p <- - sel_mod_$coefficients[["p"]]/(2*sel_mod_$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod_)$r.squared, summary(sel_mod_)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod_, update(sel_mod_, .~.-log(yield)- ( t + I(t^2) + p + I(p^2)) ) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePest$t, list(YieldWeatherUsePest$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePest$p, list(YieldWeatherUsePest$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod_$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod_$coefficients[["t"]] me_p_global <- ((2*sel_mod_$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod_$coefficients[["p"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Rapeseed --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Rapeseed" load(paste(myCrop, "YieldDataset.RData",sep="")) id <- which(YieldWeatherUse[,"iso3_o"] == "MNG") # drop outlier YieldWeatherUse <- YieldWeatherUse[-id,] sel_mod <- plm(log(yield) ~ t + I(t^2) + p + I(p^2) , data = YieldWeatherUse, index=c("iso3_o", "year"), model = "between", effect = "individual") sel_mod_coef <- summary(sel_mod)$coeff GCMeanT <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), mean)$x, na.rm = TRUE) GCMeanP <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), mean)$x, na.rm = TRUE) v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]]/(2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "P_Vertex" , "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p + I(p^2) ) ) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Rice------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Rice" load(paste(myCrop, "YieldPestDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2)+ p + I(p^2)+ pest , data = YieldWeatherUsePest, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod)$coeff v_t <- - sel_mod$coefficients[["t"]] / (2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]] / (2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "P_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2)+ p + I(p^2) ) ) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePest$t, list(YieldWeatherUsePest$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePest$p, list(YieldWeatherUsePest$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Rye------------------------------------------------------------------------------ rm(list=ls()) myCrop <- "Rye" load(paste(myCrop, "YieldIrrDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p +t:AAI + p:AAI , data = YieldWeatherUseIrr, index=c("iso3_o", "year"), model = "between", effect = "individual") sel_mod_coef <- summary(sel_mod)$coeff v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2) + p +t:AAI + p:AAI))) # sdev sdev_t <- mean(aggregate(YieldWeatherUseIrr$t, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUseIrr$p, list(YieldWeatherUseIrr$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_t_global_irr <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]]) me_p_global <- sel_mod$coefficients[["p"]]*10 me_p_global_irr <- ( sel_mod$coefficients[["p"]] + sel_mod$coefficients[["p:AAI"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_t_global_irr_sd <- me_t_global_irr * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 me_p_global_irr_sd <- me_p_global_irr * sdev_p / 10 # # # Sorghum --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Sorghum" load(paste(myCrop, "YieldDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + I(p^2) + as.factor(iso3_o):(as.numeric(year)) + as.factor(iso3_o):I((as.numeric(year))^2), data = YieldWeatherUse, index=c("iso3_o", "year"), effect = "individual", model = "within" ) sel_mod__RobCoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] CTREND_FE <- data.frame(fixef(sel_mod) ) v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]]/(2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- (t + I(t^2)+ p + I(p^2) ) ) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Soybeans --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Soybeans" load(paste(myCrop, "YieldDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + (as.numeric(year)) + I((as.numeric(year))^2) , data = YieldWeatherUse, index=c("iso3_o", "year"), model = "fd", effect = "individual") se_mod_robcoef <- coeftest(sel_mod, vcov = vcovHC(sel_mod, method = "arellano", type ="HC3", cluster = "group"))[,] v_t <- - sel_mod$coefficients[["t"]] / (2*sel_mod$coefficients[["I(t^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p )), vcov=vcovHC) # sdev sdev_t <- mean(aggregate(YieldWeatherUse$t, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUse$p, list(YieldWeatherUse$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- sel_mod$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Sugarbeet --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Sugarbeet" load(paste(myCrop, "YieldPestFertDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2)+ p +I(p^2) + pest + fert , data = YieldWeatherUsePestFert, index=c("iso3_o", "year"), model = "between", effect ="individual") sel_mod_coef <- summary(sel_mod)$coeff # # Create Data Table v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]]/(2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p , summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2)+ p +I(p^2))) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePestFert$t, list(YieldWeatherUsePestFert$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePestFert$p, list(YieldWeatherUsePestFert$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) * 10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Sunflower --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Sunflower" load(paste(myCrop, "YieldPestFertDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + pest + fert , data = YieldWeatherUsePestFert, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod)$coeff v_t <- - sel_mod$coefficients[["t"]] / (2*sel_mod$coefficients[["I(t^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p )) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePestFert$t, list(YieldWeatherUsePestFert$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePestFert$p, list(YieldWeatherUsePestFert$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- sel_mod$coefficients[["p"]]*10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # SweetPotatoes --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "SweetPotatoes" load(paste(myCrop, "YieldPestFertDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + I(p^2) + pest + fert , data = YieldWeatherUsePestFert, index=c("iso3_o", "year"), model = "between", effect = "individual") sel_mod_coef <- summary(sel_mod)$coeff v_t <- - sel_mod$coefficients[["t"]]/(2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]]/(2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "P_Vertex" , "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p + I(p^2) )) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePestFert$t, list(YieldWeatherUsePestFert$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePestFert$p, list(YieldWeatherUsePestFert$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) *10 me_t_global_sd <- me_t_global * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 # # # Wheat --------------------------------------------------------------------------------- rm(list=ls()) myCrop <- "Wheat" load(paste(myCrop, "YieldPestFertIrDataset.RData",sep="")) sel_mod <- plm(log(yield) ~ t + I(t^2) + p + I(p^2) + pest + fert +t:AAI + p:AAI, data = YieldWeatherUsePestFertIr, index=c("iso3_o", "year"), model = "between", effect ="time") sel_mod_coef <- summary(sel_mod)$coeff v_t <- - sel_mod$coefficients[["t"]] / (2*sel_mod$coefficients[["I(t^2)"]] ) v_p <- - sel_mod$coefficients[["p"]] / (2*sel_mod$coefficients[["I(p^2)"]] ) sel_mod_sum <- data.table( Vars = c("T_Vertex", "P_Vertex", "R2", "R2adj", "GCMeanT", "GCMeanP"), Value = c(v_t, v_p, summary(sel_mod)$r.squared, summary(sel_mod)$adj.r.squared, GCMeanT, GCMeanP) ) # Joint restriction Wald test waldtest(sel_mod, update(sel_mod, .~.-log(yield)- ( t + I(t^2) + p + I(p^2) +t:AAI + p:AAI )) ) # sdev sdev_t <- mean(aggregate(YieldWeatherUsePestFertIr$t, list(YieldWeatherUsePestFertIr$iso3_o), sd)$x, na.rm = TRUE) sdev_p <- mean(aggregate(YieldWeatherUsePestFertIr$p, list(YieldWeatherUsePestFertIr$iso3_o), sd)$x, na.rm = TRUE) # # Marginal effect me_t_global <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + sel_mod$coefficients[["t"]] me_t_global_irr <- (2*sel_mod$coefficients[["I(t^2)"]] ) * GCMeanT + ( sel_mod$coefficients[["t"]] + sel_mod$coefficients[["t:AAI"]]) me_p_global <- ((2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + sel_mod$coefficients[["p"]]) *10 me_p_global_irr <- ( (2*sel_mod$coefficients[["I(p^2)"]] ) * GCMeanP + ( sel_mod$coefficients[["p"]] + sel_mod$coefficients[["p:AAI"]])) * 10 me_t_global_sd <- me_t_global * sdev_t me_t_global_irr_sd <- me_t_global_irr * sdev_t me_p_global_sd <- me_p_global * sdev_p / 10 me_p_global_irr_sd <- me_p_global_irr * sdev_p / 10