library(dplyr)
library(purrr)
library(gsheet)
library(raster)
library(ncdf4)
library(lubridate)
library(readxl)
library(writexl)
library(tidyverse)
library(ggplot2)
library(INLA)
library(caret)
library(inlabru)
library(PresenceAbsence)Target spot
Packages
ma2 <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1j0WZXtJsSMN1MAbmkppnCMnr5d9LnsxV/edit?usp=sharing&ouid=112586075609758894128&rtpof=true&sd=true")ma3 = ma2 %>%
filter(year == 2025) %>%
filter(Source == "Claudia")
ma3 <- ma3 %>%
mutate(
pd_90 = planting_date + 90)
ma3$id <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)Nasapower - Claudia 90
library(janitor)
# IMPORTING DATA FROM NASAPOWER AND CREATING A NEW DATAFRAME WITH THE WEATHER DATA
library(nasapower)
box <- data.frame()
# Loop para extrair dados climáticos para cada estudo
for(i in 1:nrow(ma3)) {
# Obter os dados climáticos usando nasapower
lil_nasa <- get_power(
community = "ag",
temporal_api = "daily",
dates = c(ma3$planting_date[i], ma3$pd_90[i]), # Data de plantio e 140 dias depois
lonlat = c(ma3$longitude[i], ma3$lat[i]), # Longitude e Latitude
pars = c("T2M", "RH2M", "PRECTOTCORR", "T2M_MAX", "T2M_MIN", "T2MDEW") # Parâmetros desejados
) %>%
# Adicionar colunas de identificação ao dataframe
mutate(
id = ma3$id[i],
location = ma3$location[i], # Adiciona a localização
state = ma3$state[i],
mean_sev = ma3$mean_sev[i],# Adiciona o estado
year = format(as.Date(ma3$planting_date[i]), "%Y") # Adiciona o ano a partir da data de plantio
)
# Combina os dados atuais com os anteriores
box <- bind_rows(box, lil_nasa)
}
#write_xlsx(box, "data/weather_nasa.xlsx")
wd_f90 = box
wd_f90$wd_90 <- rep(0:90, times = 19)
wd_f90$mean_sev = wd_f90$mean_sev/100
wd_f90 = wd_f90 %>% #wd_f90
dplyr::filter(wd_90 >= 0 & wd_90 <= 40) %>%
group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M = mean(T2M),
T2M_MAX = mean(T2M_MAX),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev)
)%>%
dplyr::select(-id)Nasapower - Bio 90
ma4 = ma2 %>%
filter(Source == "Bio")
ma4 <- ma4 %>%
dplyr::mutate(
pd_90 = planting_date + 90)
ma4$id <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21)
ma4$mean_sev = ma4$mean_sev/100
library(janitor)
# IMPORTING DATA FROM NASAPOWER AND CREATING A NEW DATAFRAME WITH THE WEATHER DATA
library(nasapower)
box3 <- data.frame()
# Loop para extrair dados climáticos para cada estudo
for(i in 1:nrow(ma4)) {
# Obter os dados climáticos usando nasapower
lil_nasa <- get_power(
community = "ag",
temporal_api = "daily",
dates = c(ma4$planting_date[i], ma4$pd_90[i]), # Data de plantio e 140 dias depois
lonlat = c(ma4$longitude[i], ma4$lat[i]), # Longitude e Latitude
pars = c("T2M", "RH2M", "PRECTOTCORR", "T2M_MAX", "T2M_MIN", "T2MDEW") # Parâmetros desejados
) %>%
# Adicionar colunas de identificação ao dataframe
mutate(
id = ma4$id[i],
location = ma4$location[i], # Adiciona a localização
state = ma4$state[i], # Adiciona o estado
mean_sev = ma4$mean_sev[i],
year = format(as.Date(ma4$planting_date[i]), "%Y") # Adiciona o ano a partir da data de plantio
)
# Combina os dados atuais com os anteriores
box3 <- bind_rows(box3, lil_nasa)
}
#write_xlsx(box, "data/weather_nasa.xlsx")
wd_f90_bio = box3
wd_f90_bio$wd_90 <- rep(0:90, times = 21)
wd_f90_bio = wd_f90_bio %>% #wd_f90
dplyr::filter(wd_90 >= 0 & wd_90 <= 40) %>%
group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M = mean(T2M),
T2M_MAX = mean(T2M_MAX),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev)
)%>%
dplyr::select(-id)#wd_b90_bio = rbind(wd_f90_bio,wd_f90)
#write_xlsx(wd_b90_bio,"data/wd_f90.xlsx")
wd_b90_bio = read_xlsx("data/wd_f90.xlsx")Nasapower - All epidemics
ma2 = ma2 %>%
dplyr::mutate(
pd_90 = planting_date + 90)
ma2$id = 1:254
# IMPORTING DATA FROM NASAPOWER AND CREATING A NEW DATAFRAME WITH THE WEATHER DATA
library(nasapower)
box <- data.frame()
# Loop para extrair dados climáticos para cada estudo
for(i in 1:nrow(ma2)) {
# Obter os dados climáticos usando nasapower
lil_nasa <- get_power(
community = "ag",
temporal_api = "daily",
dates = c(ma2$planting_date[i], ma2$pd_90[i]), # Data de plantio e 140 dias depois
lonlat = c(ma2$longitude[i], ma2$lat[i]), # Longitude e Latitude
pars = c("T2M", "RH2M", "PRECTOTCORR", "T2M_MAX", "T2M_MIN", "T2MDEW") # Parâmetros desejados
) %>%
# Adicionar colunas de identificação ao dataframe
mutate(
id = ma2$id[i],
location = ma2$location[i], # Adiciona a localização
state = ma2$state[i], # Adiciona o estado
mean_sev = ma2$mean_sev[i],
planting_date = ma2$planting_date[i],
year = format(as.Date(ma2$planting_date[i]), "%Y") # Adiciona o ano a partir da data de plantio
)
# Combina os dados atuais com os anteriores
box <- bind_rows(box, lil_nasa)
}#write_xlsx(box,data/weather_total.xlsx")
box = read_xlsx("data/weather_total.xlsx")box <- box %>%
mutate(
# Criando a data permitida diretamente com base no estado
allowed_date = case_when(
state == "MT" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-09-16")),
state == "MS" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-09-16")),
state == "GO" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-09-25")),
state == "PR" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-09-11")),
state == "TO" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-10-01")),
state == "DF" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-10-01")),
state == "BA" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-10-01")),
state == "MG" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-10-01")),
state == "SP" ~ as.Date(paste0(format(YYYYMMDD, "%Y"), "-9-16"))
),
# Ajustando o ano da allowed_date para casos de janeiro e fevereiro
allowed_date = if_else(
format(YYYYMMDD, "%m") %in% c("01", "02"),
as.Date(paste0(as.numeric(format(YYYYMMDD, "%Y")) - 1, "-", format(allowed_date, "%m-%d"))),
allowed_date
),
# Calculando a diferença em dias
days_difference = as.numeric(YYYYMMDD - allowed_date)
)Inference
0 to 20
#box$mean_sev = box$mean_sev/100
wd_0_20 = box %>%
dplyr::filter(days >=0 & days <= 20)
wd_0_20 = wd_0_20 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
#ma3 <- head(ma, nrow(ma2) - 37)
#wd_0_20$year = ma2$year # ESTAVA RODANDO ESSE NO INICIO
wd_0_20 = wd_0_20 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
corr <-cor.test(x = wd_0_20$RH2M,y = wd_0_20$PRECTOTCORR, method = "spearman")
corr
Spearman's rank correlation rho
data: wd_0_20$RH2M and wd_0_20$PRECTOTCORR
S = 746631, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.7266223 #wd_0_20 <- head(wd_0_20, nrow(wd_0_20) - 4)
inla_0_20 = inla(epidemic ~T2M + RH2M+PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_0_20,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_0_20)Time used:
Pre = 0.23, Running = 0.151, Post = 0.0563, Total = 0.438
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -10.862 3.973 -18.694 -10.849 -3.107 -10.849 0
T2M 0.264 0.101 0.067 0.264 0.465 0.264 0
RH2M 0.073 0.029 0.016 0.073 0.131 0.073 0
PRECTOTCORR -0.001 0.004 -0.010 -0.001 0.007 -0.001 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 2.26 2.89 0.40 1.50 8.36 0.925
Deviance Information Criterion (DIC) ...............: 277.03
Deviance Information Criterion (DIC, saturated) ....: 276.67
Effective number of parameters .....................: 12.49
Watanabe-Akaike information criterion (WAIC) ...: 277.15
Effective number of parameters .................: 11.91
Marginal log-Likelihood: -166.23
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')20 to 40
wd_20_40 = box %>%
dplyr::filter(days >= 20 & days <= 40)
wd_20_40 = wd_20_40 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_20_40 = wd_20_40 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_20_40 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_20_40,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_20_40)Time used:
Pre = 0.143, Running = 0.147, Post = 0.0406, Total = 0.331
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -9.088 4.552 -18.023 -9.086 -0.162 -9.086 0
T2M 0.185 0.121 -0.052 0.185 0.424 0.185 0
RH2M 0.081 0.030 0.023 0.081 0.139 0.081 0
PRECTOTCORR -0.005 0.004 -0.014 -0.005 0.003 -0.005 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 3.44 5.87 0.462 1.87 14.58 1.06
Deviance Information Criterion (DIC) ...............: 281.75
Deviance Information Criterion (DIC, saturated) ....: 281.39
Effective number of parameters .....................: 11.97
Watanabe-Akaike information criterion (WAIC) ...: 282.25
Effective number of parameters .................: 11.78
Marginal log-Likelihood: -167.90
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')40 to 60
wd_40_60 = box %>%
dplyr::filter(days >= 40 & days <= 60)
wd_40_60 = wd_40_60 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_40_60 = wd_40_60 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_40_60 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_40_60,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_40_60)Time used:
Pre = 0.135, Running = 0.141, Post = 0.0466, Total = 0.323
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -9.123 5.051 -19.001 -9.133 0.813 -9.133 0
T2M 0.223 0.133 -0.038 0.223 0.484 0.223 0
RH2M 0.052 0.031 -0.008 0.052 0.112 0.052 0
PRECTOTCORR 0.004 0.004 -0.004 0.004 0.012 0.004 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 2.97 4.81 0.423 1.69 12.26 0.97
Deviance Information Criterion (DIC) ...............: 277.39
Deviance Information Criterion (DIC, saturated) ....: 277.02
Effective number of parameters .....................: 12.25
Watanabe-Akaike information criterion (WAIC) ...: 277.37
Effective number of parameters .................: 11.56
Marginal log-Likelihood: -165.72
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')60 to 80
wd_60_80 = box %>%
dplyr::filter(days >= 60 & days <= 80)
wd_60_80 = wd_60_80 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_60_80 = wd_60_80 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_60_80 = inla(epidemic ~T2M + RH2M+PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_60_80,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_60_80)Time used:
Pre = 0.129, Running = 0.148, Post = 0.0351, Total = 0.312
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -4.664 4.916 -14.289 -4.670 4.995 -4.670 0
T2M 0.056 0.125 -0.189 0.056 0.301 0.056 0
RH2M 0.062 0.033 -0.003 0.062 0.126 0.062 0
PRECTOTCORR -0.005 0.003 -0.011 -0.005 0.002 -0.005 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 3.10 4.36 0.508 1.93 12.00 1.17
Deviance Information Criterion (DIC) ...............: 285.68
Deviance Information Criterion (DIC, saturated) ....: 285.31
Effective number of parameters .....................: 12.10
Watanabe-Akaike information criterion (WAIC) ...: 286.03
Effective number of parameters .................: 11.77
Marginal log-Likelihood: -169.94
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')0 to 80
wd_0_80 = box %>%
dplyr::filter(days >= 0 & days <= 80)
wd_0_80 = wd_0_80 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_0_80 = wd_0_80 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_0_80 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_0_80,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_0_80)Time used:
Pre = 0.13, Running = 0.141, Post = 0.038, Total = 0.309
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -11.249 4.927 -20.929 -11.244 -1.600 -11.244 0
T2M 0.229 0.125 -0.016 0.229 0.475 0.229 0
RH2M 0.089 0.039 0.012 0.089 0.165 0.089 0
PRECTOTCORR -0.001 0.002 -0.004 -0.001 0.003 -0.001 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 2.76 4.06 0.433 1.67 10.89 0.997
Deviance Information Criterion (DIC) ...............: 279.06
Deviance Information Criterion (DIC, saturated) ....: 278.69
Effective number of parameters .....................: 12.13
Watanabe-Akaike information criterion (WAIC) ...: 279.36
Effective number of parameters .................: 11.76
Marginal log-Likelihood: -167.37
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')0 to 60
wd_0_60 = box %>%
dplyr::filter(days >= 0 & days <= 60)
wd_0_60 = wd_0_60 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_0_60 = wd_0_60 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_0_60 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_0_60,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_0_60)Time used:
Pre = 0.133, Running = 0.152, Post = 0.0351, Total = 0.319
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -11.345 4.738 -20.656 -11.339 -2.069 -11.339 0
T2M 0.268 0.123 0.029 0.267 0.510 0.267 0
RH2M 0.068 0.036 -0.002 0.068 0.138 0.068 0
PRECTOTCORR 0.001 0.002 -0.003 0.001 0.006 0.001 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 2.74 4.04 0.431 1.66 10.84 0.989
Deviance Information Criterion (DIC) ...............: 277.04
Deviance Information Criterion (DIC, saturated) ....: 276.67
Effective number of parameters .....................: 12.11
Watanabe-Akaike information criterion (WAIC) ...: 277.35
Effective number of parameters .................: 11.75
Marginal log-Likelihood: -166.25
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')0 to 40
wd_0_40 = box %>%
dplyr::filter(days >= 0 & days <= 40)
wd_0_40 = wd_0_40 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_0_40 = wd_0_40 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_0_40 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_0_40,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_0_40)Time used:
Pre = 0.129, Running = 0.137, Post = 0.0623, Total = 0.328
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -11.499 4.418 -20.194 -11.489 -2.862 -11.489 0
T2M 0.256 0.114 0.034 0.256 0.481 0.256 0
RH2M 0.088 0.033 0.022 0.088 0.154 0.088 0
PRECTOTCORR -0.002 0.003 -0.008 -0.002 0.004 -0.002 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 2.54 3.56 0.418 1.60 9.82 0.965
Deviance Information Criterion (DIC) ...............: 277.98
Deviance Information Criterion (DIC, saturated) ....: 277.62
Effective number of parameters .....................: 12.22
Watanabe-Akaike information criterion (WAIC) ...: 278.27
Effective number of parameters .................: 11.82
Marginal log-Likelihood: -166.80
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')20 to 60
wd_20_60 = box %>%
dplyr::filter(days >= 20 & days <= 60)
wd_20_60 = wd_20_60 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_20_60 = wd_20_60 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_20_60 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_20_60,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_20_60)Time used:
Pre = 0.131, Running = 0.142, Post = 0.0365, Total = 0.31
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -10.064 4.909 -19.694 -10.064 -0.435 -10.064 0
T2M 0.228 0.129 -0.025 0.228 0.482 0.228 0
RH2M 0.067 0.032 0.003 0.067 0.130 0.067 0
PRECTOTCORR 0.001 0.003 -0.004 0.001 0.006 0.001 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 4.35 8.33 0.497 2.08 19.30 1.13
Deviance Information Criterion (DIC) ...............: 281.28
Deviance Information Criterion (DIC, saturated) ....: 280.92
Effective number of parameters .....................: 11.73
Watanabe-Akaike information criterion (WAIC) ...: 281.71
Effective number of parameters .................: 11.51
Marginal log-Likelihood: -167.62
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')20 to 80
wd_20_80 = box %>%
dplyr::filter(days >= 20 & days <= 80)
wd_20_80 = wd_20_80 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_20_80 = wd_20_80 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_20_80 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_20_80,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_20_80)Time used:
Pre = 0.149, Running = 0.148, Post = 0.0412, Total = 0.338
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -9.514 5.077 -19.470 -9.515 0.446 -9.515 0
T2M 0.176 0.131 -0.080 0.176 0.432 0.176 0
RH2M 0.084 0.037 0.012 0.084 0.156 0.084 0
PRECTOTCORR -0.001 0.002 -0.005 -0.001 0.003 -0.001 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 3.58 6.21 0.469 1.91 15.26 1.07
Deviance Information Criterion (DIC) ...............: 282.39
Deviance Information Criterion (DIC, saturated) ....: 282.02
Effective number of parameters .....................: 11.92
Watanabe-Akaike information criterion (WAIC) ...: 282.80
Effective number of parameters .................: 11.67
Marginal log-Likelihood: -168.54
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')40 to 80
wd_40_80 = box %>%
dplyr::filter(days >= 40 & days <= 80)
wd_40_80 = wd_40_80 %>%
dplyr:: group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
year = first(YEAR)
) %>%
dplyr::select(-id)
wd_40_80 = wd_40_80 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
inla_40_80 = inla(epidemic ~T2M + RH2M+ PRECTOTCORR+
f(year, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
#f(days, model = "iid", hyper = list(prec = list(param = c(0.001,0.001)))),
data = wd_40_80,
family = "binomial",
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE))
summary(inla_40_80)Time used:
Pre = 0.138, Running = 0.14, Post = 0.0651, Total = 0.343
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -7.684 5.198 -17.860 -7.691 2.531 -7.691 0
T2M 0.142 0.134 -0.121 0.142 0.404 0.142 0
RH2M 0.069 0.036 -0.002 0.069 0.140 0.069 0
PRECTOTCORR -0.001 0.003 -0.006 -0.001 0.004 -0.001 0
Random effects:
Name Model
year IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for year 2.27 3.01 0.389 1.48 8.53 0.90
Deviance Information Criterion (DIC) ...............: 280.21
Deviance Information Criterion (DIC, saturated) ....: 279.85
Effective number of parameters .....................: 12.65
Watanabe-Akaike information criterion (WAIC) ...: 280.27
Effective number of parameters .................: 11.99
Marginal log-Likelihood: -167.83
CPO, PIT is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')PREDICTION INLAbru
box$days <- rep(0:90, times = 254)
wd_60_8022 = box %>%
dplyr::filter(days >= 0 & days <= 40)
wd_60_8022 = wd_60_8022 %>%
group_by(id) %>%
dplyr::summarise(
RH2M = mean(RH2M),
T2M_MAX = mean(T2M_MAX),
T2M = mean(T2M),
T2M_MIN = mean(T2M_MIN),
PRECTOTCORR = sum(PRECTOTCORR),
mean_sev = mean(mean_sev),
) %>%
dplyr::select(-id)
wd_60_8022 = wd_60_8022 %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))wd_0_203 = wd_60_8022 %>%
dplyr::select(-epidemic)
inla_0_20 = bru(epidemic ~T2M + RH2M,
data = wd_60_8022,
family = "binomial")
summary(inla_0_20)inlabru version: 2.12.0
INLA version: 24.12.11
Components:
T2M: main = linear(T2M), group = exchangeable(1L), replicate = iid(1L), NULL
RH2M: main = linear(RH2M), group = exchangeable(1L), replicate = iid(1L), NULL
Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
Likelihoods:
Family: 'binomial'
Tag: ''
Data class: 'tbl_df', 'tbl', 'data.frame'
Response class: 'numeric'
Predictor: epidemic ~ .
Used components: effects[T2M, RH2M, Intercept], latent[]
Time used:
Pre = 0.226, Running = 0.138, Post = 0.0222, Total = 0.386
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
T2M 0.207 0.100 0.012 0.207 0.403 0.207 0
RH2M 0.072 0.021 0.030 0.072 0.114 0.072 0
Intercept -9.616 3.836 -17.134 -9.616 -2.098 -9.616 0
Deviance Information Criterion (DIC) ...............: 289.23
Deviance Information Criterion (DIC, saturated) ....: 288.87
Effective number of parameters .....................: 2.96
Watanabe-Akaike information criterion (WAIC) ...: 289.43
Effective number of parameters .................: 3.10
Marginal log-Likelihood: -160.46
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')wd_60_8022 <- head(wd_60_8022, nrow(wd_60_8022) - 40)
wd_60_8022 %>%
filter(epidemic == 0)wd_0_203 = wd_60_8022 %>%
dplyr::select(-epidemic)
inla_0_20 = bru(epidemic ~T2M + RH2M,
data = wd_60_8022,
family = "binomial")wd_f90 = wd_b90_bio %>%
dplyr::mutate(epidemic = ifelse(mean_sev >= 0.25, 1, 0))
wd_0_203 = wd_f90 %>%
dplyr::select(-epidemic)#set.seed(1234)
t = predict(object= inla_0_20, newdata = wd_0_203, formula = ~ c(pred = Intercept + T2M + RH2M), n.samples = 1000)
#write_xlsx(t,"data/t.xlsx")t = read_xlsx("data/t.xlsx")
wd_f90$mean = plogis(t$mean)
wd_60_8 = wd_f90 %>%
dplyr::select(epidemic,mean) %>%
mutate(
ID = nrow(epidemic)
)
wd_60_8 = as.data.frame(wd_60_8)
dat_lasso <- data.frame(1, wd_f90$epidemic, plogis(t$mean))
optimal.thresholds(dat_lasso)auc.roc.plot(dat_lasso)
cm_rf_052 = confusionMatrix(data = as.factor(as.numeric(wd_f90$mean > 0.78)), mode= "everything", reference = as.factor(wd_f90$epidemic))
cm_rf_052Confusion Matrix and Statistics
Reference
Prediction 0 1
0 13 5
1 3 19
Accuracy : 0.8
95% CI : (0.6435, 0.9095)
No Information Rate : 0.6
P-Value [Acc > NIR] : 0.006065
Kappa : 0.5918
Mcnemar's Test P-Value : 0.723674
Sensitivity : 0.8125
Specificity : 0.7917
Pos Pred Value : 0.7222
Neg Pred Value : 0.8636
Precision : 0.7222
Recall : 0.8125
F1 : 0.7647
Prevalence : 0.4000
Detection Rate : 0.3250
Detection Prevalence : 0.4500
Balanced Accuracy : 0.8021
'Positive' Class : 0