assignment_6

Author

Sam Kuhn 

## Load data
#| echo: false
#| message: false
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(haven)
library(gt)
library(modelsummary)
library(rddensity)
library(rdrobust)

meyersson_df <- haven::read_dta(here::here("data", "Meyersson2014.dta"))
superfund_df <- haven::read_dta(here::here("data", "superfund.dta"))

Part 1

1. Create a scatterplot that shows a relationship between Y and X with a linear fitted line.

meyersson_df |>
  ggplot(aes(X, Y)) +
  geom_point() +
  geom_smooth(method = "lm", linetype = "dashed") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

2. RDD plot

rdplot(meyersson_df$Y, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

3. Estimation using OLS

ols_mod <- lm(Y ~ T + X + X * T, data = meyersson_df)
modelsummary(ols_mod)
 (1)
(Intercept) 16.199
(0.424)
T −0.554
(0.856)
X −0.012
(0.011)
T × X −0.144
(0.051)
Num.Obs. 2629
R2 0.012
R2 Adj. 0.011
AIC 19322.1
BIC 19351.5
Log.Lik. −9656.041
F 10.416
RMSE 9.53

4. Estimation using local polynomials

Uniform kernel

poly_mod <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "uniform", p = 1, h = 20)
summary(poly_mod)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                      Manual
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             608          280
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  20.000       20.000
BW bias (b)                  20.000       20.000
rho (h/b)                     1.000        1.000
Unique Obs.                    2314          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     2.927     1.235     2.371     0.018     [0.507 , 5.347]     
        Robust         -         -     1.636     0.102    [-0.582 , 6.471]     
=============================================================================

Triangular kernel

triangle_mod <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "triangular", p = 1, h = 20)
summary(triangle_mod)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                      Manual
Kernel                   Triangular
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             608          280
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  20.000       20.000
BW bias (b)                  20.000       20.000
rho (h/b)                     1.000        1.000
Unique Obs.                    2314          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     2.937     1.343     2.187     0.029     [0.305 , 5.569]     
        Robust         -         -     1.379     0.168    [-1.117 , 6.414]     
=============================================================================
triangle_mod_2 <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "triangular", p = 2, h = 20)
summary(triangle_mod_2)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                      Manual
Kernel                   Triangular
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             608          280
Order est. (p)                    2            2
Order bias  (q)                   3            3
BW est. (h)                  20.000       20.000
BW bias (b)                  20.000       20.000
rho (h/b)                     1.000        1.000
Unique Obs.                    2314          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     2.649     1.921     1.379     0.168    [-1.117 , 6.414]     
        Robust         -         -     0.420     0.674    [-3.969 , 6.135]     
=============================================================================

5. Automatic bandwith selection

automatic_bandwidth_mod <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "uniform", bwselect = "mserd")
summary(automatic_bandwidth_mod)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             474          251
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  15.449       15.449
BW bias (b)                  28.309       28.309
rho (h/b)                     0.546        0.546
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     3.202     1.357     2.360     0.018     [0.542 , 5.861]     
        Robust         -         -     2.044     0.041     [0.134 , 6.339]     
=============================================================================
  1. MSE-optimal bandwidth and triangular kernel
automatic_bandwidth_tri_mod <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "triangular", bwselect = "mserd")
summary(automatic_bandwidth_tri_mod)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             529          266
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  17.240       17.240
BW bias (b)                  28.576       28.576
rho (h/b)                     0.603        0.603
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     3.020     1.427     2.116     0.034     [0.223 , 5.816]     
        Robust         -         -     1.776     0.076    [-0.309 , 6.276]     
=============================================================================
  1. MSE-optimal bandwidth, triangular kernel, and p = 2.
automatic_bandwidth_tri_mod_2 <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "triangular", bwselect = "mserd", p = 2)
summary(automatic_bandwidth_tri_mod_2)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             702          291
Order est. (p)                    2            2
Order bias  (q)                   3            3
BW est. (h)                  23.121       23.121
BW bias (b)                  35.191       35.191
rho (h/b)                     0.657        0.657
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     2.772     1.808     1.533     0.125    [-0.772 , 6.315]     
        Robust         -         -     1.325     0.185    [-1.276 , 6.600]     
=============================================================================
  1. Two different MSE-optimal bandwidth (msetwo), uniform kernel, and p = 1
automatic_bandwidth_uniform <- rdrobust(meyersson_df$Y, meyersson_df$X, kernel = "uniform", bwselect = "msetwo", p = 1)
summary(automatic_bandwidth_uniform)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                      msetwo
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             533          242
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  17.370       14.890
BW bias (b)                  36.306       26.804
rho (h/b)                     0.478        0.555
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     2.912     1.355     2.149     0.032     [0.256 , 5.567]     
        Robust         -         -     1.971     0.049     [0.016 , 6.071]     
=============================================================================
  1. MSE-optimal bandwidth, uniform kernel, p=1, and with controls of vshr_islam1994 partycount lpop1994 merkezi merkezp subbuyuk buyuk
covariates_mod <- rdrobust(
  meyersson_df$Y,
  meyersson_df$X,
  kernel = "uniform",
  p = 1,
  bwselect = "msetwo", 
  covs = meyersson_df$vshr_islam1994 + meyersson_df$lpop1994 + meyersson_df$merkezi + meyersson_df$merkezp +
    meyersson_df$subbuyuk + meyersson_df$buyuk
)
summary(covariates_mod)
Covariate-adjusted Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                      msetwo
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             534          248
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  17.442       15.017
BW bias (b)                  32.838       26.905
rho (h/b)                     0.531        0.558
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     3.228     1.293     2.496     0.013     [0.694 , 5.762]     
        Robust         -         -     2.181     0.029     [0.331 , 6.209]     
=============================================================================
  1. The same model as iv) with the clustered standard error at the province level
covariates_mod <- rdrobust(
  meyersson_df$Y,
  meyersson_df$X,
  kernel = "uniform",
  p = 1,
  covs = meyersson_df$vshr_islam1994 + meyersson_df$lpop1994 + meyersson_df$merkezi + meyersson_df$merkezp +
    meyersson_df$subbuyuk + meyersson_df$buyuk,
  cluster = meyersson_df$prov_num
)
summary(covariates_mod)
Covariate-adjusted Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             529          266
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  17.230       17.230
BW bias (b)                  37.748       37.748
rho (h/b)                     0.456        0.456
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     3.068     1.436     2.136     0.033     [0.252 , 5.883]     
        Robust         -         -     2.132     0.033     [0.269 , 6.377]     
=============================================================================

Validation tests

RD plots

rdplot(meyersson_df$vshr_islam1994, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

rdplot(meyersson_df$partycount, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

rdplot(meyersson_df$lpop1994, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

rdplot(meyersson_df$merkezi, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

rdplot(meyersson_df$merkezp, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

rdplot(meyersson_df$subbuyuk, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

rdplot(meyersson_df$buyuk, meyersson_df$X, c = 0, p = 4, nbins = c(20, 20))

RD tests

validation_1 <- rdrobust(meyersson_df$vshr_islam1994, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_1)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             404          233
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  13.428       13.428
BW bias (b)                  24.567       24.567
rho (h/b)                     0.547        0.547
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.171     1.390     0.123     0.902    [-2.553 , 2.896]     
        Robust         -         -     0.117     0.907    [-2.994 , 3.376]     
=============================================================================
validation_2 <- rdrobust(meyersson_df$partycount, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_2)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             346          209
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  11.113       11.113
BW bias (b)                  21.590       21.590
rho (h/b)                     0.515        0.515
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.008     0.451    -0.018     0.986    [-0.893 , 0.877]     
        Robust         -         -    -0.212     0.832    [-1.130 , 0.909]     
=============================================================================
validation_3 <- rdrobust(meyersson_df$lpop1994, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_3)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             451          242
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  14.574       14.574
BW bias (b)                  26.531       26.531
rho (h/b)                     0.549        0.549
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.064     0.246     0.261     0.794    [-0.418 , 0.546]     
        Robust         -         -     0.266     0.790    [-0.482 , 0.634]     
=============================================================================
validation_4 <- rdrobust(meyersson_df$merkezi, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_4)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             584          276
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  19.018       19.018
BW bias (b)                  39.245       39.245
rho (h/b)                     0.485        0.485
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.045     0.069    -0.650     0.516    [-0.180 , 0.090]     
        Robust         -         -    -0.654     0.513    [-0.207 , 0.103]     
=============================================================================
validation_5 <- rdrobust(meyersson_df$merkezp, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_5)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             295          192
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   9.934        9.934
BW bias (b)                  20.838       20.838
rho (h/b)                     0.477        0.477
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.038     0.036     1.047     0.295    [-0.033 , 0.109]     
        Robust         -         -     1.120     0.263    [-0.034 , 0.126]     
=============================================================================
validation_6 <- rdrobust(meyersson_df$subbuyuk, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_6)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             368          220
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  11.855       11.855
BW bias (b)                  23.709       23.709
rho (h/b)                     0.500        0.500
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.006     0.035     0.156     0.876    [-0.064 , 0.075]     
        Robust         -         -    -0.102     0.919    [-0.081 , 0.073]     
=============================================================================
validation_7 <- rdrobust(meyersson_df$buyuk, meyersson_df$X, kernel = "uniform", bwselect = "mserd", p = 1)
summary(validation_7)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 2629
BW type                       mserd
Kernel                      Uniform
VCE method                       NN

Number of Obs.                 2314          315
Eff. Number of Obs.             630          284
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  20.791       20.791
BW bias (b)                  44.189       44.189
rho (h/b)                     0.470        0.470
Unique Obs.                    2311          315

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.006     0.018     0.347     0.729    [-0.029 , 0.041]     
        Robust         -         -     0.205     0.837    [-0.037 , 0.045]     
=============================================================================

Part 2

  1. Estimate the effect of the NPL status by 2000 on the house price using the OLS
vars <- c(
  "lnmeanhs8", "tothsun8", "occupied80", "ownocc8", "bedrms0_80occ", "bedrms1_80occ", "bedrms2_80occ",
  "bedrms3_80occ", "bedrms4_80occ", "bedrms5_80occ", "detach80occ", "attach80occ", "mobile80occ", "blt0_1yrs80occ", "blt2_5yrs80occ",
  "blt6_10yrs80occ", "blt10_20yrs80occ", "blt20_30yrs80occ", "blt30_40yrs80occ", "blt40_yrs80occ", "nofullkitchen80",
  "firestoveheat80", "noaircond80", "zerofullbath80", "avhhin8", "povrat8", "unemprt8", "welfare8", "pop_den8", "shrblk8",
  "shrhsp8", "child8", "old8", "shrfor8", "ffh8", "smhse8", "hsdrop8", "no_hs_diploma8", "ba_or_better8"
)

ols_mod_part2 <- lm(reformulate(vars, "lnmdvalhs0"), data = superfund_df)
modelsummary(ols_mod_part2, stars = TRUE)
 (1)
(Intercept) −1086459.718
(945807.621)
lnmeanhs8 0.762***
(0.054)
tothsun8 0.000
(0.000)
occupied80 0.804*
(0.346)
ownocc8 0.000*
(0.000)
bedrms0_80occ 776748.135
(682871.038)
bedrms1_80occ 776753.732
(682871.126)
bedrms2_80occ 776753.574
(682871.140)
bedrms3_80occ 776753.277
(682871.135)
bedrms4_80occ 776753.292
(682871.154)
bedrms5_80occ 776753.577
(682871.114)
detach80occ 309708.734
(709363.826)
attach80occ 309708.404
(709363.827)
mobile80occ 309708.834
(709363.826)
blt0_1yrs80occ 0.306
(0.511)
blt2_5yrs80occ −0.270
(0.213)
blt6_10yrs80occ 0.089
(0.197)
blt10_20yrs80occ −0.126
(0.125)
blt20_30yrs80occ −0.291**
(0.112)
blt30_40yrs80occ −0.454*
(0.202)
nofullkitchen80 −0.577
(0.827)
firestoveheat80 0.151
(0.186)
noaircond80 0.211***
(0.060)
zerofullbath80 1.064
(0.692)
avhhin8 0.000***
(0.000)
povrat8 −0.143
(0.336)
unemprt8 −0.966**
(0.354)
welfare8 0.722*
(0.334)
pop_den8 0.000
(0.000)
shrblk8 −0.387***
(0.101)
shrhsp8 −0.192
(0.185)
child8 −0.772*
(0.386)
old8 −1.068**
(0.377)
shrfor8 0.761**
(0.276)
ffh8 0.748***
(0.177)
smhse8 0.423**
(0.142)
hsdrop8 0.123
(0.155)
no_hs_diploma8 −0.291
(0.193)
ba_or_better8 −0.067
(0.284)
Num.Obs. 487
R2 0.779
R2 Adj. 0.760
AIC 42.7
BIC 210.2
Log.Lik. 18.655
RMSE 0.23
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
  1. Estimate the effect of the NPL status by 2000 on the house price using the instrumental variable approach
iv_mod <- rdrobust(superfund_df$npl2000, superfund_df$hrs_82,
  c = 28.5,
  covs = superfund_df$lnmeanhs8 + superfund_df$tothsun8 + superfund_df$occupied80 + superfund_df$ownocc8 +
    superfund_df$bedrms0_80occ + superfund_df$bedrms1_80occ + superfund_df$bedrms2_80occ + superfund_df$bedrms3_80occ +
    superfund_df$bedrms4_80occ + superfund_df$bedrms5_80occ + superfund_df$detach80occ + superfund_df$attach80occ +
    superfund_df$mobile80occ + superfund_df$blt0_1yrs80occ + superfund_df$blt2_5yrs80occ + superfund_df$blt6_10yrs80occ +
    superfund_df$blt10_20yrs80occ + superfund_df$blt20_30yrs80occ + superfund_df$blt30_40yrs80occ + superfund_df$blt40_yrs80occ +
    superfund_df$nofullkitchen80 + superfund_df$firestoveheat80 + superfund_df$noaircond80 + superfund_df$zerofullbath80 +
    superfund_df$avhhin8 + superfund_df$povrat8 + superfund_df$unemprt8 + superfund_df$welfare8 +
    superfund_df$pop_den8 + superfund_df$shrblk8 + superfund_df$shrhsp8 + superfund_df$child8 +
    superfund_df$old8 + superfund_df$shrfor8 + superfund_df$ffh8 + superfund_df$smhse8 +
    superfund_df$hsdrop8 + superfund_df$no_hs_diploma8 + superfund_df$ba_or_better8 
)
summary(iv_mod)
Covariate-adjusted Sharp RD estimates using local polynomial regression.

Number of Obs.                  487
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                  181          306
Eff. Number of Obs.              74          115
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  10.243       10.243
BW bias (b)                  16.211       16.211
rho (h/b)                     0.632        0.632
Unique Obs.                     168          263

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional     0.604     0.130     4.662     0.000     [0.350 , 0.858]     
        Robust         -         -     3.780     0.000     [0.284 , 0.896]     
=============================================================================
superfund_df <- superfund_df |> 
  filter(abs(hrs_82 <= 40.5) & abs(hrs_82 >= 16.5))

controls <- superfund_df |> 
  select(any_of(vars), statefips) |> 
  mutate(statefips = factor(statefips))

controls_matrix <- model.matrix(~., data = controls)

5. Create the following figure with the probability of placement on NPL by 2000 on the y-axis and the 1982 HRS score on the x-axis

superfund_df |> 
  group_by(hrsgroup) |> 
  mutate(mean_npl = mean(npl2000)) |> 
  ggplot(aes(hrsgroup, mean_npl)) +
  geom_point() +
  theme_minimal()

6. bar chart

superfund_df |> 
  ggplot(aes(hrsgroup, npl2000)) +
  geom_bar(stat = "identity") +
  theme_minimal()

superfund_df |> 
  group_by(hrsgroup) |> 
  mutate(mean_value = mean(lnmdvalhs0)) |> 
  ggplot(aes(hrsgroup, mean_value)) +
  geom_point() +
  theme_minimal()