ПРОВЕРКА ПРЕДПОСЫЛКИ МОДЕЛИ О ПАРАЛЛЕЛЬНОСТИ ТРЕНДОВ
region_means <- data %>%
filter(region %in% c(12, 137, 129), id_w >= 14 & id_w <= 25) %>%
group_by(region, id_w) %>%
summarise(mean_cons = mean(cons, na.rm = TRUE))
## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
ggplot(region_means, aes(x = id_w, y = mean_cons, color = factor(region))) +
geom_line() +
geom_point() +
labs(
title = "Среднее потребление электроэнергии для регионов",
subtitle = "Волны с 14 по 25",
x = "Волна (id_w)",
y = "Среднее потребление электроэнергии",
color = "Регион"
) +
scale_x_continuous(breaks = 14:25) +
theme_minimal() +
theme(legend.position = "bottom")
ПОСТРОЕНИЕ БАЗОВОЙ И РАСШИРЕННОЙ МОДЕЛЕЙ РЕГРЕССИЙ
DIFF-IN-DIFFS
did_data <- data %>%
filter(region %in% c(129, 137))
did_data <- did_data %>%
mutate(
treatment = ifelse(region == 137, 1, 0),
post = ifelse(id_w > 22, 1, 0), # 22 bc id_w=22 values represent pre and post periods
treatment_effect = treatment * post)
# 1. Basic DiD model without controls using lm instead of feols
basic_model <- lm(cons ~ treatment + post + treatment_effect, data = did_data)
# 2. DiD model with controls using lm
control_vars <- c("c71", "c77", "c78", "c91a","e99a", "e123", "f14")
control_formula <- paste(control_vars, collapse = " + ")
full_formula <- paste("cons ~ treatment + post + treatment_effect +", control_formula)
full_model <- lm(as.formula(full_formula), data = did_data)
summary(basic_model)
##
## Call:
## lm(formula = cons ~ treatment + post + treatment_effect, data = did_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -63.93 -26.36 -13.35 11.40 933.07
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.9281 1.1799 54.182 <2e-16 ***
## treatment -0.3315 1.7628 -0.188 0.8508
## post 5.2901 2.3803 2.222 0.0263 *
## treatment_effect -2.3774 3.5039 -0.679 0.4975
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 44.91 on 3505 degrees of freedom
## Multiple R-squared: 0.001866, Adjusted R-squared: 0.001012
## F-statistic: 2.184 on 3 and 3505 DF, p-value: 0.0878
summary(full_model)
##
## Call:
## lm(formula = as.formula(full_formula), data = did_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -129.70 -20.59 -7.53 10.82 695.67
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.781e+01 4.159e+00 11.494 < 2e-16 ***
## treatment 5.020e+00 1.587e+00 3.163 0.00157 **
## post 1.094e+01 2.050e+00 5.337 1.00e-07 ***
## treatment_effect -5.273e+00 2.926e+00 -1.802 0.07167 .
## c71 -4.097e+00 1.524e+00 -2.689 0.00721 **
## c77 -2.106e-01 1.960e+00 -0.107 0.91445
## c78 -2.047e+01 3.477e+00 -5.888 4.28e-09 ***
## c91a 6.650e+00 1.313e+00 5.064 4.33e-07 ***
## e99a NA NA NA NA
## e123 2.525e-01 6.670e-03 37.864 < 2e-16 ***
## f14 -4.026e-04 2.045e-05 -19.690 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.27 on 3498 degrees of freedom
## (1 пропущенное наблюдение удалено)
## Multiple R-squared: 0.3135, Adjusted R-squared: 0.3117
## F-statistic: 177.5 on 9 and 3498 DF, p-value: < 2.2e-16
# Визуализация результатов DiD
# Рассчитаем средние значения потребления для групп до и после
plot_data <- did_data %>%
group_by(treatment, post) %>%
summarise(mean_cons = mean(cons, na.rm = TRUE), .groups = 'drop') %>%
mutate(
group = ifelse(treatment == 1, "137 (Целевой регион)", "129 (Контрольный регион)"),
period = ifelse(post == 1, "После", "До")
)
# Создание графика DiD
ggplot(plot_data, aes(x = period, y = mean_cons, group = group, color = group)) +
geom_line(size = 1) +
geom_point(size = 3) +
theme_minimal() +
labs(
title = "Влияние введения социальной нормы на потребление электроэнергии",
x = "Период",
y = "Среднее потребление",
color = "Регион"
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Получение значения коэффициента взаимодействия и его p-значения
interaction_coef <- coef(summary(full_model))["treatment_effect", ]
cat("Коэффициент взаимодействия (DiD эффект):", interaction_coef[1], "\n")
## Коэффициент взаимодействия (DiD эффект): -5.272621
cat("Стандартная ошибка:", interaction_coef[2], "\n")
## Стандартная ошибка: 2.926371
cat("t-значение:", interaction_coef[3], "\n")
## t-значение: -1.801761
cat("p-значение:", interaction_coef[4], "\n")
## p-значение: 0.0716691
Вывод: введение соц нормы в 2014 году было связано со снижением
потребления электроэнергии на одного члена семьи в Ростовской об-ти на
5.02 киловатт
ВИЗУАЛИЗАЦИЯ
# 1. График с доверительными интервалами
did_ci_data <- did_data %>%
group_by(treatment, post) %>%
summarise(
mean_cons = mean(cons, na.rm = TRUE),
sd_cons = sd(cons, na.rm = TRUE),
n = n(),
se_cons = sd_cons / sqrt(n),
lower_ci = mean_cons - 1.96 * se_cons,
upper_ci = mean_cons + 1.96 * se_cons,
.groups = 'drop'
) %>%
mutate(
group = ifelse(treatment == 1, "137 (Целевой регион)", "129 (Контрольный регион)"),
period = ifelse(post == 1, "После", "До")
)
# График с доверительными интервалами
did_ci_plot <- ggplot(did_ci_data, aes(x = period, y = mean_cons, group = group, color = group)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
theme_minimal() +
labs(
title = "Влияние социальной нормы на потребление электроэнергии",
subtitle = "С 95% доверительными интервалами",
x = "Период",
y = "Среднее потребление (кВт·ч)",
color = "Регион"
) +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
print(did_ci_plot)
# 2. Барплот для сравнения изменений до и после
change_data <- did_data %>%
group_by(treatment, post) %>%
summarise(mean_cons = mean(cons, na.rm = TRUE), .groups = 'drop') %>%
pivot_wider(names_from = post, values_from = mean_cons) %>%
mutate(
change = `1` - `0`,
percent_change = (change / `0`) * 100,
group = ifelse(treatment == 1, "137 (Целевой регион)", "129 (Контрольный регион)")
)
bar_plot <- ggplot(change_data, aes(x = group, y = change, fill = group)) +
geom_col(width = 0.6) +
geom_text(aes(label = sprintf("%.2f", change)), vjust = -0.5) +
theme_minimal() +
labs(
title = "Абсолютное изменение в потреблении электроэнергии",
subtitle = "Сравнение целевого и контрольного регионов",
x = "",
y = "Изменение в потреблении (кВт·ч)",
fill = "Регион"
) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold")
)
percent_bar_plot <- ggplot(change_data, aes(x = group, y = percent_change, fill = group)) +
geom_col(width = 0.6) +
geom_text(aes(label = sprintf("%.1f%%", percent_change)), vjust = -0.5) +
theme_minimal() +
labs(
title = "Процентное изменение в потреблении электроэнергии",
subtitle = "Сравнение целевого и контрольного регионов",
x = "",
y = "Изменение в % от исходного уровня",
fill = "Регион"
) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold")
)
# 3. Распределение потребления по группам и периодам (плотность распределения)
density_plot <- ggplot(did_data, aes(x = cons, fill = interaction(treatment, post))) +
geom_density(alpha = 0.5) +
scale_fill_discrete(
name = "Группа и период",
labels = c("129 (Контрольный), До", "129 (Контрольный), После",
"137 (Целевой), До", "137 (Целевой), После")
) +
theme_minimal() +
labs(
title = "Распределение потребления электроэнергии",
subtitle = "По регионам и периодам до/после",
x = "Потребление (кВт·ч)",
y = "Плотность"
) +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold")
)
# 4. Boxplot для сравнения распределений
boxplot <- ggplot(did_data, aes(x = interaction(treatment, post), y = cons, fill = factor(treatment))) +
geom_boxplot() +
scale_x_discrete(labels = c("Контрольный, До", "Целевой, До", "Контрольный, После", "Целевой, После")) +
scale_fill_discrete(
name = "Регион",
labels = c("129 (Контрольный)", "137 (Целевой)")
) +
theme_minimal() +
labs(
title = "Распределение потребления электроэнергии",
subtitle = "Боксплот по регионам и периодам",
x = "",
y = "Потребление (кВт·ч)"
) +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1)
)
# 5. Визуализация эффекта DiD в виде arrows plot
did_arrow_data <- did_data %>%
group_by(treatment, post) %>%
summarise(mean_cons = mean(cons, na.rm = TRUE), .groups = 'drop') %>%
mutate(
group = ifelse(treatment == 1, "Целевой", "Контрольный"),
period = ifelse(post == 1, "После", "До"),
x_pos = ifelse(post == 0, 1, 2),
group_offset = ifelse(treatment == 1, 0.15, -0.15)
)
arrow_plot <- ggplot(did_arrow_data, aes(x = x_pos + group_offset, y = mean_cons, color = group)) +
geom_point(size = 3) +
geom_line(aes(group = group)) +
geom_text(aes(label = sprintf("%.1f", mean_cons)), vjust = -1) +
scale_x_continuous(
breaks = c(1, 2),
labels = c("До", "После"),
limits = c(0.7, 2.3)
) +
theme_minimal() +
labs(
title = "Эффект разности разностей (DiD)",
subtitle = "Изменение среднего потребления электроэнергии",
x = "",
y = "Среднее потребление (кВт·ч)",
color = "Регион"
) +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold")
)
# Вывод всех графиков
print(did_ci_plot)
print(bar_plot)
print(percent_bar_plot)
print(density_plot)
print(boxplot)
print(arrow_plot)
ggsave("did_ci_plot.png", did_ci_plot, width = 8, height = 6, dpi = 300)
ggsave("did_bar_plot.png", bar_plot, width = 8, height = 6, dpi = 300)
ggsave("did_percent_bar_plot.png", percent_bar_plot, width = 8, height = 6, dpi = 300)
ggsave("did_density_plot.png", density_plot, width = 8, height = 6, dpi = 300)
ggsave("did_boxplot.png", boxplot, width = 8, height = 6, dpi = 300)
ggsave("did_arrow_plot.png", arrow_plot, width = 8, height = 6, dpi = 300)