Data preparation EVS, ESS and ISSP

In this script, I select, recode and create variables of the European Values Survey, European Social Survey and International Social Survey Program.

rm(list = ls())
library(plyr)
library(dplyr)
library(foreign)
library(tidyverse)
library(labelled)
library(questionr)
library(psych)

European Values Survey

#Load in EVS trend file
evs <- foreign::read.spss("/Users/anuschka/Documents/gesis_dir/evs/ZA7503_v2-0-0.sav", use.value.labels = T,  to.data.frame = T)

# So apparently, in the Netherlands in 1981 and in 2017, the question for the dependent var was not asked. This means that I cannot use those years from the EVS. A pity. 

evs$S003_n <- as.numeric(evs$S003)
evs$S020 <- as.numeric(as.character(evs$S020))

evssel <- evs %>% filter(S003_n == 31)  %>%
                filter (S020 >1981 & S020 <2017) %>%
            select(stdyno_w, S001, S002EVS, S006, S017, S018,S020, B001, X001, X002, X003, X003R, X003R2, X011, X023, X023R, X025, X025A, X025CSEVS, X007, G027A, X028, X047CS, X047_EVS, X047R_EVS, A006, E033)

table(evssel$stdyno_w, evssel$X025CSEVS, useNA = "always")

freq(evssel$X025CSEVS)

attributes(evs$B001)
table(evssel$B001, useNA = "always")

#Recode it, because now a higher number means lower willingness to make sacrifices. I want everything in the direction that a higher score equals more positive attitudes towards the climate. 

evssel$climate <- as.numeric(evssel$B001)
evssel$religious <- as.numeric(evssel$A006)
evssel <- evssel %>% 
     mutate_at(c("climate", "religious"), funs(recode(., `1` = 4, `2` = 3, `3` = 2, `4` = 1)))
evssel$climate5 <- (evssel$climate-1)*(4/3)+1
evssel$religious[is.na(evssel$religious)] <- 88 #11 missings get value 88
evssel$religious <- as.factor(evssel$religious)
evssel$religious <- revalue(evssel$religious, c("1"="Not at all important", "2"="Not very important", "3"="Rather important", "4" = "Very important", "88"="Missing"))

#Independent variables and missings
table(evssel$X047_EVS, useNA = "always") #Scale of household incomes from 1 to 10, with 540 missings (gross)
# Fix income in quartiles 

# For the EVS I only have one variable for income, but of course the quantiles and mean can be different per year. 
# So I have to take that into account. Following code is definitely not the most efficient way to do so, but it works
evssel_w1 <- as.data.frame(evssel[evssel$S002EVS == "1990-1993", ])
evssel_w1$income <- as.numeric(evssel_w1$X047_EVS)
mean(evssel_w1$income, na.rm=T)
evssel_w1$income[is.na(evssel_w1$income)] <- 5.614
evssel_w1$income_quart <- with(evssel_w1, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
evssel_w1$income_quart <- as.numeric(evssel_w1$income_quart)

evssel_w2 <- as.data.frame(evssel[evssel$S002EVS == "1999-2001", ])
evssel_w2$income <- as.numeric(evssel_w2$X047_EVS)
mean(evssel_w2$income, na.rm=T)
evssel_w2$income[is.na(evssel_w2$income)] <- 6.403
evssel_w2$income_quart <- with(evssel_w2, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
evssel_w2$income_quart <- as.numeric(evssel_w2$income_quart)

evssel_w3 <- as.data.frame(evssel[evssel$S002EVS == "2008-2010", ])
evssel_w3$income <- as.numeric(evssel_w3$X047_EVS)
mean(evssel_w3$income, na.rm=T)
evssel_w3$income[is.na(evssel_w3$income)] <- 5.277
evssel_w3$income_quart <- with(evssel_w3, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
evssel_w3$income_quart <- as.numeric(evssel_w3$income_quart)

evssel_w1 <- rbind(evssel_w1, evssel_w2)
evssel_w1 <- rbind(evssel_w1, evssel_w3)
evssel_w1 <- evssel_w1 %>% select(income_quart)
evssel <- cbind(evssel, evssel_w1)


evssel$sex <- revalue(evssel$X001, c("Male"="1", "Female"="2"))
evssel$sex <- evssel$sex - 1
table(evssel$sex)
table(evssel$E033, useNA = "always")
evssel$lrscale <- as.numeric(evssel$E033)

evssel$X007 <- as.numeric(evssel$X007)
evssel$marstat[evssel$X007 <= 2] <- 1 #Living together
evssel$marstat[evssel$X007 > 2] <- 2 #Not living together
evssel$marstat <- as.factor(evssel$marstat)

#You generally start school at 6
evssel$edu_yrs <- as.numeric(as.character(evssel$X023))
evssel$edu_yrs <- evssel$edu_yrs - 6
table(evssel$edu_yrs, useNA = "always")
mean(evssel$edu_yrs, na.rm=T)
evssel$edu_yrs[is.na(evssel$edu_yrs)] <- 13.518
#Education as isced
evssel$isced[evssel$edu_yrs <=4] <- 0
evssel$isced[evssel$edu_yrs > 4 & evssel$edu_yrs <= 6] <- 1
evssel$isced[evssel$edu_yrs > 6 & evssel$edu_yrs <= 10] <- 2
evssel$isced[evssel$edu_yrs > 10 & evssel$edu_yrs <= 13] <- 3
evssel$isced[evssel$edu_yrs > 13 & evssel$edu_yrs <= 15] <- 4
evssel$isced[evssel$edu_yrs > 15 & evssel$edu_yrs <= 18] <- 5
evssel$isced[evssel$edu_yrs > 18] <- 6
table(evssel$isced, useNA = "always")

table(evssel$X003, useNA = "always")
evssel$X003 <- revalue(evssel$X003, c("82 and more (in EVS 2017)" = "82"))
evssel$age <- as.numeric(as.character(evssel$X003))
table(evssel$age, useNA = "always")
# 3 missings on age
mean(evssel$age, na.rm=T)
evssel$age[is.na(evssel$age)] <- 49.122

# Missings can be done faster
mis_vars <- c("lrscale", "edu_yrs") 

for (var in mis_vars) {
  evssel[is.na(evssel[,var]), var] <- mean(evssel[,var], na.rm = TRUE)
}


#Let's first save the prepped data. 
save(evssel, file="./data/final_data/evssel.Rdata")

evssel <- evssel %>% select(stdyno_w, S001, S002EVS, S006, S017, S018, S020, climate5, sex, X002, age, edu_yrs, isced,  religious, lrplace, income, income_quart, marstat, weightvec, X023, X023R, X047CS, X047_EVS, X047R_EVS) %>% 
        dplyr::rename(surveyyear = S020)

save(evssel, file="./data/final_data/evssel.Rdata")

European Social Survey

# Let's do the same with the ESS, for which I will use the 2016 and 2020 wave. 
ess2016 <- foreign::read.spss("/Users/anuschka/Documents/gesis_dir/ess/ESS8e02_2.sav", use.value.labels = T,  to.data.frame = T)

ess16sel <- ess2016 %>% filter(cntry=="Netherlands")  %>%
            select(name, essround, proddate, idno, ccnthum, ccrdprs, wrclmch, edulvlb, edlvenl, eduyrs, mainact, gndr, agea, yrbrn, lrscale, domicil, hinctnta, rlgatnd, marsts, brncntr)


#Climate change caused by natural processes, human activity, or both. 4 people who do not believe are missing
attributes(ess16sel$ccnthum)
table(ess16sel$cause, useNA = "always")
ess16sel$cause <- as.numeric(ess16sel$ccnthum)
ess16sel$cause[ess16sel$cause==6] <- NA

#Personal responsibility to reduce climate change
attributes(ess16sel$ccrdprs)
table(ess16sel$pers_resp, useNA = "always")
ess16sel$pers_resp <- (as.numeric(ess16sel$ccrdprs) - 1)* (4/10) + 1

#How worried about climate change
table(ess16sel$wrclmch, useNA = "always")
ess16sel$worry <- as.numeric(ess16sel$wrclmch)

save(ess16sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/ess16sel.Rdata")
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/ess16sel.Rdata")


#Independent vars
table(ess16sel$eduyrs, useNA = "always")
ess16sel$eduyrs <- as.numeric(as.character(ess16sel$eduyrs))
table(ess16sel$edulvlb, useNA = "always")
ess16sel$edulvlb <- as.numeric(ess16sel$edulvlb)
#Education as isced
ess16sel$isced[ess16sel$edulvlb <=1] <- 0
ess16sel$isced[ess16sel$edulvlb ==2] <- 1
ess16sel$isced[ess16sel$edulvlb > 2 & ess16sel$edulvlb <= 5] <- 2
ess16sel$isced[ess16sel$edulvlb > 5 & ess16sel$edulvlb <= 15] <- 3
ess16sel$isced[ess16sel$edulvlb > 15 & ess16sel$edulvlb <= 20] <- 4
ess16sel$isced[ess16sel$edulvlb > 20 & ess16sel$edulvlb <= 26] <- 5
ess16sel$isced[ess16sel$edulvlb ==27] <- 6
ess16sel$isced[ess16sel$edulvlb ==28] <- NA
table(ess16sel$isced, useNA = "always")

table(ess16sel$gndr, useNA = "always")
ess16sel$sex <- revalue(ess16sel$gndr, c("Male"="1", "Female"="2"))
table(ess16sel$lrscale, useNA = "always")
ess16sel$lrscale <- as.numeric(ess16sel$lrscale)
ess16sel$lrscale <- (ess16sel$lrscale - 1) * (9/10) + 1
attributes(ess16sel$urban)
ess16sel$urban <- revalue(ess16sel$domicil, c("Farm or home in countryside"= "Low urbanity", "Country village"="Low urbanity", "Town or small city"="Medium urbanity", "Suburbs or outskirts of big city"="High urbanity", "A big city"="High urbanity"))
ess16sel$urban <- factor(ess16sel$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)
ess16sel$ch_attend <- as.numeric(ess16sel$rlgatnd)
ess16sel <- ess16sel %>% 
     mutate_at(c("ch_attend"), funs(recode(., `1` = 7, `2` = 6, `3` = 5, `4` = 4, `5` = 3, `6` = 2, `7` =1)))
table(ess16sel$ch_attend, useNA = "always")
ess16sel$income <- as.numeric(ess16sel$hinctnta) #On a scale from 1 to 10, household income

mean(ess16sel$income, na.rm=T)
ess16sel$income[is.na(ess16sel$income)] <- 5.993
ess16sel$income_quart <- with(ess16sel, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
ess16sel$income_quart <- as.numeric(ess16sel$income_quart)

table(ess16sel$marsts, useNA = "always")                    
ess16sel$marsts <- as.numeric(ess16sel$marsts)
ess16sel$marstat[ess16sel$marsts <= 2] <- 1 #Living together
ess16sel$marstat[ess16sel$marsts > 2] <- 2 #Not living together
ess16sel$marstat[is.na(ess16sel$marsts)] <- 3 #821 missings on this variable...
ess16sel$marstat <- as.factor(ess16sel$marstat)
ess16sel$age <- as.numeric(as.character(ess16sel$agea))

#Substitution of missings
lapply(ess16sel, table, useNA = "always")

mis_vars <- c("lrscale") 

for (var in mis_vars) {
  ess16sel[is.na(ess16sel[,var]), var] <- mean(ess16sel[,var], na.rm = TRUE)
}

save(ess16sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/ess16sel.Rdata")


ess2020 <- foreign::read.spss("/Users/anuschka/Documents/gesis_dir/ess/ESS10.sav", use.value.labels = T,  to.data.frame = T)

ess20sel <- ess2020 %>% filter(cntry=="Netherlands")  %>%
            select(name, essround, proddate, idno, ccnthum, ccrdprs, wrclmch, edulvlb, edlvenl, eduyrs,gndr, agea, yrbrn, prtvthnl, lrscale, domicil, hinctnta, rlgatnd, marsts, brncntr)

#Recode the dependent vars just like the previous wave
attributes(ess20sel$ccnthum)
table(ess20sel$cause, useNA = "always")
ess20sel$cause <- as.numeric(ess20sel$ccnthum)
ess20sel$cause[ess20sel$cause==6] <- NA

#Personal responsibility to reduce climate change
attributes(ess20sel$ccrdprs)
table(ess20sel$pers_resp, useNA = "always")
ess20sel$pers_resp <- (as.numeric(ess20sel$ccrdprs) - 1)* (4/10) + 1

#How worried about climate change
table(ess20sel$wrclmch, useNA = "always")
ess20sel$worry <- as.numeric(ess20sel$wrclmch)

save(ess20sel, file="/Users/anuschka/Documents/climatechange/climatechange/ess20sel.Rdata")

# Independent vars 
table(ess20sel$eduyrs, useNA = "always")
ess20sel$eduyrs <- as.numeric(as.character(ess20sel$eduyrs))
table(as.numeric(ess20sel$edulvlb), useNA = "always")
ess20sel$edulvlb <- as.numeric(ess20sel$edulvlb)
ess20sel$isced[ess20sel$edulvlb <=1] <- 0
ess20sel$isced[ess20sel$edulvlb ==2] <- 1
ess20sel$isced[ess20sel$edulvlb > 2 & ess20sel$edulvlb <= 5] <- 2
ess20sel$isced[ess20sel$edulvlb > 5 & ess20sel$edulvlb <= 15] <- 3
ess20sel$isced[ess20sel$edulvlb > 15 & ess20sel$edulvlb <= 20] <- 4
ess20sel$isced[ess20sel$edulvlb > 20 & ess20sel$edulvlb <= 26] <- 5
ess20sel$isced[ess20sel$edulvlb ==27] <- 6
ess20sel$isced[ess20sel$edulvlb ==28] <- NA
table(ess20sel$isced, useNA = "always")

table(ess20sel$gndr, useNA = "always")
ess20sel$sex <- revalue(ess20sel$gndr, c("Male"="1", "Female"="2"))
table(ess20sel$lrscale, useNA = "always")
ess20sel$lrscale <- as.numeric(ess20sel$lrscale)
ess20sel$lrscale <- (ess20sel$lrscale - 1) * (9/10) + 1
attributes(ess20sel$domicil)
ess20sel$urban <- revalue(ess20sel$domicil, c("Farm or home in countryside"= "Low urbanity", "Country village"="Low urbanity", "Town or small city"="Medium urbanity", "Suburbs or outskirts of big city"="High urbanity", "A big city"="High urbanity"))
ess20sel$urban <- factor(ess20sel$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)
ess20sel$ch_attend <- as.numeric(ess20sel$rlgatnd)
ess20sel <- ess20sel %>% 
     mutate_at(c("ch_attend"), funs(recode(., `1` = 7, `2` = 6, `3` = 5, `4` = 4, `5` = 3, `6` = 2, `7` =1)))
table(ess20sel$ch_attend, useNA = "always")
ess20sel$income <- as.numeric(ess20sel$hinctnta) 
table(ess20sel$income, useNA = "always") #152 missings on income.
mean(ess20sel$income, na.rm=T)
ess20sel$income[is.na(ess20sel$income)] <- 6.455
ess20sel$income_quart <- with(ess20sel, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
ess20sel$income_quart <- as.numeric(ess20sel$income_quart)

ess20sel$age <- as.numeric(as.character(ess20sel$agea))
table(ess20sel$marsts, useNA = "always")                    
ess20sel$marsts <- as.numeric(ess20sel$marsts)
ess20sel$marstat[ess20sel$marsts <= 2] <- 1 #Living together
ess20sel$marstat[ess20sel$marsts > 2] <- 2 #Not living together
ess20sel$marstat[is.na(ess20sel$marsts)] <- 3 #732 missings on this variable...
ess20sel$marstat <- as.factor(ess20sel$marstat)

# Missings
lapply(ess20sel, table, useNA = "always")

mis_vars <- c("lrscale", "age") 

for (var in mis_vars) {
  ess20sel[is.na(ess20sel[,var]), var] <- mean(ess20sel[,var], na.rm = TRUE)
}

save(ess20sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/ess20sel.Rdata")
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/ess20sel.Rdata")

#Merge data 2016 and 2020
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/ess16sel.Rdata")

#Select the variables that I want to keep
ess16sel <- ess16sel %>% select(name, essround, proddate, idno, worry, cause, pers_resp, sex, eduyrs, isced,urban, ch_attend, income, income_quart, yrbrn, lrscale, age, marstat, brncntr, weightvec)
                                
ess20sel <- ess20sel %>% select(name, essround, proddate, idno, worry, cause, pers_resp, sex, eduyrs, isced,urban, ch_attend, income, income_quart, yrbrn, lrscale, age, marstat, brncntr, weightvec)

#First check whether the levels/attributes of variables are the same
attributes(ess16sel$worry)
attributes(ess20sel$worry) 
# Also did this for the other variables, they're the same so I can merge

esstotal <- dplyr::bind_rows(ess16sel, ess20sel)

esstotal$surveyyear[esstotal$essround==8] <- 2016
esstotal$surveyyear[esstotal$essround==10] <- 2020
esstotal$surveyyear <- as.numeric(esstotal$surveyyear)
table(esstotal$surveyyear, useNA = "always")
esstotal$sex <- esstotal$sex - 1
table(esstotal$sex)
save(esstotal, file="./data/final_data/esstotal.Rdata")

# One time it did not correctly apply esstotal
#esstotal$urban <- as.factor(esstotal$urban)
#esstotal$urban <- revalue(esstotal$urban, c("1"="Low urbanity", "2"="Medium urbanity", "3"="High urbanity"))
#esstotal$urban <- factor(esstotal$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)

International Social Survey Program

issp1993 <- foreign::read.spss("/Users/anuschka/Documents/gesis_dir/issp/ISSP1993.sav", use.value.labels = T,  to.data.frame = T)

issp93sel <- issp1993 %>% filter(V3=="Netherlands - NL")  %>%
            select(V1, V2, V13, V14, V17, V19, V24, V25, V26, V27, V28, V54, V55, V200, V201, V202, V204, V205, V281, V253, V268, V277, V352, V411)

attributes(issp93sel$V13) #Higher score means more positive (then you think we don't worry too much)
issp93sel$worry <- as.numeric(issp93sel$V13)
attributes(issp93sel$V14) #Higher score means that you don't think that everything in modern life harms the environment. SO I want to recode this variable
attributes(issp93sel$V24) #Willing to pay higher prices in order to protect environment. Very unwilling is negative so I recode this var. Then a higher score means very willing. Same applies for v25 and v26. 
attributes(issp93sel$V28) # I do what is right for the environment, even if it costs a little bit extra time. Strongly agree is more positive attitude, so recode. 
issp93sel$lifeharm <- as.numeric(issp93sel$V14)
issp93sel$willing_price <- as.numeric(issp93sel$V24)
issp93sel$willing_tax <- as.numeric(issp93sel$V25)
issp93sel$willing_living <- as.numeric(issp93sel$V26)
issp93sel$do_right <- as.numeric(issp93sel$V28)

issp93sel <- issp93sel %>% 
     mutate_at(c("lifeharm", "willing_price", "willing_tax", "willing_living", "do_right"), funs(recode(., `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5`= 1)))

attributes(issp93sel$V17) #Higher score means that you don't think we worry too much about human progress harming the environment. A low score means you think we worry too much, so that represents more negative attitudes. No recoding.
issp93sel$progharm <- as.numeric(issp93sel$V17)
attributes(issp93sel$V19) #Difficult one. "In order to protect the environment, we need economic growth". A higher score does not necessarily mean worse or better attitudes. So I leave this one as it is. 
issp93sel$econprotect <- as.numeric(issp93sel$V19)
attributes(issp93sel$V27) #Too difficult for someone like me to do sth about the environment. This also a difficult one. If you are worried but pessimistic, you may believe that you can't solve it on your own. If you then strongly agree, it means you do care about the environment. However, it can also be interpreted as sth of an excuse, you don't have to do anything bc as an individual you can't do much anyway. I will interpret it as the last one; Strongly disagree means you're optimistic and maybe positive about climate change. No recode necessary. 
issp93sel$dodiff <- as.numeric(issp93sel$V27)
attributes(issp93sel$V54) #Score 1 means that ordinary people should decide how to tackle env. problems, score 2 means that government should pass laws. 
issp93sel$people_decide <- as.numeric(issp93sel$V54)
issp93sel$people_decide <- (issp93sel$people_decide - 1) * (4/1) + 1
attributes(issp93sel$V55) #Same as previous one, but then for businesses. I think for these statements, there is not one side that expresses more positive attitudes. 
issp93sel$bus_decide <- as.numeric(issp93sel$V55)
issp93sel$bus_decide <- (issp93sel$bus_decide - 1) * (4/1) + 1


#Independent vars
table(issp93sel$V204)
issp93sel$V204 <- revalue(issp93sel$V204 , c("No formal schooling, NAV"= "0", "1 year"="1", "GB:10 or less"="10", "GB:14 or more"="14", "AUS: 20 or more"="20", "43 years" = "43", "Still school J:high s" = "0", "Still college,uni" = "0", "Other answer" = "0"))
issp93sel$eduyrs <- as.numeric(as.character(issp93sel$V204))
issp93sel$isced[issp93sel$eduyrs <=4] <- 0
issp93sel$isced[issp93sel$eduyrs > 4 & issp93sel$eduyrs <= 6] <- 1
issp93sel$isced[issp93sel$eduyrs > 6 & issp93sel$eduyrs <= 10] <- 2
issp93sel$isced[issp93sel$eduyrs > 10 & issp93sel$eduyrs <= 13] <- 3
issp93sel$isced[issp93sel$eduyrs > 13 & issp93sel$eduyrs <= 15] <- 4
issp93sel$isced[issp93sel$eduyrs > 15 & issp93sel$eduyrs <= 18] <- 5
issp93sel$isced[issp93sel$eduyrs > 18] <- 6
table(issp93sel$isced, useNA = "always")

table(issp93sel$V200, useNA = "always")
issp93sel$sex <- revalue(issp93sel$V200, c("Male"="1", "Female"="2"))
table(issp93sel$lrscale, useNA = "always")
issp93sel$lrscale <- as.numeric(issp93sel$V281)
issp93sel$lrscale <- (issp93sel$lrscale - 1) * (9/4) + 1
issp93sel$lrscale[issp93sel$lrscale==12.25 | issp93sel$lrscale ==14.5] <- 11

table(issp93sel$urban, useNA = "always")
issp93sel$urban <- revalue(issp93sel$V352, c("Rural"= "Low urbanity", "Suburbs, city-town"="Medium urbanity", "Urban"="High urbanity"))
issp93sel$urban <- factor(issp93sel$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)
table(issp93sel$V277)
issp93sel$ch_attend <- as.numeric(issp93sel$V277)
issp93sel <- issp93sel %>% 
     mutate_at(c("ch_attend"), funs(recode(., `1` = 6, `2` = 5, `3` = 4, `4` = 3, `5` = 2, `6` = 1)))
table(issp93sel$ch_attend, useNA = "always")
table(issp93sel$V268, useNA = "always")

issp93sel$income <- as.numeric(issp93sel$V268) #Family income in categories (guilders). 526 missings
table(issp93sel$income, useNA = "always")
mean(issp93sel$income, na.rm=T) #Income is ordinal here, while in the other datasets it is interval
issp93sel$income[is.na(issp93sel$income)] <- 12.014
 
issp93sel$income_quart <- with(issp93sel, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
issp93sel$income_quart <- as.numeric(issp93sel$income_quart)

table(issp93sel$age, useNA = "always")
issp93sel$age <- as.numeric(as.character(issp93sel$V201))
issp93sel$V202 <- as.numeric(issp93sel$V202)
issp93sel$marstat[issp93sel$V202 < 2] <- 1 #Living together
issp93sel$marstat[issp93sel$V202 >= 2] <- 2 #Not living together


save(issp93sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp93sel.Rdata")
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp93sel.Rdata")

#Now 2000. 
issp2000 <- foreign::read.spss("/Users/anuschka/Documents/gesis_dir/issp/ISSP2000.sav", use.value.labels = T,  to.data.frame = T)

issp00sel <- issp2000 %>% filter(V3=="NL-Netherlands")  %>%
            select(V1, V2, V11, V12, V13, V14, V16, V17, V19, V20, V21, V22, V23, V40, V41, V24, V25, V26, V42, V200, V201, V202, V204, V205, V246, V241, V242, V243, V262, V274)


 #Dependent variables
issp00sel$lifeharm <- as.numeric(issp00sel$V12)
issp00sel$willing_price <- as.numeric(issp00sel$V19)
issp00sel$willing_tax <- as.numeric(issp00sel$V20)
issp00sel$willing_living <- as.numeric(issp00sel$V21)
issp00sel$do_right <- as.numeric(issp00sel$V23)
issp00sel$growharm <- as.numeric(issp00sel$V16)

issp00sel <- issp00sel %>% 
     mutate_at(c("lifeharm", "willing_price", "willing_tax", "willing_living", "do_right", "growharm"), funs(recode(., `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5`= 1)))

issp00sel$worry <- as.numeric(issp00sel$V11)
issp00sel$progharm <- as.numeric(issp00sel$V13)
issp00sel$econprotect <- as.numeric(issp00sel$V14)
issp00sel$dodiff <- as.numeric(issp00sel$V22)
issp00sel$people_decide <- as.numeric(issp00sel$V40)
issp00sel$people_decide <- (issp00sel$people_decide - 1) * (4/1) + 1
issp00sel$bus_decide <- as.numeric(issp00sel$V41)
issp00sel$bus_decide <- (issp00sel$bus_decide - 1) * (4/1) + 1

table(issp00sel$V24, useNA = "always") #More important things in life than protecting the environment. Higher score (strongly disagree) means that you care more about the environment than strongly disagree. 
issp00sel$moreimp <- as.numeric(issp00sel$V24)
table(issp00sel$V25, useNA = "always") #There's no point in doing what I can for the environment unless others do the same. Higher score means more positive attitudes, although again (like mentioned at the other statement) you can also interpret strongly agree as caring about the environment but not knowing what to do. 
issp00sel$othersame <- as.numeric(issp00sel$V25)
table(issp00sel$V26, useNA = "always") #Many of the claims exaggerated. Higher score means strongly disagree so more positive attitudes. 
issp00sel$exag <- as.numeric(issp00sel$V26)
table(issp00sel$V42, useNA = "always") #Higher score means that government is doing too little, which expresses worry and that we should do sth. 
issp00sel$country_effort <- as.numeric(issp00sel$V42)

#Independent vars 
table(issp00sel$V204, useNA = "always")
issp00sel$V204 <- revalue(issp00sel$V204 , c("No form school"= "0", "AUS: 0 or 1 yr"="0", "NZ:1-8 yrs"="8", "N:7-9 yrs"="9", "NZ:9-11 yrs"="11", "NZ:12-13 yrs" = "13", "NZ:14-16 yrs" = "15", "NZ:17+ yrs" = "17", "Other educ,other answer" = "0", "Still at college,uni" = "20", "Still at school,N:+uni" = "20"))
issp00sel$eduyrs <- as.numeric(as.character(issp00sel$V204))
issp00sel$isced[issp00sel$eduyrs <=4] <- 0
issp00sel$isced[issp00sel$eduyrs > 4 & issp00sel$eduyrs <= 6] <- 1
issp00sel$isced[issp00sel$eduyrs > 6 & issp00sel$eduyrs <= 10] <- 2
issp00sel$isced[issp00sel$eduyrs > 10 & issp00sel$eduyrs <= 13] <- 3
issp00sel$isced[issp00sel$eduyrs > 13 & issp00sel$eduyrs <= 15] <- 4
issp00sel$isced[issp00sel$eduyrs > 15 & issp00sel$eduyrs <= 18] <- 5
issp00sel$isced[issp00sel$eduyrs > 18] <- 6
table(issp00sel$isced, useNA = "always")

table(issp00sel$V200, useNA = "always")
issp00sel$sex <- revalue(issp00sel$V200, c("Male"="1", "Female"="2"))
table(issp00sel$lrscale, useNA = "always")
issp00sel$lrscale <- as.numeric(issp00sel$V246)
issp00sel$lrscale[issp00sel$lrscale==6 | issp00sel$lrscale ==7] <- 11
issp00sel$lrscale <- (issp00sel$lrscale - 1) * (9/4) + 1
issp00sel$lrscale[issp00sel$lrscale==23.5] <- 11
table(issp00sel$V274, useNA = "always")
issp00sel$urban <- revalue(issp00sel$V274, c("Rural,RP:total rur"= "Low urbanity", "Suburb,city,town,county seat"="Medium urbanity", "Urban,RP:total urb"="High urbanity"))
issp00sel$urban <- factor(issp00sel$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)
table(issp00sel$V243)
issp00sel$ch_attend <- as.numeric(issp00sel$V243)
issp00sel <- issp00sel %>% 
     mutate_at(c("ch_attend"), funs(recode(., `1` = 6, `2` = 5, `3` = 4, `4` = 3, `5` = 2, `6` = 1)))
table(issp00sel$ch_attend, useNA = "always")
table(issp00sel$V241, useNA = "always") 

issp00sel$income <- as.numeric(as.character(issp00sel$V241)) #Family income absolute (so different than previous wave). But with the same number of categories (17) as previous wave
table(issp00sel$income, useNA = "always")
mean(issp00sel$income, na.rm=T) 
issp00sel$income[is.na(issp00sel$income)] <- 65201.6
 
issp00sel$income_quart <- with(issp00sel, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
issp00sel$income_quart <- as.numeric(issp00sel$income_quart)

table(issp00sel$V201, useNA = "always")
issp00sel$V201 <- revalue(issp00sel$V201, c("16 years"= "16", "GB:18 yrs or above"="18"))
issp00sel$age <- as.numeric(as.character(issp00sel$V201))
table(issp00sel$age, useNA = "always")
table(issp00sel$marstat, useNA = "always")
issp00sel$V202 <- as.numeric(issp00sel$V202)
issp00sel$marstat[issp00sel$V202 < 2] <- 1 #Living together
issp00sel$marstat[issp00sel$V202 >= 2] <- 2 #Not living together


#Missings -----
lapply(issp00sel, table, useNA = "always")

mis_vars <- c("lrscale") 

for (var in mis_vars) {
  issp00sel[is.na(issp00sel[,var]), var] <- mean(issp00sel[,var], na.rm = TRUE)
}

save(issp00sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp00sel.Rdata")
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp00sel.Rdata")


#Last dataset of 2010. 
issp2010 <- foreign::read.spss("/Users/anuschka/Documents/gesis_dir/issp/ISSP2010.sav", use.value.labels = T,  to.data.frame = T)

issp10sel <- issp2010 %>% filter(country=="NL-Netherlands") %>%
            select(studyno, v23, v24, v25, v26, v27, v28, v29, v30, v31, v32, v33, v46, v47, v34, v35, v36, v48, SEX, AGE, EDUCYRS, PARTY_LR, NL_INC, ATTEND, URBRURAL, MARITAL)

#Dependent variables 
issp10sel$lifeharm <- as.numeric(issp10sel$v24)
issp10sel$willing_price <- as.numeric(issp10sel$v29)
issp10sel$willing_tax <- as.numeric(issp10sel$v30)
issp10sel$willing_living <- as.numeric(issp10sel$v31)
issp10sel$do_right <- as.numeric(issp10sel$v33)
issp10sel$growharm <- as.numeric(issp10sel$v27)

issp10sel <- issp10sel %>% 
     mutate_at(c("lifeharm", "willing_price", "willing_tax", "willing_living", "do_right", "growharm"), funs(recode(., `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5`= 1)))

issp10sel$worry <- as.numeric(issp10sel$v23)
issp10sel$progharm <- as.numeric(issp10sel$v25)
issp10sel$econprotect <- as.numeric(issp10sel$v26)
issp10sel$dodiff <- as.numeric(issp10sel$v32)
issp10sel$people_decide <- as.numeric(issp10sel$v46)
issp10sel$people_decide <- (issp10sel$people_decide - 1) * (4/1) + 1
issp10sel$bus_decide <- as.numeric(issp10sel$v47)
issp10sel$bus_decide <- (issp10sel$bus_decide - 1) * (4/1) + 1
issp10sel$moreimp <- as.numeric(issp10sel$v34)
issp10sel$othersame <- as.numeric(issp10sel$v35)
issp10sel$exag <- as.numeric(issp10sel$v36)
issp10sel$country_effort <- as.numeric(issp10sel$v48)


#Independent vars
table(issp10sel$EDUCYRS, useNA = "always")
issp10sel$EDUCYRS <- revalue(issp10sel$EDUCYRS , c("No formal schooling, no years at school"= "0", "8 years, NZ: primary"="8", "11 years, NZ: secondary school less than 3 years"="11", "13 years, NZ: secondary school more han 4 years"="13", "16 years, NZ: university less than 3 years"="16", "18 years, NZ: university more than 4 years" = "18", "Still at college, university, in vocational training" = "20"))
issp10sel$eduyrs <- as.numeric(as.character(issp10sel$EDUCYRS)) #170 missings on education
mean(issp10sel$eduyrs, na.rm = T)
issp10sel$eduyrs[is.na(issp10sel$eduyrs)] <- 13.694
issp10sel$isced[issp10sel$eduyrs <=4] <- 0
issp10sel$isced[issp10sel$eduyrs > 4 & issp10sel$eduyrs <= 6] <- 1
issp10sel$isced[issp10sel$eduyrs > 6 & issp10sel$eduyrs <= 10] <- 2
issp10sel$isced[issp10sel$eduyrs > 10 & issp10sel$eduyrs <= 13] <- 3
issp10sel$isced[issp10sel$eduyrs > 13 & issp10sel$eduyrs <= 15] <- 4
issp10sel$isced[issp10sel$eduyrs > 15 & issp10sel$eduyrs <= 18] <- 5
issp10sel$isced[issp10sel$eduyrs > 18] <- 6
table(issp10sel$isced, useNA = "always")
table(issp10sel$SEX, useNA = "always")
issp10sel$sex <- revalue(issp10sel$SEX, c("Male"="1", "Female"="2"))
table(issp10sel$lrscale, useNA = "always")
issp10sel$lrscale <- as.numeric(issp10sel$PARTY_LR)

issp10sel$lrscale <- (issp10sel$lrscale - 1) * (9/4) + 1
issp10sel$lrscale[issp10sel$lrscale==12.25] <- 11
table(issp10sel$URBRURAL, useNA = "always") #24 missings 
issp10sel$urban <- revalue(issp10sel$URBRURAL, c("A farm or home in the country"= "Low urbanity","A country village" = "Low urbanity",  "A town or a small city"="Medium urbanity", "The suburbs or outskirts of a big city" = "High urbanity", "A big city"="High urbanity"))
issp10sel$urban <- factor(issp10sel$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)
table(issp00sel$V243) #See how I can recode them to be the same as in previous waves
issp10sel$ch_attend <- as.numeric(issp10sel$ATTEND)
issp10sel <- issp10sel %>% 
     mutate_at(c("ch_attend"), funs(recode(., `1` = 6, `2` = 6, `3` = 5, `4` = 4, `5` = 3, `6` = 2, `7` = 2, `8`= 1)))
table(issp10sel$ch_attend, useNA = "always")
issp10sel$income <- as.numeric(issp10sel$NL_INC)
table(issp10sel$income, useNA = "always")
mean(issp10sel$income, na.rm=T) 
issp10sel$income[is.na(issp10sel$income)] <- 8.461
 
issp10sel$income_quart <- with(issp10sel, cut(income, 
                                breaks=quantile(income, probs=seq(0,1, by=0.25), na.rm=TRUE), 
                                include.lowest=TRUE))
issp10sel$income_quart <- as.numeric(issp10sel$income_quart)

table(issp10sel$AGE, useNA = "always")
issp10sel$age <- as.numeric(as.character(issp10sel$AGE))
table(issp10sel$age, useNA = "always")

issp10sel$MARITAL <- as.numeric(issp10sel$MARITAL)
issp10sel$marstat[issp10sel$MARITAL < 2] <- 1 #Living together
issp10sel$marstat[issp10sel$MARITAL >= 2] <- 2 #Not living together

save(issp10sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp10sel.Rdata")

#Missings later ----------------------------------
lapply(issp10sel, table, useNA = "always")

mis_vars <- c("lrscale") 

for (var in mis_vars) {
  issp10sel[is.na(issp10sel[,var]), var] <- mean(issp10sel[,var], na.rm = TRUE)
}

save(issp10sel, file="/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp10sel.Rdata")

#Now I have to make one datafile of these three datasets. 
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp93sel.Rdata")
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp00sel.Rdata")
load("/Users/anuschka/Documents/climatechange/climatechange/data/all_waves/issp10sel.Rdata")

issp93sel <- issp93sel %>% select(V1, worry, lifeharm, willing_price, willing_tax, willing_living, do_right, progharm, econprotect, people_decide, bus_decide, eduyrs, isced, sex, lrscale, urban, ch_attend, income, income_quart, age, marstat, weightvec)
                                  
issp00sel <- issp00sel %>% select(V1, worry, lifeharm, willing_price, willing_tax, willing_living, do_right, progharm, econprotect, people_decide, bus_decide, exag, moreimp, othersame, country_effort, growharm, eduyrs, isced, sex, lrscale, urban, ch_attend, income_a, income, income_quart, age, marstat, weightvec)

issp10sel <- issp10sel %>% select(studyno, worry, lifeharm, willing_price, willing_tax, willing_living, do_right, progharm, econprotect, people_decide, bus_decide, exag, moreimp, othersame, country_effort, growharm, eduyrs, isced, sex, lrscale, urban, ch_attend, income, income_quart, NL_INC, age, marstat, weightvec) 

issp93sel <- issp93sel %>% rename(studyno = V1)
issp00sel <- issp00sel %>% rename(studyno = V1)

#Check again whether variables have the same levels etc. (should be)
attributes(issp93sel$worry)
attributes(issp00sel$worry)
attributes(issp10sel$worry) #etc etc 

#Bind rows
issptotal <- plyr::rbind.fill(issp93sel, issp00sel, issp10sel)

table(issptotal$bus_decide, useNA = "always")

issptotal$surveyyear[issptotal$studyno==2450] <- 1993
issptotal$surveyyear[issptotal$studyno==3440] <- 2000
issptotal$surveyyear[issptotal$studyno=="GESIS Data Archive Study Number ZA5500"] <- 2010

# Change order of variables
issptotal <- issptotal %>% select(surveyyear, worry:growharm, income_a, studyno, NL_INC, weightvec)

table(issptotal$urban, useNA = "always")
issptotal$urban <- revalue(issptotal$urban, c("1"= "Low urbanity", "2"="Medium urbanity", "3" = "High urbanity"))
issptotal$urban <- factor(issptotal$urban, levels=c("Low urbanity", "Medium urbanity", "High urbanity"), ordered=TRUE)
issptotal$sex <- issptotal$sex - 1
save(issptotal, file="./data/final_data/issptotal.Rdata")

Merging all the data into large datasets

# For the three datasets, change the isced into isced cat
load("/Users/anuschka/Documents/climatechange/climatechange/data/final_data/evssel.Rdata")


table(evssel$isced)
evssel$isced_cat[evssel$isced <=2] <- "Basic"
evssel$isced_cat[evssel$isced == 3 | evssel$isced == 4] <- "Intermediate"
evssel$isced_cat[evssel$isced >=5] <- "Advanced"
evssel$isced_cat <- factor(evssel$isced_cat, levels=c("Basic", "Intermediate", "Advanced"), ordered=TRUE)
table(evssel$isced_cat)

save(evssel, file="/Users/anuschka/Documents/climatechange/climatechange/data/final_data/evssel.Rdata")


load("/Users/anuschka/Documents/climatechange/climatechange/data/final_data/esstotal.Rdata")


table(esstotal$isced)
esstotal$isced_cat[esstotal$isced <=2] <- "Basic"
esstotal$isced_cat[esstotal$isced == 3 | esstotal$isced == 4] <- "Intermediate"
esstotal$isced_cat[esstotal$isced >=5] <- "Advanced"
esstotal$isced_cat <- factor(esstotal$isced_cat, levels=c("Basic", "Intermediate", "Advanced"), ordered=TRUE)
table(esstotal$isced_cat)

save(esstotal, file="/Users/anuschka/Documents/climatechange/climatechange/data/final_data/esstotal.Rdata")

load("/Users/anuschka/Documents/climatechange/climatechange/data/final_data/issptotal.Rdata")


table(issptotal$isced)
issptotal$isced_cat[issptotal$isced <=2] <- "Basic"
issptotal$isced_cat[issptotal$isced == 3 | issptotal$isced == 4] <- "Intermediate"
issptotal$isced_cat[issptotal$isced >=5] <- "Advanced"
issptotal$isced_cat <- factor(issptotal$isced_cat, levels=c("Basic", "Intermediate", "Advanced"), ordered=TRUE)
table(issptotal$isced_cat)

save(issptotal, file="/Users/anuschka/Documents/climatechange/climatechange/data/final_data/issptotal.Rdata")