library(tidyverse)
library(dplyr)
library(metafor)
library(ggplot2)
library(gsheet)
library(cowplot)
library(broom)
library(writexl)
library(readxl)
library(ggthemes)
library(minpack.lm)Economic returns on fungicide use across soybean cultivars with varying tolerance to target spot
Package
Importation
ma <- gsheet2tbl("https://docs.google.com/spreadsheets/d/151DE26uSMN4WBS0PTQcjdw7BO56gXy_VilFZS3b9AAM/edit?usp=sharing")
molina = read_xlsx("data/molina.xlsx")unique(ma$cultivar) [1] "M8210 IPRO" "CD 2728 IPRO" "BMX Potência RR" "NA 6700 IPRO"
[5] "NA 5909 RG" "TMG 803" "TMG 2181 IPRO" "M 9144RR"
[9] "8579 RSF IPRO" "M 8473 IPRO" "TEC 7849 IPRO" "BMX Garra RR"
[13] "63i64RSF IPRO" "8473RSF" "TMG 7067 IPRO" "P98Y70"
[17] "CD 2827 IPRO" "M 6410 IPRO" "Juruena" "GUAIA 7487 RR"
[21] "W791 RR" "Garra IPRO" "DM 75i76 IPRO" "74178 RSF IPRO"
[25] "M8349 IPRO" "M9144 RR" "GA 76 IPRO" "Ima 84114"
[29] "NS 8338 IPRO" "BS 2606 IPRO" "TMG 2383 IPRO" "DM73i75 IPRO"
[33] "C2827 IPRO" "TMG 2381 IPRO" "RK 7518 IPARO" "CZ48B32"
[37] "M8372" "M8210" "GA76 IPRO" "CZ58B28"
[41] "CZ 37B43 IPRO" "HO MARACAÍ IPRO" "77i79RSF IPRO" "BMX DOMINIO IPRO"
[45] NA "M5947 IPRO" "80182RSF IPRO" "M220i2X"
[49] "BÔNUS" "TMG 2285 IPRO" ma |>
filter(!cultivar %in% c("63i64RSF IPRO","8579 RSF IPRO",
"CD 2728 IPRO","DM 75i76 IPRO","M8210 IPRO","TMG 2181 IPRO")) %>%
filter(!is.na(sev)) |>
ggplot(aes(sev))+
geom_histogram(fill = "black", color = "white")+
ggthemes::theme_few()+
#theme_minimal()+
#facet_wrap(~year)+
labs(x = "TS Severity (%)",
y = "Count")+
facet_wrap(~cultivar)+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
Cluster
ma1 = ma %>%
dplyr::select(study, cultivar, sev,yld,block,ai,state,year)
molina1 = molina %>%
dplyr::select(study, cultivar, sev,yield,rep,fungic,state,year)
colnames(molina1) = c("study","cultivar","sev","yld","block","ai","state","year")
ma_c = rbind(ma1,molina1)
unique(ma_c$state)[1] "MT" "GO" "MS" "TO" "BA" "PR" "DF"ma_c = ma_c %>%
filter(!is.na(yld)) %>%
filter(!is.na(sev)) %>%
filter(!is.na(cultivar))
ma_c_reg = ma_c %>%
#dplyr::select(active_ingredient, block, sev,yld) %>%
group_by(cultivar) %>%
do(tidy(lm(.$yld ~ .$sev)))
ma_c_regma_c_reg = ma_c_reg |>
filter(term %in% c("(Intercept)",".$sev"))
ma_c_reg[ma_c_reg$term== "(Intercept)",c("parameters")] <- "Intercept"
ma_c_reg[ma_c_reg$term== ".$sev",c("parameters")] <- "Slope"
b1_c = ma_c_reg |>
dplyr::filter(parameters == "Slope") |>
#filter(estimate < 0) %>%
dplyr::group_by(cultivar) |>
summarise(
Slope = estimate
)
b0_c = ma_c_reg |>
filter(parameters == "Intercept") |>
group_by(cultivar) |>
summarise(
Intercept = estimate
)
parameters_c = cbind(b1_c,b0_c)
parameters_c[,3] = NULL
parameters_c = parameters_c %>%
filter(!cultivar %in% c("BÔNUS","M9144_RR","74178 RSF IPRO"))
dc_c = parameters_c %>%
group_by(cultivar) %>%
summarise(
DC = (Slope/Intercept)*-1,
b1= Intercept/1000,
DC = DC,
DC_perc = DC*100
)
set.seed(123)
dados_dummy <- model.matrix(~ cultivar - 1, data = dc_c)
set.seed(123)
dados_preparados <- cbind(dados_dummy, DC = dc_c$DC)
set.seed(123)
resultado_kmeans <- kmeans(dados_preparados, centers = 3)library(factoextra)
p = fviz_cluster(resultado_kmeans, data = dados_preparados, label = NULL,geom = "point",pointsize = 3, shape = 20)
p2 = p + theme_few() +
scale_color_manual(values = c('#9ccb86','#cf597e','steelblue'))+
scale_fill_manual(values = c('#9ccb86','#cf597e','steelblue'))+
theme(
axis.title = element_text(size = 14, face = "bold"),
axis.text = element_text(size = 12)) +
labs(
title = "",
x = "PC1",
y = "PC2")+
theme(legend.position = "none")
p2
ggsave("fig/k-means.png", dpi = 600, width = 8, height = 6)dc_c$Cluster <- resultado_kmeans$cluster
dc_c %>%
ggplot(aes(DC))+
geom_histogram(color = "white",fill = "steelblue")+
facet_wrap(~Cluster)+
#scale_x_continuous(breaks = c(0, 0.1, 0.2, 0.3), limits = c(0, 0.3))+
theme_minimal()+
theme(text= element_text(size = 14, face = "bold"))
c_01 = dc_c %>%
filter(Cluster == 1)
c_02 = dc_c %>%
filter(Cluster == 2)
c_03 = dc_c %>%
filter(Cluster == 3)
c_01$cultivar [1] "80182RSF IPRO" "8579 RSF IPRO" "AS3730_IPRO" "BMX DOMINIO IPRO"
[5] "BMX_Ativa_RR" "BRSGO_9160_RR" "C2827 IPRO" "CZ58B28"
[9] "GA76 IPRO" "Garra IPRO" "HO MARACAÍ IPRO" "Juruena"
[13] "M8210" "M8210 IPRO" "M8210_IPRO" "M8372"
[17] "NA 6700 IPRO" "RK 7518 IPARO" "TMG 2383 IPRO" "TMG2281_IPRO"
[21] "TMG803" "W791 RR" c_02$cultivar[1] "BRSGO_8151_RR" "BRSGO_8661_RR" "GA 76 IPRO" "M 9144RR"
[5] "M8349 IPRO" "M9144 RR" "Syn1180_RR" "TMG 2285 IPRO"c_03$cultivar [1] "5G830_RR" "63i64RSF IPRO" "77i79RSF IPRO" "8473RSF"
[5] "BMX Garra RR" "BMX Potência RR" "BMX_Potencia_RR" "BS 2606 IPRO"
[9] "CD 2728 IPRO" "CD 2827 IPRO" "CZ 37B43 IPRO" "CZ48B32"
[13] "DM 75i76 IPRO" "DM73i75 IPRO" "GUAIA 7487 RR" "Ima 84114"
[17] "M 6410 IPRO" "M 8473 IPRO" "M220i2X" "M5947 IPRO"
[21] "M8336_RR" "NA 5909 RG" "NA_5909_RG" "NS 8338 IPRO"
[25] "TEC 7849 IPRO" "TMG 2181 IPRO" "TMG 2381 IPRO" "TMG 7067 IPRO"
[29] "TMG 803" "TMG1179_RR" "TMG1180_RR" Number of trials by cultivar
study_counts <- ma_c %>%
filter(!cultivar %in% c("BÔNUS","74178 RSF IPRO")) %>%
group_by(cultivar,year) %>%
summarise(n_study = n_distinct(study)) %>%
group_by(cultivar) %>%
summarise(n_study = sum(n_study))
write_xlsx(study_counts,"data/study_counts.xlsx")sum(study_counts$n_study)[1] 173Number of trials
count_study = ma_c %>%
#filter(!cultivar %in% c("BÔNUS","M9144_RR","74178 RSF IPRO")) %>%
group_by(study) %>%
summarise(n_study = n_distinct(year))sum(count_study$n_study)[1] 176Descriptive analysis
Overall
sev = ma |>
filter(!is.na(sev)) |>
ggplot(aes(sev))+
geom_histogram(fill = "black", color = "white")+
#D2691E
ggthemes::theme_few()+
#theme_minimal()+
#facet_wrap(~year)+
labs(x = "TS Severity (%)",
y = "Count")+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
#00798C
#236E96
sev_log = ma |>
filter(!is.na(sev)) |>
ggplot(aes(log10(sev)))+
geom_histogram(fill = "black", color = "white")+
#theme_minimal()+
ggthemes::theme_few()+
#facet_wrap(~year)+
labs(x = "log10(TS severity)",
y = "Count")+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
yld = ma |>
filter(!is.na(yld)) |>
ggplot(aes(yld))+
geom_histogram(fill = "black", color = "white")+
#scale_y_continuous(breaks = c(0, 200, 400, 600,800), limits = c(0, 800))+
ggthemes::theme_few()+
#theme_minimal()+
#facet_wrap(~year)+
labs(x = "Yield (Kg/ha)",
y = "Count")+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
#FF9933
plot_grid(sev,sev_log,yld, labels = c("A","B","C"), nrow = 3, align = "V")
ggsave("fig/overall.png", bg = "white", dpi = 600, height = 9, width = 10)Cultivar
Severity
sev_cultivar = ma |>
filter(ai == "_CHECK") %>%
filter(cultivar %in% c("63i64RSF IPRO","859 RSF IPRO",
"CD 2728 IPRO","DM 75i76 IPRO","M8210 IPRO","TMG 2181 IPRO","HO MARACAÍ IPRO")) %>%
filter(!is.na(sev)) |>
ggplot(aes(sev))+
geom_histogram(fill = "black", color = "white", bins = 10)+
#theme_minimal()+
ggthemes::theme_few()+
#facet_wrap(~year)+
labs(x = "TS severity",
y = "Count")+
#facet_wrap(~condition)+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
sev_cultivar
ggsave("fig/cultivar_severity.png", bg = "white", dpi = 600, height = 9, width = 10)Yield
yld_cultivar = ma |>
filter(cultivar %in% c("63i64RSF IPRO","879 RSF IPRO",
"CD 2728 IPRO","DM 75i76 IPRO","M8210 IPRO","TMG 2181 IPRO","HO MARACAÍ IPRO")) %>%
filter(!is.na(yld)) |>
ggplot(aes(yld))+
geom_histogram(fill = "black", color = "white", bins = 10)+
#scale_y_continuous(breaks = c(0, 200, 400, 600,800), limits = c(0, 800))+
ggthemes::theme_few()+
#theme_minimal()+
#facet_wrap(~year)+
labs(x = "Yield (Kg/ha)",
y = "Count")+
#facet_wrap(~condition)+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
yld_cultivar
ggsave("fig/cultivar_yield.png", bg = "white", dpi = 600, height = 9, width = 10)Active ingredients
Severity
sev_ai = ma |>
filter(ai %in% c("BIX+PROT+TRIFL","AZOX+TEBU+MANC","FLUX+PYRA","_CHECK")) %>%
filter(!is.na(sev)) |>
ggplot(aes(sev))+
geom_histogram(fill = "black", color = "white", bins = 10)+
#theme_minimal()+
ggthemes::theme_few()+
#facet_wrap(~year)+
labs(x = "TS severity",
y = "Count")+
facet_wrap(~ai)+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
sev_ai
ggsave("fig/ai_severity.png", bg = "white", dpi = 600, height = 9, width = 10)Yield
yld_ai = ma |>
filter(ai %in% c("BIX+PROT+TRIFL",
"AZOX+TEBU+MANC","FLUX+PYRA","_CHECK")) %>%
filter(!is.na(yld)) |>
ggplot(aes(yld))+
geom_histogram(fill = "black", color = "white", bins = 15)+
#scale_y_continuous(breaks = c(0, 200, 400, 600,800), limits = c(0, 800))+
ggthemes::theme_few()+
#theme_minimal()+
#facet_wrap(~year)+
labs(x = "Yield (Kg/ha)",
y = "Count")+
facet_wrap(~ai)+
theme(axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))
yld_ai
ggsave("fig/ai_yield.png", bg = "white", dpi = 600, height = 9, width = 10)Cultivar x Active ingredients
Severity
sev_cultivar_ai = ma |>
filter(!is.na(sev)) %>%
filter(cultivar %in% c("63i64RSF IPRO","8579 RSF IPRO","CD 2728 IPRO", "DM 75i76 IPRO","M8210 IPRO","TMG 2181 IPRO")) %>%
filter(ai %in% c("BIX+PROT+TRIFL","AZOX+TEBU+MANC","FLUX+PYRA")) %>%
ggplot(aes(reorder(ai,+sev),sev))+
#geom_jitter(color = "black",alpha = .2, width = 0.09, size = 2)+
geom_boxplot(fill = NA, color = "steelblue",outlier.colour = NA,position = position_dodge(width = 0.95))+
theme_minimal()+
coord_flip()+
#facet_wrap(~condition)+
labs(x = "",
y = "TS Severity (%)")+
theme(strip.text = element_text(size = 12, face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))Yield
yld_cultivar_ai = ma |>
filter(!is.na(sev)) %>%
filter(cultivar %in% c("63i64RSF IPRO","8579 RSF IPRO","CD 2728 IPRO", "DM 75i76 IPRO","M8210 IPRO","TMG 2181 IPRO")) %>%
filter(ai %in% c("BIX+PROT+TRIFL","AZOX+TEBU+MANC","FLUX+PYRA")) %>%
ggplot(aes(reorder(ai,+yld),yld))+
#geom_jitter(color = "black",alpha = .2, width = 0.09, size = 2)+
geom_boxplot(fill = NA, color = "steelblue",outlier.colour = NA,position = position_dodge(width = 0.95))+
theme_minimal()+
coord_flip()+
#facet_wrap(~condition)+
labs(x = "",
y = "Yield (kg/ha)")+
theme(strip.text = element_text(size = 12, face = "bold"),
axis.title.x = element_text(face = "bold"),
axis.title.y = element_text(face = "bold"))Framework simulation
gera.norm.bid.geral<-function(tamanho.amostra,correlacao,m1,m2,sigma1,sigma2)
{
ro<-correlacao
n<-tamanho.amostra
x<-matrix(0,n,2)
for (i in 1:n)
{x[i,1]<-rnorm(1,m1,sigma1)
x[i,2]<-rnorm(1,m2+ro*sigma1/sigma2*(x[i,1]-m1),sigma2*(sqrt(1-ro^2)))
}
return(x)
}
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 2487387 132.9 4366044 233.2 4366044 233.2
Vcells 4919307 37.6 10146329 77.5 8387619 64.0HIGH-TOLERANCE
res = ma_c %>%
filter(cultivar %in% c_03$cultivar)
unique(res$cultivar) [1] "CD 2728 IPRO" "BMX Potência RR" "NA 5909 RG" "TMG 803"
[5] "TMG 2181 IPRO" "M 8473 IPRO" "TEC 7849 IPRO" "BMX Garra RR"
[9] "63i64RSF IPRO" "8473RSF" "TMG 7067 IPRO" "CD 2827 IPRO"
[13] "M 6410 IPRO" "GUAIA 7487 RR" "DM 75i76 IPRO" "Ima 84114"
[17] "NS 8338 IPRO" "BS 2606 IPRO" "DM73i75 IPRO" "TMG 2381 IPRO"
[21] "CZ48B32" "CZ 37B43 IPRO" "77i79RSF IPRO" "M5947 IPRO"
[25] "M220i2X" "5G830_RR" "BMX_Potencia_RR" "M8336_RR"
[29] "NA_5909_RG" "TMG1179_RR" "TMG1180_RR" res_regression = res %>%
group_by(study) %>%
do(tidy(lm(.$yld ~ .$sev)))
res_regression = as.data.frame(res_regression)
res_regressionres_severity = res |>
filter(ai == "_CHECK") %>%
ggplot(aes(sev))+
geom_histogram(color = "white",fill = "steelblue", bins = 10)+
scale_x_continuous(limits = c(0, 100),
breaks = seq(0, 1000, by = 10))+
theme_bw()+
labs(x = "Severity (%)",
y = "Density")
res_yield = res |>
filter(ai == "_CHECK") %>%
ggplot(aes(yld))+
geom_histogram(color = "white",fill = "brown", bins = 10)+
scale_x_continuous(limits = c(2000, 7000),
breaks = seq(2000, 7000, by = 500))+
theme_bw()+
labs(x = "Yield (kg/ha)",
y = "")
plot_grid(res_severity,res_yield)
res_regression = res_regression |>
filter(term %in% c("(Intercept)",".$sev"))
res_regression[res_regression$term== "(Intercept)",c("parameters")] <- "Intercept"
res_regression[res_regression$term== ".$sev",c("parameters")] <- "Slope"
res_b0 = res_regression |>
filter(parameters == "Intercept")
res_b1= res_regression %>%
filter(parameters == "Slope") %>%
filter(!estimate > 0)
mean(res_b0$estimate)[1] 3902.37sd(res_b0$estimate)[1] 662.0473mean(res_b1$estimate)[1] -14.06671sd(res_b1$estimate)[1] 10.05459Intercept
intercept_tss_res = res_regression |>
filter(parameters == "Intercept")
mean_intercept_res = mean(intercept_tss_res$estimate)
sd_intercept_res = sd(intercept_tss_res$estimate)
plot(ecdf(intercept_tss_res$estimate))
curve(pnorm(x, mean_intercept_res,sd_intercept_res), add = T, col = "red")
Kolmogorov-Smirnov Test
Fx_res_b0 = environment(ecdf(intercept_tss_res$estimate))$y
x_res_b0= environment(ecdf(intercept_tss_res$estimate))$x
ks.test(Fx_res_b0, pnorm(x_res_b0, mean(intercept_tss_res$estimate), sd(intercept_tss_res$estimate)))
Exact two-sample Kolmogorov-Smirnov test
data: Fx_res_b0 and pnorm(x_res_b0, mean(intercept_tss_res$estimate), sd(intercept_tss_res$estimate))
D = 0.060606, p-value = 1
alternative hypothesis: two-sidedplot(Fx_res_b0, pnorm(x_res_b0, mean(intercept_tss_res$estimate), sd(intercept_tss_res$estimate)))
abline(a=0,b=1)
intercep_plot_res = intercept_tss_res %>%
ggplot(aes(estimate))+
geom_histogram(aes(y = ..density..), bins = 8,color = "white", fill = "#FF4917")+
#stat_function(fun=function(x) dnorm(x, mean_intercept_res, sd_intercept_res), color ="black", size =1.2)+
theme_bw()+
labs(x="Intercept (kg/ha)", y = "Density")
intercep_plot_res
Slope
slope_tss_res = res_regression |>
filter(parameters == "Slope") %>%
filter(!estimate > 0)
mean_slope = mean(slope_tss_res$estimate)
sd_slope = sd(slope_tss_res$estimate)
plot(ecdf(slope_tss_res$estimate))
curve(pnorm(x, mean_slope, sd_slope), add = T, col = "red")
Kolmogorov-Smirnov Test
library(minpack.lm)
Fx =environment(ecdf(-slope_tss_res$estimate))$y
x = environment(ecdf(-slope_tss_res$estimate))$x
slope_reg = nlsLM(Fx ~ pgamma(x, shape, rate,log = FALSE) ,
start = c(shape = 2.5, rate = 0.13),
control = nls.lm.control(maxiter = 1024))
summary(slope_reg)
Formula: Fx ~ pgamma(x, shape, rate, log = FALSE)
Parameters:
Estimate Std. Error t value Pr(>|t|)
shape 1.59411 0.18559 8.589 6.28e-09 ***
rate 0.10944 0.01482 7.382 9.85e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.05392 on 25 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 1.49e-08shape_res = summary(slope_reg)$coef[1]
rate_res = summary(slope_reg)$coef[2]
Fx =environment(ecdf(-slope_tss_res$estimate))$y
x = environment(ecdf(-slope_tss_res$estimate))$x
ks.test(Fx, pgamma(x, shape_res, rate_res))
Exact two-sample Kolmogorov-Smirnov test
data: Fx and pgamma(x, shape_res, rate_res)
D = 0.14815, p-value = 0.9357
alternative hypothesis: two-sidedshape_res = summary(slope_reg)$coef[1]
rate_res = summary(slope_reg)$coef[2]
slope_plot_res = slope_tss_res %>%
ggplot(aes(estimate))+
geom_histogram(aes(y = ..density..),bins = 20, color = "white", fill = "#159EE6")+
scale_x_continuous(breaks = c(-200,-150,-100,-50,0),
limits = c(-200,0))+
theme_bw()+
labs(x="Slope (kg/p.p.)", y = "Density")
slope_plot_res
Correlation
slope_coefficient_res = res_regression |>
dplyr::filter(parameters == "Slope") |>
#filter(estimate < 0) %>%
dplyr::group_by(study) |>
summarise(
Slope = estimate
)
slope_coefficient_res[,1] = NULL
slope_coefficient_res |>
filter(!Slope == "NA") |>
summarise(
mean = mean(Slope))intercept_coefficient_res = res_regression |>
filter(parameters == "Intercept") |>
group_by(study) |>
summarise(
Intercept = estimate
)
mean(intercept_coefficient_res$Intercept)[1] 3902.37sd(intercept_coefficient_res$Intercept)[1] 662.0473regression_res = cbind(slope_coefficient_res,intercept_coefficient_res)
regression_res = regression_res %>%
filter(Slope <0)
correlation <- function(x, y) {
if (length(x) != length(y)) {
stop("Os vetores 'x' e 'y' devem ter o mesmo comprimento")
}
cor_result <- cor(x, y, use = "complete.obs")
return(cor_result)
}
correlation(regression_res$Slope, regression_res$Intercept)[1] -0.518309corr_actual_plot_res = regression_res %>%
mutate(alfa =(Slope/Intercept)*100) %>%
#filter(alfa > -3 & alfa < 0) %>%
ggplot(aes(Intercept,Slope))+
geom_point( size =2, color = "black")+
# geom_density_2d(color = "black")+
geom_smooth(method = lm, color = "red", se = F, size = 1.2)+
#scale_x_continuous(breaks = c(3500,4000,4500,5000,5500), limits = c(3500, 5500))+
scale_y_continuous(breaks = c(-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-160, 0))+
ggthemes::theme_few()+
labs(y= "Slope",
x = "Intercept",
title = "Empirical: r = -0.518")+
theme(text = element_text(face = "bold"))
# coord_cartesian(
# xlim = c(0,10000))
corr_actual_plot_res
j_res<-gera.norm.bid.geral(10000,0.51,0,0,1,1)
plot(j_res[,1],j_res[,2])
Severity
Kolmogorov-Smirnov Test
sev_res = ma_c %>%
filter(cultivar %in% c_03$cultivar) %>%
filter(ai %in% c("ZZ_CHECK","_CHECK"))
sev_r = sev_res$sev
Fx_res = environment(ecdf(sev_r))$y
x_res = environment(ecdf(sev_r))$x/100
summary(nlsLM(Fx_res ~ pbeta(x_res, shape1, shape2, log = FALSE) ,
start = c(shape1 = 1, shape2 = 1),
control = nls.lm.control(maxiter = 100000)))
Formula: Fx_res ~ pbeta(x_res, shape1, shape2, log = FALSE)
Parameters:
Estimate Std. Error t value Pr(>|t|)
shape1 3.42200 0.09256 36.97 <2e-16 ***
shape2 4.75064 0.13386 35.49 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03186 on 172 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 1.49e-08ks.test(Fx_res,pbeta(x_res,3.42200,4.75064))
Asymptotic two-sample Kolmogorov-Smirnov test
data: Fx_res and pbeta(x_res, 3.422, 4.75064)
D = 0.063218, p-value = 0.8777
alternative hypothesis: two-sidedPlot
plot(x_res,Fx_res)
curve(pbeta(x, 3.42200,4.75064),0,1, add = T)
sev_dist_plot_res = sev_res %>%
ggplot(aes(sev/100))+
geom_histogram(aes(y = ..density..),bins = 8, color = "white", fill = "darkgreen")+
scale_x_continuous(breaks = c(1.00,0.75,0.50,0.25,0.00),
limits = c(0.00,1.00))+
theme_bw()+
labs(x="Severity (proportion)", y = "Density")
sev_dist_plot_res
Relative yield loss
empiric_ryl_res= regression_res %>%
mutate(cc = (Slope/Intercept)*100)
# filter( cc > -3 & cc <0 )
head(empiric_ryl_res)stat_emp_ryl_res =empiric_ryl_res%>%
summarise(data = "empirical",
mean = mean(cc),
median = median(cc),
variance = var(cc))real_RYL_res= empiric_ryl_res %>%
ggplot(aes(cc))+
geom_histogram(aes(y = ..density..), bins = 12, color = "white", fill = "black")+
# scale_y_continuous(limits = c(0,1.4),breaks = seq(0, 1.4,by = 0.4))+
geom_vline(data = stat_emp_ryl_res, aes(xintercept = mean, color = "Mean"),size =1)+
geom_vline(data = stat_emp_ryl_res, aes(xintercept = median, color = "Median"),size =1)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
scale_color_calc()+
labs(x = "",
y = "Frequency",
color = "",
title = "Empirical | High")+
# xlim(-3,0.2)+
theme(legend.position = "none",
text = element_text(face = "bold"))
real_RYL_res
mean(res_b0$estimate)[1] 3902.37sd(res_b0$estimate)[1] 662.0473b0_res = pnorm(j_res[,2])
b0_t_res = qnorm(b0_res, 3902.37, 662*sqrt(33))
hist(b0_t_res, prob = T)
curve(dnorm(x, mean_intercept_res, sd_intercept_res), add = T)
mean(res_b0$estimate)[1] 3902.37sd(res_b0$estimate)[1] 662.0473b1_res = pnorm(j_res[,1])
b1_t_res = qgamma(b1_res, shape = shape_res, rate = rate_res)*-1
hist(b1_t_res, prob = T)
curve(dgamma(-x,shape=shape_res, rate = rate_res), add = T)
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")
correlation(b1_t_res, b0_t_res)[1] -0.4762788cor.test(b1_t_res, b0_t_res, method = "pearson")
Pearson's product-moment correlation
data: b1_t_res and b0_t_res
t = -54.161, df = 9998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4912926 -0.4609820
sample estimates:
cor
-0.4762788 actual_cc_res = regression_res %>%
mutate(cc = (Slope/Intercept)*100)
simulated_cc_res = data.frame(b0_t_res, b1_t_res,
cc = (b1_t_res/b0_t_res)*100)%>%
filter( cc > -3 & cc <0 )
ks.test(actual_cc_res$cc, simulated_cc_res$cc)
Asymptotic two-sample Kolmogorov-Smirnov test
data: actual_cc_res$cc and simulated_cc_res$cc
D = 0.14465, p-value = 0.6265
alternative hypothesis: two-sidedfx_actual_res = environment(ecdf(actual_cc_res$cc))$y
fx_simu_res = environment(ecdf(simulated_cc_res$cc))$y
ks.test((fx_actual_res),(fx_simu_res))
Asymptotic two-sample Kolmogorov-Smirnov test
data: (fx_actual_res) and (fx_simu_res)
D = 0.037033, p-value = 1
alternative hypothesis: two-sidedecdf_damage_res = ggplot()+
stat_ecdf(aes(simulated_cc_res$cc,color = "Simulated"), size=1.2,geom = "step")+
stat_ecdf(aes(actual_cc_res$cc, color = "Empirical"), size=1.2,geom = "step")+
ggthemes::theme_few()+
scale_color_manual(values = c("black","orange"))+
labs(y = "Probability",
x = "Relative percent yield loss",
color ="",
title = "CDF | High")+
theme(legend.position = "top",
legend.background = element_blank(),
text = element_text(face = "bold"))
ecdf_damage_res
corr_sim_plot_res = data.frame(b0_t_res,b1_t_res) %>%
mutate(alfa =(b1_t_res/b0_t_res)*100 ) %>%
#filter(alfa > -3 & alfa < 0) %>%
ggplot(aes(b0_t_res,b1_t_res))+
geom_point( size =2, color = "orange", alpha =0.3)+
# geom_density_2d(color = "black")+
geom_smooth(method = lm, color = "red", se = F, size = 1.2)+
scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000), limits = c(0, 10000))+
scale_y_continuous(breaks = c(-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-160, 0))+
ggthemes::theme_few()+
labs(y= "Slope",
x = "Intercept",
title = "Simulated: r = -0.476")+
coord_cartesian(
xlim = c(0,7000),
ylim = c(-160,0)
)+
theme(text = element_text(face = "bold"))
corr_sim_plot_res
MEDIUM-TOLERANCE
m_sus = ma_c %>%
filter(cultivar %in% c_01$cultivar)
unique(m_sus$cultivar) [1] "M8210 IPRO" "NA 6700 IPRO" "8579 RSF IPRO" "Juruena"
[5] "W791 RR" "Garra IPRO" "TMG 2383 IPRO" "C2827 IPRO"
[9] "RK 7518 IPARO" "M8372" "M8210" "GA76 IPRO"
[13] "CZ58B28" "HO MARACAÍ IPRO" "BMX DOMINIO IPRO" "80182RSF IPRO"
[17] "TMG803" "BRSGO_9160_RR" "BMX_Ativa_RR" "AS3730_IPRO"
[21] "M8210_IPRO" "TMG2281_IPRO" m_sus_regression = m_sus %>%
group_by(study) %>%
do(tidy(lm(.$yld ~ .$sev)))
m_sus_regression = as.data.frame(m_sus_regression)
m_sus_regressionm_sus_severity = m_sus |>
filter(ai == "_CHECK") %>%
ggplot(aes(sev))+
geom_histogram(color = "white",fill = "steelblue", bins = 10)+
scale_x_continuous(limits = c(0, 100),
breaks = seq(0, 1000, by = 10))+
theme_bw()+
labs(x = "Severity (%)",
y = "Density")
m_sus_yield = m_sus |>
filter(ai == "_CHECK") %>%
ggplot(aes(yld))+
geom_histogram(color = "white",fill = "brown", bins = 10)+
scale_x_continuous(limits = c(2000, 7000),
breaks = seq(2000, 7000, by = 500))+
theme_bw()+
labs(x = "Yield (kg/ha)",
y = "")
plot_grid(m_sus_severity,m_sus_yield)
m_sus_regression = m_sus_regression |>
filter(term %in% c("(Intercept)",".$sev"))
m_sus_regression[m_sus_regression$term== "(Intercept)",c("parameters")] <- "Intercept"
m_sus_regression[m_sus_regression$term== ".$sev",c("parameters")] <- "Slope"
m_sus_b0 = m_sus_regression |>
filter(parameters == "Intercept")
m_sus_b1 = m_sus_regression |>
filter(parameters == "Slope") %>%
filter(!estimate >0)
mean(m_sus_b0$estimate)[1] 4198.677sd(m_sus_b0$estimate)[1] 619.4013mean(m_sus_b1$estimate)[1] -32.01062sd(m_sus_b1$estimate)[1] 19.07879Intercept
intercept_tss_m_sus = m_sus_regression |>
filter(parameters == "Intercept")
mean_intercept_m_sus = mean(intercept_tss_m_sus$estimate)
sd_intercept_m_sus = sd(intercept_tss_m_sus$estimate)
plot(ecdf(intercept_tss_m_sus$estimate))
curve(pnorm(x, mean_intercept_m_sus,
sd_intercept_m_sus), add = T, col = "red")
Kolmogorov-Smirnov Test
Fx_m_sus_b0 = environment(ecdf(intercept_tss_m_sus$estimate))$y
x_m_sus_b0= environment(ecdf(intercept_tss_m_sus$estimate))$x
ks.test(Fx_m_sus_b0, pnorm(x_m_sus_b0, mean(intercept_tss_m_sus$estimate), sd(intercept_tss_m_sus$estimate)))
Exact two-sample Kolmogorov-Smirnov test
data: Fx_m_sus_b0 and pnorm(x_m_sus_b0, mean(intercept_tss_m_sus$estimate), sd(intercept_tss_m_sus$estimate))
D = 0.11111, p-value = 0.9974
alternative hypothesis: two-sidedplot(Fx_m_sus_b0, pnorm(x_m_sus_b0, mean(intercept_tss_m_sus$estimate), sd(intercept_tss_m_sus$estimate)))
abline(a=0,b=1)
intercep_plot_m_sus = intercept_tss_m_sus %>%
ggplot(aes(estimate))+
geom_histogram(aes(y = ..density..),bins = 8, color = "white", fill = "#FF4917")+
#stat_function(fun=function(x) dnorm(x, mean_intercept_m_sus, sd_intercept_m_sus), color ="black", size =1.2)+
theme_bw()+
labs(x="Intercept (kg/ha)", y = "Density")
intercep_plot_m_sus
Slope
slope_tss_m_sus = m_sus_regression |>
filter(parameters == "Slope") %>%
filter(!estimate > 0) %>%
filter(!estimate > -10)
mean_slope_m_sus = mean(slope_tss_m_sus$estimate)
sd_slope_m_sus = sd(slope_tss_m_sus$estimate)
plot(ecdf(slope_tss_m_sus$estimate))
curve(pnorm(x, mean_slope_m_sus, sd_slope_m_sus), add = T, col = "red")
Kolmogorov-Smirnov Test
Fx_m_sus_b1 =environment(ecdf(-slope_tss_m_sus$estimate))$y
x_m_sus_b1 = environment(ecdf(-slope_tss_m_sus$estimate))$x
slope_m_sus = nlsLM(Fx_m_sus_b1 ~ pgamma(x_m_sus_b1,
shape, rate,log = FALSE) ,
start = c(shape = 2.5, rate = 0.13),
control = nls.lm.control(maxiter = 1024))
summary(slope_m_sus)
Formula: Fx_m_sus_b1 ~ pgamma(x_m_sus_b1, shape, rate, log = FALSE)
Parameters:
Estimate Std. Error t value Pr(>|t|)
shape 3.959797 0.186074 21.28 3.29e-15 ***
rate 0.117930 0.006016 19.60 1.58e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.02042 on 20 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 1.49e-08shape_m_sus = summary(slope_m_sus)$coef[1]
rate_m_sus = summary(slope_m_sus)$coef[2]
Fx_m_sus_b1 =environment(ecdf(-slope_tss_m_sus$estimate))$y
x_m_sus_b1 = environment(ecdf(-slope_tss_m_sus$estimate))$x
ks.test(Fx_m_sus_b1, pgamma(x_m_sus_b1, shape_m_sus, rate_m_sus))
Exact two-sample Kolmogorov-Smirnov test
data: Fx_m_sus_b1 and pgamma(x_m_sus_b1, shape_m_sus, rate_m_sus)
D = 0.045455, p-value = 1
alternative hypothesis: two-sidedslope_plot_m_sus = slope_tss_m_sus %>%
ggplot(aes(estimate))+
geom_histogram(aes(y = ..density..),bins = 8, color = "white", fill = "#159EE6")+
theme_bw()+
scale_x_continuous(breaks = c(-200,-150,-100,-50,0),
limits = c(-200,0))+
labs(x="Slope (kg/p.p.)", y = "Density")
slope_plot_m_sus
Correlation
slope_coefficient_m_sus = m_sus_regression |>
filter(parameters == "Slope") |>
#filter(estimate < 0) %>%
group_by(study) |>
summarise(
Slope = estimate
)
slope_coefficient_m_sus[,1] = NULL
slope_coefficient_m_sus |>
filter(!Slope == "NA") |>
summarise(
mean = mean(Slope))intercept_coefficient_m_sus = m_sus_regression |>
filter(parameters == "Intercept") |>
group_by(study) |>
summarise(
Intercept = estimate
)
mean(intercept_coefficient_m_sus$Intercept)[1] 4198.677sd(intercept_coefficient_m_sus$Intercept)[1] 619.4013regression_m_sus = cbind(slope_coefficient_m_sus,intercept_coefficient_m_sus)
regression_m_sus = regression_m_sus %>%
filter(Slope <0)
correlation(regression_m_sus$Slope, regression_m_sus$Intercept)[1] -0.7341185corr_actual_plot_m_sus = regression_m_sus %>%
mutate(alfa =(Slope/Intercept)*100) %>%
#filter(alfa > -3 & alfa < 0) %>%
ggplot(aes(Intercept,Slope))+
geom_point( size =2, color = "black")+
# geom_density_2d(color = "black")+
geom_smooth(method = lm, color = "red", se = F, size = 1.2)+
# scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000), limits = c(0, 10000))+
scale_y_continuous(breaks = c(-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-160, 0))+
ggthemes::theme_few()+
labs(y= "Slope",
x = "Intercept",
title = "Empirical: r = -0.734")+
theme(text = element_text(face = "bold"))
# coord_cartesian(
# xlim = c(0,10000))
corr_actual_plot_m_sus
plot_grid(corr_actual_plot_res,corr_actual_plot_m_sus,ncol = 2, labels = c("AUTO"))
j_m_sus<-gera.norm.bid.geral(10000,0.73,0,0,1,1)
plot(j_m_sus[,1],j_m_sus[,2])
Severity
Kolmogorov-Smirnov Test
sev_m_sus = ma_c %>%
filter(cultivar %in% c_01$cultivar) %>%
filter(ai %in% c("ZZ_CHECK","_CHECK"))
sev_m_s = sev_m_sus$sev
Fx_m_sus= environment(ecdf(sev_m_s))$y
x_m_sus = environment(ecdf(sev_m_s))$x/100
summary(nlsLM(Fx_m_sus ~ pbeta(x_m_sus, shape1, shape2, log = FALSE) ,
start = c(shape1 = 1, shape2 = 1),
control = nls.lm.control(maxiter = 100000)))
Formula: Fx_m_sus ~ pbeta(x_m_sus, shape1, shape2, log = FALSE)
Parameters:
Estimate Std. Error t value Pr(>|t|)
shape1 4.4564 0.1497 29.77 <2e-16 ***
shape2 7.8184 0.2690 29.07 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03163 on 97 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 1.49e-08ks.test(Fx_m_sus,pbeta(x_m_sus,4.4564,7.8184))
Exact two-sample Kolmogorov-Smirnov test
data: Fx_m_sus and pbeta(x_m_sus, 4.4564, 7.8184)
D = 0.080808, p-value = 0.9053
alternative hypothesis: two-sidedPlot
plot(x_m_sus,Fx_m_sus)
curve(pbeta(x, 4.4564,7.8184),0,1, add = T)
sev_dist_plot_m_sus = sev_m_sus %>%
ggplot(aes(sev/100))+
geom_histogram(aes(y = ..density..), bins = 8, color = "white", fill = "darkgreen")+
scale_x_continuous(breaks = c(1.00,0.75,0.50,0.25,0.00),
limits = c(0.00,1.00))+
theme_bw()+
labs(x="Severity (proportion)", y = "Density")
sev_dist_plot_m_sus
Relative yield loss
empiric_ryl_m_sus= regression_m_sus %>%
mutate(cc = (Slope/Intercept)*100)
# filter( cc > -3 & cc <0 )
head(empiric_ryl_m_sus)stat_emp_ryl_m_sus =empiric_ryl_m_sus%>%
summarise(data = "empirical",
mean = mean(cc),
median = median(cc),
variance = var(cc))real_RYL_m_sus= empiric_ryl_m_sus %>%
ggplot(aes(cc))+
geom_histogram(aes(y = ..density..), bins = 10, color = "white", fill = "black")+
# scale_y_continuous(limits = c(0,1.4),breaks = seq(0, 1.4,by = 0.4))+
geom_vline(data = stat_emp_ryl_m_sus, aes(xintercept = mean, color = "Mean"),size =1)+
geom_vline(data = stat_emp_ryl_m_sus, aes(xintercept = median, color = "Median"),size =1)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
scale_color_calc()+
labs(x = "",
y = "",
color = "",
title = "Empirical | Medium")+
# xlim(-3,0.2)+
theme(legend.position = "none",
text = element_text(face = "bold"))
real_RYL_m_sus
mean(m_sus_b0$estimate)[1] 4198.677sd(m_sus_b0$estimate)[1] 619.4013b0_m_sus = pnorm(j_m_sus[,2])
b0_t_m_sus = qnorm(b0_m_sus, 4198.677, 619*sqrt(27))
hist(b0_t_m_sus, prob = T)
curve(dnorm(x, mean_intercept_m_sus, sd_intercept_m_sus), add = T)
b1_m_sus = pnorm(j_m_sus[,1])
b1_t_m_sus = qgamma(b1_m_sus, shape = shape_m_sus, rate = rate_m_sus)*-1
hist(b1_t_m_sus, prob = T)
curve(dgamma(-x,shape=shape_m_sus, rate = rate_m_sus), add = T)
correlation(b1_t_m_sus, b0_t_m_sus)[1] -0.7118311cor.test(b1_t_m_sus, b0_t_m_sus, method = "pearson")
Pearson's product-moment correlation
data: b1_t_m_sus and b0_t_m_sus
t = -101.34, df = 9998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.7213667 -0.7020257
sample estimates:
cor
-0.7118311 actual_cc_m_sus = regression_m_sus %>%
mutate(cc = (Slope/Intercept)*100)
simulated_cc_m_sus = data.frame(b0_t_m_sus, b1_t_m_sus,
cc = (b1_t_m_sus/b0_t_m_sus)*100)%>%
filter( cc > -3 & cc <0 )
ks.test(actual_cc_m_sus$cc, simulated_cc_m_sus$cc)
Asymptotic two-sample Kolmogorov-Smirnov test
data: actual_cc_m_sus$cc and simulated_cc_m_sus$cc
D = 0.16419, p-value = 0.5123
alternative hypothesis: two-sidedfx_actual_m_sus = environment(ecdf(actual_cc_m_sus$cc))$y
fx_simu_m_sus = environment(ecdf(simulated_cc_m_sus$cc))$y
ks.test((fx_actual_m_sus),(fx_simu_m_sus))
Asymptotic two-sample Kolmogorov-Smirnov test
data: (fx_actual_m_sus) and (fx_simu_m_sus)
D = 0.039995, p-value = 1
alternative hypothesis: two-sidedecdf_damage_m_sus = ggplot()+
stat_ecdf(aes(simulated_cc_m_sus$cc,color = "Simulated"), size=1.2,geom = "step")+
stat_ecdf(aes(actual_cc_m_sus$cc, color = "Empirical"), size=1.2,geom = "step")+
ggthemes::theme_few()+
scale_color_manual(values = c("black","orange"))+
labs(y = "",
x = "Relative percent yield loss",
color ="",
title = "CDF | Medium")+
theme(legend.position = "top",
legend.background = element_blank(),
text = element_text(face = "bold"))
ecdf_damage_m_sus
corr_sim_plot_m_sus = data.frame(b0_t_m_sus,b1_t_m_sus) %>%
mutate(alfa =(b0_t_m_sus/b1_t_m_sus)*100 ) %>%
#filter(alfa > -3 & alfa < 0) %>%
ggplot(aes(b0_t_m_sus,b1_t_m_sus))+
geom_point( size =2, color = "orange", alpha =0.3)+
# geom_density_2d(color = "black")+
geom_smooth(method = lm, color = "red", se = F, size = 1.2)+
scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000), limits = c(0, 10000))+
scale_y_continuous(breaks = c(-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-160, 0))+
ggthemes::theme_few()+
labs(y= "Slope",
x = "Intercept",
title = "Simulated: r = -0.711")+
coord_cartesian(
xlim = c(0,7000),
ylim = c(-160,0))+
theme(text = element_text(face = "bold"))
corr_sim_plot_m_sus
plot_grid(corr_actual_plot_res,corr_actual_plot_m_sus,corr_sim_plot_res,corr_sim_plot_m_sus, ncol = 2, labels = c("AUTO"))
LOW-TOLERANCE
susceptible = ma_c %>%
filter(cultivar %in% c_02$cultivar)
sus_regression = susceptible%>%
group_by(study,block) %>%
do(tidy(lm(.$yld ~ .$sev)))
sus_regression = as.data.frame(sus_regression)
sus_regressionsus_severity = susceptible |>
# filter(ai == "_CHECK") %>%
ggplot(aes(sev))+
geom_histogram(color = "white",fill = "steelblue", bins = 15)+
scale_x_continuous(limits = c(0, 100),
breaks = seq(0, 1000, by = 10))+
theme_bw()+
labs(x = "Severity (%)",
y = "Density")
sus_yield = susceptible |>
#filter(ai == "_CHECK") %>%
ggplot(aes(yld))+
geom_histogram(color = "white",fill = "brown", bins = 10)+
scale_x_continuous(limits = c(2000, 7000),
breaks = seq(2000, 7000, by = 500))+
theme_bw()+
labs(x = "Yield (kg/ha)",
y = "")
plot_grid(sus_severity,sus_yield)
sus_regression = sus_regression |>
filter(term %in% c("(Intercept)",".$sev"))
sus_regression[sus_regression$term== "(Intercept)",c("parameters")] <- "Intercept"
sus_regression[sus_regression$term== ".$sev",c("parameters")] <- "Slope"
sus_b0 = sus_regression |>
filter(parameters == "Intercept") %>%
filter(!estimate < 3000)
sus_b1 = sus_regression |>
filter(parameters == "Slope") %>%
filter(!estimate >0)
mean(sus_b0$estimate)[1] 4818.456sd(sus_b0$estimate)[1] 426.985mean(sus_b1$estimate)[1] -76.35969sd(sus_b1$estimate)[1] 40.67328Intercept
intercept_tss_sus = sus_regression |>
filter(parameters == "Intercept") %>%
filter(!estimate < 3000)
mean_intercept_sus = mean(intercept_tss_sus$estimate)
sd_intercept_sus = sd(intercept_tss_sus$estimate)
plot(ecdf(intercept_tss_sus$estimate))
curve(pnorm(x, mean_intercept_sus,sd_intercept_sus), add = T, col = "red")
Kolmogorov-Smirnov Test
Fx_sus_b0 = environment(ecdf(intercept_tss_sus$estimate))$y
x_sus_b0= environment(ecdf(intercept_tss_sus$estimate))$x
ks.test(Fx_sus_b0, pnorm(x_sus_b0, mean(intercept_tss_sus$estimate), sd(intercept_tss_sus$estimate)))
Exact two-sample Kolmogorov-Smirnov test
data: Fx_sus_b0 and pnorm(x_sus_b0, mean(intercept_tss_sus$estimate), sd(intercept_tss_sus$estimate))
D = 0.2, p-value = 0.832
alternative hypothesis: two-sidedplot(Fx_sus_b0, pnorm(x_sus_b0, mean(intercept_tss_sus$estimate), sd(intercept_tss_sus$estimate)))
abline(a=0,b=1)
intercep_plot_sus = intercept_tss_sus %>%
ggplot(aes(estimate))+
geom_histogram(aes(y = ..density..),bins = 8, color = "white", fill = "#FF4917")+
# stat_function(fun=function(x) dnorm(x, mean_intercept_sus, sd_intercept_sus), color ="black", size =1.2)+
theme_bw()+
labs(x="Intercept (kg/ha)", y = "Density")
intercep_plot_sus
Slope
slope_tss_sus = sus_regression |>
filter(parameters == "Slope") %>%
filter(!estimate > 0)
mean_slope = mean(slope_tss_sus$estimate)
sd_slope = sd(slope_tss_sus$estimate)
plot(ecdf(slope_tss_sus$estimate))
curve(pnorm(x, mean_slope, sd_slope), add = T, col = "red")
Kolmogorov-Smirnov Test
Fx =environment(ecdf(-slope_tss_sus$estimate))$y
x = environment(ecdf(-slope_tss_sus$estimate))$x
slope_sus = nlsLM(Fx ~ pgamma(x, shape, rate,log = FALSE) ,
start = c(shape = 2.5, rate = 0.13),
control = nls.lm.control(maxiter = 1024))
summary(slope_sus)
Formula: Fx ~ pgamma(x, shape, rate, log = FALSE)
Parameters:
Estimate Std. Error t value Pr(>|t|)
shape 3.595815 0.174872 20.56 7.46e-16 ***
rate 0.048475 0.002561 18.93 4.21e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.02137 on 22 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 1.49e-08shape_sus = summary(slope_sus)$coef[1]
rate_sus = summary(slope_sus)$coef[2]
Fx =environment(ecdf(-slope_tss_sus$estimate))$y
x = environment(ecdf(-slope_tss_sus$estimate))$x
ks.test(Fx, pgamma(x, shape_sus, rate_sus))
Exact two-sample Kolmogorov-Smirnov test
data: Fx and pgamma(x, shape_sus, rate_sus)
D = 0.041667, p-value = 1
alternative hypothesis: two-sidedshape_sus = summary(slope_sus)$coef[1]
rate_sus = summary(slope_sus)$coef[2]
slope_plot_sus = slope_tss_sus %>%
ggplot(aes(estimate))+
geom_histogram(aes(y = ..density..),bins = 8, color = "white", fill = "#159EE6")+
scale_x_continuous(breaks = c(-200,-150,-100,-50,0), limits = c(-200,0))+
theme_bw()+
labs(x="Slope (kg/p.p.)", y = "Density")
slope_plot_sus
Correlation
slope_coefficient_sus = sus_regression |>
filter(parameters == "Slope") |>
#filter(estimate < 0) %>%
group_by(block) |>
summarise(
Slope = estimate
)
slope_coefficient_sus[,1] = NULL
slope_coefficient_sus |>
filter(!Slope == "NA") |>
summarise(
mean = mean(Slope))intercept_coefficient_sus = sus_regression |>
filter(parameters == "Intercept") |>
group_by(block) |>
summarise(
Intercept = estimate
)
mean(intercept_coefficient_sus$Intercept)[1] 4472.231sd(intercept_coefficient_sus$Intercept)[1] 882.7419regression_sus = cbind(slope_coefficient_sus,intercept_coefficient_sus)
regression_sus = regression_sus %>%
filter(Slope <0)
correlation(regression_sus$Slope, regression_sus$Intercept)[1] -0.5292833corr_actual_plot_sus = regression_sus %>%
mutate(alfa =(Slope/Intercept)*100) %>%
#filter(alfa > -3 & alfa < 0) %>%
ggplot(aes(Intercept,Slope))+
geom_point( size =2, color = "black")+
# geom_density_2d(color = "black")+
geom_smooth(method = lm, color = "red", se = F, size = 1.2)+
# scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000), limits = c(0, 10000))+
scale_y_continuous(breaks = c(-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-160, 0))+
ggthemes::theme_few()+
labs(y= "Slope",
x = "Intercept",
title = "Empirical: r = -0.529")+
theme(text = element_text(face = "bold"))
# coord_cartesian(
# xlim = c(0,10000))
corr_actual_plot_sus
j_sus<-gera.norm.bid.geral(10000,0.529,0,0,1,1)
plot(j_sus[,1],j_sus[,2])
Severity
Kolmogorov-Smirnov Test
sev_sus = ma_c %>%
filter(cultivar %in% c_02$cultivar) %>%
filter(ai %in% c("ZZ_CHECK","_CHECK"))
sev_s = sev_sus$sev
Fx_sus= environment(ecdf(sev_s))$y
x_sus = environment(ecdf(sev_s))$x/100
summary(nlsLM(Fx_sus ~ pbeta(x_sus, shape1, shape2, log = FALSE) ,
start = c(shape1 = 1, shape2 = 1),
control = nls.lm.control(maxiter = 100000)))
Formula: Fx_sus ~ pbeta(x_sus, shape1, shape2, log = FALSE)
Parameters:
Estimate Std. Error t value Pr(>|t|)
shape1 1.4193 0.1828 7.766 5.32e-08 ***
shape2 5.9858 0.9731 6.151 2.35e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07092 on 24 degrees of freedom
Number of iterations to convergence: 7
Achieved convergence tolerance: 1.49e-08ks.test(Fx_sus,pbeta(x_sus,1.4193,5.9858))
Exact two-sample Kolmogorov-Smirnov test
data: Fx_sus and pbeta(x_sus, 1.4193, 5.9858)
D = 0.15385, p-value = 0.926
alternative hypothesis: two-sidedPlot
plot(x_sus,Fx_sus)
curve(pbeta(x, 1.4193,5.9858),0,1, add = T)
sev_dist_plot_sus = sev_sus %>%
ggplot(aes(sev/100))+
geom_histogram(aes(y = ..density..), bins = 12, color = "white", fill = "darkgreen")+
scale_x_continuous(breaks = c(1.00,0.75,0.50,0.25,0.00),
limits = c(0.00,1.00))+
theme_bw()+
labs(x="Severity (proportion)", y = "Density")
sev_dist_plot_sus
Relative yield loss
empiric_ryl_sus= regression_sus %>%
mutate(cc = (Slope/Intercept)*100)
# filter( cc > -3 & cc <0 )
head(empiric_ryl_sus)stat_emp_ryl_sus =empiric_ryl_sus%>%
summarise(data = "empirical",
mean = mean(cc),
median = median(cc),
variance = var(cc))real_RYL_sus= empiric_ryl_sus %>%
ggplot(aes(cc))+
geom_histogram(aes(y = ..density..), bins = 10, color = "white", fill = "black")+
# scale_y_continuous(limits = c(0,1.4),breaks = seq(0, 1.4,by = 0.4))+
geom_vline(data = stat_emp_ryl_sus, aes(xintercept = mean, color = "Mean"),size =1)+
geom_vline(data = stat_emp_ryl_sus, aes(xintercept = median, color = "Median"),size =1)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
scale_color_calc()+
labs(x = "",
y = "",
color = "",
title = "Empirical | Low")+
# xlim(-3,0.2)+
theme(legend.position = "none",
text = element_text(face = "bold"))
real_RYL_sus
mean(sus_b0$estimate)[1] 4818.456sd(sus_b0$estimate)[1] 426.985b0_sus = pnorm(j_sus[,2])
b0_t_sus = qnorm(b0_sus, 4818.456, 426*sqrt(24))
hist(b0_t_sus, prob = T)
curve(dnorm(x, mean_intercept_sus, sd_intercept_sus), add = T)
b1_sus = pnorm(j_sus[,1])
b1_t_sus = qgamma(b1_sus, shape = shape_sus, rate = rate_sus)*-1
hist(b1_t_sus, prob = T)
curve(dgamma(-x,shape=shape_sus, rate = rate_sus), add = T)
correlation(b1_t_sus, b0_t_sus)[1] -0.5110051cor.test(b1_t_sus, b0_t_sus, method = "pearson")
Pearson's product-moment correlation
data: b1_t_sus and b0_t_sus
t = -59.442, df = 9998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5253435 -0.4963766
sample estimates:
cor
-0.5110051 actual_cc_sus = regression_sus %>%
mutate(cc = (Slope/Intercept)*100)
simulated_cc_sus = data.frame(b0_t_sus, b1_t_sus, cc = (b1_t_sus/b0_t_sus)*100)%>%
filter( cc > -3 & cc <0 )
ks.test(actual_cc_sus$cc, simulated_cc_sus$cc)
Asymptotic two-sample Kolmogorov-Smirnov test
data: actual_cc_sus$cc and simulated_cc_sus$cc
D = 0.14025, p-value = 0.7341
alternative hypothesis: two-sidedfx_actual_sus = environment(ecdf(actual_cc_sus$cc))$y
fx_simu_sus = environment(ecdf(simulated_cc_sus$cc))$y
ks.test((fx_actual_sus),(fx_simu_sus))
Asymptotic two-sample Kolmogorov-Smirnov test
data: (fx_actual_sus) and (fx_simu_sus)
D = 0.041657, p-value = 1
alternative hypothesis: two-sidedecdf_damage_sus = ggplot()+
stat_ecdf(aes(simulated_cc_sus$cc,color = "Simulated"), size=1.2,geom = "step")+
stat_ecdf(aes(actual_cc_sus$cc, color = "Empirical"), size=1.2,geom = "step")+
ggthemes::theme_few()+
scale_color_manual(values = c("black","orange"))+
labs(y = "",
x = "Relative percent yield loss",
color ="",
title = "CDF | Low")+
theme(legend.position = "top",
legend.background = element_blank(),
text = element_text(face = "bold"))
ecdf_damage_sus
corr_sim_plot_sus = data.frame(b0_t_sus,b1_t_sus) %>%
mutate(alfa =(b1_t_sus/b0_t_sus)*100 ) %>%
#filter(alfa > -3 & alfa < 0) %>%
ggplot(aes(b0_t_sus,b1_t_sus))+
geom_point( size =2, color = "orange", alpha =0.3)+
# geom_density_2d(color = "black")+
geom_smooth(method = lm, color = "red", se = F, size = 1.2)+
scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000), limits = c(0, 10000))+
scale_y_continuous(breaks = c(-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-160, 0))+
ggthemes::theme_few()+
labs(y= "Slope",
x = "Intercept",
title = "Simulated: r = -0.511")+
coord_cartesian(
xlim = c(0,7000),
ylim = c(-160,0))+
theme(text = element_text(face = "bold"))
corr_sim_plot_sus
plot_grid(corr_actual_plot_sus,corr_sim_plot_sus, ncol = 2, labels = c("AUTO"))
Joining
plot_grid(intercep_plot_res,intercep_plot_m_sus,intercep_plot_sus,slope_plot_res,slope_plot_m_sus,slope_plot_sus,sev_dist_plot_res,sev_dist_plot_m_sus,sev_dist_plot_sus, ncol = 3, labels = c("AUTO"))
ggsave("fig/empirical_distribution.png", bg = "white",
dpi = 600, width = 10, height = 8)Soybean price
soybean = gsheet2tbl("https://docs.google.com/spreadsheets/d/1-jQ9OgWdLQCb0iB0FqbrhuVi7LiNhqxvf9QU4-iuc3o/edit#gid=1085329359") sbr_price = soybean %>%
filter(year>=2018) %>%
mutate(price = (price/60)/4,
national_price = (national_price/60)/4)
sbr_pricemean(sbr_price$price)[1] 0.2932792sd(sbr_price$price)[1] 0.02167031Observed distribution
sbr_price %>%
ggplot(aes(price))+
geom_histogram(bins = 10, fill = "steelblue", color = "white")+
ggthemes::theme_few()+
labs(x = "Soybean prince",
y = "Frequency")+
scale_x_continuous(breaks = seq(0,1,by=0.025))+
theme(text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 14, face = "bold"),
legend.position = "none")
hist((sbr_price$price), prob = T)
curve(dnorm(x, mean(sbr_price$price), sd(sbr_price$price)),0.15,0.35, add = T)
plot(ecdf(sbr_price$price))
curve(pnorm(x, mean(sbr_price$price), sd(sbr_price$price)),0.2,0.35, add = T)
mean(sbr_price$price)[1] 0.2932792median(sbr_price$price)[1] 0.2935692Shapiro
shapiro.test(sbr_price$price)
Shapiro-Wilk normality test
data: sbr_price$price
W = 0.99021, p-value = 0.6164Kolmogorov-Smirnov Test
Fx = environment(ecdf(sbr_price$price))$y
x= environment(ecdf(sbr_price$price))$x
ks.test(Fx, pnorm(x, mean(sbr_price$price), sd(sbr_price$price)))
Asymptotic two-sample Kolmogorov-Smirnov test
data: Fx and pnorm(x, mean(sbr_price$price), sd(sbr_price$price))
D = 0.054545, p-value = 0.9967
alternative hypothesis: two-sidedplot(Fx, pnorm(x, mean(sbr_price$price), sd(sbr_price$price)))
Vizualization
price_plot = sbr_price %>%
ggplot(aes(price))+
geom_histogram(aes(y = ..density..),bins = 10, color = "white", fill = "black")+ #"#1C8C20"
#stat_function(fun=function(x) dnorm(x, mean(sbr_price$price), sd(sbr_price$price)), color= "black", size = 1.2)+
ggthemes::theme_few()+
labs(x="Soybean price (USD/kg)", y = "Frequency")+
scale_x_continuous(breaks = seq(0,1,by=0.025))+
theme(text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 14, face = "bold"),
legend.position = "none")
price_plot
# ggsave("fig/soybean_price.png", bg = "white",
# dpi = 600, height = 5, width = 6)Simulated variables
#Soybean
mean_price = mean(sbr_price$price)
sd_price = sd(sbr_price$price)
# High
mean_intercept_res = mean(intercept_coefficient_res$Intercept)
sd_intercept_res = sd(intercept_coefficient_res$Intercept)
slope_coefficient_res = slope_coefficient_res %>%
filter(!Slope < -40) %>%
filter(!Slope > 0)
mean_res = mean(slope_coefficient_res$Slope)
sd_res = sd(slope_coefficient_res$Slope)
s1_sev_res = 3.42200
s2_sev_res = 4.75064
# Medium
mean_intercept_m_sus = mean(intercept_coefficient_m_sus$Intercept)
sd_intercept_m_sus = sd(intercept_coefficient_m_sus$Intercept)
slope_coefficient_m_sus = slope_coefficient_m_sus %>%
filter(!Slope > 0) %>%
filter(!Slope > -10)
mean_m_sus = mean(slope_coefficient_m_sus$Slope)
sd_m_sus = sd(slope_coefficient_m_sus$Slope)
s1_sev_m_sus = 4.4564
s2_sev_m_sus = 7.8184
# Low
mean_intercept_sus = mean(intercept_coefficient_sus$Intercept)
sd_intercept_sus = sd(intercept_coefficient_sus$Intercept)
slope_coefficient_sus = slope_coefficient_sus
mean_sus = mean(slope_coefficient_sus$Slope)
sd_sus = sd(slope_coefficient_sus$Slope)
s1_sev_sus = 1.4193
s2_sev_sus = 5.9858Monte Carlo simulation
set.seed(1)
n=40000
lambda = seq(0,1, by=0.05)
fun_price = seq(-10, 260, by=15)
n_aplication = 1
operational_cost = 10
comb_matrix = as.matrix(data.table::CJ(lambda,fun_price))
colnames(comb_matrix) = c("lambda","fun_price")
comb_matrix = cbind(comb_matrix,operational_cost, n_aplication)
C = comb_matrix[,"n_aplication"]*(comb_matrix[,"operational_cost"]+comb_matrix[,"fun_price"] )
comb_matrix = cbind(comb_matrix,C)
N = length(comb_matrix[,1])*n
big_one = matrix(0, ncol = 24, nrow =N)
big_one[,1] = rep(comb_matrix[,1],n)
big_one[,2] = rep(comb_matrix[,2],n)
big_one[,3] = rep(comb_matrix[,3],n)
big_one[,4] = rep(comb_matrix[,4],n)
big_one[,5] = rep(comb_matrix[,5],n)
set.seed(1)
#Moderate Resistant (MR)
sn_res = rbeta(N, s1_sev_res, s2_sev_res)
sf_res = sn_res*(1-big_one[,1])
# Moderate Susceptible (MS)
sn_m_sus = rbeta(N, s1_sev_m_sus, s2_sev_m_sus)
sf_m_sus = sn_m_sus*(1-big_one[,1])
# Susceptible (S)
sn_sus = rbeta(N, s1_sev_sus, s2_sev_sus)
sf_sus = sn_sus*(1-big_one[,1])
# simulating the coeficientes
## Moderate Resistant (MR)
set.seed(1)
normal_correlated_res<-gera.norm.bid.geral(N,0.51,0,0,1,1)
b0_n_res = pnorm(normal_correlated_res[,2])
b1_n_res = pnorm(normal_correlated_res[,1])
b0_res = qnorm(b0_n_res, mean_intercept_res,sd_intercept_res)
b1_res = -qgamma(b1_n_res, shape_res, rate_res,)
rm(b0_n_res,b1_n_res,normal_correlated_res)
## Moderate Susceptible (MS)
set.seed(1)
normal_correlated_m_sus<-gera.norm.bid.geral(N,0.73,0,0,1,1)
b0_n_m_sus = pnorm(normal_correlated_m_sus[,2])
b1_n_m_sus = pnorm(normal_correlated_m_sus[,1])
b0_m_sus = qnorm(b0_n_m_sus, mean_intercept_m_sus,sd_intercept_m_sus)
b1_m_sus = -qgamma(b1_n_m_sus, shape_m_sus, rate_m_sus,)
rm(b0_n_m_sus,b1_n_m_sus,normal_correlated_m_sus)
## Susceptible (S)
set.seed(1)
normal_correlated_sus<-gera.norm.bid.geral(N,0.52,0,0,1,1)
b0_n_sus = pnorm(normal_correlated_sus[,2])
b1_n_sus = pnorm(normal_correlated_sus[,1])
b0_sus = qnorm(b0_n_sus, mean_intercept_sus,sd_intercept_sus)
b1_sus = -qgamma(b1_n_sus, shape_sus, rate_sus,)
rm(b0_n_sus,b1_n_sus,normal_correlated_sus)
# Calculating the alha coeficient
## Moderate Resistant (MR)
alfa_res = (b1_res/b0_res)*100
## Moderate Susceptible (MS)
alfa_m_sus = (b1_m_sus/b0_m_sus)*100
## Susceptible (S)
alfa_sus = (b1_sus/b0_sus)*100
# Calculating yield gain
## Moderate Resistant (MR)
yn_res = b0_res - (-alfa_res*b0_res*sn_res)
yf_res = b0_res - (-alfa_res*b0_res*sf_res)
## Moderate Susceptible (MS)
yn_m_sus = b0_m_sus - (-alfa_m_sus*b0_m_sus*sn_m_sus)
yf_m_sus = b0_m_sus - (-alfa_m_sus*b0_m_sus*sf_m_sus)
## Susceptible (S)
yn_sus = b0_sus - (-alfa_sus*b0_sus*sn_sus)
yf_sus = b0_sus - (-alfa_sus*b0_sus*sf_sus)
# Simulating soybean price
set.seed(1)
soy_price = rnorm(N, mean(sbr_price$price),sd(sbr_price$price))
#income = yield_gain*soy_price # calculating the income
big_one[,6] = yn_res
big_one[,7] = yf_res
big_one[,8] = yn_m_sus
big_one[,9] = yf_m_sus
big_one[,10] = yn_sus
big_one[,11] = yf_sus
big_one[,12] = soy_price
big_one[,13] = b1_res
big_one[,14] = b1_m_sus
big_one[,15] = b1_sus
big_one[,16] = alfa_res
big_one[,17] = alfa_m_sus
big_one[,18] = alfa_sus
big_one[,19] = b0_res
big_one[,20] = b0_m_sus
big_one[,21] = b0_sus
big_one[,22] = sn_res
big_one[,23] = sn_m_sus
big_one[,24] = sn_sus
colnames(big_one) = c("lambda","fun_price","operational_cost","n_aplication","C","yn_res","yf_res","yn_m_sus","yf_m_sus","yn_sus","yf_sus","soy_price","b1_res","b1_m_sus","b1_sus","alfa_res","alfa_m_sus","alfa_sus","b0_res","b0_m_sus","b0_sus","sn_res","sn_m_sus" ,"sn_sus")big_one_df = as.data.frame(big_one) %>%
filter(b0_res>=0) %>%
filter(b0_sus>=0) %>%
filter(yn_res>0) %>%
filter(yn_sus>0) %>%
filter(alfa_res > -3 & alfa_res < 0) %>%
filter(alfa_sus > -3 & alfa_sus < 0) %>%
mutate(yield_gain_res = yf_res-yn_res,
yield_gain_m_sus = yf_m_sus-yn_m_sus,
yield_gain_sus = yf_sus-yn_sus,
yield_gain_perc_res = ((yf_res/yn_res)-1)*100,
yield_gain_perc_m_sus = ((yf_m_sus/yn_m_sus)-1)*100,
yield_gain_perc_sus = ((yf_sus/yn_sus)-1)*100,
income_res = yield_gain_res*soy_price,
income_m_sus = yield_gain_m_sus*soy_price,
income_sus = yield_gain_sus*soy_price,
CP = C/soy_price,
profit_res = (income_res>=C)*1,
profit_m_sus = (income_m_sus>=C)*1,
profit_sus = (income_sus>=C)*1) %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60)
gc()#write_csv(big_one_df,"data/big_one_df.csv")
big_one_df = read_csv("data/big_one_df.csv")Plotting variables
sn_res_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(sn_res))+
geom_histogram(color = "white", fill = "steelblue", bins = 10, alpha= .5)+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "Frequency",
x = "Target spot severity (%)")
gc()
sn_m_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(sn_m_sus))+
geom_histogram(color = "white", fill = "darkgreen", bins = 10, alpha= .5)+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Target spot severity (%)")
gc()
sn_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(sn_sus))+
geom_histogram(color = "white", fill = "brown", bins = 10, alpha= .5)+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Target spot severity (%)")
gc()b0_res_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(b0_res))+
geom_histogram(color = "white", fill = "steelblue", bins = 10, alpha= .5)+
scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000), limits = c(0, 7000))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "Frequency",
x = "Intercept (kg/ha)")
gc()
b0_m_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(b0_m_sus))+
geom_histogram(color = "white", fill = "darkgreen", bins = 10, alpha= .5)+
scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000), limits = c(0, 7000))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Intercept (kg/ha)")
gc()
b0_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(b0_sus))+
geom_histogram(color = "white", fill = "brown", bins = 10, alpha= .5)+
scale_x_continuous(breaks = c(0,1000,2000,3000,4000,5000,6000,7000), limits = c(0, 7000))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Intercept (kg/ha)")
gc()b1_res_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(b1_res))+
geom_histogram(color = "white", fill = "steelblue", bins = 10, alpha= .5)+
scale_x_continuous(breaks = c(-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-100,0))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "Frequency",
x = "Slope (kg/p.p.)")
gc()
b1_m_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(b1_m_sus))+
geom_histogram(color = "white", fill = "darkgreen", bins = 10, alpha= .5)+
scale_x_continuous(breaks = c(-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0), limits = c(-100,0))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Slope (kg/p.p.)")
gc()
b1_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(b1_sus))+
geom_histogram(color = "white", fill = "brown", bins = 10, alpha= .5)+
scale_x_continuous(breaks = c(-200,-180,-160,
-140,-120, -100,-80,-60,-40,-20,0), limits = c(-200,0))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Slope (kg/p.p.)")
gc()alpha_res_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(alfa_res))+
geom_histogram(color = "white", fill = "black", bins = 10)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "Frequency",
x = "Yield loss (%/p.p.)")
gc()
alpha_m_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(alfa_m_sus))+
geom_histogram(color = "white", fill = "black", bins = 10)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Yield loss (%/p.p.)")
gc()
alpha_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(alfa_sus))+
geom_histogram(color = "white", fill = "black", bins = 10)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
theme(text = element_text(size = 12),
axis.title = element_text(size = 14,face = "bold"))+
labs(y = "",
x = "Yield loss (%/p.p.)")
gc()#alpha_res_graphic,alpha_m_sus_graphic,alpha_sus_graphic,
plot_grid(sn_res_graphic,sn_m_sus_graphic,sn_sus_graphic,b0_res_graphic,b0_m_sus_graphic,b0_sus_graphic,b1_res_graphic,b1_m_sus_graphic,b1_sus_graphic, ncol = 3, labels = c("(a)", "(b)", "(c)", "(d)","(e)","(f)","(g)","(h)","(i)"),label_x = -0.03)
ggsave("fig/simulated_variables.png", dpi = 1000, bg = "white",
height = 8, width = 12)
big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(soy_price))+
geom_histogram(color = "white", fill = "black", bins = 20)+
ggthemes::theme_few()+
labs(x="Soybean price (USD/kg)", y = "Frequency")+
scale_x_continuous(breaks = seq(0,1,by=0.025))+
theme(text = element_text(size = 12, face = "bold"),
axis.title = element_text(size = 14, face = "bold"),
legend.position = "none")
ggsave("fig/soybean_price.png", bg = "white",
dpi = 600, height = 5, width = 6)Plotting RYL
stat_simu_ryl_res = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
summarise(data = "Simulated",
mean = mean(alfa_res),
median = median(alfa_res),
variance = var(alfa_res))
alpha_res_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(alfa_res))+
geom_histogram(color = "white", fill = "black", bins = 10)+
geom_vline(data = stat_simu_ryl_res, aes(xintercept = mean, color = "Mean"),size =1)+
geom_vline(data = stat_simu_ryl_res, aes(xintercept = median, color = "Median"),size =1)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
scale_color_calc()+
labs(x = "Relative yield loss (%/p.p.)",
y = "Frequency",
color = "",
title = "Simulated | High")+
# xlim(-3,0.2)+
theme(legend.position = "none",
text = element_text(face = "bold"))
gc()
stat_simu_ryl_m_sus = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
summarise(data = "Simulated",
mean = mean(alfa_m_sus),
median = median(alfa_m_sus),
variance = var(alfa_m_sus))
alpha_m_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(alfa_m_sus))+
geom_histogram(color = "white", fill = "black", bins = 10)+
geom_vline(data = stat_simu_ryl_m_sus, aes(xintercept = mean, color = "Mean"),size =1)+
geom_vline(data = stat_simu_ryl_m_sus, aes(xintercept = median, color = "Median"),size =1)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
scale_color_calc()+
labs(x = "Relative yield loss (%/p.p.)",
y = "",
color = "",
title = "Simulated | Medium")+
# xlim(-3,0.2)+
theme(legend.position = "none",
text = element_text(face = "bold"))
gc()
stat_simu_ryl_sus = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
summarise(data = "Simulated",
mean = mean(alfa_sus),
median = median(alfa_sus),
variance = var(alfa_sus))
alpha_sus_graphic = big_one_df %>%
filter(C <= 180) %>%
filter(lambda >= 0.4) %>%
filter(lambda <= 0.8) %>%
filter(C >= 60) %>%
ggplot(aes(alfa_sus))+
geom_histogram(color = "white", fill = "black", bins = 10)+
geom_vline(data = stat_simu_ryl_sus, aes(xintercept = mean, color = "Mean"),size =1)+
geom_vline(data = stat_simu_ryl_sus, aes(xintercept = median, color = "Median"),size =1)+
scale_x_continuous(breaks = c(-2.0,-1.75,-1.50,-1.25,-1.00,-0.75,-0.50,-0.25,0.0), limits = c(-2.0,0.0))+
ggthemes::theme_few()+
scale_color_calc()+
labs(x = "Relative yield loss (%/p.p.)",
y = "",
color = "",
title = "Simulated | Low")+
# xlim(-3,0.2)+
theme(legend.position = "none",
text = element_text(face = "bold"))
gc()
plot_grid(corr_actual_plot_res,corr_actual_plot_m_sus,corr_actual_plot_sus,corr_sim_plot_res,corr_sim_plot_m_sus,corr_sim_plot_sus,
real_RYL_res,real_RYL_m_sus,real_RYL_sus,alpha_res_graphic,alpha_m_sus_graphic,alpha_sus_graphic,
ecdf_damage_res,ecdf_damage_m_sus,ecdf_damage_sus,
ncol = 3, labels = c("(a)", "(b)", "(c)", "(d)","(e)","(f)","(g)","(h)","(i)", "(j)","(k)","(l)","(m)","(n)","(o)"), label_x = -0.03,label_size = 24)
ggsave("fig/RYL.png", dpi = 600, height = 20, width = 18)