All installed packages were loaded to be used in the R Program.
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
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ 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 conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(plm)
##
## Attaching package: 'plm'
##
## The following objects are masked from 'package:dplyr':
##
## between, lag, lead
test.dat <- load("Ozone_Drought_Final.RData")
test.dat <- combinedAir.final
regions <- read.csv("region_code.csv")
test.dat$State.Code <- as.numeric(test.dat$State.Code)
test.dat$USDM.categorical <- factor(test.dat$USDM.categorical, levels = c("NoDrought", "ModerateDrought", "SevereDrought"))
test.dat$USDM.categorical <- relevel(test.dat$USDM.categorical, ref = "NoDrought")
test.dat$month = as.numeric(test.dat$month)
## Filter for summer period only
test.dat1 <- test.dat %>%
filter(dplyr::between(month, 5, 9))
md2_random<- plm(Max.Ozone ~ USDM.categorical +elevation+ Longitude + Latitude,
data=test.dat1,
index = c("GEOID", "Year","month"),
model = ("random"))
## 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")
md2_fix<- plm(Max.Ozone ~ USDM.categorical +elevation+ Longitude + Latitude,
data=test.dat1,
index = c("GEOID", "Year","month"),
effect = "twoway")
## 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(md2_random, md2_fix, type = "text", title = "All Regions", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## All Regions
## ==================================================================================
## Dependent Variable: Ozone Max
## --------------------------------------------------
## Random Effect Model Fixed Effect Model
## ----------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 2.313*** 2.036***
## (0.019) (0.019)
##
## USDM.categoricalSevereDrought 2.608*** 2.652***
## (0.027) (0.028)
##
## elevation 0.004*** 0.004***
## (0.0001) (0.0001)
##
## Longitude -0.009 -0.412***
## (0.013) (0.044)
##
## Latitude 0.299*** 1.157***
## (0.034) (0.061)
##
## Constant 31.497***
## (1.746)
##
## ----------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 3,437,966 3,437,966
## R2 0.010 0.007
## Adjusted R2 0.010 0.007
## F Statistic 25,216.030*** 4,731.343*** (df = 5; 3437021)
## ==================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
# Adding regions column
test.dat2 <- test.dat1 %>%
merge(regions, by = "State.Code")
ne <- test.dat2 %>%
filter(noaa_region == "northeast")
ne_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = ne,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
ne_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = ne,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(ne_random, ne_fix, type = "text", title = "Northeast", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Northeast
## ===============================================================================
## Dependent Variable: Ozone Max
## -----------------------------------------------
## Random Effect Model Fixed Effect Model
## -------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 2.142*** -0.407***
## (0.069) (0.074)
##
## USDM.categoricalSevereDrought 1.346*** -0.740**
## (0.298) (0.298)
##
## elevation 0.012*** 0.013***
## (0.0003) (0.0003)
##
## Longitude 0.654*** 3.458***
## (0.073) (0.195)
##
## Latitude -3.257*** -5.019***
## (0.123) (0.257)
##
## Constant 226.366***
## (9.237)
##
## -------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 519,759 519,759
## R2 0.013 0.005
## Adjusted R2 0.013 0.005
## F Statistic 3,283.838*** 526.668*** (df = 5; 519583)
## ===============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
nr <- test.dat2 %>%
filter(noaa_region == "northern_rockies")
nr_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = nr,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
nr_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = nr,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(nr_random, nr_fix, type = "text", title = "Northern Rockies", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Northern Rockies
## ==============================================================================
## Dependent Variable: Ozone Max
## ----------------------------------------------
## Random Effect Model Fixed Effect Model
## ------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 1.107*** 0.227***
## (0.070) (0.076)
##
## USDM.categoricalSevereDrought 1.594*** -0.343***
## (0.105) (0.116)
##
## elevation 0.002*** 0.002***
## (0.0003) (0.0003)
##
## Longitude -0.321*** -1.184***
## (0.064) (0.136)
##
## Latitude -0.825*** 0.149
## (0.114) (0.206)
##
## Constant 45.117***
## (7.443)
##
## ------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 101,196 101,196
## R2 0.040 0.002
## Adjusted R2 0.040 0.001
## F Statistic 647.878*** 36.421*** (df = 5; 101122)
## ==============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
nw <- test.dat2 %>%
filter(noaa_region == "northwest")
nw_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = nw,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
nw_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = nw,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(nw_random, nw_fix, type = "text", title = "Northwest", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Northwest
## ==============================================================================
## Dependent Variable: Ozone Max
## ----------------------------------------------
## Random Effect Model Fixed Effect Model
## ------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 0.375*** -0.328***
## (0.090) (0.118)
##
## USDM.categoricalSevereDrought 0.049 -2.239***
## (0.217) (0.247)
##
## elevation 0.006*** 0.007***
## (0.0002) (0.0003)
##
## Longitude 0.122 -0.942***
## (0.135) (0.225)
##
## Latitude -2.372*** -5.575***
## (0.227) (0.406)
##
## Constant 159.847***
## (17.309)
##
## ------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 79,052 79,052
## R2 0.040 0.020
## Adjusted R2 0.040 0.019
## F Statistic 1,431.794*** 317.315*** (df = 5; 78999)
## ==============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
ohv <- test.dat2 %>%
filter(noaa_region == "ohio_valley")
ohv_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = ohv,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
ohv_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = ohv,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(ohv_random, ohv_fix, type = "text", title = "Ohio Valley", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Ohio Valley
## ==============================================================================
## Dependent Variable: Ozone Max
## ----------------------------------------------
## Random Effect Model Fixed Effect Model
## ------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 3.845*** 0.745***
## (0.055) (0.057)
##
## USDM.categoricalSevereDrought 6.411*** 2.089***
## (0.124) (0.127)
##
## elevation -0.001*** -0.001***
## (0.0004) (0.0004)
##
## Longitude 0.258*** 2.509***
## (0.048) (0.239)
##
## Latitude -0.564*** -1.091***
## (0.086) (0.287)
##
## Constant 92.040***
## (5.664)
##
## ------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 628,182 628,182
## R2 0.029 0.001
## Adjusted R2 0.029 0.0005
## F Statistic 7,291.707*** 98.720*** (df = 5; 627976)
## ==============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
sth <- test.dat2 %>%
filter(noaa_region == "south")
sth_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = sth,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
sth_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = sth,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(sth_random, sth_fix, type = "text", title = "South", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## South
## ===============================================================================
## Dependent Variable: Ozone Max
## -----------------------------------------------
## Random Effect Model Fixed Effect Model
## -------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 2.715*** 1.467***
## (0.057) (0.061)
##
## USDM.categoricalSevereDrought 3.977*** 2.156***
## (0.078) (0.091)
##
## elevation 0.016*** 0.016***
## (0.001) (0.001)
##
## Longitude 0.199** -0.836***
## (0.101) (0.267)
##
## Latitude 1.331*** 6.426***
## (0.102) (0.311)
##
## Constant 15.543
## (10.119)
##
## -------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 413,667 413,667
## R2 0.024 0.005
## Adjusted R2 0.024 0.004
## F Statistic 4,634.372*** 376.171*** (df = 5; 413516)
## ===============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
se <- test.dat2 %>%
filter(noaa_region == "southeast")
se_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = se,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
se_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = se,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(se_random, se_fix, type = "text", title = "Southeast", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Southeast
## =================================================================================
## Dependent Variable: Ozone Max
## -------------------------------------------------
## Random Effect Model Fixed Effect Model
## ---------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 5.511*** 2.985***
## (0.049) (0.054)
##
## USDM.categoricalSevereDrought 8.146*** 4.472***
## (0.069) (0.076)
##
## elevation 0.006*** 0.010***
## (0.001) (0.001)
##
## Longitude -0.178** 0.683***
## (0.076) (0.236)
##
## Latitude 1.381*** 2.223***
## (0.075) (0.281)
##
## Constant -19.468**
## (7.611)
##
## ---------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 530,095 530,095
## R2 0.051 0.010
## Adjusted R2 0.051 0.010
## F Statistic 23,555.350*** 1,055.552*** (df = 5; 529907)
## =================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
sw <- test.dat2 %>%
filter(noaa_region == "southwest")
sw_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = sw,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
sw_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = sw,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(sw_random, sw_fix, type = "text", title = "Southwest", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Southwest
## ===============================================================================
## Dependent Variable: Ozone Max
## -----------------------------------------------
## Random Effect Model Fixed Effect Model
## -------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 0.885*** 0.629***
## (0.041) (0.050)
##
## USDM.categoricalSevereDrought 2.394*** 1.293***
## (0.047) (0.067)
##
## elevation 0.003*** 0.003***
## (0.0001) (0.0001)
##
## Longitude 1.121*** 2.690***
## (0.068) (0.092)
##
## Latitude 1.130*** 2.758***
## (0.074) (0.113)
##
## Constant 125.974***
## (8.046)
##
## -------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 325,309 325,309
## R2 0.044 0.012
## Adjusted R2 0.044 0.012
## F Statistic 4,669.218*** 776.754*** (df = 5; 325220)
## ===============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
um <- test.dat2 %>%
filter(noaa_region == "upper_midwest")
um_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = um,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
um_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = um,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(um_random, um_fix, type = "text", title = "Upper Midwest", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## Upper Midwest
## ==============================================================================
## Dependent Variable: Ozone Max
## ----------------------------------------------
## Random Effect Model Fixed Effect Model
## ------------------------------------------------------------------------------
## USDM.categoricalModerateDrought 1.182*** -0.545***
## (0.080) (0.083)
##
## USDM.categoricalSevereDrought 2.240*** -1.409***
## (0.248) (0.250)
##
## elevation 0.015*** 0.006**
## (0.002) (0.003)
##
## Longitude 0.589*** 2.650***
## (0.063) (0.430)
##
## Latitude -0.202 3.463***
## (0.126) (0.264)
##
## Constant 100.356***
## (7.172)
##
## ------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 265,885 265,885
## R2 0.010 0.001
## Adjusted R2 0.010 0.001
## F Statistic 406.934*** 66.289*** (df = 5; 265776)
## ==============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
wst <- test.dat2 %>%
filter(noaa_region == "west")
wst_random <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = wst,
index = c("GEOID", "Year", "month"),
model = ("random"))
## 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")
wst_fix <- plm(Max.Ozone ~ USDM.categorical + elevation + Longitude + Latitude,
data = wst,
index = c("GEOID", "Year", "month"),
effect = "twoway")
## 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(wst_random, wst_fix, type = "text", title = "West", align = TRUE,
add.lines = list(c("Model Name", "Random Effect Model", "Fixed Effect Model")),
column.labels = c("Random Effect Model", "Fixed Effect Model"),
dep.var.caption = "Dependent Variable: Ozone Max",
dep.var.labels.include = FALSE,
omit.table.layout = "#")
##
## West
## ===============================================================================
## Dependent Variable: Ozone Max
## -----------------------------------------------
## Random Effect Model Fixed Effect Model
## -------------------------------------------------------------------------------
## USDM.categoricalModerateDrought -0.185*** 0.476***
## (0.039) (0.059)
##
## USDM.categoricalSevereDrought -1.847*** -0.200**
## (0.047) (0.083)
##
## elevation 0.004*** 0.004***
## (0.0001) (0.0001)
##
## Longitude -1.036*** -1.098***
## (0.056) (0.056)
##
## Latitude 0.407*** 0.666***
## (0.090) (0.090)
##
## Constant -91.856***
## (7.236)
##
## -------------------------------------------------------------------------------
## Model Name Random Effect Model Fixed Effect Model
## Observations 566,959 566,959
## R2 0.011 0.008
## Adjusted R2 0.011 0.008
## F Statistic 5,934.351*** 940.701*** (df = 5; 566877)
## ===============================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01