## library
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.1.8
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(dplyr)
library(ggplot2)
library(plm)
##
## Attaching package: 'plm'
##
## The following objects are masked from 'package:dplyr':
##
## between, lag, lead
library(fixest)
library(broom)
library(forestplot)
## Loading required package: grid
## Loading required package: checkmate
## Loading required package: abind
library(stargazer)
##
## Please cite as:
##
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
## upload data
df <- readRDS("ozone_daily_drought_dat.rds")
## check urban
## understand data
urban_table <- table(df$Urban.2013.bin)
prop.table(urban_table) * 100
##
## Metro Nonmetro
## 74.70993 25.29007
# it appears that metro has 75%.
### re regression
# benchmark model
df$wday <- lubridate::wday(df$date)
df$omax2 <- df$Ozone.Max*1000
fe_model <- plm(
omax2 ~ USDM.categorical +Urban.2013.bin+ splines::ns(tmean, 4) + +precip + factor(wday) +
mean_elev + factor(holiday) + I(long^2) + I(lat^2) + I(long*lat)+long + lat,
data=df,
index = c("GEOID", "date"),
model = "within",
effect = "twoways",
# vcov. = vcovHC, type = "HC1"
)
stargazer(fe_model,type = "text")
##
## ===============================================================
## Dependent variable:
## -------------------------------
## omax2
## ---------------------------------------------------------------
## USDM.categoricalModerateDrought 1.609***
## (0.021)
##
## USDM.categoricalSevereDrought 2.752***
## (0.034)
##
## splines::ns(tmean, 4)1 15.434***
## (0.248)
##
## splines::ns(tmean, 4)2 28.246***
## (0.164)
##
## splines::ns(tmean, 4)3 22.753***
## (0.530)
##
## splines::ns(tmean, 4)4 46.087***
## (0.244)
##
## precip -0.293***
## (0.001)
##
## ---------------------------------------------------------------
## Observations 1,971,247
## R2 0.180
## Adjusted R2 0.179
## F Statistic 61,732.420*** (df = 7; 1967573)
## ===============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
As you can see that if we use county fixed effect model, we don’t
see
urban.2013.bin, now let us try state fiexed effect model.
fe_model <- plm(
omax2 ~ USDM.categorical +Urban.2013.bin +
mean_elev+ splines::ns(tmean, 4) + +precip + factor(wday) + factor(holiday) + I(long^2) + I(lat^2) + I(long*lat) +long + lat,
data=df,
index = c("state_fips", "date"),
model = "within",
effect = "twoways",
# vcov. = vcovHC, type = "HC1"
)
## Warning in pdata.frame(data, index): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
stargazer(fe_model,type = "text")
##
## ================================================================
## Dependent variable:
## --------------------------------
## omax2
## ----------------------------------------------------------------
## USDM.categoricalModerateDrought 1.644***
## (0.022)
##
## USDM.categoricalSevereDrought 2.920***
## (0.034)
##
## Urban.2013.binNonmetro -0.558***
## (0.020)
##
## mean_elev 0.011***
## (0.00004)
##
## splines::ns(tmean, 4)1 14.878***
## (0.254)
##
## splines::ns(tmean, 4)2 29.245***
## (0.167)
##
## splines::ns(tmean, 4)3 20.958***
## (0.543)
##
## splines::ns(tmean, 4)4 45.104***
## (0.237)
##
## precip -0.293***
## (0.001)
##
## I(long2) -0.021***
## (0.0002)
##
## I(lat2) -0.118***
## (0.001)
##
## I(long * lat) -0.007***
## (0.0004)
##
## long -3.309***
## (0.035)
##
## lat 9.063***
## (0.062)
##
## ----------------------------------------------------------------
## Observations 1,971,247
## R2 0.225
## Adjusted R2 0.223
## F Statistic 40,708.360*** (df = 14; 1968431)
## ================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
As you can see that if we use state fixed effect model, we can see
both
urban.2013.bin and elevation. elevation has positive association
while rural has negative association. which is quite understanable.