library(tidyr)
library(dplyr)
library(ggplot2)
library(scales)
# install.packages("patchwork")
# install.packages("cowplot")
library(grid)
library(patchwork)
library(cowplot)


# ======================================================
# 1. Sector information: unchanged from the mid-adoption code
# ======================================================

Sector_list <- c(
  "Cereals", "Fruit & Veg", "Fish & Seafood", "Meat & Meat products",
  "Veg Oil & Fats", "Dairy", "Sugar & Sugar products",
  "Bev & Tobacco", "Processed food"
)

sector_size <- data.frame(
  Sector = Sector_list,
  Share = c(
    0.007422, 0.008407, 0.002751, 0.008285, 0.005419,
    0.004928, 0.002099, 0.008368, 0.021580
  )
)

sector_order_key <- sector_size %>%
  arrange(Share) %>%
  pull(Sector)

sector_colors <- setNames(
  scales::hue_pal()(nrow(sector_size)),
  sector_size$Sector
)

make_2030_df <- function(values) {
  stopifnot(length(values) == length(Sector_list))
  data.frame(
    Sector = Sector_list,
    check.names = FALSE,
    `2030` = values
  )
}

# ======================================================
# 2. LOW-ADOPTION PATHWAY: 2030 values
# ======================================================

# ------------------------------------------------------
# 2.1 Current clinical coverage: U.S.-only shock
# ------------------------------------------------------

Low_Current_qo_us <- make_2030_df(c(
  -1.67, -2.05, -1.97, -3.27, -1.31,
  -3.60, -4.50, -4.42, -3.68
))

Low_Current_po_us <- make_2030_df(c(
  -0.28, -0.50, -0.97, -0.22, -0.26,
  -0.13, -0.16, -0.12, -0.16
))

Low_Current_qmw_us <- make_2030_df(c(
  -3.66, -3.96, -3.30, -4.64, -3.47,
  -4.60, -4.85, -4.63, -4.32
))

Low_Current_qxw_us <- make_2030_df(c(
  1.17, 1.22, 1.70, 1.74, 1.30,
  1.04, 0.85, 0.29, 0.53
))

# ------------------------------------------------------
# 2.2 Current clinical coverage: North America-wide shock
# U.S. results from the Scenario 2 qo/po/qmw/qxw*US columns
# ------------------------------------------------------

Low_Current_qo_na <- make_2030_df(c(
  -1.72, -2.14, -2.03, -3.30, -1.35,
  -3.62, -4.52, -4.45, -3.73
))

Low_Current_po_na <- make_2030_df(c(
  -0.28, -0.52, -0.99, -0.23, -0.27,
  -0.13, -0.16, -0.12, -0.17
))

Low_Current_qmw_na <- make_2030_df(c(
  -3.65, -3.96, -3.28, -4.59, -3.49,
  -4.61, -4.81, -4.63, -4.32
))

Low_Current_qxw_na <- make_2030_df(c(
  1.04, 0.94, 1.43, 1.49, 1.24,
  0.76, 0.45, -0.09, 0.10
))

# ------------------------------------------------------
# 2.3 Expanded weight-management access: U.S.-only shock
# ------------------------------------------------------

Low_Expanded_qo_us <- make_2030_df(c(
  -2.97, -3.67, -3.54, -5.83, -2.31,
  -6.42, -8.02, -7.88, -6.56
))

Low_Expanded_po_us <- make_2030_df(c(
  -0.51, -0.90, -1.74, -0.41, -0.47,
  -0.24, -0.29, -0.22, -0.30
))

Low_Expanded_qmw_us <- make_2030_df(c(
  -6.48, -7.06, -5.87, -8.22, -6.14,
  -8.17, -8.61, -8.24, -7.68
))

Low_Expanded_qxw_us <- make_2030_df(c(
  2.09, 2.19, 3.07, 3.15, 2.33,
  1.89, 1.53, 0.52, 0.96
))

# ------------------------------------------------------
# 2.4 Expanded weight-management access: North America-wide shock
# U.S. results from the Scenario 2 qo/po/qmw/qxw*US columns
# ------------------------------------------------------

Low_Expanded_qo_na <- make_2030_df(c(
  -3.05, -3.81, -3.62, -5.87, -2.37,
  -6.44, -8.04, -7.91, -6.62
))

Low_Expanded_po_na <- make_2030_df(c(
  -0.52, -0.93, -1.77, -0.41, -0.48,
  -0.24, -0.30, -0.23, -0.31
))

Low_Expanded_qmw_na <- make_2030_df(c(
  -6.48, -7.06, -5.84, -8.15, -6.17,
  -8.18, -8.56, -8.24, -7.68
))

Low_Expanded_qxw_na <- make_2030_df(c(
  1.89, 1.80, 2.69, 2.77, 2.24,
  1.46, 0.95, -0.02, 0.34
))

# ======================================================
# 3. HIGH-ADOPTION PATHWAY: 2030 values
# ======================================================

# ------------------------------------------------------
# 3.1 Current clinical coverage: U.S.-only shock
# ------------------------------------------------------

High_Current_qo_us <- make_2030_df(c(
  -1.66, -2.05, -1.97, -3.28, -1.30,
  -3.60, -4.50, -4.42, -3.68
))

High_Current_po_us <- make_2030_df(c(
  -0.28, -0.50, -0.97, -0.22, -0.26,
  -0.13, -0.16, -0.12, -0.16
))

High_Current_qmw_us <- make_2030_df(c(
  -3.65, -3.96, -3.30, -4.65, -3.46,
  -4.60, -4.84, -4.62, -4.31
))

High_Current_qxw_us <- make_2030_df(c(
  1.17, 1.22, 1.70, 1.74, 1.30,
  1.04, 0.85, 0.29, 0.53
))

# ------------------------------------------------------
# 3.2 Current clinical coverage: North America-wide shock
# U.S. results from the Scenario 2 qo/po/qmw/qxw*US columns
# ------------------------------------------------------

High_Current_qo_na <- make_2030_df(c(
  -1.73, -2.15, -2.05, -3.31, -1.35,
  -3.61, -4.51, -4.44, -3.73
))

High_Current_po_na <- make_2030_df(c(
  -0.28, -0.52, -1.00, -0.23, -0.27,
  -0.13, -0.16, -0.12, -0.17
))

High_Current_qmw_na <- make_2030_df(c(
  -3.65, -3.96, -3.27, -4.59, -3.49,
  -4.60, -4.80, -4.62, -4.31
))

High_Current_qxw_na <- make_2030_df(c(
  1.02, 0.91, 1.34, 1.45, 1.23,
  0.71, 0.40, -0.14, 0.04
))

# ------------------------------------------------------
# 3.3 Expanded weight-management access: U.S.-only shock
# ------------------------------------------------------

High_Expanded_qo_us <- make_2030_df(c(
  -2.97, -3.67, -3.54, -5.84, -2.32,
  -6.42, -8.03, -7.88, -6.57
))

High_Expanded_po_us <- make_2030_df(c(
  -0.51, -0.91, -1.74, -0.41, -0.47,
  -0.24, -0.29, -0.22, -0.30
))

High_Expanded_qmw_us <- make_2030_df(c(
  -6.49, -7.06, -5.87, -8.24, -6.14,
  -8.18, -8.62, -8.25, -7.69
))

High_Expanded_qxw_us <- make_2030_df(c(
  2.09, 2.19, 3.07, 3.15, 2.33,
  1.90, 1.53, 0.52, 0.97
))

# ------------------------------------------------------
# 3.4 Expanded weight-management access: North America-wide shock
# U.S. results from the Scenario 2 qo/po/qmw/qxw*US columns
# ------------------------------------------------------

High_Expanded_qo_na <- make_2030_df(c(
  -3.05, -3.81, -3.63, -5.88, -2.38,
  -6.45, -8.05, -7.92, -6.63
))

High_Expanded_po_na <- make_2030_df(c(
  -0.52, -0.93, -1.78, -0.41, -0.48,
  -0.24, -0.30, -0.23, -0.31
))

High_Expanded_qmw_na <- make_2030_df(c(
  -6.48, -7.06, -5.85, -8.17, -6.17,
  -8.19, -8.56, -8.25, -7.69
))

High_Expanded_qxw_na <- make_2030_df(c(
  1.88, 1.79, 2.68, 2.76, 2.24,
  1.45, 0.94, -0.03, 0.33
))

# ======================================================
# 4. Scenario definitions and combine panel data
# ======================================================
scenario_levels <- c("Current: U.S. only","Current: North America","Expanded: U.S. only","Expanded: North America")
scenario_shapes <- c("Current: U.S. only"=21,"Current: North America"=21,"Expanded: U.S. only"=23,"Expanded: North America"=23)
scenario_offsets <- c("Current: U.S. only"=-0.24,"Current: North America"=-0.08,"Expanded: U.S. only"=0.08,"Expanded: North America"=0.24)
make_panel_df <- function(current_us,current_na,expanded_us,expanded_na) {
  current_us %>%
    rename(Current_US=`2030`) %>%
    left_join(current_na %>% rename(Current_NA=`2030`),by="Sector") %>%
    left_join(expanded_us %>% rename(Expanded_US=`2030`),by="Sector") %>%
    left_join(expanded_na %>% rename(Expanded_NA=`2030`),by="Sector") %>%
    left_join(sector_size,by="Sector") %>%
    mutate(Sector=factor(Sector,levels=sector_order_key),x_id=as.numeric(Sector)) %>%
    pivot_longer(cols=c(Current_US,Current_NA,Expanded_US,Expanded_NA),names_to="Scenario_code",values_to="Impact") %>%
    mutate(
      Scenario=case_when(
        Scenario_code=="Current_US"~"Current: U.S. only",
        Scenario_code=="Current_NA"~"Current: North America",
        Scenario_code=="Expanded_US"~"Expanded: U.S. only",
        Scenario_code=="Expanded_NA"~"Expanded: North America"
      ),
      Scenario=factor(Scenario,levels=scenario_levels),
      x_plot=x_id+unname(scenario_offsets[as.character(Scenario)]),
      PointFill=case_when(
        Scenario_code %in% c("Current_US","Expanded_US")~as.character(Sector),
        Scenario_code %in% c("Current_NA","Expanded_NA")~"Hollow"
      )
    )
}
Low_Po_rng <- make_panel_df(Low_Current_po_us,Low_Current_po_na,Low_Expanded_po_us,Low_Expanded_po_na)
Low_Qo_rng <- make_panel_df(Low_Current_qo_us,Low_Current_qo_na,Low_Expanded_qo_us,Low_Expanded_qo_na)
Low_Qmw_rng <- make_panel_df(Low_Current_qmw_us,Low_Current_qmw_na,Low_Expanded_qmw_us,Low_Expanded_qmw_na)
Low_Qxw_rng <- make_panel_df(Low_Current_qxw_us,Low_Current_qxw_na,Low_Expanded_qxw_us,Low_Expanded_qxw_na)
High_Po_rng <- make_panel_df(High_Current_po_us,High_Current_po_na,High_Expanded_po_us,High_Expanded_po_na)
High_Qo_rng <- make_panel_df(High_Current_qo_us,High_Current_qo_na,High_Expanded_qo_us,High_Expanded_qo_na)
High_Qmw_rng <- make_panel_df(High_Current_qmw_us,High_Current_qmw_na,High_Expanded_qmw_us,High_Expanded_qmw_na)
High_Qxw_rng <- make_panel_df(High_Current_qxw_us,High_Current_qxw_na,High_Expanded_qxw_us,High_Expanded_qxw_na)
# ======================================================
# 5. Plot function
# ======================================================
make_plot_fourtype <- function(df_rng,y_min,y_max) {
  connector_df <- df_rng %>%
    select(Sector,x_id,Scenario_code,Impact) %>%
    pivot_wider(names_from=Scenario_code,values_from=Impact)
  ggplot() +
    geom_hline(yintercept=0,linewidth=0.9,color="black") +
    geom_segment(data=connector_df,aes(x=x_id-0.24,xend=x_id-0.08,y=Current_US,yend=Current_NA),color="grey50",linewidth=0.6,alpha=0.80) +
    geom_segment(data=connector_df,aes(x=x_id+0.08,xend=x_id+0.24,y=Expanded_US,yend=Expanded_NA),color="grey50",linewidth=0.6,alpha=0.80) +
    geom_point(data=df_rng,aes(x=x_plot,y=Impact,size=Share,color=Sector,fill=PointFill,shape=Scenario),alpha=0.90,stroke=1.7) +
    scale_x_continuous(breaks=seq_along(sector_order_key),labels=NULL,expand=expansion(mult=c(0.04,0.04))) +
    coord_cartesian(ylim=c(y_min,y_max),clip="off") +
    scale_shape_manual(
      values=scenario_shapes,
      breaks=scenario_levels,
      name="GLP-1 adoption scenario",
      guide=guide_legend(
        order=1,nrow=2,byrow=TRUE,
        override.aes=list(size=8,color="black",fill=c("black",NA,"black",NA),stroke=1.7,alpha=0.90)
      )
    ) +
    scale_color_manual(
      values=sector_colors,
      breaks=sector_order_key,
      name="Sector",
      guide=guide_legend(
        order=2,nrow=2,byrow=TRUE,
        override.aes=list(size=10,shape=21,fill=unname(sector_colors[sector_order_key]),alpha=0.90,stroke=1)
      )
    ) +
    scale_fill_manual(values=c(sector_colors,"Hollow"=NA),guide="none") +
    scale_size_area(
      max_size=23,
      breaks=c(0.0025,0.0050,0.0100,0.0150),
      labels=c("0.25%","0.50%","1.00%","1.50%"),
      name="Sector share (2030)",
      guide=guide_legend(
        order=3,nrow=1,
        override.aes=list(shape=21,color="grey30",fill="grey30",alpha=0.90,stroke=1)
      )
    ) +
    labs(x=NULL,y=NULL) +
    theme_minimal(base_size=46) +
    theme(
      aspect.ratio=1,
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.text.y=element_text(size=32,color="black"),
      panel.grid.major.x=element_line(linewidth=0.5,color="grey88"),
      panel.grid.minor.x=element_blank(),
      legend.key.width=unit(1.7,"lines"),
      legend.key.height=unit(1.7,"lines"),
      legend.text=element_text(size=32,color="black"),
      legend.title=element_text(size=35,face="bold",hjust=0,color="black"),
      legend.title.align=0,
      legend.spacing.y=unit(0.25,"cm"),
      plot.margin=margin(10,10,10,10)
    )
}
# ======================================================
# 6. Panel-label helper
# ======================================================
panel <- function(p,title) {
  label_strip <- ggplot() +
    annotate("rect",xmin=-Inf,xmax=Inf,ymin=-Inf,ymax=Inf,fill="grey90",color=NA) +
    annotate("text",x=0.5,y=0.5,label=title,size=12,fontface="bold") +
    scale_x_continuous(limits=c(0,1),expand=c(0,0)) +
    scale_y_continuous(limits=c(0,1),expand=c(0,0)) +
    theme_void() +
    theme(plot.margin=margin(0,0,0,0))
  (p/label_strip)+plot_layout(heights=c(1,0.10))
}
# ======================================================
# 7. Final four-panel figure functions
# ======================================================
make_final_figure_four_panels <- function(pA,pB,pC,pD,y_label_text) {
  panel_grid <- ((pA|pB)/(pC|pD)) +
    plot_layout(guides="collect") &
    theme(legend.position="bottom",legend.box="vertical",legend.box.spacing=unit(0.6,"cm"))
  y_label <- ggplot() +
    annotate("text",x=0.5,y=0.5,label=y_label_text,angle=90,size=16,fontface="bold") +
    scale_x_continuous(limits=c(0,1),expand=c(0,0)) +
    scale_y_continuous(limits=c(0,1),expand=c(0,0)) +
    theme_void()
  plot_grid(y_label,panel_grid,ncol=2,rel_widths=c(0.05,0.95))
}
build_adoption_figure <- function(Po_rng,Qo_rng,Qmw_rng,Qxw_rng,y_label_text) {
  pA <- panel(make_plot_fourtype(Po_rng,-10,5),"Panel A: Farmgate price")
  pB <- panel(make_plot_fourtype(Qo_rng,-10,5),"Panel B: Domestic production volume")
  pC <- panel(make_plot_fourtype(Qmw_rng,-10,5),"Panel C: Import volume")
  pD <- panel(make_plot_fourtype(Qxw_rng,-10,5),"Panel D: Export volume")
  make_final_figure_four_panels(pA,pB,pC,pD,y_label_text)
}
# ======================================================
# 8. Low- and high-adoption figures
# ======================================================
Figure_low_adoption <- build_adoption_figure(
  Low_Po_rng,Low_Qo_rng,Low_Qmw_rng,Low_Qxw_rng,
  "Sectoral Impact—Low Pathway (%)"
)
Figure_high_adoption <- build_adoption_figure(
  High_Po_rng,High_Qo_rng,High_Qmw_rng,High_Qxw_rng,
  "Sectoral Impact—High Pathway (%)"
)
Figure_low_adoption
Figure_high_adoption












#############################################
# Figure: Land-use change
#
# Bubble center = Current clinical coverage, Shock US
# Bubble size = same sector share as Figure 1
# Open circle = Current clinical coverage, Shock NA
# Black short bar = Expanded weight-management, Shock US
# Gray short bar = Expanded weight-management, Shock NA
# Vertical line links the two expanded weight-management bars
#
# Updated for:
#   - Low pathway
#   - High pathway
#############################################

# ======================================================
# IMPORTANT
# This block assumes the following already exist
# from your Figure 1 code:
#   Sector_list
#   sector_size
#   sector_order_key
#   sector_colors
# ======================================================


# ======================================================
# 1. Land-use-change data: LOW pathway
# ======================================================

# Low pathway - Current clinical coverage - US
Land_low_current_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.12, -1.49, -1.65, -0.96, -0.62,
    -1.12, -1.78, -1.52, -1.24
  )
)

# Low pathway - Current clinical coverage - NA
Land_low_current_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.15, -1.56, -1.70, -0.95, -0.64,
    -1.11, -1.77, -1.51, -1.24
  )
)

# Low pathway - Expanded weight-management access - US
Land_low_expanded_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.99, -2.68, -2.96, -1.72, -1.10,
    -2.01, -3.22, -2.74, -2.23
  )
)

# Low pathway - Expanded weight-management access - NA
Land_low_expanded_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -2.05, -2.79, -3.04, -1.72, -1.13,
    -2.00, -3.20, -2.73, -2.24
  )
)


# ======================================================
# 2. Land-use-change data: HIGH pathway
# ======================================================

# High pathway - Current clinical coverage - US
Land_high_current_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.11, -1.49, -1.65, -0.96, -0.62,
    -1.12, -1.78, -1.51, -1.23
  )
)

# High pathway - Current clinical coverage - NA
Land_high_current_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.14, -1.57, -1.71, -0.95, -0.64,
    -1.11, -1.77, -1.51, -1.24
  )
)

# High pathway - Expanded weight-management access - US
Land_high_expanded_us <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -1.99, -2.68, -2.96, -1.73, -1.12,
    -2.02, -3.22, -2.74, -2.23
  )
)

# High pathway - Expanded weight-management access - NA
Land_high_expanded_na <- data.frame(
  Sector = Sector_list,
  check.names = FALSE,
  `2030` = c(
    -2.05, -2.79, -3.03, -1.72, -1.13,
    -2.01, -3.21, -2.73, -2.24
  )
)


# ======================================================
# 3. Scenario definitions
# ======================================================
scenario_levels <- c("Current: U.S. only","Current: North America","Expanded: U.S. only","Expanded: North America")
scenario_shapes <- c("Current: U.S. only"=21,"Current: North America"=21,"Expanded: U.S. only"=23,"Expanded: North America"=23)
scenario_offsets <- c("Current: U.S. only"=-0.24,"Current: North America"=-0.08,"Expanded: U.S. only"=0.08,"Expanded: North America"=0.24)
# ======================================================
# 5. Combine land-use data
# ======================================================
make_land_rng <- function(current_us,current_na,expanded_us,expanded_na) {
  current_us %>%
    rename(Current_US=`2030`) %>%
    left_join(current_na %>% rename(Current_NA=`2030`),by="Sector") %>%
    left_join(expanded_us %>% rename(Expanded_US=`2030`),by="Sector") %>%
    left_join(expanded_na %>% rename(Expanded_NA=`2030`),by="Sector") %>%
    left_join(sector_size,by="Sector") %>%
    mutate(Sector=factor(Sector,levels=sector_order_key),x_id=as.numeric(Sector)) %>%
    pivot_longer(cols=c(Current_US,Current_NA,Expanded_US,Expanded_NA),names_to="Scenario_code",values_to="Impact") %>%
    mutate(
      Scenario=case_when(
        Scenario_code=="Current_US"~"Current: U.S. only",
        Scenario_code=="Current_NA"~"Current: North America",
        Scenario_code=="Expanded_US"~"Expanded: U.S. only",
        Scenario_code=="Expanded_NA"~"Expanded: North America"
      ),
      Scenario=factor(Scenario,levels=scenario_levels),
      x_plot=x_id+unname(scenario_offsets[as.character(Scenario)]),
      PointFill=case_when(
        Scenario_code %in% c("Current_US","Expanded_US")~as.character(Sector),
        Scenario_code %in% c("Current_NA","Expanded_NA")~"Hollow"
      )
    )
}
Land_rng_low <- make_land_rng(Land_low_current_us,Land_low_current_na,Land_low_expanded_us,Land_low_expanded_na)
Land_rng_high <- make_land_rng(Land_high_current_us,Land_high_current_na,Land_high_expanded_us,Land_high_expanded_na)
# ======================================================
# 6. Land-use plot function
# Circle = current coverage
# Diamond = expanded access
# Filled = U.S. only
# Hollow = North America
# Size = sector share
# ======================================================
make_land_plot <- function(df_rng,y_min,y_max,y_label_text) {
  connector_df <- df_rng %>%
    select(Sector,x_id,Scenario_code,Impact) %>%
    pivot_wider(names_from=Scenario_code,values_from=Impact)
  ggplot() +
    geom_hline(yintercept=0,linewidth=0.9,color="black") +
    geom_segment(data=connector_df,aes(x=x_id-0.24,xend=x_id-0.08,y=Current_US,yend=Current_NA),color="grey55",linewidth=0.6,alpha=0.80) +
    geom_segment(data=connector_df,aes(x=x_id+0.08,xend=x_id+0.24,y=Expanded_US,yend=Expanded_NA),color="grey55",linewidth=0.6,alpha=0.80) +
    geom_point(data=df_rng,aes(x=x_plot,y=Impact,size=Share,color=Sector,fill=PointFill,shape=Scenario),alpha=0.90,stroke=1.7) +
    scale_x_continuous(breaks=seq_along(sector_order_key),labels=NULL,expand=expansion(mult=c(0.04,0.04))) +
    coord_cartesian(ylim=c(y_min,y_max),clip="off") +
    scale_shape_manual(
      values=scenario_shapes,
      breaks=scenario_levels,
      name="GLP-1 adoption scenario",
      guide=guide_legend(
        order=1,
        nrow=2,
        byrow=TRUE,
        override.aes=list(size=7,color="black",fill=c("black",NA,"black",NA),stroke=1.5,alpha=0.90)
      )
    ) +
    scale_color_manual(
      values=sector_colors,
      breaks=sector_order_key,
      name="Sector",
      guide=guide_legend(
        order=2,
        nrow=3,
        byrow=TRUE,
        override.aes=list(size=7.5,shape=21,fill=unname(sector_colors[sector_order_key]),alpha=0.90,stroke=1)
      )
    ) +
    scale_fill_manual(values=c(sector_colors,"Hollow"=NA),guide="none") +
    scale_size_area(
      max_size=23,
      breaks=c(0.0025,0.0050,0.0100,0.0150),
      labels=c("0.25%","0.50%","1.00%","1.50%"),
      name="Sector share (2030)",
      guide=guide_legend(
        order=3,
        nrow=1,
        override.aes=list(shape=21,color="grey35",fill="grey35",alpha=0.90,stroke=1)
      )
    ) +
    labs(x=NULL,y=y_label_text) +
    theme_minimal(base_size=40) +
    theme(
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.x=element_blank(),
      axis.text.y=element_text(size=28,color="black"),
      axis.title.y=element_text(size=30,face="bold",margin=margin(r=10),color="black"),
      panel.grid.major.x=element_line(linewidth=0.5,color="grey88"),
      panel.grid.minor.x=element_blank(),
      legend.key.width=unit(1.15,"lines"),
      legend.key.height=unit(1.10,"lines"),
      legend.text=element_text(size=24,color="black"),
      legend.title=element_text(size=26,face="bold",hjust=0,color="black"),
      legend.title.align=0,
      legend.spacing.x=unit(0.10,"cm"),
      legend.spacing.y=unit(0.05,"cm"),
      legend.margin=margin(0,0,0,0),
      plot.margin=margin(5,8,0,8)
    )
}
# ======================================================
# 7. Low- and high-pathway plots
# ======================================================
p_land_low <- make_land_plot(
  Land_rng_low,
  y_min=-4,
  y_max=1,
  y_label_text="Land-use Change—Low Pathway (%)"
)
p_land_high <- make_land_plot(
  Land_rng_high,
  y_min=-4,
  y_max=1,
  y_label_text="Land-use Change—High Pathway (%)"
)
# ======================================================
# 8. Separate figures if needed
# ======================================================
Figure_land_use_low <- p_land_low + theme(legend.position="bottom",legend.box="vertical")
Figure_land_use_high <- p_land_high + theme(legend.position="bottom",legend.box="vertical")
# ======================================================
# 9. Combined figure with ONE shared legend
# ======================================================
Figure_land_use_low_high <- (
  p_land_low |
    p_land_high
) +
  plot_layout(guides="collect") &
  theme(
    legend.position="bottom",
    legend.box="vertical",
    legend.box.spacing=unit(0.05,"cm")
  )
# ======================================================
# 10. Display
# ======================================================
Figure_land_use_low
Figure_land_use_high
Figure_land_use_low_high




#############################################
# Appendix Figure: Decomposed U.S. Trade Differences
# Difference = Shock NA - Shock US
# Low and high adoption pathways
#
# Colored bubble = Current clinical coverage
# Bubble size = 2030 sector share
# Horizontal line = Expanded weight-management access
#############################################

# ======================================================
# 0. Required packages
# ======================================================
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(patchwork)
library(cowplot)
library(grid)

# ======================================================
# 1. Sector information
# Same sector shares as the original code
# ======================================================

Sector_list <- c(
  "Cereals",
  "Fruit & Veg",
  "Fish & Seafood",
  "Meat & Meat products",
  "Veg Oil & Fats",
  "Dairy",
  "Sugar & Sugar products",
  "Bev & Tobacco",
  "Processed food"
)

sector_size <- data.frame(
  Sector = Sector_list,
  Share = c(
    0.007422,
    0.008407,
    0.002751,
    0.008285,
    0.005419,
    0.004928,
    0.002099,
    0.008368,
    0.021580
  )
)

sector_order_key <- sector_size %>%
  arrange(Share) %>%
  pull(Sector)

sector_colors <- setNames(
  scales::hue_pal()(nrow(sector_size)),
  sector_size$Sector
)

# ======================================================
# 2. LOW ADOPTION PATHWAY
# Difference = Shock NA - Shock US
# ======================================================

# ------------------------------------------------------
# 2.1 Current clinical coverage: low pathway
# Colored bubbles
# ------------------------------------------------------

# Export from U.S. to CA, MX, and ROW
Low_current_export_diff <- data.frame(
  Sector = Sector_list,
  CA  = c(-0.65, -0.90, -0.75, -0.81, -0.54, -1.28, -1.32, -1.33, -1.23),
  MX  = c(-0.69, -0.58, -0.73, -0.94, -0.65, -0.78, -0.86, -0.76, -0.71),
  ROW = c( 0.03,  0.05,  0.06,  0.03,  0.03,  0.02,  0.02,  0.01,  0.02),
  Type = "Export",
  check.names = FALSE
)

# Import to U.S. from CA, MX, and ROW
Low_current_import_diff <- data.frame(
  Sector = Sector_list,
  CA  = c( 0.09,  0.00,  0.20,  0.19,  0.07,  0.17,  0.06,  0.04,  0.07),
  MX  = c( 0.47,  0.10,  0.30,  0.39,  0.12,  0.19,  0.27,  0.05,  0.11),
  ROW = c(-0.08, -0.11, -0.14, -0.11, -0.08, -0.03, -0.07, -0.01, -0.04),
  Type = "Import",
  check.names = FALSE
)

# ------------------------------------------------------
# 2.2 Expanded weight-management access: low pathway
# Horizontal line markers
# ------------------------------------------------------

# Export from U.S. to CA, MX, and ROW
Low_expanded_export_diff <- data.frame(
  Sector = Sector_list,
  CA  = c(-0.93, -1.27, -1.06, -1.16, -0.77, -1.81, -1.88, -1.89, -1.74),
  MX  = c(-1.06, -0.91, -1.16, -1.46, -1.01, -1.22, -1.34, -1.17, -1.10),
  ROW = c( 0.04,  0.08,  0.09,  0.06,  0.04,  0.04,  0.03,  0.01,  0.02),
  Type = "Export",
  check.names = FALSE
)

# Import to U.S. from CA, MX, and ROW
Low_expanded_import_diff <- data.frame(
  Sector = Sector_list,
  CA  = c( 0.12, -0.01,  0.28,  0.25,  0.08,  0.23,  0.07,  0.05,  0.09),
  MX  = c( 0.73,  0.15,  0.46,  0.60,  0.18,  0.29,  0.41,  0.08,  0.17),
  ROW = c(-0.13, -0.15, -0.21, -0.16, -0.12, -0.05, -0.11, -0.02, -0.05),
  Type = "Import",
  check.names = FALSE
)

# ======================================================
# 3. HIGH ADOPTION PATHWAY
# Difference = Shock NA - Shock US
# ======================================================

# ------------------------------------------------------
# 3.1 Current clinical coverage: high pathway
# Colored bubbles
# ------------------------------------------------------

# Export from U.S. to CA, MX, and ROW
High_current_export_diff <- data.frame(
  Sector = Sector_list,
  CA  = c(-0.71, -1.00, -1.00, -0.96, -0.58, -1.43, -1.48, -1.49, -1.38),
  MX  = c(-0.80, -0.67, -0.84, -1.11, -0.77, -0.93, -1.02, -0.90, -0.84),
  ROW = c( 0.03,  0.07,  0.08,  0.04,  0.03,  0.02,  0.02,  0.01,  0.02),
  Type = "Export",
  check.names = FALSE
)

# Import to U.S. from CA, MX, and ROW
High_current_import_diff <- data.frame(
  Sector = Sector_list,
  CA  = c( 0.09, -0.01,  0.26,  0.22,  0.08,  0.19,  0.06,  0.04,  0.07),
  MX  = c( 0.56,  0.11,  0.32,  0.47,  0.13,  0.23,  0.32,  0.06,  0.13),
  ROW = c(-0.10, -0.12, -0.19, -0.13, -0.09, -0.03, -0.09, -0.02, -0.05),
  Type = "Import",
  check.names = FALSE
)

# ------------------------------------------------------
# 3.2 Expanded weight-management access: high pathway
# Horizontal line markers
# ------------------------------------------------------

# Export from U.S. to CA, MX, and ROW
High_expanded_export_diff <- data.frame(
  Sector = Sector_list,
  CA  = c(-0.94, -1.29, -1.08, -1.18, -0.77, -1.85, -1.91, -1.92, -1.77),
  MX  = c(-1.09, -0.92, -1.17, -1.51, -1.04, -1.27, -1.40, -1.21, -1.14),
  ROW = c( 0.05,  0.09,  0.09,  0.05,  0.04,  0.04,  0.03,  0.01,  0.02),
  Type = "Export",
  check.names = FALSE
)

# Import to U.S. from CA, MX, and ROW
High_expanded_import_diff <- data.frame(
  Sector = Sector_list,
  CA  = c( 0.12, -0.02,  0.28,  0.25,  0.09,  0.23,  0.07,  0.05,  0.09),
  MX  = c( 0.75,  0.15,  0.46,  0.62,  0.18,  0.31,  0.43,  0.08,  0.17),
  ROW = c(-0.13, -0.16, -0.20, -0.16, -0.12, -0.04, -0.11, -0.02, -0.06),
  Type = "Import",
  check.names = FALSE
)

# ======================================================
# 4. Functions for long format and data combination
# ======================================================

to_long_current <- function(df) {
  df %>%
    pivot_longer(
      cols = c(CA, MX, ROW),
      names_to = "Partner",
      values_to = "Current"
    ) %>%
    select(Sector, Type, Partner, Current)
}

to_long_expanded <- function(df) {
  df %>%
    pivot_longer(
      cols = c(CA, MX, ROW),
      names_to = "Partner",
      values_to = "Expanded"
    ) %>%
    select(Sector, Type, Partner, Expanded)
}

make_trade_pathway_df <- function(
    current_export,
    current_import,
    expanded_export,
    expanded_import
) {
  current_long <- bind_rows(
    to_long_current(current_export),
    to_long_current(current_import)
  )
  
  expanded_long <- bind_rows(
    to_long_expanded(expanded_export),
    to_long_expanded(expanded_import)
  )
  
  current_long %>%
    left_join(
      expanded_long,
      by = c("Sector", "Type", "Partner")
    ) %>%
    left_join(
      sector_size,
      by = "Sector"
    ) %>%
    mutate(
      Sector = factor(Sector, levels = sector_order_key),
      Partner = factor(Partner, levels = c("CA", "MX", "ROW")),
      Type = factor(Type, levels = c("Import", "Export"))
    ) %>%
    arrange(Type, Partner, Sector)
}

trade_low <- make_trade_pathway_df(
  Low_current_export_diff,
  Low_current_import_diff,
  Low_expanded_export_diff,
  Low_expanded_import_diff
)

trade_high <- make_trade_pathway_df(
  High_current_export_diff,
  High_current_import_diff,
  High_expanded_export_diff,
  High_expanded_import_diff
)

# ======================================================
# 5. Plot function
# Bubble center = current clinical coverage
# Bubble size = relative sector size
# Horizontal line = expanded weight-management access
# ======================================================

scenario_levels <- c("Current clinical coverage","Expanded weight-management access")
scenario_shapes <- c("Current clinical coverage"=21,"Expanded weight-management access"=23)
scenario_offsets <- c("Current clinical coverage"=-0.14,"Expanded weight-management access"=0.14)
make_trade_pathway_plot <- function(df,partner,type,y_min,y_max) {
  df2 <- df %>%
    filter(Partner==partner,Type==type) %>%
    mutate(Sector=factor(Sector,levels=sector_order_key),x_id=as.numeric(Sector)) %>%
    pivot_longer(cols=c(Current,Expanded),names_to="Scenario_code",values_to="Impact") %>%
    mutate(
      Scenario=case_when(
        Scenario_code=="Current"~"Current clinical coverage",
        Scenario_code=="Expanded"~"Expanded weight-management access"
      ),
      Scenario=factor(Scenario,levels=scenario_levels),
      x_plot=x_id+unname(scenario_offsets[as.character(Scenario)])
    )
  ggplot(df2) +
    geom_hline(yintercept=0,linewidth=1,color="black") +
    geom_point(
      aes(x=x_plot,y=Impact,size=Share,color=Sector,shape=Scenario),
      fill=NA,alpha=0.90,stroke=1.7
    ) +
    scale_x_continuous(
      breaks=seq_along(sector_order_key),
      labels=sector_order_key,
      expand=expansion(mult=c(0.05,0.05))
    ) +
    coord_cartesian(ylim=c(y_min,y_max),clip="off") +
    scale_shape_manual(
      values=scenario_shapes,
      breaks=scenario_levels,
      name="GLP-1 adoption scenario",
      guide=guide_legend(
        order=1,nrow=1,
        override.aes=list(size=9,color="black",fill=NA,stroke=1.7,alpha=0.90)
      )
    ) +
    scale_size_area(
      max_size=23,
      breaks=c(0.0025,0.0050,0.0100,0.0150),
      labels=c("0.25%","0.50%","1.00%","1.50%"),
      name="Sector share (2030)",
      guide=guide_legend(
        order=3,nrow=1,
        override.aes=list(shape=21,color="grey30",fill="grey30",alpha=0.90,stroke=1)
      )
    ) +
    scale_color_manual(
      values=sector_colors,
      breaks=sector_order_key,
      name="Sector",
      guide=guide_legend(
        order=2,nrow=2,byrow=TRUE,
        override.aes=list(
          size=10,
          shape=21,
          fill=unname(sector_colors[sector_order_key]),
          alpha=0.90,
          stroke=1
        )
      )
    ) +
    labs(x=NULL,y=NULL) +
    theme_minimal(base_size=40) +
    theme(
      aspect.ratio=1,
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.x=element_blank(),
      legend.key.size=unit(2,"lines"),
      legend.text=element_text(size=32),
      legend.title=element_text(size=34,face="bold",hjust=0),
      legend.title.align=0,
      plot.margin=margin(10,10,10,10)
    )
}

# ======================================================
# 6. Panel wrapper with grey label bar
# ======================================================

panel_diff <- function(plot_obj, label_text) {
  
  label_plot <- ggplot(
    data.frame(x = 0.5, y = 0.5),
    aes(x, y)
  ) +
    geom_rect(
      aes(
        xmin = 0,
        xmax = 1,
        ymin = 0,
        ymax = 1
      ),
      fill = "grey90",
      color = NA
    ) +
    geom_text(
      label = label_text,
      size = 12,
      fontface = "bold",
      hjust = 0.5,
      vjust = 0.5
    ) +
    scale_x_continuous(
      limits = c(0, 1),
      expand = c(0, 0)
    ) +
    scale_y_continuous(
      limits = c(0, 1),
      expand = c(0, 0)
    ) +
    theme_void()
  
  plot_obj / label_plot +
    plot_layout(heights = c(12, 2))
}

# ======================================================
# 7. Build one complete 3-by-2 figure for a pathway
# Rows: CA, MX, ROW
# Columns: Import and Export
# Each pathway receives its own y-axis label
# ======================================================

build_trade_pathway_figure <- function(df, y_label_text) {
  
  pI_CA <- panel_diff(
    make_trade_pathway_plot(
      df,
      partner = "CA",
      type = "Import",
      y_min = -5,
      y_max = 5
    ),
    "Import from Canada"
  )
  
  pE_CA <- panel_diff(
    make_trade_pathway_plot(
      df,
      partner = "CA",
      type = "Export",
      y_min = -5,
      y_max = 5
    ),
    "Export to Canada"
  )
  
  pI_MX <- panel_diff(
    make_trade_pathway_plot(
      df,
      partner = "MX",
      type = "Import",
      y_min = -5,
      y_max = 5
    ),
    "Import from Mexico"
  )
  
  pE_MX <- panel_diff(
    make_trade_pathway_plot(
      df,
      partner = "MX",
      type = "Export",
      y_min = -5,
      y_max = 5
    ),
    "Export to Mexico"
  )
  
  pI_ROW <- panel_diff(
    make_trade_pathway_plot(
      df,
      partner = "ROW",
      type = "Import",
      y_min = -10,
      y_max = 5
    ),
    "Import from Rest of World"
  )
  
  pE_ROW <- panel_diff(
    make_trade_pathway_plot(
      df,
      partner = "ROW",
      type = "Export",
      y_min = -10,
      y_max = 5
    ),
    "Export to Rest of World"
  )
  
  panel_grid <- (
    (pI_CA | pE_CA) /
      (pI_MX | pE_MX) /
      (pI_ROW | pE_ROW)
  ) +
    plot_layout(guides = "collect") &
    theme(
      legend.position = "bottom",
      legend.box = "vertical",
      legend.box.spacing = unit(1, "cm")
    )
  
  pathway_y_label <- ggplot() +
    annotate(
      "text",
      x = 0.5,
      y = 0.5,
      label = y_label_text,
      angle = 90,
      size = 14,
      fontface = "bold"
    ) +
    theme_void()
  
  plot_grid(
    pathway_y_label,
    panel_grid,
    ncol = 2,
    rel_widths = c(0.04, 0.96)
  )
}

# ======================================================
# 8. Create separate low- and high-pathway figures
# ======================================================

figure_trade_low <- build_trade_pathway_figure(
  trade_low,
  "U.S. Trade Response Differences—Low Pathway (percentage points)"
)

figure_trade_high <- build_trade_pathway_figure(
  trade_high,
  "U.S. Trade Response Differences—High Pathway (percentage points)"
)

# ======================================================
# 9. Display the two figures separately
# ======================================================

figure_trade_low
figure_trade_high






#############################################
# Appendix Figures: GHG and CO2 Emissions by Activities
# Low and High adoption pathways plotted separately
#
# Bubble center = Current Clinical Coverage, Shock US
# Bubble size   = sector share
# Open circle   = Current Clinical Coverage, Shock NA
# Short bars    = Expanded Weight-Management Access, Shock US and Shock NA
# Vertical line links the two expanded-access bars
#############################################

# ======================================================
# 0. Required packages
# ======================================================
library(dplyr)
library(ggplot2)
library(scales)
library(patchwork)
library(cowplot)
library(grid)

# ======================================================
# 1. Sector information
# Keep the original sector shares unchanged
# ======================================================

Sector_list <- c(
  "Cereals",
  "Fruit & Veg",
  "Fish & Seafood",
  "Meat & Meat products",
  "Veg Oil & Fats",
  "Dairy",
  "Sugar & Sugar products",
  "Bev & Tobacco",
  "Processed food"
)

sector_size <- data.frame(
  Sector = Sector_list,
  
  # CO2 sectoral shares from the original R code
  CO2_share = c(
    0.002572671,  # Cereal
    0.000902421,  # Veg_Fru
    0.000404898,  # Fish
    0.002338355,  # Meat
    0.002401265,  # Vegoil_F
    0.001744601,  # Dairy
    0.001547423,  # Sugar
    0.001911849,  # Beve_Tob
    0.007876525   # Oth_Fod
  ),
  
  # Total GHG sectoral shares from the original R code
  GHG_share = c(
    0.018127851,  # Cereal
    0.004815414,  # Veg_Fru
    0.000309525,  # Fish
    0.001740481,  # Meat
    0.006855848,  # Vegoil_F
    0.001310898,  # Dairy
    0.001148547,  # Sugar
    0.001473797,  # Beve_Tob
    0.008712234   # Oth_Fod
  )
)

sector_legend_order <- sector_size %>%
  arrange(GHG_share) %>%
  pull(Sector)

sector_colors <- setNames(
  scales::hue_pal()(length(Sector_list)),
  Sector_list
)

Year_used <- "2030"

# The supplied Excel screenshots are used at 2030, matching the appendix figure.
if (Year_used != "2030") {
  stop("This script currently contains the supplied 2030 Excel values only.")
}

# ======================================================
# 2. Low-adoption pathway: 2030 Excel values
# ======================================================

# ------------------------------
# Low pathway: Total GHG emissions
# ------------------------------

GHG_low_current_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.71, -2.13, -2.26, -3.27, -1.34, -3.61, -3.71, -4.43, -3.55),
  check.names = FALSE
)

GHG_low_current_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.76, -2.23, -2.32, -3.30, -1.38, -3.62, -3.72, -4.45, -3.59),
  check.names = FALSE
)

GHG_low_expanded_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.04, -3.82, -4.04, -5.83, -2.38, -6.43, -6.63, -7.89, -6.33),
  check.names = FALSE
)

GHG_low_expanded_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.12, -3.97, -4.13, -5.88, -2.44, -6.45, -6.64, -7.93, -6.39),
  check.names = FALSE
)

# ------------------------------
# Low pathway: Total CO2 emissions
# ------------------------------

CO2_low_current_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.81, -2.25, -2.27, -3.28, -1.47, -3.62, -4.60, -4.44, -3.73),
  check.names = FALSE
)

CO2_low_current_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.87, -2.36, -2.33, -3.30, -1.51, -3.63, -4.61, -4.46, -3.77),
  check.names = FALSE
)

CO2_low_expanded_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.22, -4.04, -4.05, -5.84, -2.60, -6.44, -8.19, -7.91, -6.65),
  check.names = FALSE
)

CO2_low_expanded_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.31, -4.19, -4.15, -5.89, -2.67, -6.46, -8.21, -7.94, -6.71),
  check.names = FALSE
)

# ======================================================
# 3. High-adoption pathway: 2030 Excel values
# ======================================================

# ------------------------------
# High pathway: Total GHG emissions
# ------------------------------

GHG_high_current_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.70, -2.13, -2.26, -3.28, -1.34, -3.60, -3.71, -4.42, -3.54),
  check.names = FALSE
)

GHG_high_current_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.77, -2.24, -2.34, -3.31, -1.39, -3.62, -3.72, -4.45, -3.59),
  check.names = FALSE
)

GHG_high_expanded_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.04, -3.83, -4.04, -5.84, -2.38, -6.44, -6.64, -7.90, -6.34),
  check.names = FALSE
)

GHG_high_expanded_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.13, -3.97, -4.14, -5.89, -2.45, -6.46, -6.65, -7.93, -6.39),
  check.names = FALSE
)

# ------------------------------
# High pathway: Total CO2 emissions
# ------------------------------

CO2_high_current_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.81, -2.25, -2.27, -3.28, -1.47, -3.61, -4.59, -4.43, -3.72),
  check.names = FALSE
)

CO2_high_current_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-1.88, -2.37, -2.35, -3.32, -1.52, -3.62, -4.61, -4.46, -3.77),
  check.names = FALSE
)

CO2_high_expanded_us <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.23, -4.04, -4.06, -5.86, -2.61, -6.45, -8.20, -7.92, -6.65),
  check.names = FALSE
)

CO2_high_expanded_na <- data.frame(
  Sector = Sector_list,
  `2030` = c(-3.32, -4.19, -4.15, -5.90, -2.67, -6.47, -8.22, -7.95, -6.72),
  check.names = FALSE
)

# ======================================================
# 4. Combine data for an emission panel
# ======================================================
scenario_levels <- c("Current: U.S. only","Current: North America","Expanded: U.S. only","Expanded: North America")
scenario_shapes <- c("Current: U.S. only"=21,"Current: North America"=21,"Expanded: U.S. only"=23,"Expanded: North America"=23)
scenario_offsets <- c("Current: U.S. only"=-0.24,"Current: North America"=-0.08,"Expanded: U.S. only"=0.08,"Expanded: North America"=0.24)
make_emission_panel_df <- function(current_us,current_na,expanded_us,expanded_na,emission_type,share_column) {
  current_us %>%
    select(Sector,all_of(Year_used)) %>% rename(Current_US=all_of(Year_used)) %>%
    left_join(current_na %>% select(Sector,all_of(Year_used)) %>% rename(Current_NA=all_of(Year_used)),by="Sector") %>%
    left_join(expanded_us %>% select(Sector,all_of(Year_used)) %>% rename(Expanded_US=all_of(Year_used)),by="Sector") %>%
    left_join(expanded_na %>% select(Sector,all_of(Year_used)) %>% rename(Expanded_NA=all_of(Year_used)),by="Sector") %>%
    left_join(sector_size,by="Sector") %>%
    mutate(Share=.data[[share_column]],Emission=emission_type) %>%
    arrange(Share) %>%
    mutate(Sector=factor(Sector,levels=Sector),x_id=row_number()) %>%
    pivot_longer(cols=c(Current_US,Current_NA,Expanded_US,Expanded_NA),names_to="Scenario_code",values_to="Impact") %>%
    mutate(
      Scenario=case_when(
        Scenario_code=="Current_US"~"Current: U.S. only",
        Scenario_code=="Current_NA"~"Current: North America",
        Scenario_code=="Expanded_US"~"Expanded: U.S. only",
        Scenario_code=="Expanded_NA"~"Expanded: North America"
      ),
      Scenario=factor(Scenario,levels=scenario_levels),
      x_plot=x_id+unname(scenario_offsets[as.character(Scenario)]),
      PointFill=case_when(
        Scenario_code %in% c("Current_US","Expanded_US")~as.character(Sector),
        Scenario_code %in% c("Current_NA","Expanded_NA")~"Hollow"
      )
    )
}
GHG_low_rng <- make_emission_panel_df(GHG_low_current_us,GHG_low_current_na,GHG_low_expanded_us,GHG_low_expanded_na,"Total GHG emissions","GHG_share")
CO2_low_rng <- make_emission_panel_df(CO2_low_current_us,CO2_low_current_na,CO2_low_expanded_us,CO2_low_expanded_na,"Total CO2 emissions","CO2_share")
GHG_high_rng <- make_emission_panel_df(GHG_high_current_us,GHG_high_current_na,GHG_high_expanded_us,GHG_high_expanded_na,"Total GHG emissions","GHG_share")
CO2_high_rng <- make_emission_panel_df(CO2_high_current_us,CO2_high_current_na,CO2_high_expanded_us,CO2_high_expanded_na,"Total CO2 emissions","CO2_share")
# ======================================================
# 5. Plot function
# Circle = Current clinical coverage
# Diamond = Expanded weight-management access
# Filled = U.S.-only
# Hollow = North America-wide
# Marker size = sector emissions share
# ======================================================
make_emission_plot <- function(df_rng,y_min,y_max) {
  sector_order_key <- levels(df_rng$Sector)
  ggplot(df_rng) +
    geom_hline(yintercept=0,linewidth=1,color="black") +
    geom_point(
      aes(x=x_plot,y=Impact,size=Share,color=Sector,fill=PointFill,shape=Scenario),
      alpha=0.90,stroke=1.7
    ) +
    scale_x_continuous(
      breaks=seq_along(sector_order_key),
      labels=sector_order_key,
      expand=expansion(mult=c(0.05,0.05))
    ) +
    coord_cartesian(ylim=c(y_min,y_max),clip="off") +
    scale_shape_manual(
      values=scenario_shapes,
      breaks=scenario_levels,
      name="GLP-1 adoption scenario",
      guide=guide_legend(
        order=1,nrow=2,byrow=TRUE,
        override.aes=list(
          size=9,
          color="black",
          fill=c("black",NA,"black",NA),
          stroke=1.7,
          alpha=0.90
        )
      )
    ) +
    scale_fill_manual(
      values=c(sector_colors,"Hollow"=NA),
      guide="none"
    ) +
    scale_color_manual(
      values=sector_colors,
      breaks=sector_legend_order,
      limits=sector_legend_order,
      name="Sector",
      guide=guide_legend(
        order=2,nrow=2,byrow=TRUE,
        override.aes=list(
          size=10,
          shape=21,
          fill=unname(sector_colors[sector_legend_order]),
          alpha=0.90,
          stroke=1
        )
      )
    ) +
    scale_size_area(
      max_size=26,
      limits=c(0,0.020),
      breaks=c(0.0050,0.0100,0.0150,0.0200),
      labels=c("0.50%","1.00%","1.50%","2.00%"),
      name="Sector share (2030)",
      guide=guide_legend(
        order=3,nrow=1,
        override.aes=list(
          shape=21,
          color="grey30",
          fill="grey30",
          alpha=0.90,
          stroke=1
        )
      )
    ) +
    labs(x=NULL,y=NULL) +
    theme_minimal(base_size=40) +
    theme(
      aspect.ratio=1,
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.x=element_blank(),
      legend.key.size=unit(2,"lines"),
      legend.text=element_text(size=32),
      legend.title=element_text(size=34,face="bold",hjust=0),
      legend.title.align=0,
      plot.margin=margin(10,10,10,10)
    )
}
# ======================================================
# 6. Panel wrapper with grey label bar
# ======================================================

panel_emission <- function(plot_obj, label_text) {
  label_plot <- ggplot(
    data.frame(x = 0.5, y = 0.5),
    aes(x, y)
  ) +
    geom_rect(
      aes(
        xmin = 0,
        xmax = 1,
        ymin = 0,
        ymax = 1
      ),
      fill = "grey90",
      color = NA
    ) +
    geom_text(
      label = label_text,
      size = 12,
      fontface = "bold",
      hjust = 0.5,
      vjust = 0.5
    ) +
    scale_x_continuous(
      limits = c(0, 1),
      expand = c(0, 0)
    ) +
    scale_y_continuous(
      limits = c(0, 1),
      expand = c(0, 0)
    ) +
    theme_void()
  
  plot_obj / label_plot +
    plot_layout(heights = c(12, 1))
}

# ======================================================
# 7. Build one separate figure for each pathway
# ======================================================

build_emission_figure <- function(
    GHG_rng,
    CO2_rng,
    y_label_text,
    y_min = -10,
    y_max = 2
) {
  p_GHG <- panel_emission(
    make_emission_plot(
      GHG_rng,
      y_min = y_min,
      y_max = y_max
    ),
    paste0("Total GHG Emissions")
  )
  
  p_CO2 <- panel_emission(
    make_emission_plot(
      CO2_rng,
      y_min = y_min,
      y_max = y_max
    ),
    paste0("Total CO2 Emissions")
  )
  
  panel_grid <- (
    p_GHG | p_CO2
  ) +
    plot_layout(guides = "collect") &
    theme(
      legend.position = "bottom",
      legend.box = "vertical",
      legend.box.spacing = unit(1, "cm")
    )
  
  y_label <- ggplot() +
    annotate(
      "text",
      x = 0.5,
      y = 0.5,
      label = y_label_text,
      angle = 90,
      size = 14,
      fontface = "bold"
    ) +
    theme_void()
  
  plot_grid(
    y_label,
    panel_grid,
    ncol = 2,
    rel_widths = c(0.04, 0.96)
  )
}

# ======================================================
# 8. Separate final figures
# ======================================================

figure_emissions_low <- build_emission_figure(
  GHG_low_rng,
  CO2_low_rng,
  y_label_text = "Sectoral Emission Impact—Low Pathway (%)"
)

figure_emissions_high <- build_emission_figure(
  GHG_high_rng,
  CO2_high_rng,
  y_label_text = "Sectoral Emission Impact—High Pathway (%)"
)

# Display the two figures separately
figure_emissions_low
figure_emissions_high