Create a gist now

Instantly share code, notes, and snippets.

@briatte /ess.r
Last active Aug 29, 2015

What would you like to do?
library(dplyr)
library(ggplot2)
library(readr)
library(survey)
ess5 = read_csv("ESS5_FR_SDDF.csv") %>%
left_join(read_csv("ESS5FR.csv"), by = c("cntry", "idno")) %>%
ess_prepare
ess5 = svydesign(
ids = ~psu ,
strata = ~stratify ,
probs = ~prob ,
data = ess5
)
svymean(~ pbldmn, ess5, na.rm = TRUE)
svymean(~ catho, ess5, na.rm = TRUE)
m5 = svyglm(pbldmn == "Yes" ~ I(agea^2) + (gndr == "Female") +
catho + eduyrs + lrscale + region,
family = quasibinomial(link = "logit"), ess5)
# m5 = svyglm(pbldmn == "Yes" ~ region,
# family = quasibinomial(link = "logit"), ess5)
summary(m5)
m5 = ess_coefs(m5)
m5$wave = "ESS5 (2010)"
ess6 = read_csv("ESS6_FR_SDDF.csv") %>%
left_join(read_csv("ESS6FR.csv"), by = c("cntry", "idno")) %>%
ess_prepare
ess6 = svydesign(
ids = ~psu,
strata = ~stratify,
probs = ~prob,
data = ess6
)
svymean(~ pbldmn, ess6, na.rm = TRUE)
svymean(~ catho, ess6, na.rm = TRUE)
m6 = svyglm(pbldmn == "Yes" ~ I(agea^2) + (gndr == "Female") +
catho + eduyrs + lrscale + region,
family = quasibinomial(link = "logit"), ess6)
# m6 = svyglm(pbldmn == "Yes" ~ region,
# family = quasibinomial(link = "logit"), ess6)
summary(m6)
m6 = ess_coefs(m6)
m6$wave = "ESS6 (2012)"
ggplot(rbind(m5, m6) %>% filter(panel != "Contrôles"),
aes(y = b, ymin = lo, x = reorder(coef, b), ymax = hi, group = wave, color = wave)) +
geom_point(position = position_dodge(.5)) +
geom_linerange(position = position_dodge(.5)) +
geom_hline(x = 0, lty = "dashed") +
scale_color_discrete("") +
coord_flip() +
labs(y = NULL, x = NULL) +
theme_bw(14) +
theme(panel.grid = element_blank())
ggsave("coef_regions.png", height = 9, width = 7, dpi = 72)
ggplot(rbind(m5, m6) %>% filter(panel == "Contrôles"),
aes(y = b, ymin = lo, x = reorder(coef, b), ymax = hi, group = wave, color = wave)) +
geom_point(position = position_dodge(.5)) +
geom_linerange(position = position_dodge(.5)) +
geom_hline(x = 0, lty = "dashed") +
scale_color_discrete("") +
coord_flip() +
labs(y = NULL, x = NULL) +
theme_bw(14) +
theme(panel.grid = element_blank())
ggsave("coef_controles.png", height = 4, width = 7, dpi = 72)
ess_prepare <- function(x) {
x$pbldmn[ x$pbldmn %in% c("Don't know", "Refusal") ] = NA
x$rlgblg[ x$rlgblg %in% c("Don't know", "Refusal") ] = NA
x$eduyrs[ x$eduyrs %in% c("Don't know", "Refusal") ] = NA
x$eduyrs = as.numeric(x$eduyrs)
x$lrscale[ x$lrscale == "Left" ] = "0"
x$lrscale[ x$lrscale == "Right" ] = "10"
x$lrscale[ x$lrscale %in% c("Don't know", "Refusal") ] = NA
x$lrscale = as.numeric(x$lrscale)
x$catho = x$rlgdnm == "Roman Catholic"
return(x)
}
ess_coefs <- function(m) {
m = data_frame(x = names(coef(m)), b = coef(m), lo = confint(m)[, 1], hi = confint(m)[, 2])
m$x = c("(Intercept)" = "Intercept",
"I(agea^2)" = "Âge (quadratique)",
"gndr == \"Female\"TRUE" = "Femme",
"cathoTRUE" = "Religion catholique",
"eduyrs" = "Années d'éducation",
"lrscale" = "Positionnement gauche-droite",
"regionFR10" = "Île de France",
"regionFR21" = "Champagne-Ardenne",
"regionFR22" = "Picardie",
"regionFR23" = "Haute-Normandie",
"regionFR24" = "Centre",
"regionFR25" = "Basse-Normandie",
"regionFR26" = "Bourgogne",
"regionFR30" = "Nord - Pas-de-Calais",
"regionFR41" = "Lorraine",
"regionFR42" = "Alsace",
"regionFR43" = "Franche-Comté",
"regionFR51" = "Pays de la Loire",
"regionFR52" = "Bretagne",
"regionFR53" = "Poitou-Charentes",
"regionFR61" = "Aquitaine",
"regionFR62" = "Midi-Pyrénées",
"regionFR63" = "Limousin",
"regionFR71" = "Rhône-Alpes",
"regionFR72" = "Auvergne",
"regionFR81" = "Languedoc-Roussillon",
"regionFR82" = "Provence-Alpes-Côte d'Azur")[ m$x ]
m$panel = "Contrôles"
# "classification Todd"
m$panel[ m$x %in% c("Alsace", "Lorraine", "Rhône-Alpes",
# Pyrénées-Atlantiques
"Midi-Pyrénées", "Aquitaine",
# Grand Ouest
"Bretagne", "Basse-Normandie", "Pays de la Loire",
"Poitou-Charentes", "Haute-Normandie",
# Massif central méridional
"Franche-Comté") ] = "Région périphérique"
m$panel[ m$x %in% c("Picardie", "Nord - Pas-de-Calais", "Auvergne",
"Provence-Alpes-Côte d'Azur", "Champagne-Ardenne",
"Centre", "Île de France", "Limousin", "Bourgogne",
"Languedoc-Roussillon") ] = "Région centrale"
return(m)
}
library(readr)
library(dplyr)
library(survey)
library(rgdal)
library(maptools)
library(ggplot2)
library(gridExtra)
# map theme
theme_mapped = theme_bw(14) +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
strip.text = element_text(size = rel(1)),
strip.background = element_rect(fill = "grey90"),
legend.title = element_text(size = rel(1)),
legend.text = element_text(size = rel(1)),
legend.position = "bottom")
ess5 = read_csv("/Users/fr/Documents/Data/ESS/ESS_FR/ESS5FR.csv") %>%
full_join(read_csv("/Users/fr/Documents/Data/ESS/ESS_FR/ESS5_FR_SDDF.csv"), by = c("idno", "cntry"))
ess5$pbldmn[ ess5$pbldmn == "Don't know" ] = NA
ess5$pbldmn = as.integer(ess5$pbldmn == "Yes")
ess6 = read_csv("~/Documents/Data/ESS/ESS_FR/ESS6FR.csv") %>%
full_join(read_csv("~/Documents/Data/ESS/ESS_FR/ESS6_FR_SDDF.csv"), by = c("idno", "cntry"))
ess6$pbldmn[ ess6$pbldmn == "Don't know" ] = NA
ess6$pbldmn = as.integer(ess6$pbldmn == "Yes")
ess5.design = svydesign(ids = ~psu, strata = ~stratify, probs = ~prob, data = ess5, nest = TRUE)
ess6.design = svydesign(ids = ~psu, strata = ~stratify, probs = ~prob, data = ess6, nest = TRUE)
ess5 = svyby(~ pbldmn, ~ region, ess5.design, svymean, na.rm = TRUE)
ess6 = svyby(~ pbldmn, ~ region, ess6.design, svymean, na.rm = TRUE)
# EN: (c) EuroGeographics for the administrative boundaries
# http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units
eu = readOGR(dsn = "/Users/fr/Documents/Data/Maps/NUTS_2010_03M_SH/NUTS_2010_03M_SH/Data",
layer = "NUTS_RG_03M_2010")
eu = subset(eu, NUTS_ID %in% unique(c( ess5$region, ess6$region )))
map5 = fortify(eu, region = "NUTS_ID") %>%
filter(id %in% unique(ess5$region)) %>%
left_join(select(ess5, id = region, pbldmn), by = "id") %>%
mutate(pbldmn = 100 * pbldmn, panel = "ESS5 (2010)")
map5 = ggplot(map5, aes(map_id = id)) +
geom_map(aes(fill = pbldmn), map = map5, color = "white", size = 1) +
expand_limits(x = map5$long, y = map5$lat) +
scale_fill_gradient2("%", na.value = "grey80",
limits = range(100 * ess5$pbldmn),
low = "#046C9A", mid = "grey90",
midpoint = mean(100 * ess5$pbldmn),
high = "#B40F20",
guide = guide_colorbar(title.vjust = 0.75, barwidth = 15)) +
facet_grid(. ~ panel) +
coord_map("gilbert") +
theme_mapped +
labs(y = NULL, x = NULL)
map6 = fortify(eu, region = "NUTS_ID") %>%
filter(id %in% unique(ess6$region)) %>%
left_join(select(ess6, id = region, pbldmn), by = "id") %>%
mutate(pbldmn = 100 * pbldmn, panel = "ESS6 (2012)")
map6 = ggplot(map6, aes(map_id = id)) +
geom_map(aes(fill = pbldmn), map = map6, color = "white", size = 1) +
expand_limits(x = map6$long, y = map6$lat) +
scale_fill_gradient2("%", na.value = "grey80",
limits = range(100 * ess6$pbldmn),
low = "#046C9A", mid = "grey90",
midpoint = mean(100 * ess6$pbldmn),
high = "#B40F20",
guide = guide_colorbar(title.vjust = 0.75, barwidth = 15)) +
facet_grid(. ~ panel) +
coord_map("gilbert") +
theme_mapped +
labs(y = NULL, x = NULL)
png("ess_pbldmn.png", width = 1000, height = 500, units = "px", res = 72)
grid.arrange(map5, map6, ncol = 2)
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment