library(readr)
library(fixest)

# ============================================================
# 1. Reproduce the first column of Table 1
#    Model includes PSI as the focal explanatory variable,
#    together with county and year fixed effects.
# ============================================================

m1 <- feols(
  ShannonBD ~ PI | county + year,
  data = d1,
  cluster = ~county
)

summary(m1)


# ============================================================
# 2. Reproduce the second column of Table 1
#    Add meteorological, socioeconomic, observation-effort,
#    and carbon-emission covariates.
# ============================================================

m2 <- feols(
  ShannonBD ~ PI + Temp + Wind + Pop + Duration + Carbon |
    county + year,
  data = d1,
  cluster = ~county
)

summary(m2)


# ============================================================
# 3. Reproduce the third column of Table 1
#    Full model including all explanatory covariates.
#    County and year fixed effects are included.
# ============================================================

m3 <- feols(
  ShannonBD ~ PI + Temp + Wind + Pop + Duration + Carbon +
    Water + Green + Farm + Grass |
    county + year,
  data = d1,
  cluster = ~county
)

summary(m3)


# ============================================================
# 4. Estimate the Within R2 associated with Duration alone
#    This model includes Duration as the only time-varying
#    explanatory variable, in addition to the fixed effects.
# ============================================================

m4 <- feols(
  ShannonBD ~ Duration | county + year,
  data = d1,
  cluster = ~county
)

summary(m4)


# ============================================================
# 5. Calculate the unique contribution of PI
#
#    The unique contribution of PI is calculated as the
#    difference between the Within R2 of the full model and
#    that of a reduced model in which PI is removed:
#
#    Unique contribution of PI =
#    Within R2(full model) - Within R2(model without PI)
# ============================================================

# Full model containing PI and all other covariates
fullmodel <- feols(
  ShannonBD ~ PI + Temp + Wind + Pop + Duration + Carbon +
    Water + Green + Farm + Grass |
    county + year,
  data = d1,
  cluster = ~county
)

# Reduced model excluding PI
reducePImodel <- feols(
  ShannonBD ~ Temp + Wind + Pop + Duration + Carbon +
    Water + Green + Farm + Grass |
    county + year,
  data = d1,
  cluster = ~county
)

# Unique contribution of PI to the Within R2
unique_PI <- fitstat(
  fullmodel, "wr2", simplify = TRUE
) -
  fitstat(
    reducePImodel, "wr2", simplify = TRUE
  )

unique_PI


# ============================================================
# 6. Calculate the unique contribution of Duration
#
#    Unique contribution of Duration =
#    Within R2(full model) - Within R2(model without Duration)
# ============================================================

# Reduced model excluding Duration
reduceDurationmodel <- feols(
  ShannonBD ~ PI + Temp + Wind + Pop + Carbon +
    Water + Green + Farm + Grass |
    county + year,
  data = d1,
  cluster = ~county
)

# Unique contribution of Duration to the Within R2
unique_Duration <- fitstat(
  fullmodel, "wr2", simplify = TRUE
) -
  fitstat(
    reduceDurationmodel, "wr2", simplify = TRUE
  )

unique_Duration


# ============================================================
# 7. Compare the unique contributions of Duration and PI
#
#    Ratio =
#    Unique contribution of Duration /
#    Unique contribution of PI
#
#    A ratio of 83.84 indicates that the unique contribution
#    of Duration is approximately 84 times that of PI.
# ============================================================

0.038567 / 0.00046