---
title: "Shared vision for a decarbonized future energy system in the United States"
author: "Deidra Miniard"
date: "February 26, 2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### Load Libraries
```{r loadLibraries, warning=F}
library(tidyverse)
library(ggplot2)
library(kableExtra)
library(psych)
```

## Load Data
```{r loadData}
#load data 
EnergyMix <- read.csv("C:/EnergyMix.csv", fileEncoding = "UTF-8-BOM")
```

## Remove Bots
``` {r RemoveBots}
#eliminate participants who failed the attention check
EnergyMix <- subset(EnergyMix, aboutSense2 != "0")     
```


##Energy Mix estimates that did not add up to 100
```{r}
#Sum up current energy estimates
EnergyMix$curSum <- (EnergyMix$c_biomass + EnergyMix$c_coal + EnergyMix$c_oil + EnergyMix$c_ng + EnergyMix$c_nuclear + EnergyMix$c_solar + EnergyMix$c_wind + EnergyMix$c_hydro + EnergyMix$c_geotherm)

#sum up future estimates
EnergyMix$futSum <- (EnergyMix$f_biomass + EnergyMix$f_coal + EnergyMix$f_oil + EnergyMix$f_ng + EnergyMix$f_nuclear + EnergyMix$f_solar + EnergyMix$f_wind + EnergyMix$f_hydro + EnergyMix$f_geotherm)

#look at descriptive statistics
describe(EnergyMix$curSum)
describe(EnergyMix$futSum)

#Subset observations with current estimate sum not equal to 100
NotEqCur <- subset(EnergyMix, curSum != 100)
NotEqCur
#n = 818
818/2429
#33.3%

#Subset observations with future estimate sum not equal to 100
NotEqFut <- subset(EnergyMix, futSum != 100)
NotEqFut
#n = 656
656/2429
#27.0%
```
## standardize Current Energy Mix Percentage Estimates
``` {r normalizeCurrent}

#Normalize current energy mix estimates
EnergyMix$normc_oil <- (((as.numeric(as.character(EnergyMix$c_oil)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_ng <- (((as.numeric(as.character(EnergyMix$c_ng)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_coal <- (((as.numeric(as.character(EnergyMix$c_coal)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_nuclear <- (((as.numeric(as.character(EnergyMix$c_nuclear)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_geotherm <- (((as.numeric(as.character(EnergyMix$c_geotherm)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_solar <- (((as.numeric(as.character(EnergyMix$c_solar)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_wind <- (((as.numeric(as.character(EnergyMix$c_wind)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_biomass <- (((as.numeric(as.character(EnergyMix$c_biomass)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)

EnergyMix$normc_hydro <- (((as.numeric(as.character(EnergyMix$c_hydro)))/
                          (as.numeric(as.character(EnergyMix$c_oil)) + 
                             as.numeric(as.character(EnergyMix$c_ng)) + 
                             as.numeric(as.character(EnergyMix$c_coal)) + 
                             as.numeric(as.character(EnergyMix$c_nuclear)) + 
                             as.numeric(as.character(EnergyMix$c_geotherm)) + 
                             as.numeric(as.character(EnergyMix$c_solar)) + 
                             as.numeric(as.character(EnergyMix$c_wind)) + 
                             as.numeric(as.character(EnergyMix$c_biomass)) + 
                             as.numeric(as.character(EnergyMix$c_hydro)) ) ) * 100)
```

## Normalize Future Energy Mix  Percentage Estimates
``` {r}
#Normalize Future energy mix estimates
EnergyMix$normf_oil <- (((as.numeric(as.character(EnergyMix$f_oil)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_ng <- (((as.numeric(as.character(EnergyMix$f_ng)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_coal <- (((as.numeric(as.character(EnergyMix$f_coal)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_nuclear <- (((as.numeric(as.character(EnergyMix$f_nuclear)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_geotherm <- (((as.numeric(as.character(EnergyMix$f_geotherm)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_solar <- (((as.numeric(as.character(EnergyMix$f_solar)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_wind <- (((as.numeric(as.character(EnergyMix$f_wind)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_biomass <- (((as.numeric(as.character(EnergyMix$f_biomass)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)

EnergyMix$normf_hydro <- (((as.numeric(as.character(EnergyMix$f_hydro)))/
                          (as.numeric(as.character(EnergyMix$f_oil)) + 
                             as.numeric(as.character(EnergyMix$f_ng)) + 
                             as.numeric(as.character(EnergyMix$f_coal)) + 
                             as.numeric(as.character(EnergyMix$f_nuclear)) + 
                             as.numeric(as.character(EnergyMix$f_geotherm)) + 
                             as.numeric(as.character(EnergyMix$f_solar)) + 
                             as.numeric(as.character(EnergyMix$f_wind)) + 
                             as.numeric(as.character(EnergyMix$f_biomass)) + 
                             as.numeric(as.character(EnergyMix$f_hydro)) ) ) * 100)
head(EnergyMix)
```

##Policy Score
# Higher value means greater support for decarbonization policies and greater opposition to anti-decarbonization policies
```{r}
#Create political ideology groupings
EnergyMix$polIdeoGroup = dplyr::recode( EnergyMix$politicalBelief, "1" = "Liberal" , "2" = "Liberal" , "3" = "Liberal" , "4" = "Moderate" , "5" = "Conservative" , "6" = "Conservative" , "7" = "Conservative")

##Reverse code anti-decarbonization policies so that 1 means strongly support and 5 means strongly oppose
EnergyMix$Repp_fuelEcon <- dplyr::recode( EnergyMix$pp_fuelEcon, "1" = 5 , "2" = 4 , "3" = 3 , "4" = 2 , "5" = 1)
EnergyMix$Repp_lowerSubs <- dplyr::recode( EnergyMix$pp_lowerSubs, "1" = 5 , "2" = 4 , "3" = 3 , "4" = 2 , "5" = 1)
EnergyMix$Repp_constructPipe <- dplyr::recode( EnergyMix$pp_constructPipe, "1" = 5 , "2" = 4 , "3" = 3 , "4" = 2 , "5" = 1)
EnergyMix$Repp_tariffs <- dplyr::recode( EnergyMix$pp_tariffs, "1" = 5 , "2" = 4 , "3" = 3 , "4" = 2 , "5" = 1)
EnergyMix$Repp_relaxReg <- dplyr::recode( EnergyMix$pp_relaxReg, "1" = 5 , "2" = 4 , "3" = 3 , "4" = 2 , "5" = 1)
EnergyMix$Repp_coal <- dplyr::recode( EnergyMix$pp_coal, "1" = 5 , "2" = 4 , "3" = 3 , "4" = 2 , "5" = 1)

#Create Policy score: 1-5: Low value means against decarbonization policies and high value means support for decarbonization
#Nuclear included as decarbonization policy
EnergyMix$polScoreNuc <- ((EnergyMix$pp_nucCoal + 
                             EnergyMix$pp_carbonTax + 
                             EnergyMix$pp_fundRenew + 
                             EnergyMix$pp_subsidCars + 
                             EnergyMix$pp_haltPipe + 
                             EnergyMix$pp_fundBattery + 
                             EnergyMix$Repp_fuelEcon + 
                             EnergyMix$Repp_lowerSubs + 
                             EnergyMix$Repp_constructPipe + 
                             EnergyMix$Repp_tariffs + 
                             EnergyMix$Repp_relaxReg + 
                             EnergyMix$Repp_coal) / 12 )

#Policy score wtihout nuclear
EnergyMix$polScoreNoNuc <- ((EnergyMix$pp_carbonTax + 
                               EnergyMix$pp_fundRenew + 
                               EnergyMix$pp_subsidCars + 
                               EnergyMix$pp_haltPipe + 
                               EnergyMix$pp_fundBattery + 
                               EnergyMix$Repp_fuelEcon + 
                               EnergyMix$Repp_lowerSubs + 
                               EnergyMix$Repp_constructPipe + 
                               EnergyMix$Repp_tariffs + 
                               EnergyMix$Repp_relaxReg +
                               EnergyMix$Repp_coal) / 11 )

#check distribution and summarys stats
describe(EnergyMix$polScoreNuc)
describe(EnergyMix$polScoreNoNuc)
hist(EnergyMix$polScoreNoNuc)
hist(EnergyMix$polScoreNuc)
```


#Decarbonization score
```{r}
#Sum of future decarbonized energy resource estimates
EnergyMix$DecarbFut <- (EnergyMix$normf_biomass + 
                          EnergyMix$normf_geotherm + 
                          EnergyMix$normf_hydro + 
                          EnergyMix$normf_solar + 
                          EnergyMix$normf_wind + 
                          EnergyMix$normf_nuclear)

#Sum of current decarbonized energy resource estimates
EnergyMix$DecarbCur <- (EnergyMix$normc_biomass + 
                          EnergyMix$normc_geotherm + 
                          EnergyMix$normc_hydro + 
                          EnergyMix$normc_solar + 
                          EnergyMix$normc_wind + 
                          EnergyMix$normc_nuclear)

#Future Decarbonization score - current decarbonization score
##Variable for regression model 
EnergyMix$DecarbFutCur <- (EnergyMix$DecarbFut - EnergyMix$DecarbCur)

#Look at summary statistics for Decarb Score
describe(EnergyMix$DecarbFutCur)

#Check for NA in DecarbFutCur calculation
EnergyMix[!complete.cases(EnergyMix$DecarbFutCur),]
#One person put all 0s for future energy estimate 
```


#Climate Change Attitude, Importance, and Behavior variables
```{r}
#climate change attitude score; higher value means more belief in cc and more importance for climate change
EnergyMix$ccAttitudeBase <- ( ( EnergyMix$ccBelief + EnergyMix$ccImportant ) / 2 )

#Create data frame with climate change variables to look at alpha of ccAttitudeBase measure
ccdf <- data.frame(EnergyMix$ccBelief , 
                   EnergyMix$ccImportant )
alpha(ccdf)

#create climate change issue importance score relative to other issues
EnergyMix$ccImpUSRel <- (EnergyMix$US_cc - ( (EnergyMix$US_econ + EnergyMix$US_health + EnergyMix$US_imm) / 3 ) )

#Behavior measure for climate change
EnergyMix$ccBehavior <- ( ( EnergyMix$volunt_cc + EnergyMix$donate_cc + EnergyMix$govt_cc ) / 3 )

#Look at correlation between climate change variables
#Create data frame for cc variables
ccVar <- data.frame(EnergyMix$ccImpUSRel , EnergyMix$ccBehavior, EnergyMix$ccBelief , EnergyMix$ccImportant ,  EnergyMix$ccAttitudeBase)
corr<- cor(ccVar , use = "complete.obs")
round(corr , digits=2) 
```


##Demographics 
```{r}
#Income summary stats
describe(EnergyMix$income)

#Recode income by 1,000's 
EnergyMix$Re_income <- (EnergyMix$income / 1000)

#mean and median for age
describe(EnergyMix$age)
#range is 18 - 1992

#Remove 3 age outliers for regression model and summary statistics for age
EnergyMixRegData <- subset(EnergyMix , age < 1992)
range(EnergyMixRegData$age)
#range is 18 - 1983
EnergyMixRegData <- subset(EnergyMix , age < 1983)
range(EnergyMixRegData$age)
#range is 18-465
EnergyMixRegData <- subset(EnergyMix , age < 465)
range(EnergyMixRegData$age)
#range 18 - 84
#summary statistics for age
describe(EnergyMixRegData$age)

#Gender summary stats
table(EnergyMix$gender)
#Percentages for gender distribution
#male
1227/2429
#female
1194/2429
#other
8/2429

#Education summary stats
table(EnergyMix$education)
#1 = Some high school
9/2429
#2 = High school diploma or GED
265/2429
#3 = Some college
657/2429
#4 = College degree
1067/2429
#5 = Some graduate school
96/2429
#6 = Graduate degree
335/2429
#hold a college degree
(1067+96+335) / 2429

#Look at distribution of political party
table(EnergyMix$politicalParty)
#Get percentages of Political party
#Republican
574/2429
#Democrat
1110/2429
#Independent
745/2429

#Look at distribution of political ideology
table(EnergyMix$polIdeoGroup)
#Conservative
707/2429
#Moderate
494/2429
#Liberal
1228/2429

#correlation between policies
plot.plot <- pairs.panels(EnergyMix[48:59])
ggsave("plots/plot.pdf" , plot.plot, height = 15 , width = 20)
```

#Check for imagination condition effects on constructed variables or variables of interest
#Check for block order effects on constructed variables or variables of interest
```{r}
##Check for differences in variables of interest by question block order
lmPol <- lm(EnergyMix$polScoreNuc ~ EnergyMix$blockOrder)
anova(lmPol)

lmBehCC <- lm(EnergyMix$ccBehavior ~ EnergyMix$blockOrder)
anova(lmBehCC)

lmIssUSToday <- lm(EnergyMix$issImpUS ~ EnergyMix$blockOrder)
anova(lmIssUSToday)

lmIssGlobal <- lm(EnergyMix$issImpGlobal ~ EnergyMix$blockOrder)
anova(lmIssGlobal)

lmIssUSFut <- lm(EnergyMix$issImpUSFut ~ EnergyMix$blockOrder)
anova(lmIssUSFut)

lmIssGlFut <- lm(EnergyMix$issImpGlFut ~ EnergyMix$blockOrder)
anova(lmIssGlFut)

lmCC <- lm(EnergyMix$ccAttitudeBase ~ EnergyMix$blockOrder)
anova(lmCC)

##Check for differences in variables of interest from imagination experimental group
lmPol <- lm(EnergyMix$polScoreNoNuc ~ EnergyMix$imagGroup)
anova(lmPol)

lmBehCC <- lm(EnergyMix$ccBehavior ~ EnergyMix$imagGroup)
anova(lmBehCC)

lmIssUSToday <- lm(EnergyMix$issImpUS ~ EnergyMix$imagGroup)
anova(lmIssUSToday)

lmIssGlobal <- lm(EnergyMix$issImpGlobal ~ EnergyMix$imagGroup)
anova(lmIssGlobal)

lmIssUSFut <- lm(EnergyMix$issImpUSFut ~ EnergyMix$imagGroup)
anova(lmIssUSFut)

lmIssGlFut <- lm(EnergyMix$issImpGlFut ~ EnergyMix$imagGroup)
anova(lmIssGlFut)

lmCC <- lm(EnergyMix$ccAttitudeBase ~ EnergyMix$imagGroup)
anova(lmCC)

##Check for differences in variables of interest for an interaction between imagination group and block order
lmPol <- lm(EnergyMix$polScoreNuc ~ EnergyMix$imagGroup*EnergyMix$blockOrder)
anova(lmPol)

lmIssUSToday <- lm(EnergyMix$issImpUS ~ EnergyMix$imagGroup*EnergyMix$blockOrder)
anova(lmIssUSToday)

lmIssGlobal <- lm(EnergyMix$issImpGlobal ~ EnergyMix$imagGroup*EnergyMix$blockOrder)
anova(lmIssGlobal)

lmIssUSFut <- lm(EnergyMix$issImpUSFut ~ EnergyMix$imagGroup*EnergyMix$blockOrder)
anova(lmIssUSFut)

lmIssGlFut <- lm(EnergyMix$issImpGlFut ~ EnergyMix$imagGroup*EnergyMix$blockOrder)
anova(lmIssGlFut)

lmCC <- lm(EnergyMix$ccAttitudeBase ~ EnergyMix$imagGroup*EnergyMix$blockOrder)
anova(lmCC)

lmPol <- lm(EnergyMix$polScoreNuc ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmPol)

lmBehCC <- lm(EnergyMix$ccBehavior ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmBehCC)

lmIssUSToday <- lm(EnergyMix$issImpUS ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmIssUSToday)

lmIssGlobal <- lm(EnergyMix$issImpGlobal ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmIssGlobal)

lmIssUSFut <- lm(EnergyMix$issImpUSFut ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmIssUSFut)

lmIssGlFut <- lm(EnergyMix$issImpGlFut ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmIssGlFut)

lmCC <- lm(EnergyMix$ccAttitudeBase ~ EnergyMix$imagGroup + EnergyMix$blockOrder)
anova(lmCC)
```

#Energy Mix Analysis Results
#Calculating percent error
```{r}
#Calculating percent error for current energy mix estimates
( ( (mean(EnergyMix$normc_ng)) - 30.6 ) / 30.6 ) * 100
( ( (mean(EnergyMix$normc_oil)) - 36.5 ) / 36.5 ) * 100
( ( (mean(EnergyMix$normc_coal)) - 13.1 ) / 13.1 ) * 100
( ( (mean(EnergyMix$normc_nuclear)) -  8.4) / 8.4 ) * 100
( ( (mean(EnergyMix$normc_solar)) - 0.9 ) / 0.9 ) * 100
( ( (mean(EnergyMix$normc_wind)) - 2.5 ) / 2.5 ) * 100
( ( (mean(EnergyMix$normc_biomass)) - 5.1 ) / 5.1 ) * 100
( ( (mean(EnergyMix$normc_hydro)) - 2.7 ) / 2.7 ) * 100
( ( (mean(EnergyMix$normc_geotherm)) - 0.2 ) / 0.2 ) * 100

#calculating percent error for all renewables
EnergyMix$Renew <- (EnergyMix$normc_solar + EnergyMix$normc_wind + EnergyMix$normc_biomass + EnergyMix$normc_hydro + EnergyMix$normc_geotherm + EnergyMix$normc_nuclear)
0.9 + 2.5 + 5.1 + 2.7 + 0.2 + 8.4
( ( (mean(EnergyMix$Renew)) - 19.8) / 19.8 ) * 100
mean(EnergyMix$Renew)

#Percent error for all fossil fuels
EnergyMix$ff <- (EnergyMix$normc_ng + EnergyMix$normc_oil + EnergyMix$normc_coal)
30.6 + 36.5 + 13.1
( ( (mean(EnergyMix$ff)) - 80.2 ) / 80.2 ) * 100
```

#Energy Mix Analysis
#Slopes for change from current to future estimates
```{r}
#Calculate slope for each energy source by future estimate minus current estimate
mean(EnergyMix$normf_oil, na.rm = TRUE) - mean(EnergyMix$normc_oil)
mean(EnergyMix$normf_ng, na.rm = TRUE) - mean(EnergyMix$normc_ng)
mean(EnergyMix$normf_coal, na.rm = TRUE) - mean(EnergyMix$normc_coal)
mean(EnergyMix$normf_nuclear, na.rm = TRUE) - mean(EnergyMix$normc_nuclear)
mean(EnergyMix$normf_solar, na.rm = TRUE) - mean(EnergyMix$normc_solar)
mean(EnergyMix$normf_wind, na.rm = TRUE) - mean(EnergyMix$normc_wind)
mean(EnergyMix$normf_biomass, na.rm = TRUE) - mean(EnergyMix$normc_biomass)
mean(EnergyMix$normf_hydro, na.rm = TRUE) - mean(EnergyMix$normc_hydro)
mean(EnergyMix$normf_geotherm, na.rm = TRUE) - mean(EnergyMix$normc_geotherm)
```

#Energy Mix Analysis
#Percent participants who want a fully decarbonized future
```{r}
#Subset individuals who have a decarbonized future sum equal to 100 (fully decarbonized)
FullDecarb <- subset(EnergyMix, DecarbFut == 100)
#703 participants

#look at participants by political ideology
table(FullDecarb$polIdeoGroup)
table(EnergyMix$polIdeoGroup)
#Conservatives
127/707
#moderates
141/494
#liberals
435/1228
```

#Energy Mix Analysis Results
#Analyze differences between current and future energy estimates using ANOVAS and TTESTs
```{r}
#Alpha level for ANOVAs after applying bonferroni correction
.05/9
#0.00556

#Alpha level for t-tests after applying bonferroni correction
.05/27
#0.00185

#ANOVAs for ENERGY SOURCE BY POL IDEO GROUP
biomass <- aov(EnergyMix$normc_biomass ~ EnergyMix$polIdeoGroup)
summary(biomass)

coal <- aov(EnergyMix$normc_coal ~ EnergyMix$polIdeoGroup)
summary(coal)

geo <- aov(EnergyMix$normc_geotherm ~ EnergyMix$polIdeoGroup)
summary(geo)

hydro <- aov(EnergyMix$normc_hydro ~ EnergyMix$polIdeoGroup)
summary(hydro)

ng <- aov(EnergyMix$normc_ng ~ EnergyMix$polIdeoGroup)
summary(ng)

nuclear <- aov(EnergyMix$normc_nuclear ~ EnergyMix$polIdeoGroup)
summary(nuclear)

oil <- aov(EnergyMix$normc_oil ~ EnergyMix$polIdeoGroup)
summary(oil)

solar <- aov(EnergyMix$normc_solar ~ EnergyMix$polIdeoGroup)
summary(solar)

wind <- aov(EnergyMix$normc_wind ~ EnergyMix$polIdeoGroup)
summary(wind)

#ANOVA FUTURE ENERGY MIX BY POL IDEO GROUP
biomass <- aov(EnergyMix$normf_biomass ~ EnergyMix$polIdeoGroup)
summary(biomass)

coal <- aov(EnergyMix$normf_coal ~ EnergyMix$polIdeoGroup)
summary(coal)

geo <- aov(EnergyMix$normf_geotherm ~ EnergyMix$polIdeoGroup)
summary(geo)

hydro <- aov(EnergyMix$normf_hydro ~ EnergyMix$polIdeoGroup)
summary(hydro)

ng <- aov(EnergyMix$normf_ng ~ EnergyMix$polIdeoGroup)
summary(ng)

nuclear <- aov(EnergyMix$normf_nuclear ~ EnergyMix$polIdeoGroup)
summary(nuclear)

oil <- aov(EnergyMix$normf_oil ~ EnergyMix$polIdeoGroup)
summary(oil)

solar <- aov(EnergyMix$normf_solar ~ EnergyMix$polIdeoGroup)
summary(solar)

wind <- aov(EnergyMix$normf_wind ~ EnergyMix$polIdeoGroup)
summary(wind)


## T-Tests
EnMix3Lib <- subset(EnergyMix, polIdeoGroup %in% c("Liberal"))
EnMix3Con <- subset(EnergyMix, polIdeoGroup %in% c("Conservative"))
EnMix3Mod <- subset(EnergyMix, polIdeoGroup %in% c("Moderate"))

#Hydro
t.test(EnMix3Lib$normc_hydro , EnMix3Con$normc_hydro )
t.test(EnMix3Lib$normc_hydro , EnMix3Mod$normc_hydro )
t.test(EnMix3Con$normc_hydro , EnMix3Mod$normc_hydro )

#NUCLEAR
t.test(EnMix3Lib$normc_nuclear , EnMix3Mod$normc_nuclear )
t.test(EnMix3Lib$normc_nuclear , EnMix3Con$normc_nuclear )
t.test(EnMix3Con$normc_nuclear , EnMix3Mod$normc_nuclear )

#OIL
t.test(EnMix3Lib$normc_oil , EnMix3Con$normc_oil )
t.test(EnMix3Lib$normc_oil , EnMix3Mod$normc_oil )
t.test(EnMix3Mod$normc_oil , EnMix3Con$normc_oil )

#WIND
t.test(EnMix3Lib$normc_wind , EnMix3Con$normc_wind )
t.test(EnMix3Lib$normc_wind , EnMix3Mod$normc_wind )
t.test(EnMix3Mod$normc_wind , EnMix3Con$normc_wind )


######T tests for FUTURE SIG DIFS#######
#NATURAL GAS
t.test(EnMix3Lib$normf_ng , EnMix3Con$normf_ng )
t.test(EnMix3Lib$normf_ng , EnMix3Mod$normf_ng )
t.test(EnMix3Con$normf_ng , EnMix3Mod$normf_ng )

#NUCLEAR
t.test(EnMix3Lib$normf_nuclear , EnMix3Mod$normf_nuclear )
t.test(EnMix3Lib$normf_nuclear , EnMix3Con$normf_nuclear )
t.test(EnMix3Con$normf_nuclear , EnMix3Mod$normf_nuclear )

#OIL
t.test(EnMix3Lib$normf_oil , EnMix3Con$normf_oil )
t.test(EnMix3Lib$normf_oil , EnMix3Mod$normf_oil )
t.test(EnMix3Mod$normf_oil , EnMix3Con$normf_oil )

#WIND
t.test(EnMix3Lib$normf_wind , EnMix3Con$normf_wind )
t.test(EnMix3Lib$normf_wind , EnMix3Mod$normf_wind )
t.test(EnMix3Mod$normf_wind , EnMix3Con$normf_wind )

#COAL
t.test(EnMix3Lib$normf_coal , EnMix3Con$normf_coal )
t.test(EnMix3Lib$normf_coal , EnMix3Mod$normf_coal )
t.test(EnMix3Con$normf_coal , EnMix3Mod$normf_coal )

#SOLAR
t.test(EnMix3Lib$normf_solar , EnMix3Con$normf_solar )
t.test(EnMix3Lib$normf_solar , EnMix3Mod$normf_solar )
t.test(EnMix3Mod$normf_solar , EnMix3Con$normf_solar )
```

#Preferences for Policy Support Analysis
```{r}
#Create decarb policy support score
EnergyMix$DecarbPolScore <- ((EnergyMix$pp_carbonTax + EnergyMix$pp_haltPipe + EnergyMix$pp_nucCoal + EnergyMix$pp_fundRenew + EnergyMix$pp_subsidCars + EnergyMix$pp_fundBattery) / 6)

#Create anti-decarb policy support score
EnergyMix$AntiDecarbPolScore <- ( (EnergyMix$pp_relaxReg + EnergyMix$pp_constructPipe + EnergyMix$pp_coal + EnergyMix$pp_lowerSubs + EnergyMix$pp_fuelEcon + EnergyMix$pp_tariffs) / 6 )

#look at alpha for Decarb score
DecarbPolAlpha <- data.frame(EnergyMix$pp_carbonTax , 
                       EnergyMix$pp_fundBattery , 
                       EnergyMix$pp_fundRenew , 
                       EnergyMix$pp_haltPipe , 
                       EnergyMix$pp_subsidCars,
                       EnergyMix$pp_nucCoal)
alpha(DecarbPolAlpha)

#look at alpha for decarb score without nuclear
DecarbPolAlphaNoNuc <- data.frame(EnergyMix$pp_carbonTax , 
                       EnergyMix$pp_fundBattery , 
                       EnergyMix$pp_fundRenew , 
                       EnergyMix$pp_haltPipe , 
                       EnergyMix$pp_subsidCars)
alpha(DecarbPolAlphaNoNuc)

AntiDecarbAlpha <- data.frame(EnergyMix$pp_relaxReg ,
                              EnergyMix$pp_constructPipe ,
                              EnergyMix$pp_coal ,
                              EnergyMix$pp_lowerSubs ,
                              EnergyMix$pp_fuelEcon , 
                              EnergyMix$pp_tariffs)
alpha(AntiDecarbAlpha)

#ttest on antidecarb vs decarb
t.test(EnergyMix$AntiDecarbPolScore , EnergyMix$DecarbPolScore)

#Anova on decarb and anti-decarb scores for political ideological groups
DecarbPol <- aov(EnergyMix$DecarbPolScore ~ EnergyMix$polIdeoGroup)
summary(DecarbPol)

AntiDecarbPol <- aov(EnergyMix$AntiDecarbPolScore ~ EnergyMix$polIdeoGroup)
summary(AntiDecarbPol)

## T-Tests for political ideological differences
EnMix3Lib <- subset(EnergyMix, polIdeoGroup %in% c("Liberal"))
EnMix3Con <- subset(EnergyMix, polIdeoGroup %in% c("Conservative"))
EnMix3Mod <- subset(EnergyMix, polIdeoGroup %in% c("Moderate"))
       
##Decarb
t.test(EnMix3Lib$DecarbPolScore , EnMix3Con$DecarbPolScore)
t.test(EnMix3Lib$DecarbPolScore , EnMix3Mod$DecarbPolScore)
t.test(EnMix3Mod$DecarbPolScore , EnMix3Con$DecarbPolScore)

#AntiDecarb
t.test(EnMix3Lib$AntiDecarbPolScore , EnMix3Con$AntiDecarbPolScore)
t.test(EnMix3Lib$AntiDecarbPolScore , EnMix3Mod$AntiDecarbPolScore)
t.test(EnMix3Mod$AntiDecarbPolScore , EnMix3Con$AntiDecarbPolScore)
```


#Issue importance
#Percentages across all participants
```{r}
#convert count to percentage
#Convert data to long format
EnergyMix4 <- EnergyMix %>%
  dplyr::select(c(ID,issImpUS,issImpGlobal,issImpUSFut,issImpGlFut,polIdeoGroup))  %>%
  group_by(ID) %>%
  gather(key = "factorLevel" , value = "issue" , -c(ID,polIdeoGroup))
table(EnergyMix$ID)

#Recode issue to category name and set order
EnergyMix4$issue <- dplyr::recode( EnergyMix4$issue, "1" = "Economy and jobs" , "2" = "Climate change" , "3" = "Immigration" , "4" = "Access to quality health care")
EnergyMix4$issue <- factor(EnergyMix4$issue , levels = c("Climate change" , "Economy and jobs" , "Access to quality health care" , "Immigration"))

EnergyMix4$focalPoint <- dplyr::recode(EnergyMix4$factorLevel , "issImpUS" = "U.S. today" , "issImpGlobal" = "world today" , "issImpUSFut" = "U.S. future" , "issImpGlFut" = "world future")

#obtain count of participants
FreqIss <- xtabs(~focalPoint + issue , data = EnergyMix4)
FreqIss

#turn frequency count into data frame
FreqIssDF <- data.frame(FreqIss)
head(FreqIssDF)

#convert data into wide format
FreqIssDF <- spread(FreqIssDF , issue , Freq)

#convert count to percentage
FreqIssDF$cc <- (FreqIssDF$`Climate change` / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100
FreqIssDF$econ <- (FreqIssDF$`Economy and jobs` / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100
FreqIssDF$health <- (FreqIssDF$`Access to quality health care` / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100
FreqIssDF$imm <- (FreqIssDF$Immigration / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100

#convert back into long format to get data frame with percentages
FreqIssLong <- FreqIssDF %>%
  dplyr::select(c(focalPoint, cc, econ, health, imm)) %>%
  gather(key = "percent" , value = "percentScore" , -c(focalPoint))
```

##Regression Analysis and Recoding
```{r}
##Recode gender and income
#Recode gender to male = 1, else = 0
EnergyMixRegData$Re_gender = dplyr::recode( EnergyMixRegData$gender, "1" = "1" , "2" = "0" , "3" = "0")

#Alpha for decarbonization policy support dependent variable
#look at alpha for policy score
policydf <- data.frame(EnergyMixRegData$Repp_fuelEcon , 
                       EnergyMixRegData$Repp_coal , 
                       EnergyMixRegData$Repp_constructPipe , 
                       EnergyMixRegData$Repp_lowerSubs , 
                       EnergyMixRegData$Repp_tariffs , 
                       EnergyMixRegData$Repp_relaxReg , 
                       EnergyMixRegData$pp_carbonTax , 
                       EnergyMixRegData$pp_fundBattery , 
                       EnergyMixRegData$pp_fundRenew , 
                       EnergyMixRegData$pp_haltPipe , 
                       EnergyMixRegData$pp_subsidCars , 
                       EnergyMixRegData$pp_nucCoal)
alpha(policydf)

#simple model
lmDecarbPolSimple <- lm(polScoreNuc ~ DecarbFutCur , data = EnergyMixRegData)
summary(lmDecarbPolSimple)

#Include demographic control variables
lmDecarbPolDemo <- lm(polScoreNuc ~ DecarbFutCur + 
                    politicalBelief + 
                    Re_gender +
                    Re_income +
                    age +
                    education, 
                  data = EnergyMixRegData)
summary(lmDecarbPolDemo)

#Full model with climate change variables
lmDecarbPolFull <- lm(polScoreNuc ~ DecarbFutCur +
                        politicalBelief + 
                        Re_gender +
                        Re_income +
                        age +
                        education +
                        ccAttitudeBase + 
                        ccImpUSRel , 
                      data = EnergyMixRegData )
summary(lmDecarbPolFull)

#Look at regression residuals for Final model
par(mfrow = c(2, 2))  # Split the plotting panel into a 2 x 2 grid
plot(lmDecarbPolFull)
```


##Figures
#Energy Estimate Figure
```{r}
head(EnergyMix)
#Convert data into long format
EnergyMix2 <- EnergyMix %>%
  dplyr::select(c(ID:imagGroup, normc_oil:normf_hydro, Repp_fuelEcon:polIdeoGroup)) %>%
  group_by(ID) %>%
  gather(key = "energySource", value = "energyEstimate" , -c(ID:imagGroup, Repp_fuelEcon:polIdeoGroup))

#seperate energy source to name of source only by seperating string text
EnergyMix2$energySource1 <- (substr(EnergyMix2$energySource, 7 , 20))

#separate time (c or f) by seperate text string
EnergyMix2$time <- (substr(EnergyMix2$energySource , 5 , 5))

#convert c or f to current or future respectively
EnergyMix2$timeEst <- dplyr::recode( EnergyMix2$time, "c" = "current" , "f" = "future" )

#check that variables were created correctly
head(EnergyMix2)

#Reorder energy sources for facet wrap
EnergyMix2$energySource1 <- factor(EnergyMix2$energySource1 , levels = c("oil" , "ng" , "coal" , "nuclear" , "biomass" , "hydro" , "wind" , "solar" , "geotherm") )

#Reorder political ideological groups
EnergyMix2$polIdeoGroup <- factor(EnergyMix2$polIdeoGroup , levels = c("Conservative" , "Moderate" , "Liberal"))

##http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
library(Rmisc)
#Use summary statistics to create figure by making summary statistics into a data frame
enMixSum <- summarySE(EnergyMix2 , measurevar="energyEstimate" , na.rm = TRUE, groupvars = c("energySource1" , "polIdeoGroup" , "timeEst"))
enMixSum

#Rename labels of energy sources for facet wrap command
facLabel <- c(biomass = "biomass" , coal = "coal" , geotherm = "geothermal" , hydro = "hydroelectricity" , ng = "natural gas" , nuclear = "nuclear" , oil = "oil" , solar = "solar" , wind = "wind")

#Energy Estimate Figure
enEst.plot <- ggplot(enMixSum , aes(x= timeEst , y = energyEstimate , group=polIdeoGroup , colour = polIdeoGroup )) +
  scale_colour_manual(values = c("Conservative"="red", "Moderate" = "darkgray" , "Liberal" = "blue" , "Actual" = "black" ) , name = "Key") + ##set colors for polIdeoGroup 
  geom_line(size = 1.5 , position = position_dodge(width = .25) ) + #Create line to connect points and generate random distance to space out lines
  geom_point(size = 1.5 , position = position_dodge(width = .25)) + #change size and type of point with randome distance
  geom_errorbar(aes(ymin=energyEstimate-ci, ymax=energyEstimate+ci) , width=.20 , position = position_dodge(width=.25)) + #error bars
  facet_wrap( ~energySource1 , labeller = labeller(energySource1 = facLabel) , ncol = 3 ) + #seperate by energy source and indicate order
  ggtitle(label = "Current and future energy estimates by energy source and political ideology") + #title & subtitle
  xlab("Time") + #x label
  ylab("Percentage estimate") +  #y label
  theme_bw() + #plot theme
  theme(plot.title = element_text(size = 18 , hjust = 0.5 , margin = margin(t = 15 , r = 5 , b = 20 , l = 5))) + #plot title settings
  theme(axis.text.x = element_text( size = 14 , hjust = 0.5)) + #x axis settings
  theme(axis.text.y = element_text( size = 14 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=14)) + #legend settings
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 16 , hjust = 0.5 )) + #legend title settings
  theme(strip.text.x = element_text(size = 16 )) + #plot subtitle settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(legend.position = "bottom") +
  ylim(0,40) + #set y axis limit
    #plot actual energy mix values onto figure
  geom_point(data = subset(EnergyMix2, energySource1 == "oil") , aes(x="current", y=36.5 , fill = "Actual Contribution" ), colour="black" , shape = 16 , size = 3.5 ) +
  geom_point(data = subset(EnergyMix2, energySource1 == "ng") , aes(x="current", y=30.6  ), colour="black" , shape = 16 , size = 3.5) +
  geom_point(data = subset(EnergyMix2, energySource1 == "coal") , aes(x="current", y=13.1  ), colour="black" , shape = 16 , size = 3.5) +
  geom_point(data = subset(EnergyMix2, energySource1 == "nuclear") , aes(x="current", y=8.4  ), colour="black" , shape = 16 , size = 3.5) + 
  geom_point(data = subset(EnergyMix2, energySource1 == "biomass") , aes(x="current", y=5.1  ), colour="black" , shape = 16 , size = 3.5) +
  geom_point(data = subset(EnergyMix2, energySource1 == "hydro") , aes(x="current", y=2.7  ), colour="black" , shape = 16 , size = 3.5) +
  geom_point(data = subset(EnergyMix2, energySource1 == "wind") , aes(x="current", y=2.5  ), colour="black" , shape = 16 , size = 3.5) +
  geom_point(data = subset(EnergyMix2, energySource1 == "solar") , aes(x="current", y=0.9  ), colour="black" , shape = 16 , size = 3.5) +
  geom_point(data = subset(EnergyMix2, energySource1 == "geotherm") , aes(x="current", y=0.2  ), colour="black" , shape = 16 , size = 3.5) +
  theme(legend.title = element_blank()) +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm'))) 

#save plot as PDF
ggsave("C:/Users/Owner/Desktop/enEst.pdf" , enEst.plot , height = 11, width = 10)  
```


##Policy Figure
```{r}
#Convert data into long format
head(EnergyMix)
EnergyMix3 <- EnergyMix %>%
  dplyr::select(c(ID:issImpGlFut, pp_nucCoal:pp_coal, volunt_cc:AntiDecarbPolScore)) %>%
  group_by(ID) %>%
  gather(key = "policy" , value = "policyScore" , -c(ID:issImpGlFut, volunt_cc:AntiDecarbPolScore))

#Create a frequency data table using xtabs command
Freq <- xtabs(~policyScore + polIdeoGroup  + policy, data = EnergyMix3 )

#Turn Frequency table into a data frame
Freq.df <- data.frame(Freq)
#Check the data frame

#Convert data back into wide format
FreqWide <- spread(Freq.df , policyScore , Freq)

#Create participant percentages
FreqWide$StronglyOppose <- ( FreqWide$`1` / ( FreqWide$`1` + FreqWide$`2` + FreqWide$`3` + FreqWide$`4` + FreqWide$`5`) ) * 100
FreqWide$Oppose <- ( FreqWide$`2` / ( FreqWide$`1` + FreqWide$`2` + FreqWide$`3` + FreqWide$`4` + FreqWide$`5`) ) * 100
FreqWide$Neutral <- ( FreqWide$`3` / ( FreqWide$`1` + FreqWide$`2` + FreqWide$`3` + FreqWide$`4` + FreqWide$`5`) ) * 100
FreqWide$Support <- ( FreqWide$`4` / ( FreqWide$`1` + FreqWide$`2` + FreqWide$`3` + FreqWide$`4` + FreqWide$`5`) ) * 100
FreqWide$StronglySupport <- ( FreqWide$`5` / ( FreqWide$`1` + FreqWide$`2` + FreqWide$`3` + FreqWide$`4` + FreqWide$`5`) ) * 100

#convert data back into long format
FreqLong <- FreqWide %>%
  dplyr::select(c(polIdeoGroup:policy , StronglyOppose:StronglySupport)) %>%
  gather(key = "percent" , value = "percentScore" , -c(polIdeoGroup:policy))

#specify order for policy support
FreqLong$percent <- factor(FreqLong$percent , levels = c("StronglyOppose" , "Oppose" , "Neutral" , "Support" , "StronglySupport"))

#Rename policy labels
facLabel <- c(pp_nucCoal = "replace coal with nuclear" , pp_carbonTax = "carbon tax" , pp_fundRenew = "fund renewable resources" , pp_subsidCars = "subsidize electric cars" , pp_haltPipe = "halt oil & NG pipelines" , pp_fundBattery = "fund battery storage research" , pp_fuelEcon = "lower fuel economy standards" , pp_lowerSubs = "lower renewable energy subsidies" , pp_constructPipe = "construct pipelines" , pp_tariffs = "tax imported solar panels" , pp_relaxReg = "relax drilling regulations" , pp_coal = "invest in new coal plants")

#Indicate policy order
FreqLong$policy <- factor(FreqLong$policy , levels = c("pp_carbonTax" , "pp_haltPipe" , "pp_nucCoal" , "pp_fundRenew" , "pp_subsidCars" , "pp_fundBattery" , "pp_relaxReg" , "pp_constructPipe" , "pp_coal" , "pp_lowerSubs" , "pp_fuelEcon" , "pp_tariffs"))

FreqLong$polIdeoGroup <- factor(FreqLong$polIdeoGroup , levels = c("Conservative" , "Moderate" , "Liberal"))

#Make x axis 2 lineshea
xLabel <- c("strongly\noppose",
            "oppose",
            "neutral",
            "support",
            "strongly\nsupport")

#Create 2 sperate figures for decarbonization policies and anti-decarbonization policies
library(grid)
library(gridBase)
library(gridGraphics)
library(gridExtra)

#Seperate policies into anti-decarbonization and decarbonization dataframes
FreqLongDecarb <- subset(FreqLong, policy %in% c("pp_carbonTax" , "pp_haltPipe" , "pp_nucCoal" , "pp_fundRenew" , "pp_subsidCars" , "pp_fundBattery"))
FreqLongAntiDecarb <- subset(FreqLong, policy %in% c("pp_relaxReg" , "pp_constructPipe" , "pp_coal", "pp_lowerSubs" , "pp_fuelEcon" , "pp_tariffs"))

#Specify order for policy support
FreqLongDecarb$percent <- factor(FreqLongDecarb$percent , levels = c("StronglyOppose" , "Oppose" , "Neutral" , "Support" , "StronglySupport"))
FreqLongAntiDecarb$percent <- factor(FreqLongAntiDecarb$percent , levels = c("StronglyOppose" , "Oppose" , "Neutral" , "Support" , "StronglySupport"))

#Specify order for policies
FreqLongDecarb$policy <- factor(FreqLongDecarb$policy , levels = c("pp_carbonTax" , "pp_haltPipe" , "pp_nucCoal" , "pp_fundRenew" , "pp_subsidCars" , "pp_fundBattery"))
FreqLongAntiDecarb$policy <- factor(FreqLongAntiDecarb$policy , levels = c("pp_relaxReg" , "pp_constructPipe" , "pp_coal" , "pp_lowerSubs" , "pp_fuelEcon" , "pp_tariffs"))

#Rename Decarbonization policies
facLabelDecarb <- c(pp_nucCoal = "replace coal with nuclear" , pp_carbonTax = "carbon tax" , pp_fundRenew = "fund renewable resources" , pp_subsidCars = "subsidize electric cars" , pp_haltPipe = "halt oil pipelines" , pp_fundBattery = "fund battery storage research")

#Rename anti-decarbonization policies
facLabelAntiDecarb <- c(pp_fuelEcon = "lower fuel economy standards" , pp_lowerSubs = "lower renewable energy subsidies" , pp_constructPipe = "construct oil pipelines" , pp_tariffs = "tax imported solar panels" , pp_relaxReg = "relax drilling regulations" , pp_coal = "invest in new coal plants")

#plot decarbonization policy support
Freq.plot.Decarb <- ggplot(FreqLongDecarb , aes( x = percent, y = percentScore, group = polIdeoGroup, colour = polIdeoGroup)) + 
  scale_colour_manual(values = c("Conservative"="red", "Moderate" = "darkgray" , "Liberal" = "blue" ) , name = "Key") +
  scale_x_discrete(labels = xLabel) +
  geom_line(size = 1.5  )  + #Create line to connect points 
  geom_point(size = 1.5 ) + #change size and type of point + 
  facet_wrap( ~ policy , labeller = labeller(policy = facLabelDecarb), ncol = 3) +
  theme_bw() +
  theme(legend.position = "none") +
  ggtitle(label = "Decarbonization policies") + 
  theme(axis.title.x = element_blank())+
  ylab("Percentage of participants") +
  theme(plot.title = element_text(size = 18 , hjust = 0.5 , margin = margin(t = 10 , r = 5 , b = 20 , l = 5) )) +
  theme(axis.text.x = element_text( size = 12 , hjust = 0.5 ) )  +
  theme(axis.text.y = element_text( size = 14 , hjust = 0.5)) +
  theme(legend.text=element_text(size=14)) + 
  theme(legend.title = element_text(size = 14 , margin = margin( r = 2 , l = 0.01))) + 
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 , colour = "white" )) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 16 , hjust = 0.5 )) + #legend title settings
  theme(axis.title.x = element_text(margin = margin(t = 30 , r = 0.5 , b = 0 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(strip.text = element_text(size = 16) ) +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm') )) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) 

#Plot anitdecarb policy support
Freq.plot.AntiDecarb <- ggplot(FreqLongAntiDecarb , aes( x = percent, y = percentScore, group = polIdeoGroup, colour = polIdeoGroup)) + 
  scale_colour_manual(values = c("Conservative"="red", "Moderate" = "darkgray" , "Liberal" = "blue" ) , name = "Key") +
  scale_x_discrete(labels = xLabel) +
  geom_line(size = 1.5  )  + #Create line to connect points and generate random distance to space out lines
  geom_point(size = 1.5 ) + #change size and type of point with randome distance + 
  facet_wrap( ~ policy , labeller = labeller(policy = facLabelAntiDecarb), ncol = 3) +
  theme_bw() +
  theme(legend.position = "bottom") +
  ggtitle(label = "Anti-decarbonization policies") + 
  xlab("Support or opposition") +
  ylab("Percentage of participants") +
  theme(plot.title = element_text(size = 18 , hjust = 0.5 , margin = margin(t = 0 , r = 5 , b = 20 , l = 5) )) +
  theme(axis.text.x = element_text( size = 12 , hjust = 0.5 ) )  +
  theme(axis.text.y = element_text( size = 14 , hjust = 0.5)) +
  theme(legend.text=element_text(size=14)) + 
  theme(legend.title = element_text(size = 14 , margin = margin( r = 2 , l = 0.01))) + 
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 16 , hjust = 0.5 )) + #legend title settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(strip.text = element_text(size = 16) ) +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm') )) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) 

#Combine Decarb and AntiDecarb figures
library(ggpubr)
figure <- ggarrange(Freq.plot.Decarb , Freq.plot.AntiDecarb,
                    ncol = 1 ) %>%
  ggexport(filename = "C:/Users/Owner/Desktop/figure.pdf" , height = 22 , width = 18 , unit = "cm")
```


#Issue Importance All participants
```{r}
#Convert data to long format
EnergyMix4 <- EnergyMix %>%
  dplyr::select(c(ID,issImpUS,issImpGlobal,issImpUSFut,issImpGlFut,polIdeoGroup))  %>%
  group_by(ID) %>%
  gather(key = "factorLevel" , value = "issue" , -c(ID,polIdeoGroup))

#Recode issue to category name and set order
EnergyMix4$issue <- dplyr::recode( EnergyMix4$issue, "1" = "Economy and jobs" , "2" = "Climate change" , "3" = "Immigration" , "4" = "Access to quality health care")
EnergyMix4$issue <- factor(EnergyMix4$issue , levels = c("Climate change" , "Economy and jobs" , "Access to quality health care" , "Immigration"))

EnergyMix4$focalPoint <- dplyr::recode(EnergyMix4$factorLevel , "issImpUS" = "U.S. today" , "issImpGlobal" = "world today" , "issImpUSFut" = "U.S. future" , "issImpGlFut" = "world future")

#obtain count of participants
FreqIss <- xtabs(~focalPoint + issue , data = EnergyMix4)
FreqIss

#turn frequency count into data frame
FreqIssDF <- data.frame(FreqIss)
head(FreqIssDF)

#convert data into wide format
FreqIssDF <- spread(FreqIssDF , issue , Freq)

#convert count to percentage
FreqIssDF$cc <- (FreqIssDF$`Climate change` / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100
FreqIssDF$econ <- (FreqIssDF$`Economy and jobs` / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100
FreqIssDF$health <- (FreqIssDF$`Access to quality health care` / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100
FreqIssDF$imm <- (FreqIssDF$Immigration / (FreqIssDF$`Climate change` + FreqIssDF$`Economy and jobs` + FreqIssDF$`Access to quality health care` + FreqIssDF$Immigration)) * 100

#convert back into long format
FreqIssLong <- FreqIssDF %>%
  dplyr::select(c(focalPoint, cc, econ, health, imm)) %>%
  gather(key = "percent" , value = "percentScore" , -c(focalPoint))

#Rename issues
FreqIssLong$percent <- dplyr::recode( FreqIssLong$percent, "cc" = "Climate change" , "econ" = "Economy and jobs" , "health" = "Access to quality health care" , "imm" = "Immigration")

#set order
FreqIssLong$percent <- factor(FreqIssLong$percent , levels = c("Climate change" , "Economy and jobs" , "Access to quality health care" , "Immigration"))
FreqIssLong$focalPoint <- factor(FreqIssLong$focalPoint , levels = c("U.S. today" , "world today" , "U.S. future" , "world future"))

library(viridis)
all_plot_Per.plot <- ggplot(FreqIssLong , aes(x= focalPoint , y = percentScore, colour = percent, fill = percent)) +
  geom_col(position = "dodge") +
  scale_colour_viridis(discrete = TRUE , name = "Key") + ##set colors for Ideology
    scale_fill_viridis(discrete = TRUE , name = "Key") + ## remove option = "C" for basic viridis
  theme_bw() +
  theme(legend.position="bottom") +
  ggtitle(label = "Most important issue across spatial and temporal scales") +
    xlab("Focal point") + #x label
  ylab("Percentage of participants") + #y label
  theme(plot.title = element_text(size = 18, hjust = 0.5 , margin = margin(t = 15 , r = 5 , b = 20 , l = 5))) +
  theme(axis.text.x = element_text( size = 16 , hjust = 0.5)) + #x axis settings
  theme(axis.text.y = element_text( size = 16 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=14.5)) + #legend settings
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 14 , hjust = 0.5 )) + #legend title settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(plot.subtitle = element_text(margin = margin(b = 10))) + #plot subtitle margin
  theme(legend.background = element_rect(color = NA)) +
  theme(legend.text = element_text(margin = margin(r = 0.3 , l = 0.01 , unit = 'cm'))) +
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) +
  theme(legend.key.width = unit(.3, "cm")) +
  theme(legend.key.height = unit(.2,"cm"))
 # theme(legend.justification = "right")

ggsave("C:/Users/Owner/Desktop/all_plot_Per.pdf" , all_plot_Per.plot , height = 10, width = 10)
```


#Issue importance by political ideology
```{r}
#Frequency table by issue and political ideology
FreqIssPol <- xtabs(~focalPoint + polIdeoGroup + issue, data = EnergyMix4)
FreqIssPol

#convert to data frame
FreqIsspolDF <- data.frame(FreqIssPol)

#convert data frame to wide format
FreqIsspolWide <- spread(FreqIsspolDF , issue , Freq)

#convert count to percentage
FreqIsspolWide$cc <- (FreqIsspolWide$`Climate change` / (FreqIsspolWide$`Climate change` + FreqIsspolWide$`Economy and jobs` + FreqIsspolWide$`Access to quality health care` + FreqIsspolWide$Immigration)) * 100
FreqIsspolWide$econ <- (FreqIsspolWide$`Economy and jobs` / (FreqIsspolWide$`Climate change` + FreqIsspolWide$`Economy and jobs` + FreqIsspolWide$`Access to quality health care` + FreqIsspolWide$Immigration)) * 100
FreqIsspolWide$health <- (FreqIsspolWide$`Access to quality health care` / (FreqIsspolWide$`Climate change` + FreqIsspolWide$`Economy and jobs` + FreqIsspolWide$`Access to quality health care` + FreqIsspolWide$Immigration)) * 100
FreqIsspolWide$imm <- (FreqIsspolWide$Immigration / (FreqIsspolWide$`Climate change` + FreqIsspolWide$`Economy and jobs` + FreqIsspolWide$`Access to quality health care` + FreqIsspolWide$Immigration)) * 100

#Convert back into long format
FreqIssPolLong <- FreqIsspolWide %>%
  dplyr::select(c(focalPoint, polIdeoGroup , cc, econ, health, imm)) %>%
  gather(key = "issue" , "percent" , -c(focalPoint, polIdeoGroup))

#Rename issues
FreqIssPolLong$issue <- dplyr::recode( FreqIssPolLong$issue, "cc" = "climate change" , "econ" = "economy and jobs" , "health" = "access to quality health care" , "imm" = "immigration")
head(FreqIssPolLong)
xLabel <- c("U.S.\ntoday",
            "world\ntoday",
            "U.S.\nfuture",
            "world\nfuture")

#set order
FreqIssPolLong$issue <- factor(FreqIssPolLong$issue , levels = c("climate change" , "economy and jobs" , "access to quality health care" , "immigration"))
FreqIssPolLong$focalPoint <- factor(FreqIssPolLong$focalPoint , levels = c("U.S. today" , "world today" , "U.S. future" , "world future"))
FreqIssPolLong$polIdeoGroup <- factor(FreqIssPolLong$polIdeoGroup , levels = c("Conservative" , "Moderate", "Liberal"))

dotplot.plot <- ggplot(FreqIssPolLong , aes(x= focalPoint , y = percent , group= polIdeoGroup , colour = polIdeoGroup, fill = polIdeoGroup)) +
  scale_colour_manual(values = c("red", "darkgray", "blue") , name = "Key") + 
  scale_fill_manual(values = c("red", "darkgray","blue") , name = "Key") +
  geom_point(size = 3 , position = position_dodge(width = .5)) +
  #geom_errorbar(aes(ymin=Pct-SE_P, ymax=Pct+SE_P) , width=.20 , position = position_dodge(width=.5)) +
  facet_wrap( ~FreqIssPolLong$issue , ncol = 2 ) + 
  scale_x_discrete(labels = xLabel) +
  theme_bw() +
  theme(legend.position="bottom") +
  ggtitle(label = "Most important issue by political ideology across spatial and temporal scales") +
   ylim(0,80) + 
    xlab("Focal point") + 
  ylab("Percentage of participants") +
  theme(plot.title = element_text(size = 18 , hjust = 0.5 , margin = margin(t = 15 , r = 5 , b = 20 , l = 5))) + #plot title settings
  theme(axis.text.x = element_text( size = 14 , hjust = 0.5)) + #x axis settings
  theme(axis.text.y = element_text( size = 14 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=14)) + #legend settings
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 14 , hjust = 0.5 )) + #legend title settings
  #theme(strip.text.x = element_text(size = 16 )) + #plot subtitle settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(plot.title = element_text(margin = margin(t = 10 , r = 5 , b = 25 , l = 5))) + #Plot title margin
  theme(strip.text = element_text(size = 16) ) +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm')))+
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) 
  
ggsave("C:/Users/Owner/Desktop/dotplot.pdf" , dotplot.plot , height = 10, width = 10)
```


##Supplemental Text
#Descriptive statistics, Political ideology and political party and policy score
```{r}
#Descriptive statistics on continuous variables and constructed variables
library(psych)
describe(EnergyMix$polScoreNuc)
describe(EnergyMix$ccAttitudeBase)
describe(EnergyMixRegData$age)
describe(EnergyMix$Re_income)
describe(EnergyMix$ccBelief)
describe(EnergyMix$ccImportant)
describe(EnergyMix$ccBehavior)
describe(EnergyMix$ccImpUSRel)
describe(EnergyMix$DecarbFutCur)
describe(EnergyMixRegData$polScoreNuc)

library(gmodels)
#tabulate political ideoogy by political party affiliation
CrossTable(EnergyMix$politicalBelief, EnergyMix$politicalParty)

#Tabulate political ideology for each policy: Examples of code that can be used
CrossTable(EnergyMix$politicalBelief , EnergyMix$pp_carbonTax)
xtabs(EnergyMix$pp_carbonTax ~ EnergyMix$politicalBelief, EnergyMix)

#Look at mean policy score by policy by political ideology
tapply(EnergyMix$pp_carbonTax , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_coal , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_constructPipe , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_fuelEcon , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_fundBattery , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_fundRenew , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_haltPipe , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_lowerSubs , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_nucCoal , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_relaxReg, EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_subsidCars , EnergyMix$polIdeoGroup , FUN=mean)
tapply(EnergyMix$pp_tariffs , EnergyMix$polIdeoGroup , FUN=mean)
```

#SI text Table S5
```{r}
#Summary statistics for current energy source percent estimates
mean(EnergyMix$normc_oil)
mean(EnergyMix$normc_ng)
mean(EnergyMix$normc_coal)
mean(EnergyMix$normc_nuclear)
mean(EnergyMix$normc_biomass)
mean(EnergyMix$normc_hydro)
mean(EnergyMix$normc_wind)
mean(EnergyMix$normc_solar)
mean(EnergyMix$normc_geotherm)

#mean statistics for future energy source percent estimates
mean(EnergyMix$normf_oil , na.rm = TRUE)
mean(EnergyMix$normf_ng , na.rm = TRUE)
mean(EnergyMix$normf_coal , na.rm = TRUE)
mean(EnergyMix$normf_nuclear, na.rm = TRUE ) 
mean(EnergyMix$normf_biomass , na.rm = TRUE)
mean(EnergyMix$normf_hydro , na.rm = TRUE)
mean(EnergyMix$normf_wind, na.rm = TRUE)
mean(EnergyMix$normf_solar, na.rm = TRUE)
mean(EnergyMix$normf_geotherm, na.rm = TRUE)

##Mean Estimates for current energy mix estimates by political ideology
##Create different data frames for liberals, moderates, and conservatives

#Liberal dataframe
EnergyMixLiberal <- subset(EnergyMix, polIdeoGroup %in% c("Liberal"))

#Conservative dataframe
EnergyMixConservative <- subset(EnergyMix, polIdeoGroup %in% c("Conservative"))

#Moderate dataframe
EnergyMixModerate <- subset(EnergyMix, polIdeoGroup %in% c("Moderate"))

#mean estimates for normalized current mix for liberals
mean(EnergyMixLiberal$normc_oil)
mean(EnergyMixLiberal$normc_coal)
mean(EnergyMixLiberal$normc_ng)
mean(EnergyMixLiberal$normc_nuclear)
mean(EnergyMixLiberal$normc_biomass)
mean(EnergyMixLiberal$normc_wind)
mean(EnergyMixLiberal$normc_hydro)
mean(EnergyMixLiberal$normc_solar)
mean(EnergyMixLiberal$normc_geotherm)

#mean estimates for normalized current mix for moderates
mean(EnergyMixModerate$normc_oil)
mean(EnergyMixModerate$normc_ng)
mean(EnergyMixModerate$normc_coal)
mean(EnergyMixModerate$normc_nuclear)
mean(EnergyMixModerate$normc_biomass)
mean(EnergyMixModerate$normc_hydro)
mean(EnergyMixModerate$normc_wind)
mean(EnergyMixModerate$normc_solar)
mean(EnergyMixModerate$normc_geotherm)

#mean estimates for normalized current mix for conservatives
mean(EnergyMixConservative$normc_oil)
mean(EnergyMixConservative$normc_ng)
mean(EnergyMixConservative$normc_coal)
mean(EnergyMixConservative$normc_nuclear)
mean(EnergyMixConservative$normc_biomass)
mean(EnergyMixConservative$normc_hydro)
mean(EnergyMixConservative$normc_wind)
mean(EnergyMixConservative$normc_solar)
mean(EnergyMixConservative$normc_geotherm)

## Mean Estimates for future energy mix estimates by political ideology
#mean estimates for normalized future mix for liberals
mean(EnergyMixLiberal$normf_oil , na.rm = TRUE)
mean(EnergyMixLiberal$normf_ng , na.rm = TRUE)
mean(EnergyMixLiberal$normf_coal , na.rm = TRUE)
mean(EnergyMixLiberal$normf_nuclear , na.rm = TRUE)
mean(EnergyMixLiberal$normf_biomass , na.rm = TRUE)
mean(EnergyMixLiberal$normf_hydro , na.rm = TRUE)
mean(EnergyMixLiberal$normf_wind , na.rm = TRUE)
mean(EnergyMixLiberal$normf_solar , na.rm = TRUE)
mean(EnergyMixLiberal$normf_geotherm , na.rm = TRUE)

#mean estimates for normalized future mix for moderates
mean(EnergyMixModerate$normf_oil)
mean(EnergyMixModerate$normf_ng)
mean(EnergyMixModerate$normf_coal)
mean(EnergyMixModerate$normf_nuclear)
mean(EnergyMixModerate$normf_biomass)
mean(EnergyMixModerate$normf_hydro)
mean(EnergyMixModerate$normf_wind)
mean(EnergyMixModerate$normf_solar)
mean(EnergyMixModerate$normf_geotherm)

#mean estimates for normalized future mix for conservatives
mean(EnergyMixConservative$normf_oil)
mean(EnergyMixConservative$normf_ng)
mean(EnergyMixConservative$normf_coal)
mean(EnergyMixConservative$normf_nuclear)
mean(EnergyMixConservative$normf_biomass)
mean(EnergyMixConservative$normf_hydro)
mean(EnergyMixConservative$normf_wind)
mean(EnergyMixConservative$normf_solar)
mean(EnergyMixConservative$normf_geotherm)
```

#SI Table S6
```{r}
#Current energy mix estimate confidence intervals for all participants
ci(EnergyMix$normc_biomass , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_oil , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_ng , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_coal , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_solar , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_wind , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_hydro , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_geotherm , confidence = 0.95 , alpha = 0.05)
ci(EnergyMix$normc_nuclear , confidence = 0.95 , alpha = 0.05)

#Confidence intervals for current energy mix estimates by political ideology
#Create data frame of confidence intervals for current estimates only
enMixCIdf <- data.frame(enMixSum)
enMixCI <- enMixCIdf %>%
  dplyr::select(c(energySource1 , polIdeoGroup , timeEst , energyEstimate, ci))

#select only current energy mix estimates
enMixCI <- subset(enMixCI , timeEst != "future")

#set lower and upper bounds as minus and plus CI
enMixCI$ciLower <- (enMixCI$energyEstimate - enMixCI$ci)
enMixCI$ciUpper <- (enMixCI$energyEstimate + enMixCI$ci)

#View data frame with confidence intervals
enMixCI <- enMixCI %>%
  dplyr::select(c(energySource1 , polIdeoGroup , energyEstimate , ciLower , ciUpper))
enMixCI
```

#SI Table S7
```{r}
#Differences between future and current energy mix estiamtes by political ideology
#all participants
mean(EnergyMix$normf_oil, na.rm = TRUE) - mean(EnergyMix$normc_oil)
mean(EnergyMix$normf_ng, na.rm = TRUE) - mean(EnergyMix$normc_ng)
mean(EnergyMix$normf_coal, na.rm = TRUE) - mean(EnergyMix$normc_coal)
mean(EnergyMix$normf_nuclear, na.rm = TRUE) - mean(EnergyMix$normc_nuclear)
mean(EnergyMix$normf_biomass, na.rm = TRUE) - mean(EnergyMix$normc_biomass)
mean(EnergyMix$normf_hydro, na.rm = TRUE) - mean(EnergyMix$normc_hydro)
mean(EnergyMix$normf_wind, na.rm = TRUE) - mean(EnergyMix$normc_wind)
mean(EnergyMix$normf_solar, na.rm = TRUE) - mean(EnergyMix$normc_solar)
mean(EnergyMix$normf_geotherm, na.rm = TRUE) - mean(EnergyMix$normc_geotherm)

#conservatives
mean(EnMix3Con$normf_oil , na.rm = TRUE) - mean(EnMix3Con$normc_oil )
mean(EnMix3Con$normf_ng , na.rm = TRUE) - mean(EnMix3Con$normc_ng )
mean(EnMix3Con$normf_coal , na.rm = TRUE) - mean(EnMix3Con$normc_coal )
mean(EnMix3Con$normf_nuclear , na.rm = TRUE) - mean(EnMix3Con$normc_nuclear )
mean(EnMix3Con$normf_biomass , na.rm = TRUE) - mean(EnMix3Con$normc_biomass )
mean(EnMix3Con$normf_hydro , na.rm = TRUE) - mean(EnMix3Con$normc_hydro )
mean(EnMix3Con$normf_wind , na.rm = TRUE) - mean(EnMix3Con$normc_wind )
mean(EnMix3Con$normf_solar , na.rm = TRUE) - mean(EnMix3Con$normc_solar )
mean(EnMix3Con$normf_geotherm , na.rm = TRUE) - mean(EnMix3Con$normc_geotherm )

#moderates
mean(EnMix3Mod$normf_oil , na.rm = TRUE) - mean(EnMix3Mod$normc_oil )
mean(EnMix3Mod$normf_ng , na.rm = TRUE) - mean(EnMix3Mod$normc_ng )
mean(EnMix3Mod$normf_coal , na.rm = TRUE) - mean(EnMix3Mod$normc_coal )
mean(EnMix3Mod$normf_nuclear , na.rm = TRUE) - mean(EnMix3Mod$normc_nuclear )
mean(EnMix3Mod$normf_biomass , na.rm = TRUE) - mean(EnMix3Mod$normc_biomass )
mean(EnMix3Mod$normf_hydro , na.rm = TRUE) - mean(EnMix3Mod$normc_hydro )
mean(EnMix3Mod$normf_wind , na.rm = TRUE) - mean(EnMix3Mod$normc_wind )
mean(EnMix3Mod$normf_solar , na.rm = TRUE) - mean(EnMix3Mod$normc_solar )
mean(EnMix3Mod$normf_geotherm , na.rm = TRUE) - mean(EnMix3Mod$normc_geotherm )

#liberals
mean(EnMix3Lib$normf_oil , na.rm = TRUE) - mean(EnMix3Lib$normc_oil )
mean(EnMix3Lib$normf_ng , na.rm = TRUE) - mean(EnMix3Lib$normc_ng )
mean(EnMix3Lib$normf_coal , na.rm = TRUE) - mean(EnMix3Lib$normc_coal )
mean(EnMix3Lib$normf_nuclear , na.rm = TRUE) - mean(EnMix3Lib$normc_nuclear )
mean(EnMix3Lib$normf_biomass , na.rm = TRUE) - mean(EnMix3Lib$normc_biomass )
mean(EnMix3Lib$normf_hydro , na.rm = TRUE) - mean(EnMix3Lib$normc_hydro )
mean(EnMix3Lib$normf_wind , na.rm = TRUE) - mean(EnMix3Lib$normc_wind )
mean(EnMix3Lib$normf_solar , na.rm = TRUE) - mean(EnMix3Lib$normc_solar )
mean(EnMix3Lib$normf_geotherm , na.rm = TRUE) - mean(EnMix3Lib$normc_geotherm )
```

#SI text S8
#Mean policy score and standard error by political ideology for each policy 
```{r}
describe(EnMix3Lib$pp_carbonTax)
describe(EnMix3Mod$pp_carbonTax)
describe(EnMix3Con$pp_carbonTax)

describe(EnMix3Lib$pp_fundBattery)
describe(EnMix3Mod$pp_fundBattery)
describe(EnMix3Con$pp_fundBattery)

describe(EnMix3Lib$pp_fundRenew)
describe(EnMix3Mod$pp_fundRenew)
describe(EnMix3Con$pp_fundRenew)

describe(EnMix3Lib$pp_haltPipe)
describe(EnMix3Mod$pp_haltPipe)
describe(EnMix3Con$pp_haltPipe)

describe(EnMix3Lib$pp_subsidCars)
describe(EnMix3Mod$pp_subsidCars)
describe(EnMix3Con$pp_subsidCars)

describe(EnMix3Lib$pp_nucCoal)
describe(EnMix3Mod$pp_nucCoal)
describe(EnMix3Con$pp_nucCoal)

describe(EnMix3Lib$pp_lowerSubs)
describe(EnMix3Mod$pp_lowerSubs)
describe(EnMix3Con$pp_lowerSubs)

describe(EnMix3Lib$pp_relaxReg)
describe(EnMix3Mod$pp_relaxReg)
describe(EnMix3Con$pp_relaxReg)

describe(EnMix3Lib$pp_tariffs)
describe(EnMix3Mod$pp_tariffs)
describe(EnMix3Con$pp_tariffs)

describe(EnMix3Lib$pp_coal)
describe(EnMix3Mod$pp_coal)
describe(EnMix3Con$pp_coal)

describe(EnMix3Lib$pp_constructPipe)
describe(EnMix3Mod$pp_constructPipe)
describe(EnMix3Con$pp_constructPipe)

describe(EnMix3Lib$pp_fuelEcon)
describe(EnMix3Mod$pp_fuelEcon)
describe(EnMix3Con$pp_fuelEcon)
```
#SI text Section 4
#Anovas for climate change importance ratings
```{r}
#alpha level for anova
0.05/4
#0.0125

#alpha level for t tests
0.05/12
#0.004

#US today
usToday <- aov(EnergyMix$US_cc ~ EnergyMix$polIdeoGroup)
summary(usToday)
t.test(EnMix3Con$US_cc , EnMix3Lib$US_cc)
t.test(EnMix3Lib$US_cc , EnMix3Mod$US_cc)
t.test(EnMix3Con$US_cc , EnMix3Mod$US_cc)

#World today
worldToday <- aov(EnergyMix$global_cc ~ EnergyMix$polIdeoGroup)
summary(worldToday)
t.test(EnMix3Con$global_cc , EnMix3Lib$global_cc)
t.test(EnMix3Lib$global_cc , EnMix3Mod$global_cc)
t.test(EnMix3Con$global_cc , EnMix3Mod$global_cc)

#US future
usFut <- aov(EnergyMix$USFut_cc ~ EnergyMix$polIdeoGroup)
summary(usFut)
t.test(EnMix3Con$USFut_cc , EnMix3Lib$USFut_cc)
t.test(EnMix3Lib$USFut_cc , EnMix3Mod$USFut_cc)
t.test(EnMix3Con$USFut_cc , EnMix3Mod$USFut_cc)

#world future
worldFut <- aov(EnergyMix$glFut_cc ~ EnergyMix$polIdeoGroup)
summary(worldFut)
t.test(EnMix3Con$glFut_cc , EnMix3Lib$glFut_cc)
t.test(EnMix3Lib$glFut_cc , EnMix3Mod$glFut_cc)
t.test(EnMix3Con$glFut_cc , EnMix3Mod$glFut_cc)
```

#SI text section 5
#Regression Model Diagnostics: For SI text REMOVE CC BEHAVIOR
```{r}
#Checking Regression Model Assumptions
lmDecarbPolFull <- lm(polScoreNuc ~ DecarbFutCur +
                        politicalBelief + 
                        Re_gender +
                        Re_income +
                        age +
                        education +
                        ccAttitudeBase + 
                        ccImpUSRel ,
                      data = EnergyMixRegData )
summary(lmDecarbPolFull)
library(gvlma)
summary(lmDecarbPolFull)

#Normality
plot(lmDecarbPolFull$fitted.values , lmDecarbPolFull$residuals)
plot(lmDecarbPolFull$residuals)
hist(lmDecarbPolFull$residuals)
skew(lmDecarbPolFull$residuals)
describe(lmDecarbPolFull$residuals)

#Linearity
plot(lmDecarbPolFull$fitted.values , lmDecarbPolFull$residuals)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$DecarbFutCur)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$politicalBelief)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$Re_gender)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$education)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$age)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$Re_income)
plot(EnergyMixRegData$polScoreNuc , EnergyMixRegData$ccImpUSRel)
modelVar <- data.frame(EnergyMixRegData$polScoreNuc,
                       EnergyMixRegData$DecarbFutCur,
                       EnergyMixRegData$politicalBelief,
                       EnergyMixRegData$education,
                       EnergyMixRegData$age,
                       EnergyMixRegData$Re_income,
                       EnergyMixRegData$ccAttitudeBase,
                       EnergyMixRegData$ccImpUSRel)
corModelVar<-cor(modelVar , use = "complete.obs")
round(corModelVar , digits = 2)
EnergyMixRegData$Re_gender_numer <- as.numeric(EnergyMixRegData$Re_gender)
cor(EnergyMixRegData$polScoreNuc , EnergyMixRegData$Re_gender_numer)

##Homoscedasticity
plot(lmDecarbPolFull$fitted.values , lmDecarbPolFull$residuals)

#Indepdence
plot(lmDecarbPolFull$residuals)
library(lmtest)
dwtest(lmDecarbPolFull)

##Multicolinearity
#Bivariate distributions
IVvar <- data.frame(EnergyMixRegData$Re_gender_numer,
                       EnergyMixRegData$DecarbFutCur,
                       EnergyMixRegData$politicalBelief,
                       EnergyMixRegData$education,
                       EnergyMixRegData$age,
                       EnergyMixRegData$Re_income,
                       EnergyMixRegData$ccAttitudeBase,
                       EnergyMixRegData$ccImpUSRel)
corIVvar<-cor(IVvar , use = "complete.obs")
round(corIVvar , digits = 2)

#Multicollinearity
#Multivariate
smc(corModelVar)
1-smc(corModelVar)
1/(1-smc(corModelVar))
library(HH)
vif(lmDecarbPolFull )

#Outliers
#Univariate
EnergyMixRegData$St_polScoreNuc <- scale(EnergyMixRegData$polScoreNuc)
EnergyMixRegData$St_DecarbFutCur <- scale(EnergyMixRegData$DecarbFutCur)
EnergyMixRegData$St_politicalBelief <- scale(EnergyMixRegData$politicalBelief)
EnergyMixRegData$St_education <- scale(EnergyMixRegData$education)
EnergyMixRegData$St_age <- scale(EnergyMixRegData$age)
EnergyMixRegData$St_Re_income <- scale(EnergyMixRegData$Re_income)
EnergyMixRegData$St_ccAttitude <- scale(EnergyMixRegData$ccAttitudeBase)
EnergyMixRegData$St_ccImpRel <- scale(EnergyMixRegData$ccImpUSRel)
EnergyMixRegData$St_gender <- scale(EnergyMixRegData$Re_gender_numer)
#----
boxplot(EnergyMixRegData$St_polScoreNuc , main = "Policy Score")
boxplot(EnergyMixRegData$St_DecarbFutCur , main = "Decarb Score")
boxplot(EnergyMixRegData$St_politicalBelief , main = "Political Belief")
boxplot(EnergyMixRegData$St_education , main = "Education")
boxplot(EnergyMixRegData$St_age , main = "Age" )
boxplot(EnergyMixRegData$St_Re_income , main = "Income")
boxplot(EnergyMixRegData$St_ccAttitude , main = "Climate Change Attitude")
boxplot(EnergyMixRegData$St_ccImpRel , main = "Climate Change Importance Relative")
boxplot(EnergyMixRegData$St_gender , main = "Gender")

#Outliers
#Bivariate
library(graphics)
library(graphicsQC)
library(car)
pairs.plot <- pairs(corModelVar)

#Outliers
#Multivariate
(cutoff = 4/nrow(EnergyMixRegData))
plot(lmDecarbPolFull , which = 4)
abline(h=cutoff , col = "red" , lty = 2)
#---
EnergyMixRegData$MD_polScoreNuc <- (EnergyMixRegData$polScoreNuc)
EnergyMixRegDataRegDataMD_DecarbFutCur <- (EnergyMixRegData$DecarbFutCur)
EnergyMixRegData$MD_politicalBelief <- (EnergyMixRegData$politicalBelief)
EnergyMixRegData$MD_education <- (EnergyMixRegData$education)
EnergyMixRegData$MD_age <- (EnergyMixRegData$age)
EnergyMixRegData$MD_Re_income <- (EnergyMixRegData$Re_income)
EnergyMixRegData$MD_ccAttitude <- (EnergyMixRegData$ccAttitudeBase)
EnergyMixRegData$MD_ccImpRel <- (EnergyMixRegData$ccImpUSRel)
EnergyMixRegData$MD_gender <- (EnergyMixRegData$Re_gender_numer)
head(EnergyMixRegData)

##Other checks
summary(gvlma(lmDecarbPolFull))

#Running full regression model with log of income + 1 to check for univariate outlier effects
EnergyMixRegData$logIncome <- log(EnergyMixRegData$income + 1)
hist(EnergyMixRegData$logIncome)
range(EnergyMixRegData$logIncome)

lmFinalLogIncome <- lm( EnergyMixRegData$polScoreNuc ~ EnergyMixRegData$DecarbFutCur + 
                          EnergyMixRegData$politicalBelief + 
                          EnergyMixRegData$Re_gender + 
                          EnergyMixRegData$education + 
                          EnergyMixRegData$age + 
                          EnergyMixRegData$logIncome + 
                          EnergyMixRegData$ccAttitudeBase + 
                          EnergyMixRegData$ccBehavior + 
                          EnergyMixRegData$ccImpUSRel, na.action=na.omit)
summary(lmFinalLogIncome)
```

#SI text Section 7, Table S9
#Regression Models: Excluding Nuclear energy from full model and Predicting support for nuclear
```{r}
###Excluding nuclear from final regression model

#Sum of future renewable energy resource estimates
EnergyMixRegData$RenewFut <- (EnergyMixRegData$normf_biomass + 
                          EnergyMixRegData$normf_geotherm + 
                          EnergyMixRegData$normf_hydro + 
                          EnergyMixRegData$normf_solar + 
                          EnergyMixRegData$normf_wind)

#Sum of current renewable energy resource estimates
EnergyMixRegData$RenewCur <- (EnergyMixRegData$normc_biomass + 
                          EnergyMixRegData$normc_geotherm + 
                          EnergyMixRegData$normc_hydro + 
                          EnergyMixRegData$normc_solar + 
                          EnergyMixRegData$normc_wind)

#Renewable score
#Future Renewable - current renewable
EnergyMixRegData$RenewFutCur <- (EnergyMixRegData$RenewFut - EnergyMixRegData$RenewCur)

#Nuc only score
EnergyMixRegData$nucFutCur <- (EnergyMixRegData$normf_nuclear - EnergyMixRegData$normc_nuclear)

#Regression model predicting decarbonization policy support with Renewable score as IV
lmRenewNoNuc <- lm( EnergyMixRegData$polScoreNuc ~ EnergyMixRegData$RenewFutCur + 
                      EnergyMixRegData$politicalBelief + 
                      EnergyMixRegData$Re_gender + 
                      EnergyMixRegData$education + 
                      EnergyMixRegData$age + 
                      EnergyMixRegData$Re_income + 
                      EnergyMixRegData$ccAttitudeBase +
                      EnergyMixRegData$ccImpUSRel)
summary(lmRenewNoNuc)
```

#SI text Section 8, table S10 & S11
##Weighting on only political ideology REMOVING CC BEHAVIOR
```{r}
library(survey)
#indicate survey design
survey.design <- svydesign(ids=~1 , data = EnergyMixRegData)

#set population parameters for political ideology (from GSS)
politicalBelief.dist <- data.frame(politicalBelief = c("1" , "2" , "3" , "4" , "5" , "6" , "7" ) ,
                                   Freq = nrow(EnergyMixRegData)*c(0.054, 0.124, 0.114, 0.381, 0.126, 0.158, 0.044))

#create weights
rake_design <- rake(design = survey.design , sample.margins = list(~politicalBelief) , 
                    population.margins = list(politicalBelief.dist))
#look at weight value distribution
summary(weights(rake_design))
hist(weights(rake_design))

#trim weights
#rake_trim <- trimWeights(rake_design , lower = 0.1 , upper = 50.00 , strict = TRUE)
#distribution of trimmed weights
#summary(weights(rake_trim))
#hist(weights(rake_trim))

#Create weight variable
EnergyMixRegData$PBweight <- weights(rake_design)
summary(EnergyMixRegData$PBweight)
table(EnergyMixRegData$PBweight)

#Re-run regression with weights
lmpolBeliefWeight <- lm(polScoreNuc ~ DecarbFutCur +
                        politicalBelief + 
                        Re_gender +
                        Re_income +
                        age +
                        education +
                        ccAttitudeBase + 
                        ccImpUSRel , 
                      weights = PBweight ,
                      data = EnergyMixRegData )
summary(lmpolBeliefWeight)
```

#SI text Section S9, tables S12 and S13
```{r}
#Models predicting individualized policies
#Carbon Tax
lmCarbonTax <- lm(EnergyMixRegData$pp_carbonTax ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmCarbonTax)

#Replace Coal with Nuclear
lmNucCoal <- lm(EnergyMixRegData$pp_nucCoal ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmNucCoal)

#Fund Renewable Energy
lmFundRenew <- lm(EnergyMixRegData$pp_fundRenew ~ DecarbFutCur + politicalBelief +Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmFundRenew)

#Fund battery storage tech
lmFundBattery <- lm(EnergyMixRegData$pp_fundBattery ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmFundBattery)

#Halt pipelines
lmHaltPipe <- lm(EnergyMixRegData$pp_haltPipe ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmHaltPipe)

#Subsidize electric cars
lmSubCars <- lm(EnergyMixRegData$pp_subsidCars ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmSubCars)

#Lower fuel economy standards
lmfuelEcon <- lm(EnergyMixRegData$pp_fuelEcon ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmfuelEcon)

#Lower Renewable Subsidies
lmLowerSubs <- lm(EnergyMixRegData$pp_lowerSubs ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmLowerSubs)

#Construct oil pipelines
lmConstructPipes <- lm(EnergyMixRegData$pp_constructPipe ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel  , data = EnergyMixRegData )
summary(lmConstructPipes)

#SOlar panel tariffs
lmSolTar <- lm(EnergyMixRegData$pp_tariffs ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmSolTar)

#Relax oil/ng regulations
lmRelaxReg <- lm(EnergyMixRegData$pp_relaxReg ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmRelaxReg)

#Build more coal
lmCoal <- lm(EnergyMixRegData$pp_coal ~ DecarbFutCur + politicalBelief + Re_gender + Re_income + age + education + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData )
summary(lmCoal)
```

#SI text Section 10
#Mediation analysis
```{r}
#Predicting DV (policy score) leaving out mediating variable (decarb score)
lm1 <- lm(polScoreNuc ~ politicalBelief + education + age + Re_income + Re_gender + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData)
summary(lm1)

#Predicting Mediating variable with other IVs 
lm1a <- lm(DecarbFutCur ~ politicalBelief + education + age + Re_income + Re_gender + ccAttitudeBase + ccImpUSRel , data = EnergyMixRegData)
summary(lm1a)

#Full model with mediating variable
lm1b <- lm(polScoreNuc ~ politicalBelief + education + age + Re_income + Re_gender + ccAttitudeBase + ccImpUSRel + DecarbFutCur , data = EnergyMixRegData)
summary(lm1b)
```

##Figure S1: political support by policy and political ideolgoy
##Policy Figure
```{r}
#Convert data into long format
head(EnergyMix)
EnergyMix3 <- EnergyMix %>%
  dplyr::select(c(ID:issImpGlFut, pp_nucCoal:pp_coal, volunt_cc:AntiDecarbPolScore)) %>%
  group_by(ID) %>%
  gather(key = "policy" , value = "policyScore" , -c(ID:issImpGlFut, volunt_cc:AntiDecarbPolScore))

#Indicate policy order
EnergyMix3$policy <- factor(EnergyMix3$policy , levels = c("pp_carbonTax" , "pp_haltPipe" , "pp_nucCoal" , "pp_fundRenew" , "pp_subsidCars" , "pp_fundBattery" , "pp_relaxReg" , "pp_constructPipe" , "pp_coal" , "pp_lowerSubs" , "pp_fuelEcon" , "pp_tariffs"))

enMixPol <- summarySE(EnergyMix3 , measurevar = "policyScore" , na.rm = TRUE , groupvars = c("policy" , "polIdeoGroup"))
enMixPol

#Seperate policies into anti-decarbonization and decarbonization dataframes
policyDecarb <- subset(enMixPol, policy %in% c("pp_carbonTax" , "pp_haltPipe" , "pp_nucCoal" , "pp_fundRenew" , "pp_subsidCars" , "pp_fundBattery"))
policyAntiDecarb <- subset(enMixPol, policy %in% c("pp_relaxReg" , "pp_constructPipe" , "pp_coal", "pp_lowerSubs" , "pp_fuelEcon" , "pp_tariffs"))

#Specify order for policies
policyDecarb$policy <- factor(policyDecarb$policy , levels = c("pp_carbonTax" , "pp_haltPipe" , "pp_nucCoal" , "pp_fundRenew" , "pp_subsidCars" , "pp_fundBattery"))
policyAntiDecarb$policy <- factor(policyAntiDecarb$policy , levels = c("pp_relaxReg" , "pp_constructPipe" , "pp_coal" , "pp_lowerSubs" , "pp_fuelEcon" , "pp_tariffs"))

#Rename Decarbonization policies
facLabelDecarb <- c(pp_nucCoal = "replace coal with nuclear" , pp_carbonTax = "carbon tax" , pp_fundRenew = "fund renewable resources" , pp_subsidCars = "subsidize electric cars" , pp_haltPipe = "halt oil pipelines" , pp_fundBattery = "fund battery storage research")

#Rename anti-decarbonization policies
facLabelAntiDecarb <- c(pp_fuelEcon = "lower fuel economy standards" , pp_lowerSubs = "lower renewable energy subsidies" , pp_constructPipe = "construct oil pipelines" , pp_tariffs = "tax imported solar panels" , pp_relaxReg = "relax drilling regulations" , pp_coal = "invest in new coal plants")

policyDecarb$polIdeoGroup <- factor(policyDecarb$polIdeoGroup , levels = c("Conservative" , "Moderate" , "Liberal"))
policyAntiDecarb$polIdeoGroup <- factor(policyAntiDecarb$polIdeoGroup , levels = c("Conservative" , "Moderate" , "Liberal"))

#Create 2 sperate figures for decarbonization policies and anti-decarbonization policies
library(grid)
library(gridBase)
library(gridGraphics)
library(gridExtra)

#plot1 <- ggplot(policyDecarb , aes(x = [polIdeoGroup , y = policyScore , group = polIdeoGroup , colour = ]))

polDecarb.plot <- ggplot(policyDecarb , aes(x= polIdeoGroup , y = policyScore , group=polIdeoGroup , colour = polIdeoGroup )) +
  scale_colour_manual(values = c("Conservative"="red", "Moderate" = "darkgray" , "Liberal" = "blue" ) , name = "Key") + ##set colors for polIdeoGroup 
  geom_line(size = 3 , position = position_dodge(width = .25) ) + #Create line to connect points and generate random distance to space out lines
  geom_point(size = 3 , position = position_dodge(width = .25)) + #change size and type of point with randome distance
  #geom_errorbar(aes(ymin=policyScore-ci, ymax=policyScore+ci) , width=.20 , position = position_dodge(width=.25)) + #error bars
  scale_y_continuous(limits = c(0.5,5.5) , breaks = c(1,2,3,4,5) , labels = c("Strongly \noppose" , "Oppose" , "Neutral" , "Support" , "Strongly \nsupport")) +
  facet_wrap( ~policy , labeller = labeller(policy = facLabelDecarb) , ncol = 3 ) + #seperate by energy source and indicate order
  ggtitle(label = "Decarbonization policies") + #title & subtitle
  theme_bw() + #plot theme
  theme(legend.position = "none") +
  theme(axis.title.x = element_blank()) +
  theme(plot.title = element_text(size = 18 , hjust = 0.5 , margin = margin(t = 10 , r = 5 , b = 20 , l = 5))) + #plot title settings
  theme(axis.text.x = element_text( size = 12 , hjust = 0.5)) + #x axis settings
  theme(axis.text.y = element_text( size = 12 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=14)) + #legend settings
  theme(legend.title = element_text(size = 14 , margin = margin( r = 2 , l = 0.01))) + #legend title settings
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 , colour = "white")) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(strip.text.x = element_text(size = 16 )) + #plot subtitle settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  #xlab("Political ideology") + #x label
  ylab("Support or opposition") +  #y label 
  #ylim(0.5,5.5) +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm') )) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) 

polAntiDecarb.plot <- ggplot(policyAntiDecarb , aes(x= polIdeoGroup , y = policyScore , group=polIdeoGroup , colour = polIdeoGroup )) +
  scale_colour_manual(values = c("Conservative"="red", "Moderate" = "darkgray" , "Liberal" = "blue" ) , name = "Key") + ##set colors for polIdeoGroup 
  geom_line(size = 3 , position = position_dodge(width = .25) ) + #Create line to connect points and generate random distance to space out lines
  geom_point(size = 3 , position = position_dodge(width = .25)) + #change size and type of point with randome distance
  #geom_errorbar(aes(ymin=policyScore-ci, ymax=policyScore+ci) , width=.20 , position = position_dodge(width=.25)) + #error bars
  scale_y_continuous(limits = c(0.5,5.5) , breaks = c(1,2,3,4,5) , labels = c("Strongly \noppose" , "Oppose" , "Neutral" , "Support" , "Strongly \nsupport")) +
  facet_wrap( ~policy , labeller = labeller(policy = facLabelAntiDecarb) , ncol = 3 ) + #seperate by energy source and indicate order
  ggtitle(label = "Anti-decarbonization policies") + #title & subtitle
  theme_bw() + #plot theme
  theme(legend.position = "bottom") +
  theme(axis.title.x = element_blank()) +
  theme(plot.title = element_text(size = 18 , hjust = 0.5 , margin = margin(t = 10 , r = 5 , b = 20 , l = 5))) + #plot title settings
  theme(axis.text.x = element_text( size = 12 , hjust = 0.5 , angle = 0.4)) + #x axis settings
  theme(axis.text.y = element_text( size = 12 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=14)) + #legend settings
  theme(legend.title = element_text(size = 14 , margin = margin( r = 2 , l = 0.01))) + #legend title settings
  theme(axis.title.x = element_text( size = 16 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 16 , hjust = 0.5 )) + #y axis title settings
  theme(strip.text.x = element_text(size = 16 )) + #plot subtitle settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  xlab("Political ideology") + #x label
  ylab("Support or opposition") +  #y label
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm') )) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))

library(ggpubr)
figure <- ggarrange(polDecarb.plot , polAntiDecarb.plot,
                    ncol=1) %>%
  ggexport(filename = "C:/Users/Owner/Desktop/figure.pdf" , height = 14 , width = 12)
```

#SI Text
#Figures S2 & S3
```{r}
#Bin participants into support or oppsition groups rather than continuous
EnergyMix$Support.or.opposition <- cut(EnergyMix$polScoreNuc , breaks = c(0,1,2,3,4,5))
table(EnergyMix$Support.or.opposition)
EnergyMix$Support.or.opposition = dplyr::recode( EnergyMix$Support.or.opposition, "(0,1]" = "Strongly oppose" , "(1,2]" = "Oppose" , "(2,3]" = "Neutral" , "(3,4]" = "Support" , "(4,5]" = "Strongly support")

#Create data frame
Freq <- xtabs(~EnergyMix$Support.or.opposition)
Freq.df <- data.frame(Freq)
Freq.df$Percentage <- ((Freq.df$Freq/2429)*100)
Freq.df

Freq.df$Support.or.opposition <- factor(Freq.df$EnergyMix.Support.or.opposition , levels = c("Strongly oppose" , "Oppose" , "Neutral" , "Support" , "Strongly support"))

#Figure S2
plot3.plot <- ggplot( Freq.df , aes( x = Support.or.opposition , y = Percentage )) +
  geom_point( size = 3 ) +
  xlab("Support or opposition") +
  ylab("Percentage of participants") +
  ggtitle( label = "Decarbonization policy support across participants") + 
  theme_bw() + #plot theme
  theme(plot.title = element_text(size = 14 , hjust = 0.5 , margin = margin(t = 15 , r = 5 , b = 20 , l = 5))) + 
  theme(axis.text.x = element_text( size = 12 , hjust = 0.5)) + #x axis settings
  theme(axis.text.y = element_text( size = 12 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=12)) + #legend settings
  theme(axis.title.x = element_text( size = 14 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 14 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 14 , hjust = 0.5 )) + #legend title settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(legend.position = "bottom") +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm') )) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) 
ggsave("C:/Users/Owner/Desktop/plot3.pdf" , plot3.plot , height = 7 , width = 8)

#Figure S4:
Freq <- xtabs(~EnergyMix$Support.or.opposition + EnergyMix$polIdeoGroup)
Freq.df <- data.frame(Freq)
FreqWide <- spread( Freq.df , EnergyMix.polIdeoGroup , Freq)
FreqWide$Conservative <- ((FreqWide$Conservative/707)*100)
FreqWide$Liberal <- ((FreqWide$Liberal/1128)*100)
FreqWide$Moderate <- ((FreqWide$Moderate/494)*100)

#convert data back into long format
FreqLong <- FreqWide %>%
  dplyr::select(c(EnergyMix.Support.or.opposition, Conservative:Moderate)) %>%
  gather(key = "polIdeoGroup" , value = "Percentage" , -c(EnergyMix.Support.or.opposition))

FreqLong$EnergyMix.Support.or.opposition <- factor(FreqLong$EnergyMix.Support.or.opposition , levels = c("Strongly oppose" , "Oppose" , "Neutral" , "Support" , "Strongly support"))

plot4.plot <- ggplot( FreqLong , aes( x = EnergyMix.Support.or.opposition , y = Percentage , group = polIdeoGroup , colour = polIdeoGroup)) +
  scale_colour_manual(values = c("Conservative"="red", "Moderate" = "darkgray" , "Liberal" = "blue" ) , name = "Key" )+
  geom_line(size = 2) +
  geom_point(size = 2) +
  xlab("Support or opposition") + 
  ylab("Percentage of participants")+
  ggtitle(label = "Decarbonization policy support by political ideology") +
  theme_bw() +
  theme(plot.title = element_text( size = 14 , hjust = 0.5 , margin = margin(t = 15 , r = 5 , b = 20 , l = 5))) +
  theme(axis.text.x = element_text( size = 12 , hjust = 0.5)) + #x axis settings
  theme(axis.text.y = element_text( size = 12 , hjust = 0.5)) + #legend settings
  theme(legend.text=element_text(size=12)) + #legend settings
  theme(axis.title.x = element_text( size = 14 , hjust = 0.5 )) + #x axis title settings
  theme(axis.title.y = element_text( size = 14 , hjust = 0.5 )) + #y axis title settings
  theme(legend.title = element_text( size = 14 , hjust = 0.5 )) + #legend title settings
  theme(axis.title.x = element_text(margin = margin(t = 20 , r = 0.5 , b = 5 , l = 0.5))) + #x title margin
  theme(axis.title.y = element_text(margin = margin(t = 5 , r = 20 , b = 5, l = 5))) + #y title margin
  theme(legend.position = "bottom") +
  theme(legend.text = element_text(margin = margin(r = 0.75 , l = 0.01 , unit = 'cm') )) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) 
ggsave("C:/Users/Owner/Desktop/plot4.pdf" , plot4.plot , height = 6 , width = 8)

```