Лабораторная 6
А. С. Миндергасов
2026-03-14
Для анализа данных нам потребуются следующие пакеты:
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(psych)##
## Присоединяю пакет: 'psych'
##
## Следующие объекты скрыты от 'package:ggplot2':
##
## %+%, alphalibrary(DescTools)##
## Присоединяю пакет: 'DescTools'
##
## Следующие объекты скрыты от 'package:psych':
##
## AUC, ICC, SDlibrary(nortest)
library(ggpubr)
library(ggstatsplot)## You can cite this package as:
## Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
## Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167Задание 1. Провести комплексную проверку на нормальность переменной Sepal.Length из набора `iris.
data(iris)
str(iris$Sepal.Length)## num [1:150] 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...Проведем анализ базовых описательных статистик, также расширенный статистический анализ
summary(iris$Sepal.Length)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.300 5.100 5.800 5.843 6.400 7.900psych::describe(iris$Sepal.Length)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 150 5.84 0.83 5.8 5.81 1.04 4.3 7.9 3.6 0.31 -0.61 0.07Из полученных данных, можем отметить, что Среднее (5.843) и медиана (5.800) очень близки Разница всего ~0.043, что указывает на симметричность распределения. Mean > Median (незначительно) — очень слабая положительная асимметрия. Эксцесс (Kurtosis) составляет -0.61, Стандартная ошибка среднего (se) равна 0.07, что свидетельствует о высокой точности оценки среднего по выборке. Истинное среднее значение в генеральной совокупности
hist(iris$Sepal.Length,
col = "steelblue",
xlab = "Длина чашелистика (см)",
main = "Распределение длины чашелистика (iris)")
Видим, что длина чашелистника вариьруется от 4 до 8, наиболее часто
4,5-6,5 см, нет сильных выбросов
ggplot(iris, aes(x = Sepal.Length)) +
geom_density(fill = "lightblue", alpha = 0.7) +
geom_vline(aes(xintercept = mean(Sepal.Length)),
color = "red", linetype = "dashed", size = 1) +
labs(title = "Плотность распределения длины чашелистика",
x = "Длина чашелистика (см)",
y = "Плотность") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Результаты визуального анализа указывают на то, что также наибольшая
плотность у чашелистников средней длины, чем длиннее чашелистник, тем
меньше становится плотность.
qqnorm(iris$Sepal.Length, main = "Q-Q plot для Sepal.Length")
qqline(iris$Sepal.Length, col = "steelblue", lwd = 2)
Диаграмма рассеивания также демонстрируют, что значения распределеются
близко с прямой
boxplot(iris$Sepal.Length,
col = "steelblue",
main = "Ящичная диаграмма для Sepal.Length",
ylab = "Длина чашелистика (см)")
Далее проведем стастистический тест на нормальность
# Тест Шапиро-Уилка
shapiro.test(iris$Sepal.Length)##
## Shapiro-Wilk normality test
##
## data: iris$Sepal.Length
## W = 0.97609, p-value = 0.01018Несмотря на то, что визуально распределение выглядит симметричным (среднее 5.84 близко к медиане 5.80), статистический тест Шапиро-Уилка указывает на значимое отклонение от нормального распределения. Это означает, что длина чашелистика в наборе данных iris статистически значимо отличается от нормального распределения (при уровне значимости α = 0.05).
Рассмторим доверительные интервалы для асимметрии и эксцесса
Skew(iris$Sepal.Length, na.rm = TRUE, conf.level = 0.95)## skew lwr.ci upr.ci
## 0.30864073 0.04045809 0.53725449Kurt(iris$Sepal.Length, method = 2, conf.level = 0.95)## kurt lwr.ci upr.ci
## -0.5520640 -0.8883591 -0.1312562Оба показателя одновременно отклоняются от нуля, что объясняет, почему тест Шапиро-Уилка дал значимый результат (p = 0.010). Даже умеренные отклонения по обоим параметрам суммарно дают достаточное свидетельство против нормальности.
Задание 2. Провести одномерный анализ по переменным V12, V13, V15. Сделать двумерный анализ по региону и возрасту.
install.packages(“DataExplorer”) library(DataExplorer) introduce(df)
library(haven)
df<-read_sav("files/База_КлимРиск_2023.sav")Проведем одномерный анализ
summary(df$V12)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 2.00 2.00 10.72 3.00 99.00 14summary(df$V13)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 1.00 1.00 1.54 2.00 3.00 23summary(df$V15)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 1.000 1.000 1.516 2.000 3.000 20Рассмторим чатотные таблицы переменных
v13_table <- table(df$V13, useNA = "ifany")
print(v13_table)##
## 1 2 3 <NA>
## 575 149 166 23print(round(prop.table(v13_table) * 100, 2))##
## 1 2 3 <NA>
## 62.98 16.32 18.18 2.52v15_table <- table(df$V15, useNA = "ifany")
print(v15_table)##
## 1 2 3 <NA>
## 586 153 154 20print(round(prop.table(v15_table) * 100, 2))##
## 1 2 3 <NA>
## 64.18 16.76 16.87 2.19v12_table <- table(df$V12, useNA = "ifany")
print(v12_table)##
## 1 2 3 88 99 <NA>
## 190 437 189 19 64 14print(round(prop.table(v12_table) * 100, 2))##
## 1 2 3 88 99 <NA>
## 20.81 47.86 20.70 2.08 7.01 1.53Заметим, что у переменной v12 есть ответ 99, который следует заменить
df$V12_clean <- ifelse(df$V12 == 99, NA, df$V12)
summary(df$V12_clean)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 2.000 2.000 3.956 2.000 88.000 78Построим рафики переменныхдля более наглядного анализа.
hist(df$V12_clean, col="steelblue", main="V12", xlab="V12")
barplot(table(df$V13), col="darkgreen", main="V13")
barplot(table(df$V15), col="coral", main="V15")
Далее проведем двумерный анализ. Первоначально создадим возрастные группы
df$age_group <- cut(df$age, breaks = c(0, 30, 50, 100),
labels = c("До 30", "31-50", "50+"))Составим таблицы по регионам
table(df$V13, df$Region)##
## 1 2 3 4
## 1 106 327 142 0
## 2 84 47 18 0
## 3 76 61 29 0table(df$V15, df$Region)##
## 1 2 3 4
## 1 133 310 143 0
## 2 57 76 20 0
## 3 76 51 27 0Алтайский край отличается более равномерным распределением ответов по всем категориям. Республика Алтай и Республика Тыва демонстрируют схожую структуру с доминированием категории 1.Есть пустой регион - Монголия.
Проведем хи-квадрат тест.
chisq.test(table(df$V13, df$Region))## Warning in chisq.test(table(df$V13, df$Region)): аппроксимация на основе
## хи-квадрат может быть неправильной##
## Pearson's Chi-squared test
##
## data: table(df$V13, df$Region)
## X-squared = NaN, df = 6, p-value = NAchisq.test(table(df$V15, df$age_group))##
## Pearson's Chi-squared test
##
## data: table(df$V15, df$age_group)
## X-squared = 24.865, df = 4, p-value = 5.356e-05В таблице есть нулевая строка или столбец, поэтому тест по регионам не корректен. Тест V15 по возрасту p < 0.001 → статистически значимая связь. Возраст влияет на ответы в V15. Различия между возрастными группами не случайны
Задание 3. Проанализировать переменные с множественным выбором V14 и V16, также сделать двумерный анализ по региону. По всем видам анализа сделать таблицы и графики.
library(dplyr)
library(tidyr)
library(questionr)##
## Присоединяю пакет: 'questionr'## Следующий объект скрыт от 'package:psych':
##
## describelibrary(ggplot2)
library(forcats)для более точного анализа, очистим данныеот региона 4 - Монголия
df_clean <- df %>%
filter(!is.na(Region), Region != 4) %>%
mutate(Region = droplevels(as.factor(Region)))Создадим возрастные группы
df_clean$age_group <- cut(df_clean$age,
breaks = c(0, 30, 50, 100),
labels = c("До 30 лет", "31-50 лет", "Старше 50"))
cat("Размер очищенных данных:", nrow(df_clean), "наблюдений\n")## Размер очищенных данных: 912 наблюденийcat("Регионы:", paste(levels(df_clean$Region), collapse = ", "), "\n")## Регионы: 1, 2, 3Проведем для начала одномерный анализ переменной V14, построим таблицу
V14 <- df_clean %>%
select(V14_1, V14_2, V14_3, V14_4, V14_5, V14_6,
V14_7, V14_8, V14_9, V14_10, V14_88)
V14_clean <- as.data.frame(lapply(V14, function(x) {
if(is.numeric(x)) {
factor(x, levels = c(0, 1), labels = c("нет", "да"))
} else {
x
}
}))
V14_freq <- multi.table(V14_clean, true.codes = list("да"), freq = TRUE)
V14_results <- data.frame(
Переменная = rownames(V14_freq),
Частота = V14_freq[, 1],
Процент = round(V14_freq[, 2], 1)
) %>%
arrange(desc(Частота))
print(V14_results)## Переменная Частота Процент
## V14_2 V14_2 436 47.8
## V14_1 V14_1 428 46.9
## V14_10 V14_10 381 41.8
## V14_5 V14_5 310 34.0
## V14_4 V14_4 266 29.2
## V14_7 V14_7 211 23.1
## V14_3 V14_3 102 11.2
## V14_8 V14_8 64 7.0
## V14_6 V14_6 58 6.4
## V14_9 V14_9 42 4.6
## V14_88 V14_88 23 2.5Наиболее популярными оказались варианты V14_2 (47.8%), V14_1 (46.9%) и V14_10 (41.8%). Варианты V14_5, V14_4 и V14_7 выбрали от 23% до 34% респондентов. Остальные варианты встречаются значительно реже — менее 12%, причем вариант ‘Другое’ (V14_88) отметили лишь 2.5% опрошенных, что свидетельствует о полноте предложенного списка. Построим график.
library(ggplot2)
V14_results %>%
mutate(Переменная = factor(Переменная, levels = Переменная[order(Процент)])) %>%
ggplot(aes(x = Переменная, y = Процент)) +
geom_bar(stat = "identity", fill = "steelblue") +
geom_text(aes(label = paste0(Процент, "%")), hjust = -0.1, size = 3) +
coord_flip() +
labs(title = "V14: множественный выбор",
x = "", y = "Процент респондентов (%)") +
theme_minimal() +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
Далее проведем двумерный анализ по регионам для рассматриваемой переменной.
V14_region <- cross.multi.table(V14_clean, df_clean$Region,
true.codes = list("да"), freq = TRUE)
print(V14_region)## 1 2 3
## V14_1 63.5 43.3 31.8
## V14_2 25.5 63.0 44.3
## V14_3 17.9 9.4 5.7
## V14_4 23.4 29.6 36.5
## V14_5 22.6 39.7 37.0
## V14_6 2.2 9.0 6.2
## V14_7 17.5 26.7 22.9
## V14_8 7.7 8.1 3.6
## V14_9 1.8 6.7 3.6
## V14_10 34.7 47.3 39.1
## V14_88 0.0 4.5 1.6V14_region_df <- as.data.frame(V14_region)
V14_region_df$Переменная <- rownames(V14_region_df)
for(reg in names(V14_region_df)[1:3]) {
cat("\n---", reg, "---\n")
top3 <- V14_region_df %>%
select(Переменная, all_of(reg)) %>%
arrange(desc(.data[[reg]])) %>%
head(3)
print(top3)
}##
## --- 1 ---
## Переменная 1
## V14_1 V14_1 63.5
## V14_10 V14_10 34.7
## V14_2 V14_2 25.5
##
## --- 2 ---
## Переменная 2
## V14_2 V14_2 63.0
## V14_10 V14_10 47.3
## V14_1 V14_1 43.3
##
## --- 3 ---
## Переменная 3
## V14_2 V14_2 44.3
## V14_10 V14_10 39.1
## V14_5 V14_5 37.0также проанализируем переменную V16. Для начала проведем одномерный анализ.
V16 <- df_clean %>%
select(V16_1, V16_2, V16_3, V16_4, V16_5, V16_6,
V16_7, V16_8, V16_9, V16_88)
cat("Количество переменных в блоке V16:", ncol(V16), "\n")## Количество переменных в блоке V16: 10print(names(V16))## [1] "V16_1" "V16_2" "V16_3" "V16_4" "V16_5" "V16_6" "V16_7" "V16_8"
## [9] "V16_9" "V16_88"# Преобразуем все в факторы с правильными метками
V16_clean <- as.data.frame(lapply(V16, function(x) {
if(is.numeric(x)) {
if(all(unique(x[!is.na(x)]) %in% c(0, 1))) {
factor(x, levels = c(0, 1), labels = c("нет", "да"))
} else {
as.factor(x)
}
} else if(is.character(x)) {
as.factor(x)
} else {
x
}
}))V16_freq <- multi.table(V16_clean, true.codes = list("да"), freq = TRUE)
V16_results <- data.frame(
Переменная = rownames(V16_freq),
Частота = V16_freq[, 1],
Процент = round(V16_freq[, 2], 1)
) %>%
arrange(desc(Частота))
print(V16_results)## Переменная Частота Процент
## V16_7 V16_7 406 44.5
## V16_3 V16_3 392 43.0
## V16_1 V16_1 326 35.7
## V16_2 V16_2 290 31.8
## V16_8 V16_8 221 24.2
## V16_5 V16_5 173 19.0
## V16_9 V16_9 142 15.6
## V16_6 V16_6 80 8.8
## V16_4 V16_4 57 6.2
## V16_88 V16_88 34 3.7Построим график для данной таблицы.
p1 <- V16_results %>%
mutate(Переменная = fct_reorder(Переменная, Процент)) %>%
ggplot(aes(x = Переменная, y = Процент, fill = Процент)) +
geom_bar(stat = "identity", width = 0.7) +
geom_text(aes(label = paste0(Процент, "%")),
hjust = -0.1, size = 3.5) +
coord_flip() +
scale_fill_gradient(low = "lightcoral", high = "darkred") +
labs(title = "Распределение ответов на вопрос V16",
subtitle = "Процент респондентов, выбравших каждый вариант",
x = "", y = "Процент (%)") +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5, size = 14))
print(p1)
Анализ вопроса V16 показал, что наиболее распространенными являются
варианты V16_7 (44.5%), V16_3 (43.0%) и V16_1 (35.7%). Высокую
популярность (более 30%) также имеет вариант V16_2 (31.8%). Варианты
V16_8, V16_5 и V16_9 выбрали от 15% до 24% респондентов. Наименее
популярны варианты V16_6 (8.8%), V16_4 (6.2%) и V16_88 (3.7%). Низкий
процент выбора варианта ‘Другое’ (3.7%) свидетельствует о том, что
предложенный перечень хорошо отражает реальные практики
респондентов.
Проведем двумерный анализ по региону.
V16_region <- cross.multi.table(V16_clean, df_clean$Region,
true.codes = list("да"),
freq = TRUE)
V16_region_df <- as.data.frame(V16_region)
V16_region_df$Переменная <- rownames(V16_region_df)
print(V16_region_df)## 1 2 3 Переменная
## V16_1 23.7 45.3 30.7 V16_1
## V16_2 39.1 31.4 22.4 V16_2
## V16_3 60.2 34.1 39.1 V16_3
## V16_4 0.7 8.3 9.4 V16_4
## V16_5 6.9 27.6 16.1 V16_5
## V16_6 2.6 14.8 3.6 V16_6
## V16_7 29.2 52.2 48.4 V16_7
## V16_8 26.6 26.5 15.6 V16_8
## V16_9 20.8 15.9 7.3 V16_9
## V16_88 0.4 5.6 4.2 V16_88regions <- names(V16_region_df)[1:3]
for(reg in regions) {
cat("\n", gsub("\\.", " ", reg), ":\n")
top3 <- V16_region_df %>%
select(Переменная, all_of(reg)) %>%
arrange(desc(.data[[reg]])) %>%
head(3)
top3$Процент <- paste0(round(top3[[reg]], 1), "%")
print(top3[, c("Переменная", "Процент")])
}##
## 1 :
## Переменная Процент
## V16_3 V16_3 60.2%
## V16_2 V16_2 39.1%
## V16_7 V16_7 29.2%
##
## 2 :
## Переменная Процент
## V16_7 V16_7 52.2%
## V16_1 V16_1 45.3%
## V16_3 V16_3 34.1%
##
## 3 :
## Переменная Процент
## V16_7 V16_7 48.4%
## V16_3 V16_3 39.1%
## V16_1 V16_1 30.7%Региональный анализ V16 выявил существенные различия в структуре ответов. В Алтайском крае доминирует вариант V16_3 (60.2%), в Республике Алтай — V16_7 (52.2%) и V16_1 (45.3%), в Республике Тыва — V16_7 (48.4%) и V16_3 (39.1%). Наибольший разброс между регионами наблюдается по вариантам V16_3 (26.1%), V16_1 (21.6%) и V16_5 (20.7%), что указывает на значительную региональную специфику.
top5_vars <- V16_results %>% head(5) %>% pull(Переменная)
V16_region_long <- V16_region_df %>%
filter(Переменная %in% top5_vars) %>%
select(Переменная, `1`, `2`, `3`) %>% # обратные кавычки для цифр
pivot_longer(cols = -Переменная,
names_to = "Регион",
values_to = "Процент") %>%
mutate(Регион = case_when(
Регион == "1" ~ "Алтайский край",
Регион == "2" ~ "Республика Алтай",
Регион == "3" ~ "Республика Тыва"
))
# График
p2 <- ggplot(V16_region_long, aes(x = Переменная, y = Процент, fill = Регион)) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
geom_text(aes(label = paste0(round(Процент, 1), "%")),
position = position_dodge(width = 0.7),
vjust = -0.3, size = 3) +
scale_fill_manual(values = c("Алтайский край" = "#1f77b4",
"Республика Алтай" = "#ff7f0e",
"Республика Тыва" = "#2ca02c")) +
labs(title = "Топ-5 вариантов V16 по регионам",
x = "", y = "Процент (%)",
fill = "Регион") +
theme_minimal() +
theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 14),
legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1))
print(p2)
Данные, представленные, на диаграмме, подверждают сделанные выше
выводы.