library(haven)
df<-read_sav("C:/Users/Admin/Documents/База_КлимРиск_2023.sav")

Проанализируйте различия в оценках уверенности в нахождении заработка в случае потери актуального места работы (вопрос В21) в зависимости от образовательного уровня респондента (вопрос А6), перекодировав его в две категории - те, кто имеют высшее образование, и кто его не имеет.

library(dplyr)
## 
## Присоединяю пакет: 'dplyr'
## Следующие объекты скрыты от 'package:stats':
## 
##     filter, lag
## Следующие объекты скрыты от 'package:base':
## 
##     intersect, setdiff, setequal, union
df<-df %>%
  rowwise() %>% 
  mutate(V18_mean=mean(c_across(starts_with("V18")), na.rm = TRUE))
hist(df$V18_mean, col="steelblue")

Перекодируем переменную по образованию в 2 группы

df<-df %>% 
  mutate(V6_obr=case_when(
    V6<=3 ~ "Без высшего образования",
    V6==4 | V6==5 ~ "С высшим образованием"
  ))

Ответы по поводу уверенности в нахождении заработка

df<-df %>% 
  mutate(V21_not=case_when(
    V21==5 ~ NA,
    V21==99 ~ NA,
    .default=as.integer(V21)
  ))

Тест Манна-Уитни

wilcox.test(V21_not~V6_obr, data = df)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  V21_not by V6_obr
## W = 43706, p-value = 0.03438
## alternative hypothesis: true location shift is not equal to 0

Рассмотрим результаты теста по критерию Мана-Уитни. Они показывают, что между выборками есть статистически значимые различия. Визуализируем

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
ggbetweenstats(df, V6_obr, V21_not, xlab = "Образование", type="nonparametric")

Количество уверенных и неуверенных примерно одинаково.

Используя созданную переменную важности соблюдения традиций (V18_3), проведите анализ по уровню образования (по переменной, созданной в упр.1) и по самооценке материального положения (вопрос А10), создав три категории - «низкие доходы», «средние доходы» и «высокие доходы, обеспеченные граждане». Используйте параметрические методы там, где это уместно, и непараметрические там, где возможно применение только таких методов. Обязательно создайте визуализации к вашим результатам.

library(car)
## Загрузка требуемого пакета: carData
## 
## Присоединяю пакет: 'car'
## Следующий объект скрыт от 'package:dplyr':
## 
##     recode
df<-df %>% 
  mutate(V10_bogat=case_when(
    V10<=2 ~ "низкие доходы",
    V10==3 ~ "средние доходы",
    V10>=4 ~ "высокие доходы"
  ))

Сравним важность соблюдения традиций и образование с использованием t-тест

leveneTest(V18_3~V6_obr, data = df)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.6404 0.4238
##       840
t.test(V18_3~V6_obr, data = df)
## 
##  Welch Two Sample t-test
## 
## data:  V18_3 by V6_obr
## t = 0.78894, df = 740.15, p-value = 0.4304
## alternative hypothesis: true difference in means between group Без высшего образования and group С высшим образованием is not equal to 0
## 95 percent confidence interval:
##  -0.07070283  0.16570929
## sample estimates:
## mean in group Без высшего образования   mean in group С высшим образованием 
##                              3.458947                              3.411444

Принимается нулевая гипотеза об отсутствии различий. Сравним важность соблюдения традиций и материального положения

kruskal.test(V18_3 ~ V10_bogat, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  V18_3 by V10_bogat
## Kruskal-Wallis chi-squared = 13.604, df = 2, p-value = 0.001112
pairwise.wilcox.test(df$V18_3, df$V10_bogat,  p.adjust.method = "BH")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df$V18_3 and df$V10_bogat 
## 
##                высокие доходы низкие доходы
## низкие доходы  0.0015         -            
## средние доходы 0.0199         0.0182       
## 
## P value adjustment method: BH

Провести анализ оценки лиц, придерживающихся определенной диеты, за три промежутка времени, используя параметрический ANOVA и его непараметрический аналог с парными сравнениями. Сделать соответствующие визуализации.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.1     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ car::recode()   masks dplyr::recode()
## ✖ purrr::some()   masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rstatix)
## 
## Присоединяю пакет: 'rstatix'
## 
## Следующий объект скрыт от 'package:stats':
## 
##     filter
library(ggpubr)
library(datarium)

Загрузка данных

data("selfesteem")
head(selfesteem)
## # A tibble: 6 × 4
##      id    t1    t2    t3
##   <int> <dbl> <dbl> <dbl>
## 1     1  4.01  5.18  7.11
## 2     2  2.56  6.91  6.31
## 3     3  3.24  4.44  9.78
## 4     4  3.42  4.71  8.35
## 5     5  2.87  3.91  6.46
## 6     6  2.05  5.34  6.65

Преобразование данных в длинный формат

selfesteem_long <- selfesteem %>%
  pivot_longer(
    cols = c(t1, t2, t3),
    names_to = "time",
    values_to = "score"
  ) %>%
  convert_as_factor(id, time)

Описательные статистики

desc_stats <- selfesteem_long %>%
  group_by(time) %>%
  get_summary_stats(score, type = "mean_sd")
desc_stats
## # A tibble: 3 × 5
##   time  variable     n  mean    sd
##   <fct> <fct>    <dbl> <dbl> <dbl>
## 1 t1    score       10  3.14 0.552
## 2 t2    score       10  4.93 0.863
## 3 t3    score       10  7.64 1.14

Визуализация данных

boxplot <- ggboxplot(selfesteem_long, x = "time", y = "score", 
                     color = "time", palette = "jco",
                     add = "jitter") +
  labs(title = "Оценка самооценки по временным точкам",
       x = "Временная точка", y = "Оценка самооценки")
boxplot

Проверка нормальности распределения

normality_test <- selfesteem_long %>%
  group_by(time) %>%
  shapiro_test(score)
normality_test
## # A tibble: 3 × 4
##   time  variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 t1    score        0.967 0.859
## 2 t2    score        0.876 0.117
## 3 t3    score        0.923 0.380

Графическая проверка нормальности

qqplot <- ggqqplot(selfesteem_long, "score", facet.by = "time") +
  labs(title = "Графики Q-Q для проверки нормальности")
qqplot

Параметрический анализ (ANOVA с повторными измерениями)

Проведение ANOVA

res.aov <- anova_test(data = selfesteem_long, dv = score, 
                      wid = id, within = time)
anova_results <- get_anova_table(res.aov)
anova_results
## ANOVA Table (type III tests)
## 
##   Effect DFn DFd      F        p p<.05   ges
## 1   time   2  18 55.469 2.01e-08     * 0.829

Парные сравнения с поправкой Бонферрони

pwc_param <- selfesteem_long %>%
  pairwise_t_test(
    score ~ time, paired = TRUE,
    p.adjust.method = "bonferroni"
  )
pwc_param
## # A tibble: 3 × 10
##   .y.   group1 group2    n1    n2 statistic    df           p p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>       <dbl> <dbl> <chr>       
## 1 score t1     t2        10    10     -4.97     9 0.000772     2e-3 **          
## 2 score t1     t3        10    10    -13.2      9 0.000000334  1e-6 ****        
## 3 score t2     t3        10    10     -4.87     9 0.000886     3e-3 **

Визуализация с результатами параметрического анализа

param_plot <- boxplot +
  stat_pvalue_manual(
    pwc_param %>% add_xy_position(x = "time"),
    step.increase = 0.1
  ) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc_param)
  )
param_plot

Непараметрический анализ (тест Фридмана)

Тест Фридмана

friedman_test <- selfesteem_long %>%
  friedman_test(score ~ time | id)
friedman_test
## # A tibble: 1 × 6
##   .y.       n statistic    df        p method       
## * <chr> <int>     <dbl> <dbl>    <dbl> <chr>        
## 1 score    10      18.2     2 0.000112 Friedman test

Парные сравнения (знаковый ранговый тест Вилкоксона с поправкой)

pwc_nonparam <- selfesteem_long %>%
  pairwise_wilcox_test(
    score ~ time, paired = TRUE,
    p.adjust.method = "bonferroni"
  )
pwc_nonparam
## # A tibble: 3 × 9
##   .y.   group1 group2    n1    n2 statistic     p p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>       
## 1 score t1     t2        10    10         0 0.002 0.006 **          
## 2 score t1     t3        10    10         0 0.002 0.006 **          
## 3 score t2     t3        10    10         1 0.004 0.012 *

Визуализация с результатами непараметрического анализа

nonparam_plot <- boxplot +
  stat_pvalue_manual(
    pwc_nonparam %>% add_xy_position(x = "time"),
    step.increase = 0.1
  ) +
  labs(
    subtitle = paste0("Тест Фридмана: χ²(", friedman_test$df, ") = ", 
                      round(friedman_test$statistic, 2), 
                      ", p = ", round(friedman_test$p, 3)),
    caption = get_pwc_label(pwc_nonparam)
  )
nonparam_plot

Альтернативная визуализация с линиями соединения

lineplot <- ggline(selfesteem_long, x = "time", y = "score", 
                   group = "id", color = "gray", alpha = 0.5) +
  stat_summary(aes(group = 1), fun = mean, geom = "line", 
               color = "red", size = 1.5) +
  stat_summary(fun.data = mean_se, geom = "errorbar", 
               width = 0.1, color = "red") +
  stat_summary(fun = mean, geom = "point", 
               color = "red", size = 3) +
  labs(title = "Изменение самооценки по временным точкам",
       x = "Временная точка", y = "Оценка самооценки",
       caption = "Красная линия показывает средние значения")
## 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.
lineplot

Вывод всех графиков

gridExtra::grid.arrange(boxplot, qqplot, param_plot, nonparam_plot, lineplot, ncol = 2)

Провести анализ таблиц сопряженности по исследованию по климату (две таблицы на ваш выбор) с помощью критерия хи-квадрат. Создать мозаичную диаграмму и ассоциативный график.

library(vcd)
## Загрузка требуемого пакета: grid
options(OutDec= ",")

создадим базовую таблицу

M<- table(df$V25, df$V1)
names(dimnames(M)) = c("Ухудшилось ли Ваше здоровье в последнее время?","Пол")
labs <- round(prop.table(M,margin=2)*100, 1)
mosaic(M, pop = FALSE, shade = TRUE)
labeling_cells(text = labs, margin = 0)(M)

struct <- structable(~ V25 + V1, data = df)
assoc(struct, data = df, shade=T, labeling_args = list(set_varnames = c(V1 = "Пол", V25="Ухудшилось ли Ваше здоровье в последнее время?"))) 

library(vcd)
options(OutDec= ",") 
M<- table(df$V10, df$V1)
names(dimnames(M)) = c("Каков уровень материального достатка Вашей семьи?","Пол")
labs <- round(prop.table(M,margin=2)*100, 1)
mosaic(M, pop = FALSE, shade = TRUE)
labeling_cells(text = labs, margin = 0)(M)

struct <- structable(~ V10 + V1, data = df)
assoc(struct, data = df, shade=T, labeling_args = list(set_varnames = c(V1 = "Пол", V10="Каков уровень материального достатка Вашей семьи?")))