Spatially Distorted Signaling: How Opinions Against Wind Infrastructure Delay Our Transition to Renewable Energy
Multivariate Regression Analysis
Geospatial Analysis
MEDS
Decarbonization Transition
Author
Sofia Ingersoll
Published
December 14, 2023
Figure 1: Wind turbines on the Bishop Hill wind farm operate among the corn and soybean fields near Bishop Hill, Illinois (U.S. Department of Agriculture, 2017).
In hopes to understand potential delays in our transition to a more sustainable future this study analyzed the activity of U.S. wind power plants (2000-2016) to inspect how population density, median income, and anti-wind opinions impacted wind plant operations. An investigation was conducted to quantify propensities in relation to local resistance against renewable energy developments for the listed attributes.
Spatial Distorted Signalling (SDS) is a term that describes the mobilization of minority opinion holders to promote legislation that aligns with their beliefs and electorally push-back (Stokes, et al, 2016). The SDS phenomenon has the ability to skew the accurate representation of public opinions due to effective organization practices. Leading to a disproportionate number of minority opinion holder interactions used to influence local legislators. As well as more trends in collective voting for smaller elections. Natural occurrences of the SDS phenomenon were recently observed in rural Canadian communities. After new legislation was passed by the Green Party to provide incentive for large-scale wind power developments, minority opinion holders gathered strategically. Their intention was to vote out the Green Party in the following election out of retaliation and they succeeded. The policy was amended shortly after, halting new production of wind power infrastructure. Within the U.S. a number of wind plant developments have encountered similar outcomes. In this investigation, we will seek to understand how population density, median income, and anti-wind power opinions within 3 km of wind power plants impacted the plant’s operational activity (2000-2016) to provide insights into localized resistance and its effects on renewable energy projects.
Main Takeaways
The overarching logit model analysis did not provide significant findings, but alluded to influential covariate interactions. The multivariate interaction model was leveraged to tabulate individually significant propensities. The probabilities revealed with an increase in median income there was a decreased in active wind energy developments, regardless of population size, and wind opinion. Further exploration of the covariate interactions are encouraged with expanded datasets, as a limited amount of data was utilized for this study, restricting the sample analysis. Using what is currently available, we will interpret the socioeconomic factors in this model as they relate to real-world social and political dynamics.
Population Density
Interestingly, pro/neutral-wind power regions revealed as population density and median income increased the odds of having an operational wind power plant decreased. Areas with pro/neutral-wind power sentiments in the lowest income percentile were found to have increased odds of an active wind plant as population density increased.
Whereas in anti-wind regions, across all median income groups, areas with higher population densities had marginally greater odds of having an inactive wind plant. These anti-wind findings suggested areas with higher median income are at greater risk of experiencing SDS. Both trends revealed could be an indication of the magnitude in which socioeconomic interactions influence sustainable energy developments.
Median Income
Overall, higher-income areas had decreased odds of having an active wind plant. Regardless of a pro/neutral-wind power sentiment, as the median wealth increased the odds of operational wind plant activity decreased. This may be an indicator that wealthier areas are at greater risk of experiencing SDS due to disproportionate access to financial resources that influence energy policy decisions. Advantages such as wealth assist minority opinion holders when strategically lobbying to tailor policies to better align with their priorities. Uneven socioeconomic power-dynamics lead to minority opinion holders preventing the development of wind power infrastructure, alongside other renewable energy solutions. Overall, these findings suggest that socioeconomic factors play a heavily influential role in shaping renewable energy projects.
Anti-Wind Infrastructure Opinion
As suspected, areas with higher opposition to wind infrastructure were less likely to have operational wind plants. Conversely, the opposite was seen for pro/neutral-wind opinion holding areas. This aligns with expectations that local opposition and social mobilization impacts the establishment of wind plants. Therefore, addressing and mitigating local resistance is essential for enhancing the feasibility and acceptance of wind energy initiatives.
The data source that was utilized in this project, US Wind Data, focuses on the public stance on wind infrastructure for census tract regions within a 3 km buffer zone of a wind infrastructure project. It contains categorical variables, binary variables, continuous socioeconomic factors such as % of races, % precinct political GOP affiliated voting share, mobilization tactics, and more. For simplicity’s sake, we’re going to focus on the variables below.
Describes the project operating status. In this study, we have converted it into a binary variable: 1 is operating, 0 is not_operating.
pop_den
Tract-level 2010 census data for population density (per mi^2)
med_inc
Tract-level 2010 census data for median income ($)
is_anti_wind
Binary measure of wind opposition: 1 is against wind power developments, 0 is pro wind power developments.
U.S. Wind Power Plant Locations
Before diving in, let’s get a sense of where we’ll be investigating.
Code
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Load Libraries ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# The following libraries were selected based on their functionality and ability to optimize our data for mapping.# Loading Librarieslibrary(tidyverse) # essential r package
Warning: package 'ggplot2' was built under R version 4.4.1
Code
library(sf) # package simplifies spatial dataframeslibrary(DT)library(tmap)library(terra)library(broom)library(stars)library(sjPlot)library(naniar)library(cowplot)library(leaflet)library(maptiles) library(ggthemes)library(ggspatial)library(patchwork)library(kableExtra)set.seed(99)#| label: map_data#| message: false#| warning: false#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Read & Raster ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# Below we will use the package `sf` to convert the lat/long vector data # into a raster geometry column. In this single line, we will also be # assigning the CRS EPSG:4326 to the sf data frame.# Coordinate Reference Systems, CRS, are required in order for the data# to be projected onto a map. The CRS was selected because it # provides a relatively proportionate display of the United States.# We are open to suggestions regarding our CRS if a different project better fits our data.# reading in & storing datawind_data <-read.csv("../data/wind_data/wind_data_usa.csv") # Confirm the Data Loaded Properly#head(wind_data) # displays the first 6 rows of the data# Let's read in our datawind_sf <- wind_data %>%# creates geometry column with desired crs st_as_sf(coords =c("longitude", "latitude"), crs =4326) %>%select(plant_name, everything())# quick CRS check#glimpse(crs(wind_sf)) # output should reveal WGS84, EPSG:4326#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Check Point! ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# Let's stop and see if our outputs are what we expect.# Were the lat/long columns correctly converted into a geometry column?# setdiff() is a way to quickly determine the differences between two data sets.# Sweet! we are looking good#setdiff(colnames(wind_sf), colnames(wind_data))#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Map Wind Plants ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# Using tmap(), we can visualize the locations of wind infrastructure power plants# throughout the United States. To achieve a more granular map, we'll need to# utilize another data set to create a base layer for our map in order to # observe these wind plants with respect to state and county jurisdictions. # First visual of the U.S. wind data provided by the geometry points# make sure tmap is in interactive modetmap_mode("view")# create the interactive maptm_shape(wind_sf) +tm_dots(col ="darkgreen",size =0.05,alpha =0.5 ) +tm_basemap("OpenStreetMap")
Figure 2: Toggleable map containing the US power plant data utilized for this anaylsis. Information about the area’s activism, polticial sentiment, median income ($), population composition, and more is available when clicking on a green bubble.
Data Distribution
Before jumping into any analysis, it’s important to get a sense of how the data is distributed and if there are any underlying trends or biases. We will employ a multivariate model to describe the effect of census tract level population density on the operating status of wind power infrastructure. A combination of binary and interaction logit regression will be considered. The initial model will apply OLS regression, this is really a formality to demonstrate why OLS is not the correct approach for interpreting our relationships of interest. The following will be a model with two continuous variables. These variables focus more on regionally dependent factors that intuitively seem to have an impact on mobilization variables that we don’t have time to cover in this project. We’ll be working with a mix of discrete and continuous data, so there some wrangling will be necessary to run the regressions we’re interested in.
Code
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Inspect & Standarize Data ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# Determining Variable Assignments for OLS#unique(wind_sf$status) # displays unique values in this# Need to rename status output variables# creating two categories: operating & not_operating# We are removing 'Operating | Decommissioned' because it skews the dataunwanted_status <-"Operating | Decommissioned"replacement_status <-"Uncertain Status"wind_sf$status[wind_sf$status== unwanted_status]<-"Uncertain Status"# were we successful ?#unique(wind_sf$status) # displays unique values in this# cleaning out NAs for OLSwind_sf <- wind_sf %>%filter(is.na(status) =='FALSE') %>%filter(is.na(is_anti_wind) =='FALSE') %>%filter(is.na(pop_den) =='FALSE') %>%filter(is.na(med_inc) =='FALSE') %>%filter(is.na(median_age) =='FALSE') %>%filter(is.na(n_turbs) =='FALSE')# were we successful ?#unique(wind_sf$status) # displays unique values in this# if_else preserves the data type but replaces unwanted valueswind_us <- wind_sf %>%mutate(status =if_else( status %in%c('Cancelled', 'Out of service (temporarily)', 'Standby', 'Decommissioned', 'Uncertain Status'), 'not_operating','operating') )# are our only outputs "operating" and "not_operating"?#print(unique(wind_us$status))# status as factor and reassigned valueswind_us <- wind_us %>%mutate(status =case_when(status =="operating"~1, status =="not_operating"~0))#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Check point! ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# are our only outputs 0 or 1?#paste("The indicator column contains", {unique(wind_us$status)})#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Violin Distribution ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# create the violin plot with log scaledensity_plot <-ggplot(data = wind_us, aes(x =factor(status), y = pop_den, fill =factor(status))) +geom_violin(alpha =0.6, color ="darkblue") +geom_jitter(col ="#F84C0B",width =0,height =0.05,alpha =0.35,size =4) +labs(title ="Population Density vs Wind Power Plant Operating Status",subtitle="Logorithimic Distribution",x ="Activation Status",y =expression("Population Density (Log Scale, "~ mi^-2~")")) +# rename x-axis labels for clarityscale_x_discrete(labels =c("0"="Inactive", "1"="Active")) +# apply logarithmic scale to y-axisscale_y_log10() +theme_538() +scale_fill_manual(values =c("skyblue", "darkblue")) +# adjust title font and alignmenttheme(plot.title =element_text(size =40,family ="Georgia", face ="bold",hjust = .99,color ="#293F2C"), # adjust subtitle font and alignmentplot.subtitle =element_text(size =38,family ="Georgia",color ="#293F2C",hjust =0.5), axis.title =element_text(size =36,family ="Georgia",color ="#293F2C"),axis.text =element_text(size =34,family ="Georgia",color ="#293F2C"),# move legend to the bottomlegend.position ="top", # remove legend title if not neededlegend.title =element_blank(), # adjust legend text sizelegend.text =element_text(size =34,family ="Georgia",color ="#293F2C"), # background color for legendlegend.key =element_rect(fill ="grey94", color ="grey94"), plot.background =element_rect(color ="#FDFBF7") ) +coord_flip()#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Jitter OLS + GLM ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# jitter plot with smooth linesjitter_plot_optimized <-ggplot(data = wind_us, aes(x = pop_den, y = status)) +# assign color to legend geom_jitter(aes(color ="Data Points"),width =0,height =0.05,alpha =0.6,# adjusted size for better visibilitysize =4) +geom_smooth(method ="lm",se =FALSE,aes(color ="OLS Line"),# slightly thicker line for visibilitysize =1.2, linetype ="solid") +geom_smooth(method ="glm",se = T,aes(color ="GLM Line"),size =1.2,linetype ="dashed",method.args =list(family ="binomial")) +labs(title ="Population Density vs Wind Power Plant Operating Status",subtitle="Logorithimic Distribution Regression Comparison",y ="Activation Status",x =expression("Population Density (Log Scale, "~ mi^-2~")")) +# rename yaxis labels for clarityscale_y_continuous(breaks =c(0, 1),labels =c("0"="Inactive", "1"="Active")) +scale_x_log10() +theme_538() +# adjust title font and alignmenttheme(plot.title =element_text(size =40,family ="Georgia", face ="bold",hjust = .99,color ="#293F2C"), # adjust subtitle font and alignmentplot.subtitle =element_text(size =38,family ="Georgia",hjust =0.5,color ="#293F2C"), axis.title =element_text(size =36,family ="Georgia",color ="#293F2C"),axis.text =element_text(size =34,family ="Georgia",color ="#293F2C"),# move legend to the bottomlegend.position ="top", # remove legend title if not neededlegend.title =element_blank(), # adjust legend text sizelegend.text =element_text(size =34,family ="Georgia",color ="#293F2C"), # background color for legendlegend.key =element_rect(fill ="grey94", color ="grey94"), plot.background =element_rect(color ="#FDFBF7") ) +# assign colors + labelsscale_color_manual(name ="Legend", values =c("Data Points"="#F84C0B", "OLS Line"="blue", "GLM Line"="skyblue"),labels =c("Data Points"="Data Points", "OLS Line"="OLS Line", "GLM Line"="GLM Line"))
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Code
# combine plots horizontallycombined_plot <- density_plot + jitter_plot_optimized +# arrange plots side by sideplot_layout(ncol =2) combined_plot
Figure 3: The two visual aids above contain violin plot with jitter points (left) and a comparative regression plot using OLS and GLM (right). The x-axis is the population density and the y-axis is the operational status for the wind power plant. Binary indicator variable will be status column active is 1, and inactive will be 0. Combining two figures provides us fuller insights into both the general trend and changes in probability of the binary outcome for the population density predictor.
The visualizations display the majority of the distribution lies within the actively operating wind infrastructure plants. This is an early indication of the degree of sample bias influencing our analysis. A trend of inactive plants and lower population density is notable in both figures. Collectively they demonstrate smaller population densities contain more inactive wind infrastructure plants. This could be attributed to with weight of a singular vote in regions with smaller demographics. Local mobilization of minority opinion holders in these regions have a greater availability to push back against policymakers. However, this visual does not encapsulate all of the necessary information required to determine this with full certainty. Our data set has low availability for non-operating infrastructure and as such, in the regression figure on the right these are being treated as outliers.
Odds Ratio for a Logit Model
To interpret the model, we compute the odds ratios for each coefficient, providing insight into how each variable and its interactions affect wind plant activity.
The Summary of Coefficients Table provides a comprehensive view of how each variable and their interactions contribute to the likelihood of wind plant activity, facilitating a better understanding of the model’s results (Table 1). For those less familiar with interpreting multivariate relationships:
Estimations reveal the model’s predicted effect on the activity status from a single unit increase in a variable. Another way to interpret this is using the odds ratio to understand the individual effects of each variable on the likelihood of having an active wind plant. The interaction terms show how the relationship between each pair of variables influences the odds of operational wind plant activity.
For instance, the interaction between pop_den:is_anti_wind captures how wind plant operational activity is affected with each unit increase in population density in areas with anti-wind power sentiment. In the instance of this dataset for this interaction model, an insignificant decrease in the odds of an activate wind plant was observed as population increases in anti-wind power areas. However, this is for the overall model and invites further exploration into how different socioeconomic groups vary in propensity for an active wind plant.
Code
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Comprehensive Model ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# fit the comprehensive modelcomprehensive_model <-glm(status ~ pop_den * med_inc * is_anti_wind,data = wind_us,family ='binomial')# get summary object to extract p-values and standard errorsmodel_summary <-summary(comprehensive_model)coef_table <- model_summary$coefficients# build the summary data framecoef_summary <-data.frame(Term =rownames(coef_table),Estimate = coef_table[, "Estimate"],`Odds Ratio`=exp(coef_table[, "Estimate"]),`p-value`= coef_table[, "Pr(>|z|)"],`Std Error`= coef_table[, "Std. Error"])# display tableDT::datatable(coef_summary, caption = htmltools::tags$caption(style ="Table 1: Summary of Coefficients, Odds Ratios, and P-values"),options =list(pageLength =10, scrollX =TRUE))
Operation Activity Trends
Our analysis of wind plant activity using logistic regression models has provided valuable insights into the factors influencing the presence of operational wind plants.The probability distribution of having an active local wind plant across different population density and median income quantiles reveals some notable trends:
Code
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- GLM Prob by Cat ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# percentile grouping (0–25%, 25–50%, 50–75%, 75–100%)wind_us <- wind_us %>%mutate(income_q =ntile(med_inc, 4),density_q =ntile(pop_den, 4),income_label =factor(paste0(seq(25, 100, by =25)[income_q], "%"),levels =paste0(seq(25, 100, by =25), "%") ),density_label =factor(paste0("Percentile ", seq(25, 100, by =25)[density_q], "%"),levels =paste0("Percentile ", seq(25, 100, by =25), "%") ),anti_label =ifelse(is_anti_wind ==1, "Anti-Wind", "Pro/Neutral") )# predict probabilitieswind_us$predicted_prob <-predict(comprehensive_model, newdata = wind_us, type ="response")#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- GLM Visualation by Cat ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# average predicted probability per groupplot_data <- wind_us %>%group_by(income_label, density_label, anti_label) %>%summarise(probability =mean(predicted_prob, na.rm =TRUE), .groups ="drop") %>%mutate(income_numeric =as.numeric(factor(income_label, levels =unique(income_label))))# visualize the dataggplot(plot_data, aes(x = income_label, y = probability, fill = density_label)) +# color by density side by sidegeom_bar(stat ="identity", position ="dodge") +# split by opinionfacet_wrap(~ anti_label) +# add lm trendline to more easily interpretgeom_smooth(aes(x = income_numeric, y = probability, group = density_label),method ="lm", color ="#F84C0B", se =FALSE, size =1) +# color palette for pop denscale_fill_manual(values =c("#cae6f2", "#7abbe2", "#2a7fb8", "#084081"),name ="Pop Density Quantile" ) +# customslabs(title ="Predicted Probability of Active Wind Infrastructure Projects",subtitle ="By Median Income Percentile, Population Density Percentile, and Wind Energy Sentiment",x ="Median Income Percentile Group",y ="Predicted Probability" ) +theme_minimal() +theme(text =element_text(family ="Georgia", size =22),axis.text.x =element_text(size =12),axis.text.y =element_text(size =14),axis.title.x =element_text(size =16),axis.title.y =element_text(size =16),plot.title =element_text(size =22, face ="bold", hjust =1.3),plot.subtitle =element_text(size =18, hjust =1.1),legend.position ="top",legend.title =element_text(size =16),legend.text =element_text(size =14),legend.key.size =unit(0.5, "cm"),legend.spacing.x =unit(0.2, "cm"),legend.spacing.y =unit(0.2, "cm"),legend.background =element_rect(fill ="white", color ="grey100", linewidth =0.5) ) +guides(fill =guide_legend(title =expression("Population Density (Log Scale, "~ mi^-2~")"), title.position ="top",title.hjust =0.5))
`geom_smooth()` using formula = 'y ~ x'
Figure 4: Side-by-side bar graphs demonstrating the relationships between population density (color of bars), median income percentiles (x-axis), and the predicted probability of an active wind plant status (y-axis). The graph on the left are the anti-wind opinion areas and the graph on the right are the pro/neutral-wind opinion areas. The orange trend line in each respective graph is the logit model fit.
The Big Question
May an assumption be made that regardless of population density, areas with higher median income have a greater propensity for inactive wind power plants and are more at risk of experiencing spatially distorted signalling?
To answer this, let’s touch a bit on the social psychology at potentially at play.
Code
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# ---- Plant Acitivity Visual ----#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~# summarize grouped probabilitiesgroupwise_propensity <- wind_us %>%group_by(income_label, density_label, anti_label) %>%summarise(# model-based probability // propensity scorepropensity_score =mean(predicted_prob),# actual treatment rate // where treatment is an active plantobserved_active_status =mean(status), # number of units in groupn =n(),# p-testp_value =prop.test(round(mean(status) *n()), n())$p.value,.groups ="drop" ) %>%# clean it up for plottingrename(`Income Percentile`= income_label,`Population Density Percentile`= density_label,`Wind Opinion`= anti_label,`Estimated Probability of Active Wind Plant`= propensity_score,`Observed Wind Plant Status Rate`= observed_active_status,`No. Matching Areas`= n,`p-value`= p_value )
Warning: There were 5 warnings in `summarise()`.
The first warning was:
ℹ In argument: `p_value = prop.test(round(mean(status) * n()), n())$p.value`.
ℹ In group 5: `income_label = 25%`, `density_label = Percentile 75%`,
`anti_label = "Anti-Wind"`.
Caused by warning in `prop.test()`:
! Chi-squared approximation may be incorrect
ℹ Run `dplyr::last_dplyr_warnings()` to see the 4 remaining warnings.
Code
# we want to know our plant names when we map, so doing a minimal joingroupwise_propensity <-st_join(x = groupwise_propensity, y = wind_us[, "plant_name"]) %>%# clean up the namerename(`Plant Name`= plant_name) %>%# reorder for a more informative table select(`Plant Name`,`Estimated Probability of Active Wind Plant`,`Observed Wind Plant Status Rate`, `Income Percentile`, `Population Density Percentile`, `Wind Opinion`, `No. Matching Areas`,`p-value`)# expand multipoints into points # each point gets repeated attributesprob_points <- groupwise_propensity %>%st_cast("POINT")
Warning in st_cast.sf(., "POINT"): repeating attributes for all sub-geometries
for which they may not be constant
Code
# create color palettepal <-colorFactor(palette =c("#E1BE6A", "#40B0A6"),domain = prob_points$`Wind Opinion`)# create popup content for each pointpopup_content <-paste0("<strong>Plant Name:</strong> ", prob_points$`Plant Name`, "<br>","<strong>Income Percentile:</strong> ", prob_points$`Income Percentile`, "<br>","<strong>Population Density Percentile:</strong> ", prob_points$`Population Density Percentile`, "<br>","<strong>Wind Opinion:</strong> ", prob_points$`Wind Opinion`, "<br>","<strong>Estimated Probability of Active Wind Plant:</strong> ", round(prob_points$`Estimated Probability of Active Wind Plant`, 3), "<br>","<strong>Observed Wind Plant Status Rate:</strong> ", round(prob_points$`Observed Wind Plant Status Rate`, 3), "<br>","<strong>No. Matching Areas:</strong> ", prob_points$`No. Matching Areas`)# display our updated map!leaflet(prob_points) %>%addTiles() %>%addCircleMarkers(color =~pal(`Wind Opinion`),radius =6,stroke =FALSE,fillOpacity =0.7,popup = popup_content ) %>%addLegend("bottomright",pal = pal,values =~`Wind Opinion`,title ="Wind Opinion",opacity =1 )
