Самостоятельная работа №8. Разведочный анализ в R.
Мозер К.О.
2025-03-25
Провести комплексную проверку на нормальность переменной
Sepal.Length из набора `iris.
data(iris)library(ggplot2)
ggplot(iris, aes(x = Sepal.Length)) +
geom_density(fill="lightblue")+
ggtitle("Гистограмма и плотность распределения Sepal.Length") +
theme_bw()
На графике наблюдается достаточно плоский эксцесс. Выведем описательные статистики.
library(DescTools)
Desc(list(Sepal.Length = iris$Sepal.Length), plotit = TRUE)## ──────────────────────────────────────────────────────────────────────────────
## Describe list(Sepal.Length = iris$Sepal.Length) (list):
##
## List of 1
## 1 $ Sepal.Length: num [1:150] 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##
## ──────────────────────────────────────────────────────────────────────────────
## 1 - Sepal.Length (numeric)
##
## length n NAs unique 0s mean meanCI'
## 150 150 0 35 0 5.843 5.710
## 100.0% 0.0% 0.0% 5.977
##
## .05 .10 .25 median .75 .90 .95
## 4.600 4.800 5.100 5.800 6.400 6.900 7.255
##
## range sd vcoef mad IQR skew kurt
## 3.600 0.828 0.142 1.038 1.300 0.309 -0.606
##
## lowest : 4.3, 4.4 (3), 4.5, 4.6 (4), 4.7 (2)
## highest: 7.3, 7.4, 7.6, 7.7 (4), 7.9
##
## ' 95%-CI (classic)
Skew(iris$Sepal.Length, na.rm=TRUE, conf.level = 0.95, ci.type = "classic")## skewness lwr.ci upr.ci
## 0.3086407 -0.3804196 0.3804196Kurt(iris$Sepal.Length, na.rm=TRUE, method=2, conf.level = 0.95, ci.type = "classic")## kurtosis lwr.ci upr.ci
## -0.5520640 -0.7714087 0.7714087Skew(iris$Sepal.Length, na.rm=TRUE, conf.level = 0.95, ci.type = "classic") ## skewness lwr.ci upr.ci
## 0.3086407 -0.3804196 0.3804196Kurt(iris$Sepal.Length, na.rm=TRUE, method=2, conf.level = 0.95, ci.type = "classic")## kurtosis lwr.ci upr.ci
## -0.5520640 -0.7714087 0.7714087Тест Шапиро-Уилка
shapiro.test(iris$Sepal.Length)##
## Shapiro-Wilk normality test
##
## data: iris$Sepal.Length
## W = 0.97609, p-value = 0.01018Тест Колмогорова-Смирнова
nortest::lillie.test(iris$Sepal.Length)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: iris$Sepal.Length
## D = 0.088654, p-value = 0.005788plot(ecdf(scale(iris$Sepal.Length)), col="red", main="Эмпирическая и теоретическая кумулятивные функции")
curve(pnorm, from = -10, to = 10, add = TRUE, col="blue")
Тест Андерсона-Дарлинга
nortest::ad.test(iris$Sepal.Length)##
## Anderson-Darling normality test
##
## data: iris$Sepal.Length
## A = 0.8892, p-value = 0.02251Тест Крамера-фон-Мизеса
nortest::cvm.test(iris$Sepal.Length)##
## Cramer-von Mises normality test
##
## data: iris$Sepal.Length
## W = 0.1274, p-value = 0.04706Провести одномерный анализ по переменным V12, V13, V15. Сделать двумерный анализ по региону и возрасту.
С помощью кода подготавим данные из файла SPSS для анализа, фокусируясь на трех переменных (V12, V13, V15).
library(haven)
df<-read_sav("C:/Users/Admin/Documents/База_КлимРиск_2023.sav")
library(questionr)
library(dplyr)##
## Присоединяю пакет: 'dplyr'## Следующие объекты скрыты от 'package:stats':
##
## filter, lag## Следующие объекты скрыты от 'package:base':
##
## intersect, setdiff, setequal, unionРаспределение ответов на вопрос V12
df <- haven::as_factor(df)
prop.table(table(df$V12))*100##
## Выше Ниже Не изменилась
## 21.134594 48.609566 21.023359
## Другое Затрудняюсь ответить
## 2.113459 7.119021Распределение ответов на вопрос V13
prop.table(table(df$V13))*100##
## Стала холоднее Стала теплее Не изменилась
## 64.60674 16.74157 18.65169Распределение ответов на вопрос V15
prop.table(table(df$V15))*100##
## Стала холоднее Стала теплее Не изменилась
## 65.62150 17.13326 17.24524Распределение ответов на вопрос V12 в зависимости от региона
prop.table(table(df$Region, df$V12), margin = 2)*100##
## Выше Ниже Не изменилась Другое
## Алтайский край 35.263158 20.823799 39.153439 5.263158
## Республика Алтай 47.894737 57.665904 38.624339 68.421053
## Республика Тыва 16.842105 21.510297 22.222222 26.315789
## Монголия 0.000000 0.000000 0.000000 0.000000
##
## Затрудняюсь ответить
## Алтайский край 50.000000
## Республика Алтай 23.437500
## Республика Тыва 26.562500
## Монголия 0.000000Распределение ответов на вопрос V13 в зависимости от региона
prop.table(table(df$Region, df$V13), margin = 2)*100##
## Стала холоднее Стала теплее Не изменилась
## Алтайский край 18.43478 56.37584 45.78313
## Республика Алтай 56.86957 31.54362 36.74699
## Республика Тыва 24.69565 12.08054 17.46988
## Монголия 0.00000 0.00000 0.00000Распределение ответов на вопрос V15 в зависимости от региона
prop.table(table(df$Region, df$V15), margin = 2)*100##
## Стала холоднее Стала теплее Не изменилась
## Алтайский край 22.69625 37.25490 49.35065
## Республика Алтай 52.90102 49.67320 33.11688
## Республика Тыва 24.40273 13.07190 17.53247
## Монголия 0.00000 0.00000 0.00000Распределение ответов на вопрос V12 в зависимости от возраста
prop.table(table(df$age_cats3, df$V12), margin = 2)*100##
## Выше Ниже Не изменилась Другое
## До 30 лет 18.478261 16.197183 22.826087 5.263158
## 31-49 лет 46.195652 49.061033 46.195652 63.157895
## 50 лет и старше 35.326087 34.741784 30.978261 31.578947
##
## Затрудняюсь ответить
## До 30 лет 16.666667
## 31-49 лет 48.333333
## 50 лет и старше 35.000000Распределение ответов на вопрос V13 в зависимости от возраста
prop.table(table(df$age_cats3, df$V13), margin = 2)*100##
## Стала холоднее Стала теплее Не изменилась
## До 30 лет 16.99463 15.54054 21.51899
## 31-49 лет 48.47943 45.27027 50.00000
## 50 лет и старше 34.52594 39.18919 28.48101Распределение ответов на вопрос V15 в зависимости от возраста
prop.table(table(df$age_cats3, df$V15), margin = 2)*100##
## Стала холоднее Стала теплее Не изменилась
## До 30 лет 14.65969 31.03448 17.33333
## 31-49 лет 51.65794 41.37931 41.33333
## 50 лет и старше 33.68237 27.58621 41.33333Таблица сопряженности с множественными ответами (V14 и V16) по регионам.
library(dplyr)
library(sjlabelled)##
## Присоединяю пакет: 'sjlabelled'## Следующий объект скрыт от 'package:dplyr':
##
## as_label## Следующие объекты скрыты от 'package:haven':
##
## as_factor, read_sas, read_spss, read_stata, write_sas, zap_labels## Следующий объект скрыт от 'package:ggplot2':
##
## as_labelV14<-df %>%
select(contains("V14")) %>%
select_if(is.factor)
V14tab<-V14 %>%
multi.table(true.codes=list("да"), freq=TRUE) %>%
as.data.frame() %>%
rename(Pct=`%multi`) %>%
arrange(desc(Pct))
V14tab## n Pct
## V14_2 436 47.8
## V14_1 428 46.9
## V14_10 381 41.7
## V14_5 310 34.0
## V14_4 266 29.1
## V14_7 211 23.1
## V14_3 102 11.2
## V14_8 64 7.0
## V14_6 58 6.4
## V14_9 42 4.6
## V14_88 23 2.5library(sjlabelled)
V14labs<-df %>%
select(contains("V14")) %>%
select_if(is.factor) %>%
get_label()
V14labs## V14_1
## "Резкие перепады температуры (от холода к теплу и наоборот)"
## V14_2
## "Затяжные морозы, увеличение периодов аномального холода"
## V14_3
## "Частые оттепели"
## V14_4
## "Увеличение гололеда на дорогах"
## V14_5
## "Обильные снегопады, увеличение снежного покрова"
## V14_6
## "Сход снежных лавин"
## V14_7
## "Увеличение количества пасмурных дней, нехватка солнца"
## V14_8
## "Раннее таяние и более позднее образование речного льда"
## V14_9
## "Участились ледовые зажоры, наледи на реках"
## V14_10
## "Усилились ветра, метели и снежные наносы"
## V14_88
## "Другое"V14tab$n<-V14labs
V14tab## n Pct
## V14_2 Резкие перепады температуры (от холода к теплу и наоборот) 47.8
## V14_1 Затяжные морозы, увеличение периодов аномального холода 46.9
## V14_10 Частые оттепели 41.7
## V14_5 Увеличение гололеда на дорогах 34.0
## V14_4 Обильные снегопады, увеличение снежного покрова 29.1
## V14_7 Сход снежных лавин 23.1
## V14_3 Увеличение количества пасмурных дней, нехватка солнца 11.2
## V14_8 Раннее таяние и более позднее образование речного льда 7.0
## V14_6 Участились ледовые зажоры, наледи на реках 6.4
## V14_9 Усилились ветра, метели и снежные наносы 4.6
## V14_88 Другое 2.5V14Regiontab<-cross.multi.table(V14, df$Region, true.codes=list("да"), freq=TRUE)
V14Regiontab<-as_tibble(V14Regiontab)
V14Regiontab$n<-V14labs
V14Regiontab<-V14Regiontab %>%
relocate(n)
V14Regiontab## # A tibble: 11 × 5
## n `Алтайский край` `Республика Алтай` `Республика Тыва` Монголия
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Резкие перепа… 63.5 43.3 31.8 0
## 2 Затяжные моро… 25.5 63 44.3 0
## 3 Частые оттепе… 17.9 9.4 5.7 0
## 4 Увеличение го… 23.4 29.6 36.5 0
## 5 Обильные снег… 22.6 39.7 37 0
## 6 Сход снежных … 2.2 9 6.2 0
## 7 Увеличение ко… 17.5 26.7 22.9 0
## 8 Раннее таяние… 7.7 8.1 3.6 0
## 9 Участились ле… 1.8 6.7 3.6 0
## 10 Усилились вет… 34.7 47.3 39.1 0
## 11 Другое 0 4.5 1.6 0library(dplyr)
library(gtsummary)
library(ggplot2)
library(scales)library(forcats)##
## Присоединяю пакет: 'forcats'## Следующий объект скрыт от 'package:sjlabelled':
##
## as_factorlibrary(stringr)
V14tab %>%
mutate(n = fct_reorder(n, Pct)) %>% # переставим местами уровни
ggplot(aes(x=n, y=Pct, fill=Pct)) + # создадим типовой график
geom_bar(stat="identity", width = 0.6)+ # добавим geom
geom_text(aes(label=number(Pct,accuracy=0.1, decimal.mark = ",")), hjust=-0.3, size=4)+ # поработаем с подписями данных
expand_limits(y=c(0,80))+# увеличим лимит оси
coord_flip()+ # перевернем график в горизонтальное положение
theme_void()+ # добавим минималистичную тему
theme(axis.text.y = element_text(size = 10, hjust=1))+ # установим шрифт для подписей
theme(legend.position="none")+ # уберем легенду
scale_fill_gradient(low = "red", high = "green")+ # установим цвета градиента для столбцов
scale_x_discrete(labels = function(x) str_wrap(x, width = 30))
Распределение ответов на вопрос в зависимости от региона
V14Regiontab<-cross.multi.table(V14, df$Region, true.codes=list("да"), freq=TRUE)
V14Regiontab<-as_tibble(V14Regiontab)
V14Regiontab$n<-V14labs
V14Regiontab<-V14Regiontab %>%
relocate(n)
V14Regiontab## # A tibble: 11 × 5
## n `Алтайский край` `Республика Алтай` `Республика Тыва` Монголия
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Резкие перепа… 63.5 43.3 31.8 0
## 2 Затяжные моро… 25.5 63 44.3 0
## 3 Частые оттепе… 17.9 9.4 5.7 0
## 4 Увеличение го… 23.4 29.6 36.5 0
## 5 Обильные снег… 22.6 39.7 37 0
## 6 Сход снежных … 2.2 9 6.2 0
## 7 Увеличение ко… 17.5 26.7 22.9 0
## 8 Раннее таяние… 7.7 8.1 3.6 0
## 9 Участились ле… 1.8 6.7 3.6 0
## 10 Усилились вет… 34.7 47.3 39.1 0
## 11 Другое 0 4.5 1.6 0Анализ по переменной V16.
V16<-df %>%
select(contains("V16")) %>%
select_if(is.factor)
V16tab<-V16 %>%
multi.table(true.codes=list("да"), freq=TRUE) %>%
as.data.frame() %>%
rename(Pct=`%multi`) %>%
arrange(desc(Pct))
V16tab## n Pct
## V16_7 406 44.5
## V16_3 392 42.9
## V16_1 326 35.7
## V16_2 290 31.8
## V16_8 221 24.2
## V16_5 173 18.9
## V16_9 142 15.6
## V16_6 80 8.8
## V16_4 57 6.2
## V16_88 34 3.7V16labs<-df %>%
select(contains("V16")) %>%
select_if(is.factor) %>%
get_label()
V16labs## V16_1
## "Увеличение количества засушливых дней, без осадков"
## V16_2
## "Увеличение периодов аномальной жары"
## V16_3
## "Большое количество осадков, сильных дождей"
## V16_4
## "Увеличение количества камнепадов и оползней в горах"
## V16_5
## "Усилился паводок, разливы рек, затопление мест, которые раньше не затапливались"
## V16_6
## "Таяние вечной мерзлоты, выход грунтовых вод на поверхность"
## V16_7
## "Сильные ветры, штормы"
## V16_8
## "Рост количества насекомых, комаров, мошки"
## V16_9
## "Рост количества насекомых-вредителей, угрожающих сельскохозяйственным культурам, хвойным лесам"
## V16_88
## "Другое"V16tab$n<-V16labs
V16tab## n
## V16_7 Увеличение количества засушливых дней, без осадков
## V16_3 Увеличение периодов аномальной жары
## V16_1 Большое количество осадков, сильных дождей
## V16_2 Увеличение количества камнепадов и оползней в горах
## V16_8 Усилился паводок, разливы рек, затопление мест, которые раньше не затапливались
## V16_5 Таяние вечной мерзлоты, выход грунтовых вод на поверхность
## V16_9 Сильные ветры, штормы
## V16_6 Рост количества насекомых, комаров, мошки
## V16_4 Рост количества насекомых-вредителей, угрожающих сельскохозяйственным культурам, хвойным лесам
## V16_88 Другое
## Pct
## V16_7 44.5
## V16_3 42.9
## V16_1 35.7
## V16_2 31.8
## V16_8 24.2
## V16_5 18.9
## V16_9 15.6
## V16_6 8.8
## V16_4 6.2
## V16_88 3.7V16Regiontab<-cross.multi.table(V16, df$Region, true.codes=list("да"), freq=TRUE)
V16Regiontab<-as_tibble(V16Regiontab)
V16Regiontab$n<-V16labs
V16Regiontab<-V16Regiontab %>%
relocate(n)
V16Regiontab## # A tibble: 10 × 5
## n `Алтайский край` `Республика Алтай` `Республика Тыва` Монголия
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Увеличение ко… 23.7 45.3 30.7 0
## 2 Увеличение пе… 39.1 31.4 22.4 0
## 3 Большое колич… 60.2 34.1 39.1 0
## 4 Увеличение ко… 0.7 8.3 9.4 0
## 5 Усилился паво… 6.9 27.6 16.1 0
## 6 Таяние вечной… 2.6 14.8 3.6 0
## 7 Сильные ветры… 29.2 52.2 48.4 0
## 8 Рост количест… 26.6 26.5 15.6 0
## 9 Рост количест… 20.8 15.9 7.3 0
## 10 Другое 0.4 5.6 4.2 0V16tab %>%
mutate(n = fct_reorder(n, Pct)) %>% # переставим местами уровни
ggplot(aes(x=n, y=Pct, fill=Pct)) + # создадим типовой график
geom_bar(stat="identity", width = 0.6)+ # добавим geom
geom_text(aes(label=number(Pct,accuracy=0.1, decimal.mark = ",")), hjust=-0.3, size=4)+ # поработаем с подписями данных
expand_limits(y=c(0,80))+# увеличим лимит оси
coord_flip()+ # перевернем график в горизонтальное положение
theme_void()+ # добавим минималистичную тему
theme(axis.text.y = element_text(size = 10, hjust=1))+ # установим шрифт для подписей
theme(legend.position="none")+ # уберем легенду
scale_fill_gradient(low = "red", high = "green")+ # установим цвета градиента для столбцов
scale_x_discrete(labels = function(x) str_wrap(x, width = 30))