#NATO_14_23

#Knihovny 

library(tidyverse)
library(car)
library(plm)
library(jtools)
library(sjPlot)
library(Modeler)
library(ggrepel)

#Finální datasety jsou přiloženy a není třeba je tedy znovu vytvářet (stačí je nahrát a pracovat s existujícím kódem pod touto sekcí)

#Tvorba datasetu pro graf 1 (vychází z dat NATO)

NATO_14_23L<-NATO_14_23%>%pivot_longer(cols=3:12, names_to = "Year", values_to = "Value")

NATO_14_23W<-NATO_14_23L%>%pivot_wider(names_from = "Issue", values_from = "Value")

NATO_14_23W<-NATO_14_23W%>%mutate(Equip_tot=MILEX_tot*(Equip_pct/100))

#Dataset pro graf 2 a lineární modely (Sumarizace za dekádu a  přidání ekon proměných)
NATO_14_23_sum<-NATO_14_23W%>%group_by(Country) %>%
  summarise(GDP_pc_avg=sum(GDP_pc)/10, Eq_total = sum(Equip_tot), Milex_pct_avg = sum(MILEX_pct)/10, Eq_pct_avg = sum(Equip_pct)/10)

NATO_14_23_sum<-left_join(NATO_14_23_sum, X230707_def_exp_2023_TABLES_en_5)

NATO_14<-NATO_14_23W%>%filter(Year=="2014")
NATO_14<-NATO_14%>%subset(select = -c(2:4, 6,7))
NATO_14<-NATO_14%>%rename(GDP_pc14=GDP_pc)

#Doplnění datestu pro modely o proměnnou “HDP v r. 2014)
NATO_14_23_sum<-left_join(NATO_14_23_sum, NATO_14)

#Graf 1 Výdaje na Equipment vs. celkové výdaje

NATO_14_23W%>%group_by(Country) %>%
  summarise(
    total = sum(Equip_tot), Mtotal = sum(MILEX_tot) )%>% filter(!Country %in% c("United States"))%>%
  ggplot(aes(fct_reorder(Country, total), total/1000))+
  geom_col(fill="darkblue")+
  geom_col(aes(fct_reorder(Country, total), Mtotal/1000), fill="khaki", alpha=0.1, col="Black")+
  geom_segment(aes(x= 12, xend=12, y=200, yend=25), col="red", arrow=arrow(length = unit(0.03, "npc")))+
  geom_text(aes(label=round(total/Mtotal, digits=2), y=-20))+
  coord_flip()+
  theme_bw()+
  theme_bw(base_size = 14)+
  xlab("")+
  ylab("Výdaje mld. USD (konst. 2015)")

#Graf 2 Výdaje na Equpment na vojáka

NATO_14_23_sum%>%
  ggplot(aes(fct_reorder(Country, Eq_total/MPer23), Eq_total/MPer23))+
  geom_col(fill="darkblue")+
  geom_text(aes(label=round(Eq_total/MPer23, digits=2), y=Eq_total/MPer23+100))+
  coord_flip()+
  theme_bw()+
  theme_bw(base_size = 14)+
  xlab("")+
  ylab("Tisíce USD (konst. 2015)")

#Modely LM (a jejich export)

Mod1<-lm(log10(Eq_total/MPer23)~log10(GDP_pc_avg), data = NATO_14_23_sum)
Mod2<-lm(log10(Eq_total/MPer23)~log10(Milex_pct_avg), data = NATO_14_23_sum)
Mod3<-lm(log10(Eq_total/MPer23)~log10(Milex_pct_avg)+log10(GDP_pc_avg), data = NATO_14_23_sum)
Mod4<-lm(log10(Eq_total/MPer23)~log10(Eq_pct_avg), data = NATO_14_23_sum)
Mod5<-lm(log10(Eq_total/MPer23)~log10(Eq_pct_avg)+log10(GDP_pc_avg), data = NATO_14_23_sum)
MOD5x<-lm(log10(Eq_total/MPer23)~log10(Milex_pct_avg)*log10(Eq_pct_avg)+log10(GDP_pc_avg), data = NATO_14_23_sum)

#Přidání predikce z regresního modelu k datasetu

NATO_14_23_sumP<-add_predictions(data=NATO_14_23_sum, Mod1)

#Graf 3 Graf regrese
NATO_14_23_sumP%>%
  ggplot(aes(log10(GDP_pc_avg), log10(Eq_total/MPer23)))+
  geom_point(aes(col=if_else(Country=="Czechia", "True", "False")), show.legend = F)+
  geom_smooth(method="lm")+
  geom_text_repel(aes(label=Country))+
  theme_bw()+
  scale_colour_manual(values=c("True"="Blue", "False"="black"))

#Export tabulek regresních modelů

export_summs(Mod1, Mod2, Mod3, Mod4, Mod5, robust = "HC3", to.file = "Word", file.name = "VRTab.docx")
export_summs(Mod6, Mod7, Mod8, Mod9, Mod10, robust = "HC3")

Mod6<-lm(log10(Eq_total/MPer23)~log10(GDP_pc14), data = NATO_14_23_sum)
Mod7<-lm(log10(Eq_total/MPer23)~log10(Milex_pct_avg), data = NATO_14_23_sum)
Mod8<-lm(log10(Eq_total/MPer23)~log10(Milex_pct_avg)+log10(GDP_pc14), data = NATO_14_23_sum)
Mod9<-lm(log10(Eq_total/MPer23)~log10(Eq_pct_avg), data = NATO_14_23_sum)
Mod10<-lm(log10(Eq_total/MPer23)~log10(Eq_pct_avg)+log10(GDP_pc14), data = NATO_14_23_sum)

export_summs(Mod6, Mod7, Mod8, Mod9, Mod10, robust = "HC3", to.file = "Word", file.name = "VRTab.docx")

#Graf efekty (není v textu - zde jen jako ověření)
FE1<-effect_plot(Mod5, pred = GDP_pc_avg, interval = T, plot.points = TRUE)
FE2<-effect_plot(Mod5, pred = Eq_pct_avg, interval = T, plot.points = TRUE)
library(cowplot)

plot_grid(FE1, FE2, ncol=2)

#Diagnostika modelu (zde jen finální modely, úpravou čísla modelu lze změnit diagnostikovaný model)
plot_model(Mod5, type="diag")
plot_model(Mod10, type="diag")

plot(Mod5)