Using Spatial Random Forest Modeling to Analyze Relationships Between Socio-Economic Change, Climatic Variability, and Population Loss in Puerto Rico
Jorge R. Soldevila Irizarry
2022-12-17
Read in libraries to be used for project
#For data handling
library(tidyverse) #Data handling
library(sf) #For handling simple features (shapefiles)
library(sp) #Handling shapefiles
library(raster) #For handling raster data
library(stars) #For handling raster data
library(terra) #for handling raster data
library(tibble) #To perform data manipulation
#For accessing US Census data
library(tidycensus) #For accessing US Census Data
census_api_key("355f6d8295ad12c6779c08234877da01cd8a94b6")
#For accessing climate data
#remotes::install_github("mikejohnson51/AOI")
#remotes::install_github("mikejohnson51/climateR")
library(AOI) #Needed for downloading climate data
library(climateR) #Used to download climate data
#For generating GFR
library(ranger) #For Spatial ML
#install.packages("SpatialML")
library(SpatialML) # To run random forest analysis
#For data visualization
library(rasterVis) #For visualizing raster data
library(leaflet) #For visualizing data
library(leafem)
library(RColorBrewer)Use tidycensus package to download socio-economic data:
#Decenial Census Populaiton Count 2010 and 2020
decenial_var_list_2010 <- load_variables(2010,"pl",cache = TRUE)
pop_2010 <- get_decennial(geography = "tract",
state = "PR",
year = 2010,
variables = c(pop_2010 = "P001001"),
output = "wide")
decenial_var_list_2020 <- load_variables(2020,"pl",cache = TRUE)
pop_2020 <- get_decennial(geography = "tract",
state = "PR",
year = 2020,
variables = c(pop_2020 = "P1_001N"),
output = "wide")#ACS 5-year Estimate variables for 2011, 2013, 2015, and 2017, and 2019
#Get List of variable names from ACS 5 year.
#Note: Some variables change name (code) after 2015
acs_var_list_2011 <- load_variables(2011,"acs5",cache = TRUE)
acs_var_list_2017 <- load_variables(2017,"acs5",cache = TRUE)
#Data
data_2011 <- get_acs(geography = "tract",
state = "PR",
year = 2011,
variables = c(mhhi_2011 = "B19013_001",
unemp_2011 = "DP03_0005E"),
output = "wide")
data_2013 <- get_acs(geography = "tract",
state = "PR",
year = 2013,
variables = c(mhhi_2013 = "B19013_001",
unemp_2013 = "DP03_0005E"),
output = "wide")
data_2015 <- get_acs(geography = "tract",
state = "PR",
year = 2015,
variables = c(mhhi_2015 = "B19013_001",
unemp_2015 = "DP03_0005E"),
output = "wide")
data_2017 <- get_acs(geography = "tract",
state = "PR",
year = 2017,
variables = c(mhhi_2017 = "B19013_001",
unemp_2017 = "DP03_0005E"),
output = "wide")
data_2019 <- get_acs(geography = "tract",
state = "PR",
year = 2019,
variables = c(mhhi_2019 = "B19013_001",
unemp_2019 = "DP03_0005E"),
output = "wide")Join data to one table and shape file
#JMake list of all socio-economic data frames
data_list <-list(pop_2010,data_2011,data_2013,data_2015,data_2017,data_2019)
#Join data
data <- data_list %>% reduce(left_join, by = "GEOID")
#Convert GEOID to numeric to join with shapefile
data$GEOID <- as.numeric(data$GEOID)
pop_2020$GEOID <- as.numeric(pop_2020$GEOID)
#remove variables not used
data <- data[,-c(2,4,6,8,9,11,13,14,16,18,19,21,23,24,26,28)]
pop_2020 <- pop_2020[,-2]
#Read in shapefile of Puerto Rico with Tracts in Vieques and Culebra removed
tract_2010 <- st_read("./Clipped_shp/Tract_2019_Cliped/Tract_2019.shp")## Reading layer `Tract_2019' from data source
## `/Users/jorgesoldevila/Library/CloudStorage/OneDrive-Hunter-CUNY/R Spatial Analysis/Final Project/Final Project/Clipped_shp/Tract_2019_Cliped/Tract_2019.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 900 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -67.27162 ymin: 17.9122 xmax: -65.58908 ymax: 18.51568
## Geodetic CRS: NAD83tract_2010$GEOID <- as.numeric(tract_2010$GEOID)
tract_2020 <- st_read("./Clipped_shp/Tract_2020_Cliped/Tract_2020.shp")## Reading layer `Tract_2020' from data source
## `/Users/jorgesoldevila/Library/CloudStorage/OneDrive-Hunter-CUNY/R Spatial Analysis/Final Project/Final Project/Clipped_shp/Tract_2020_Cliped/Tract_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 938 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -67.27162 ymin: 17.9122 xmax: -65.58908 ymax: 18.51568
## Geodetic CRS: NAD83tract_2020$GEOID <- as.numeric(tract_2020$GEOID)
block_2020 <- st_read("./Clipped_shp/Block_2020_Cliped/Block_2020.shp")## Reading layer `Block_2020' from data source
## `/Users/jorgesoldevila/Library/CloudStorage/OneDrive-Hunter-CUNY/R Spatial Analysis/Final Project/Final Project/Clipped_shp/Block_2020_Cliped/Block_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 41467 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -67.27162 ymin: 17.9122 xmax: -65.58908 ymax: 18.51568
## Geodetic CRS: NAD83block_2020$GEOID20 <- as.numeric(block_2020$GEOID20)
#Join data
tract_data_2010 <- left_join(tract_2010,data,by = "GEOID")
tract_data_2010 <- tract_data_2010 %>% st_drop_geometry()
tract_data_2020 <- left_join(tract_2020,pop_2020, by = "GEOID")
#Since Census Tracts are different in 2020 we perform a population weighted
#interpolation to convert 2020 population counto to 2010 boundaries
tract_2020_to_tract_2010 <- interpolate_pw(
from = tract_data_2020,
to = tract_2010,
to_id = "GEOID",
extensive = TRUE,
weights = tract_2020,
crs = 4269
)
#st_write(tract_2020_to_tract_2010,"tract_2020_toTract_2010.shp")
#tract_2020_to_tract_2010 <- st_read("./tract_2020_to_tract_2010.shp")
#Join 2010 and ACS data to 2020 population count in 2010 boundaries
tract_data <- left_join(tract_2020_to_tract_2010, tract_data_2010, by = "GEOID")
#Calculate percent change between 2010 and 2020
tract_data <- tract_data %>% mutate(pct_chng = (pop_2020-pop_2010)/pop_2010*100)
#Calculate binary variable of population loss and population gain
tract_data <- tract_data %>% mutate(pop_chng = as.factor((ifelse(pct_chng < 0,1,0))))
tmap::qtm(tract_data,fill = "pop_chng")
Get Climate Data and Calculate Mean
- Precipitation
precip_2011 = getTerraClim(tract_2010,param = "prcp",startDate = "2011-06-01",
endDate = "2011-11-30")
precip_2013 = getTerraClim(tract_2010,param = "prcp",startDate = "2013-06-01",
endDate = "2013-11-30")
precip_2015 = getTerraClim(tract_2010,param = "prcp",startDate = "2015-06-01",
endDate = "2015-11-30")
precip_2017 = getTerraClim(tract_2010,param = "prcp",startDate = "2017-06-01",
endDate = "2017-11-30")
precip_2019 = getTerraClim(tract_2010,param = "prcp",startDate = "2019-06-01",
endDate = "2019-11-30")
mean_precip_2011 <- mean(precip_2011$terraclim_prcp)
mean_precip_2013 <- mean(precip_2013$terraclim_prcp)
mean_precip_2015 <- mean(precip_2015$terraclim_prcp)
mean_precip_2017 <- mean(precip_2017$terraclim_prcp)
mean_precip_2019 <- mean(precip_2019$terraclim_prcp)- PDSI
pdsi_2011 = getTerraClim(tract_2010,param = "palmer",startDate = "2011-02-01",
endDate = "2011-08-31")
pdsi_2013 = getTerraClim(tract_2010,param = "palmer",startDate = "2013-06-01",
endDate = "2013-08-31")
pdsi_2015 = getTerraClim(tract_2010,param = "palmer",startDate = "2015-06-01",
endDate = "2015-08-31")
pdsi_2017 = getTerraClim(tract_2010,param = "palmer",startDate = "2017-06-01",
endDate = "2017-08-31")
pdsi_2019 = getTerraClim(tract_2010,param = "palmer",startDate = "2019-06-01",
endDate = "2019-08-31")
mean_pdsi_2011 <- mean(pdsi_2011$terraclim_palmer)
mean_pdsi_2013 <- mean(pdsi_2013$terraclim_palmer)
mean_pdsi_2015 <- mean(pdsi_2015$terraclim_palmer)
mean_pdsi_2017 <- mean(pdsi_2017$terraclim_palmer)
mean_pdsi_2019 <- mean(pdsi_2019$terraclim_palmer)Max Temperature
tmax_2011 = getTerraClim(tract_2010,param = "tmax",startDate = "2011-01-01",
endDate = "2011-12-31")
tmax_2013 = getTerraClim(tract_2010,param = "tmax",startDate = "2013-01-01",
endDate = "2013-12-31")
tmax_2015 = getTerraClim(tract_2010,param = "tmax",startDate = "2015-01-01",
endDate = "2015-12-31")
mean_tmax_2011 <- mean(tmax_2011$terraclim_tmax)
mean_tmax_2013 <- mean(tmax_2013$terraclim_tmax)
mean_tmax_2015 <- mean(tmax_2015$terraclim_tmax)2017 September Precipitation (Hurricanes Irma and Maria)
hurric_2017 = getTerraClim(tract_2010,param = "prcp",startDate = "2017-09-01", endDate = "2017-09-30")
mean_hurric_2017 <- mean(hurric_2017$terraclim_prcp)View climate data
tmap::qtm(mean_pdsi_2011)
tmap::qtm(mean_tmax_2011)
tmap::qtm(mean_precip_2017)
tmap::qtm(mean_hurric_2017)
Convert environmental data to polygon
extent <- raster::extent(tract_data)
#Crop is used for cropping the raster image to the extent of the desired shapefile
#Mask is used to only show in raster values contained within that desired shapefile
mean_pdsi_2011 <- raster::crop(x = mean_pdsi_2011, y = extent)
mean_pdsi_2011 <- raster::mask(x = mean_pdsi_2011, mask = tract_data)
mean_pdsi_2013 <- raster::crop(x = mean_pdsi_2013, y = extent)
mean_pdsi_2013 <- raster::mask(x = mean_pdsi_2013, mask = tract_data)
mean_pdsi_2015 <- raster::crop(x = mean_pdsi_2015, y = extent)
mean_pdsi_2015 <- raster::mask(x = mean_pdsi_2015, mask = tract_data)
mean_pdsi_2017 <- raster::crop(x = mean_pdsi_2017, y = extent)
mean_pdsi_2017 <- raster::mask(x = mean_pdsi_2017, mask = tract_data)
mean_pdsi_2019 <- raster::crop(x = mean_pdsi_2019, y = extent)
mean_pdsi_2019 <- raster::mask(x = mean_pdsi_2019, mask = tract_data)
mean_precip_2011 <- raster::crop(x = mean_precip_2011, y = extent)
mean_precip_2011 <- raster::mask(x = mean_precip_2011, mask = tract_data)
mean_precip_2013 <- raster::crop(x = mean_precip_2013, y = extent)
mean_precip_2013 <- raster::mask(x = mean_precip_2013, mask = tract_data)
mean_precip_2015 <- raster::crop(x = mean_precip_2015, y = extent)
mean_precip_2015 <- raster::mask(x = mean_precip_2015, mask = tract_data)
mean_precip_2017 <- raster::crop(x = mean_precip_2017, y = extent)
mean_precip_2017 <- raster::mask(x = mean_precip_2017, mask = tract_data)
mean_precip_2019 <- raster::crop(x = mean_precip_2019, y = extent)
mean_precip_2019 <- raster::mask(x = mean_precip_2019, mask = tract_data)
mean_tmax_2011 <- raster::crop(x = mean_tmax_2011, y = extent)
mean_tmax_2011 <- raster::mask(x = mean_tmax_2011, mask = tract_data)
mean_tmax_2013 <- raster::crop(x = mean_tmax_2013, y = extent)
mean_tmax_2013 <- raster::mask(x = mean_tmax_2013, mask = tract_data)
mean_tmax_2015 <- raster::crop(x = mean_tmax_2015, y = extent)
mean_tmax_2015 <- raster::mask(x = mean_tmax_2015, mask = tract_data)
mean_hurric_2017 <- raster::crop(x = mean_hurric_2017, y = extent)
mean_hurric_2017 <- raster::mask(x = mean_hurric_2017, mask = tract_data)
#Convert data to raster and vector objects for raster package
vect(tract_data)
rast(mean_pdsi_2011)
rast(mean_pdsi_2013)
rast(mean_pdsi_2015)
rast(mean_pdsi_2017)
rast(mean_pdsi_2019)
rast(mean_precip_2011)
rast(mean_precip_2013)
rast(mean_precip_2015)
rast(mean_precip_2017)
rast(mean_precip_2019)
rast(mean_tmax_2011)
rast(mean_tmax_2013)
rast(mean_tmax_2015)
rast(mean_hurric_2017)#Transform tract_data to raster coordinate system
crs <- crs(mean_pdsi_2011)
st_crs(tract_data) <- crs
st_transform(tract_data)#add raster data to tract_data
tract_data$pdsi2011 <- extract(mean_pdsi_2011,tract_data,fun = mean, na.rm = TRUE)
tract_data$pdsi2013 <- extract(mean_pdsi_2013,tract_data,fun = mean, na.rm = TRUE)
tract_data$pdsi2015 <- extract(mean_pdsi_2015,tract_data,fun = mean, na.rm = TRUE)
tract_data$pdsi2017 <- extract(mean_pdsi_2017,tract_data,fun = mean, na.rm = TRUE)
tract_data$pdsi2019 <- extract(mean_pdsi_2019,tract_data,fun = mean, na.rm = TRUE)
tract_data$precip2011 <- extract(mean_precip_2011,tract_data,fun = mean, na.rm = TRUE)
tract_data$precip2013 <- extract(mean_precip_2013,tract_data,fun = mean, na.rm = TRUE)
tract_data$precip2015 <- extract(mean_precip_2015,tract_data,fun = mean, na.rm = TRUE)
tract_data$precip2017 <- extract(mean_precip_2017,tract_data,fun = mean, na.rm = TRUE)
tract_data$precip2019 <- extract(mean_precip_2019,tract_data,fun = mean, na.rm = TRUE)
tract_data$tmax2011 <- extract(mean_tmax_2011,tract_data,fun = mean, na.rm = TRUE)
tract_data$tmax2013 <- extract(mean_tmax_2013,tract_data,fun = mean, na.rm = TRUE)
tract_data$tmax2015 <- extract(mean_tmax_2015,tract_data,fun = mean, na.rm = TRUE)
tract_data$hurric_2017 <- extract(mean_hurric_2017,tract_data, fun = mean,
na.rm = TRUE)Final data preparation
tract_data <- tract_data %>% filter(pct_chng <=0)
tract_data <- na.omit(tract_data)
tract_data <- tract_data[,c("GEOID","pop_2020","pop_2010","mhhi_2011E",
"unemp_2011","mhhi_2013E","unemp_2013",
"mhhi_2015E", "unemp_2015", "mhhi_2017E","unemp_2017",
"mhhi_2019E", "unemp_2019", "pct_chng", "pdsi2011",
"pdsi2013", "pdsi2015", "pdsi2017", "pdsi2019",
"precip2011", "precip2013", "precip2015", "precip2017",
"precip2019", "tmax2011", "tmax2013","tmax2015","hurric_2017")]
tract_data <- tract_data %>%
mutate(lon = map_dbl(geometry, ~st_point_on_surface(.x)[[1]]),
lat = map_dbl(geometry, ~st_point_on_surface(.x)[[2]]))
tract_data <- st_drop_geometry(tract_data)
final_data <- as.data.frame(tract_data[,c("GEOID","pop_2020","pop_2010","mhhi_2011E"
,"unemp_2011","mhhi_2013E","unemp_2013",
"mhhi_2015E", "unemp_2015", "mhhi_2017E","unemp_2017",
"mhhi_2019E", "unemp_2019", "pct_chng", "pdsi2011",
"pdsi2013", "pdsi2015", "pdsi2017", "pdsi2019",
"precip2011", "precip2013", "precip2015", "precip2017",
"precip2019","tmax2011","tmax2013","tmax2015","hurric_2017","lon","lat")])
coord <- tract_data[,29:30]model_1 <- grf(pct_chng ~ mhhi_2011E + unemp_2011 + mhhi_2013E + unemp_2013 +
mhhi_2015E + unemp_2015 + mhhi_2017E + unemp_2017 + mhhi_2019E + unemp_2019,
dframe = final_data, kernel = "adaptive", bw = 50, coords = coord,
ntree = 500, forests = TRUE,weighted = TRUE, print.results = TRUE)## Ranger result
##
## Call:
## ranger(pct_chng ~ mhhi_2011E + unemp_2011 + mhhi_2013E + unemp_2013 + mhhi_2015E + unemp_2015 + mhhi_2017E + unemp_2017 + mhhi_2019E + unemp_2019, data = final_data, num.trees = 500, mtry = 3, importance = "impurity", num.threads = NULL)
##
## Type: Regression
## Number of trees: 500
## Sample size: 794
## Number of independent variables: 10
## Mtry: 3
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 68.27945
## R squared (OOB): 0.1385501
## mhhi_2011E unemp_2011 mhhi_2013E unemp_2013 mhhi_2015E unemp_2015 mhhi_2017E
## 6418.701 5898.794 5436.674 4792.727 5830.489 4665.421 7021.067
## unemp_2017 mhhi_2019E unemp_2019
## 5819.853 7450.925 5168.796
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -68.1636 -3.8915 0.3076 0.0972 4.1026 23.1434
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -11.32901 -0.33127 0.02076 -0.02291 0.34664 10.01526
## Min Max Mean StD
## mhhi_2011E 55.05680 1686.2982 401.2643 331.20465
## unemp_2011 45.62011 5606.0991 290.9673 539.84462
## mhhi_2013E 44.13313 6155.7544 406.5567 578.94612
## unemp_2013 36.76247 563.2976 188.0724 92.08935
## mhhi_2015E 50.13966 2013.0616 374.4276 329.81711
## unemp_2015 40.67603 1021.5672 216.5849 141.10330
## mhhi_2017E 60.12964 1983.0969 356.7437 301.08880
## unemp_2017 36.80501 1632.7242 232.2983 212.08641
## mhhi_2019E 49.78930 4540.9704 381.0616 389.22206
## unemp_2019 33.23461 623.6220 193.7760 100.73537model_2 <- grf(pct_chng ~ pdsi2011 + pdsi2013 + pdsi2015 + pdsi2017 + pdsi2019 +
precip2011 + precip2013 + precip2015 + precip2017 + precip2019 +
tmax2011 + tmax2013 + tmax2015 + hurric_2017, dframe = final_data,
kernel = "adaptive", bw = 300, coords = coord, ntree = 500,
forests = TRUE,weighted = TRUE, print.results = TRUE )## Ranger result
##
## Call:
## ranger(pct_chng ~ pdsi2011 + pdsi2013 + pdsi2015 + pdsi2017 + pdsi2019 + precip2011 + precip2013 + precip2015 + precip2017 + precip2019 + tmax2011 + tmax2013 + tmax2015 + hurric_2017, data = final_data, num.trees = 500, mtry = 4, importance = "impurity", num.threads = NULL)
##
## Type: Regression
## Number of trees: 500
## Sample size: 794
## Number of independent variables: 14
## Mtry: 4
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 73.54891
## R squared (OOB): 0.07206781
## pdsi2011 pdsi2013 pdsi2015 pdsi2017 pdsi2019 precip2011
## 2602.611 3193.649 5013.896 3029.327 3118.222 2778.165
## precip2013 precip2015 precip2017 precip2019 tmax2011 tmax2013
## 2497.355 2551.403 2725.805 4369.342 2372.978 2336.942
## tmax2015 hurric_2017
## 2426.638 3219.084
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -81.8945 -4.3438 0.5496 0.1914 5.5799 33.0128
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -65.41731 -0.78840 0.02089 0.01224 0.97270 26.26244
## Min Max Mean StD
## pdsi2011 613.6638 2183.294 1297.3773 291.3192
## pdsi2013 630.3877 2348.336 1314.4391 315.0751
## pdsi2015 568.8724 2157.663 1159.3028 426.7282
## pdsi2017 431.2770 2122.429 1072.4253 438.5086
## pdsi2019 544.5322 1795.233 968.6996 228.9808
## precip2011 548.6847 1733.005 1018.4666 300.6989
## precip2013 492.5809 1888.988 1073.4853 332.5054
## precip2015 448.7131 1734.696 999.1144 285.7203
## precip2017 496.7335 1649.667 1008.5661 276.8785
## precip2019 520.2184 2408.262 1058.3975 438.3272
## tmax2011 551.6534 1624.981 920.5630 246.4061
## tmax2013 450.8298 1466.220 896.6409 226.4486
## tmax2015 465.6939 1481.301 873.6982 218.1252
## hurric_2017 533.7712 1991.707 1114.7938 304.3689model_3 <- grf(pct_chng ~ mhhi_2011E + unemp_2011 + mhhi_2013E + unemp_2013 +
mhhi_2015E + unemp_2015 + mhhi_2017E + unemp_2017 + mhhi_2019E+ unemp_2019 +
pdsi2011 + pdsi2013 + pdsi2015 + pdsi2017 + pdsi2019 +precip2011 +
precip2013 + precip2015 + precip2017 + precip2019 + tmax2011 +
tmax2013 + tmax2015 + hurric_2017, dframe = final_data, kernel = "adaptive",
bw = 300, coords = coord, ntree = 500, forests = TRUE,
weighted = TRUE, print.results = TRUE )## Ranger result
##
## Call:
## ranger(pct_chng ~ mhhi_2011E + unemp_2011 + mhhi_2013E + unemp_2013 + mhhi_2015E + unemp_2015 + mhhi_2017E + unemp_2017 + mhhi_2019E + unemp_2019 + pdsi2011 + pdsi2013 + pdsi2015 + pdsi2017 + pdsi2019 + precip2011 + precip2013 + precip2015 + precip2017 + precip2019 + tmax2011 + tmax2013 + tmax2015 + hurric_2017, data = final_data, num.trees = 500, mtry = 8, importance = "impurity", num.threads = NULL)
##
## Type: Regression
## Number of trees: 500
## Sample size: 794
## Number of independent variables: 24
## Mtry: 8
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 56.06494
## R squared (OOB): 0.2926549
## mhhi_2011E unemp_2011 mhhi_2013E unemp_2013 mhhi_2015E unemp_2015
## 4715.581 3869.675 3705.441 2052.175 3910.160 2096.336
## mhhi_2017E unemp_2017 mhhi_2019E unemp_2019 pdsi2011 pdsi2013
## 4470.266 2309.658 5464.179 2373.345 1231.115 1561.771
## pdsi2015 pdsi2017 pdsi2019 precip2011 precip2013 precip2015
## 4494.475 1568.993 2722.453 1126.758 1019.924 1074.900
## precip2017 precip2019 tmax2011 tmax2013 tmax2015 hurric_2017
## 1152.236 2842.207 1624.731 1400.453 1384.661 1772.219
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -72.43640 -4.02217 0.39931 0.05351 4.50499 23.27514
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -19.05282 -0.30952 0.02615 -0.01115 0.35619 5.24102
## Min Max Mean StD
## mhhi_2011E 481.0269 4452.0325 2258.9430 967.6143
## unemp_2011 313.7176 7146.2668 2113.8301 2161.8070
## mhhi_2013E 527.1469 3842.5265 1763.3004 741.5167
## unemp_2013 321.5011 1790.0603 780.8321 332.3370
## mhhi_2015E 729.0500 3866.0472 1975.2798 768.1148
## unemp_2015 387.2057 1776.6897 837.2029 328.0439
## mhhi_2017E 545.0978 5424.5185 2142.8576 1201.2484
## unemp_2017 324.3331 1416.9824 793.7346 272.7082
## mhhi_2019E 522.5476 3362.7540 2155.7945 738.4731
## unemp_2019 345.5421 1621.8907 832.3009 284.9149
## pdsi2011 255.2240 1584.0151 627.9812 269.0752
## pdsi2013 283.9714 2126.0351 740.1629 371.3996
## pdsi2015 255.0205 1522.5467 601.8779 328.9413
## pdsi2017 200.2808 1398.0285 531.6237 239.0186
## pdsi2019 214.0473 1330.0435 470.0290 163.8741
## precip2011 203.7933 1076.7244 449.8690 153.6122
## precip2013 185.0954 1222.8272 481.7838 192.2166
## precip2015 195.1382 927.5143 445.1106 163.1327
## precip2017 209.7860 935.7813 466.6188 155.5405
## precip2019 243.5459 1483.0576 532.4752 277.2578
## tmax2011 215.6266 1000.6003 430.3160 145.3967
## tmax2013 221.8357 906.3878 418.9330 127.0680
## tmax2015 221.0772 869.0266 420.6154 112.2209
## hurric_2017 227.1167 1240.5800 578.2303 213.1713