1. Data Preparation

main_path <- "C:/Users/吳聿絜/Desktop/Empirical IO"
setwd(main_path)
data <- read.csv("data.csv")

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(fixest)
## Warning: package 'fixest' was built under R version 4.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
library(modelsummary)
## Warning: package 'modelsummary' was built under R version 4.5.3
library(gt)
## Warning: package 'gt' was built under R version 4.5.3
annual_data <- data %>%
  group_by(food_type, year) %>%
  summarise(
    poison = first(poison),
    price_mean = mean(price_index, na.rm = TRUE),
    price_sd   = sd(price_index, na.rm = TRUE),
    price_cv   = sd(price_index, na.rm = TRUE) / mean(price_index, na.rm = TRUE),
    price_max  = max(price_index, na.rm = TRUE),
    price_min  = min(price_index, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(food_type, year) %>%
  group_by(food_type) %>%
  mutate(lag_poison = lag(poison, 1)) %>%
  ungroup()

2. Wide Format Transformation

wide_df <- annual_data %>%
  select(year, food_type, price_mean, poison) %>%
  pivot_wider(
    names_from = food_type,
    values_from = c(price_mean, poison),
    names_sep = "_"
  ) %>%
  arrange(year) %>%
  mutate(
    lag_poison_乳類 = lag(poison_乳類),
    lag_poison_水產加工品 = lag(poison_水產加工品),
    lag_poison_水產品 = lag(poison_水產品),
    lag_poison_穀類 = lag(poison_穀類),
    lag_poison_肉類及其加工品 = lag(poison_肉類及其加工品),
    lag_poison_蔬果類 = lag(poison_蔬果類),
    lag_poison_蛋類 = lag(poison_蛋類)
  )

3. Regression Models

3.1 Overall Panel Fixed Effects Model

model_overall <- feols(
  price_mean ~ lag_poison | food_type + year,
  data = annual_data
)
## NOTE: 7 observations removed because of NA values (RHS: 7).
summary(model_overall)
## OLS estimation, Dep. Var.: price_mean
## Observations: 63
## Fixed-effects: food_type: 7,  year: 9
## Standard-errors: IID 
##            Estimate Std. Error  t value Pr(>|t|) 
## lag_poison 0.210296   0.288818 0.728124  0.47015 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 4.66869     Adj. R2: 0.707219
##                 Within R2: 0.011154

3.2 Cross-Category Substitution Models

Seafood model

model_seafood <- lm(
  price_mean_水產品 ~ lag_poison_水產品 +
    lag_poison_肉類及其加工品 +
    lag_poison_水產加工品,
  data = wide_df
)

summary(model_seafood)
## 
## Call:
## lm(formula = price_mean_水產品 ~ lag_poison_水產品 + lag_poison_肉類及其加工品 + 
##     lag_poison_水產加工品, data = wide_df)
## 
## Residuals:
##       2       3       4       5       6       7       8       9      10 
##  2.2100 -7.4335 -4.6607 -1.4983  0.4764  2.9326  3.6756  5.3726 -1.0748 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               99.879617   3.564533  28.020 1.08e-06 ***
## lag_poison_水產品         -0.841810   0.384409  -2.190   0.0801 .  
## lag_poison_肉類及其加工品  3.008674   0.806195   3.732   0.0135 *  
## lag_poison_水產加工品      0.007099   1.019316   0.007   0.9947    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.224 on 5 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7371, Adjusted R-squared:  0.5793 
## F-statistic: 4.672 on 3 and 5 DF,  p-value: 0.06505

Vegetable model

model_vegetable <- lm(
  price_mean_蔬果類 ~ lag_poison_蔬果類 +
    lag_poison_穀類,
  data = wide_df
)

summary(model_vegetable)
## 
## Call:
## lm(formula = price_mean_蔬果類 ~ lag_poison_蔬果類 + lag_poison_穀類, 
##     data = wide_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.509  -5.870  -3.331   7.411  14.241 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        91.8580     9.1053  10.088 5.51e-05 ***
## lag_poison_蔬果類   0.4021     1.8543   0.217    0.836    
## lag_poison_穀類     5.1059     3.4640   1.474    0.191    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.98 on 6 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3023, Adjusted R-squared:  0.06971 
## F-statistic:   1.3 on 2 and 6 DF,  p-value: 0.3397

Meat model

model_meat <- lm(
  price_mean_肉類及其加工品 ~ lag_poison_肉類及其加工品 +
    lag_poison_蛋類 +
    lag_poison_乳類,
  data = wide_df
)

summary(model_meat)
## 
## Call:
## lm(formula = price_mean_肉類及其加工品 ~ lag_poison_肉類及其加工品 + 
##     lag_poison_蛋類 + lag_poison_乳類, data = wide_df)
## 
## Residuals:
##        2        3        4        5        6        7        8        9 
##  -3.5498 -12.0478   1.7530   1.1361   0.3169   2.6842   5.3239   4.5419 
##       10 
##  -0.1585 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                95.3486     3.6101  26.412 1.45e-06 ***
## lag_poison_肉類及其加工品   3.8328     1.3382   2.864   0.0352 *  
## lag_poison_蛋類            -5.9394     3.0296  -1.960   0.1072    
## lag_poison_乳類             0.5009     5.0913   0.098   0.9255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.609 on 5 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.7342, Adjusted R-squared:  0.5748 
## F-statistic: 4.604 on 3 and 5 DF,  p-value: 0.06673

4. Visualization

4.1 Overall trend

ggplot(annual_data, aes(x = poison, y = price_mean)) +
  geom_point(aes(color = food_type), size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", color = "black", linetype = "dashed") +
  theme_minimal() +
  labs(
    title = "Food Price vs Food Poisoning",
    x = "Poisoning Cases",
    y = "Price Index",
    color = "Food Type"
  )
## `geom_smooth()` using formula = 'y ~ x'

4.2 By food type

ggplot(annual_data, aes(x = poison, y = price_mean, color = food_type)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(
    title = "Price Impact by Food Type",
    x = "Poisoning Cases",
    y = "Price Index"
  )
## `geom_smooth()` using formula = 'y ~ x'

5. Regression Table (GT Output)

names(model_overall$coefficients) <- c("Lag Poison (Overall)")

names(model_seafood$coefficients) <- c(
  "(Intercept)",
  "Self Shock (Seafood)",
  "Meat Substitute Shock",
  "Processed Seafood Shock"
)

names(model_vegetable$coefficients) <- c(
  "(Intercept)",
  "Self Shock (Vegetable)",
  "Grain Substitute Shock"
)

names(model_meat$coefficients) <- c(
  "(Intercept)",
  "Self Shock (Meat)",
  "Egg Shock",
  "Dairy Shock"
)

models_list <- list(
  "Overall FE" = model_overall,
  "Seafood" = model_seafood,
  "Vegetable" = model_vegetable,
  "Meat" = model_meat
)

coef_order <- c(
  "Lag Poison (Overall)",
  "Self Shock (Seafood)",
  "Self Shock (Vegetable)",
  "Self Shock (Meat)",
  "Meat Substitute Shock",
  "Processed Seafood Shock",
  "Grain Substitute Shock",
  "Egg Shock",
  "Dairy Shock",
  "(Intercept)"
)

table <- modelsummary(
  models_list,
  output = "gt",
  stars = c('*' = .1, '**' = .05, '***' = .01)
)

table
Overall FE Seafood Vegetable Meat
Lag Poison (Overall) 0.210
(0.289)
(Intercept) 99.880*** 91.858*** 95.349***
(3.565) (9.105) (3.610)
Self Shock (Seafood) -0.842*
(0.384)
Meat Substitute Shock 3.009**
(0.806)
Processed Seafood Shock 0.007
(1.019)
Self Shock (Vegetable) 0.402
(1.854)
Grain Substitute Shock 5.106
(3.464)
Self Shock (Meat) 3.833**
(1.338)
Egg Shock -5.939
(3.030)
Dairy Shock 0.501
(5.091)
Num.Obs. 63 9 9 9
R2 0.778 0.737 0.302 0.734
R2 Adj. 0.707 0.579 0.070 0.575
R2 Within 0.011
R2 Within Adj. -0.010
AIC 404.9 60.0 73.0 64.2
BIC 439.2 61.0 73.8 65.2
Log.Lik. -25.005 -32.509 -27.122
F 4.672 1.300 4.604
RMSE 4.67 3.89 8.96 4.93
FE: food_type X
FE: year X
* p < 0.1, ** p < 0.05, *** p < 0.01