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

Author

Суворов Александр Юрьевич, Вергун Мария Андреевна

1 Глава 1: Установка и загрузка пакетов в R

Пакеты — это основа для расширения возможностей R. Они содержат функции, данные и документацию для решения специфических задач, от статистических тестов до визуализации данных. В этой главе мы рассмотрим, как устанавливать и загружать пакеты.

Теперь что есть что. Пакеты, библиотеки… Что это такое и в чем различия?

Если говорить кратко:

  • Пакет (Package) — это файл (или набор файлов), который вы скачиваете и устанавливаете. Это как купить книгу и поставить ее на полку.
  • Библиотека (Library) — это папка на вашем компьютере, где хранятся все установленные пакеты. Это как сама книжная полка.
  • Функция library() — это команда, которая берет пакет с полки и делает его доступным для работы. Это как взять книгу с полки и открыть ее, чтобы читать.
Термин Аналогия Что это в R? Команда
Пакет (Package) Книга в магазине Набор кода, данных и документации для решения определенной задачи. install.packages()
Библиотека (Library) Книжная полка Папка на вашем компьютере, где хранятся установленные пакеты. .libPaths()
Функция library() Взять книгу с полки Команда, которая загружает установленный пакет в память, делая его функции доступными. library(имя_пакета)

1.1 1.1 Установка пакетов из репозитория CRAN

CRAN (The Comprehensive R Archive Network) — это официальный репозиторий пакетов для R. Для установки пакетов из CRAN используется функция install.packages().

1.1.1 Базовая установка

Чтобы установить пакет (например, ggplot2 для визуализации), выполните в консоли R:

install.packages("ggplot2")

1.1.2 Выбор зеркала (репозитория)

CRAN имеет множество зеркал по всему миру. Иногда при установке R предлагает выбрать зеркало вручную. Чтобы указать зеркало напрямую в коде (например, российское зеркало МГУ), используйте аргумент repos:

# Установка пакета с зеркала МГУ
install.packages(
  "dplyr", 
  repos = "https://cran.gis-lab.info/"
)

Совет: Вы можете установить зеркало по умолчанию для всей сессии R, чтобы не указывать его каждый раз. ::: {.cell}

options(repos = "https://cran.gis-lab.info/")
install.packages("tidyr") # Теперь зеркало будет использоваться автоматически

:::

1.2 1.2 Установка пакетов из локальных файлов

Иногда вам может понадобиться установить пакет, который не находится в CRAN. Например, это может быть старая версия пакета или пакет, разработанный вашим коллегой. В этом случае вы можете установить его из файла формата .tar.gz (для Linux/macOS) или .zip (для Windows).

Для этого используется функция install.packages() с указанием пути к файлу и аргументом type = "source" (или "binary" для Windows).

# Указываем путь к файлу на вашем компьютере
# Пример для Windows:
install.packages("C:/Users/YourUser/Downloads/some-package.zip", repos = NULL)

# Пример для Linux/macOS:
install.packages("~/Downloads/some-package.tar.gz", repos = NULL, type = "source")

Важно: При установке из исходного кода (.tar.gz) на вашем компьютере должны быть установлены необходимые инструменты разработки (например, Rtools для Windows или Xcode command line tools для macOS).

1.3 1.3 Загрузка пакетов

Установка пакета — это разовое действие. После установки пакет необходимо загрузить в каждую новую сессию R, чтобы получить доступ к его функциям. Для этого используется функция library().

# Загружаем пакет ggplot2, чтобы использовать его функции
library(ggplot2)

# Теперь можно создавать графики
qplot(mpg, wt, data = mtcars)

Если вы попытаетесь использовать функцию из пакета, который не был загружен, R выдаст ошибку.

1.3.1 Разница между library() и require()

  • library(package): Загружает пакет. Если пакет не найден, R выдаст ошибку и остановит выполнение скрипта. Это предпочтительный способ в скриптах.
  • require(package): Загружает пакет и возвращает TRUE или FALSE. Если пакет не найден, R выдаст предупреждение, но продолжит выполнение скрипта. Используется внутри функций для условной загрузки.

1.4 1.4 Полезные инструменты

1.4.1 Управление пакетами с помощью pacman

Пакет pacman (package manager) значительно упрощает управление пакетами, объединяя установку и загрузку в одной функции.

# Сначала установите pacman, если у вас его нет
install.packages("pacman")

# Загрузите pacman
library(pacman)

# Установите и загрузите пакет за один раз
p_install(ggplot2)
p_load(dplyr, tidyr) # Можно загрузить несколько пакетов сразу

1.4.2 Проверка загруженных пакетов

Чтобы увидеть, какие пакеты загружены в текущей сессии, используйте:

(.packages())

1.5 1.5 Выполнение R-скриптов с помощью source()

Функция source() — это мощный инструмент для выполнения R-кода, сохраненного во внешнем файле (скрипте с расширением .R). Это позволяет разбивать ваш анализ на логические модули, переиспользовать код и делиться им с коллегами.

1.5.1 Как работает source()?

Когда вы вызываете source("путь/к/скрипту.R"), R делает следующее: 1. Читает весь файл скрипт.R. 2. Выполняет каждую строку кода последовательно, как если бы вы вводили ее прямо в консоль. 3. По умолчанию все созданные в скрипте объекты (функции, переменные, данные) загружаются в вашу текущую рабочую среду.

Это особенно полезно для загрузки кастомных функций, которые вы написали для своего проекта.

1.5.2 Пример: Загрузка скрипта с GitHub

Вы можете выполнять скрипты, которые хранятся в удаленных репозиториях, например, на GitHub. Для этого нужно указать прямой URL на “сырой” (raw) файл.

Давайте загрузим и выполним скрипт с набором полезных функций для медицинской статистики, который вы предоставили.

# Указываем URL на raw-файл на GitHub
url_script <- "https://raw.githubusercontent.com/aysuvorov/medstats/refs/heads/master/R_scripts/medstats.R"

# Выполняем скрипт с помощью source()
source(url_script)

После выполнения этой команды все функции, определенные в файле medstats.R, станут доступны в вашей сессии R. Вы можете проверить это, например, посмотрев на список функций в среде или вызвав одну из них.

# Проверим, что функция desc_stats теперь доступна
# Для этого нам понадобятся данные, например, встроенный набор данных iris
summary_all(iris$Sepal.Length)

1.5.3 Аргументы функции source()

Функция source() имеет несколько полезных аргументов для управления процессом выполнения:

  • echo = TRUE: По умолчанию source() выполняет код “молча”. Если вы хотите видеть в консоли каждую строку кода из скрипта по мере ее выполнения, установите echo = TRUE.
  • local = TRUE: По умолчанию (local = FALSE) все объекты создаются в вашей глобальной среде. Если вы хотите, чтобы объекты из скрипта были созданы во временной, изолированной среде (и не “засоряли” вашу рабочую область), используйте local = TRUE.
  • encoding = "UTF-8": Если ваш скрипт содержит символы на русском или других языках, может потребоваться явно указать кодировку.

Пример с аргументами:

# Выполнить скрипт, выводя код в консоль и работая в изолированной среде
# source(url_script, echo = TRUE, local = TRUE)

1.5.4 Важное замечание о безопасности

Использование source() для загрузки кода из интернета очень удобно, но сопряжено с риском. Выполняя код из неизвестного источника, вы потенциально запускаете на своем компьютере вредоносную программу.

Всегда проверяйте код, который вы собираетесь выполнить! Перед тем как использовать source() с URL, откройте эту ссылку в браузере и хотя бы бегло просмотрите содержимое файла. Загружайте скрипты только из доверенных источников (официальные репозитории, проверенные коллеги).

Отличное дополнение! Это одна из самых частых причин ошибок у начинающих. Внесу этот раздел в начало главы.

2 Глава 2: Загрузка и сохранение данных

Любой анализ начинается с данных. В этой главе мы рассмотрим основные способы загрузки табличных данных в R из различных форматов файлов, с которыми вы столкнетесь в биомедицинских исследованиях.

2.1 2.1 Путь к файлу: Как правильно указать, где лежат данные

Прежде чем загружать файл, R нужно точно указать, где его найти. Этот “адрес” называется путем к файлу (file path). Ошибки в пути — одна из самых частых проблем при работе с данными.

2.1.1 Относительные и абсолютные пути

  • Абсолютный путь — это полный адрес файла, начиная от корневого каталога диска. Он выглядит по-разному в зависимости от вашей операционной системы.
  • Относительный путь — это путь к файлу относительно вашей рабочей директории в R. Это более гибкий и предпочтительный способ для совместной работы.

Как узнать рабочую директорию? Используйте команду getwd().

# Показывает текущую рабочую директорию
getwd()

Как изменить рабочую директорию? Используйте setwd(), указав путь к папке с вашим проектом.

# Устанавливает новую рабочую директорию
# setwd("C:/Users/YourUser/Documents/MyProject")

Лучшая практика: Всегда начинайте свой R-проект с установки рабочей директории в папку проекта. Тогда все пути к файлам можно будет указывать относительно нее.

2.1.2 Синтаксис путей для разных операционных систем

Главная проблема — обратный слэш (\), который используется в Windows как разделитель пути. В R и многих других языках программирования \ является специальным символом (escape-символом).

2.1.2.1 Windows

Неправильно: C:\Users\YourUser\Documents\data.csv (R поймет \U, \D как спецсимволы и выдаст ошибку).

Правильные способы:

  1. Использовать двойные обратные слэши: C:\\Users\\YourUser\\Documents\\data.csv
  2. Использовать прямые слэши (рекомендуется): C:/Users/YourUser/Documents/data.csv Этот вариант работает во всех ОС и является самым надежным.

2.1.2.2 macOS и Linux

Здесь все проще. В этих системах по умолчанию используется прямой слэш (/).

Правильно: /home/youruser/documents/data.csv

2.1.3 Универсальное решение: Функция file.path()

Чтобы не думать о том, какая операционная система у вас или у вашего коллеги, используйте функцию file.path(). Она автоматически соберет путь, используя правильный разделитель для текущей системы.

# Создадим путь к файлу 'data.csv' в подпапке 'data'
# Этот код будет работать и на Windows, и на Mac, и на Linux
correct_path <- file.path("data", "clinical_trial_data.xlsx")

print(correct_path)
# Вывод на Windows будет: "data\\clinical_trial_data.xlsx"
# Вывод на Mac/Linux будет: "data/clinical_trial_data.xlsx"

Использование file.path() делает ваш код более надежным и переносимым.

2.2 2.2 Работа с файлами Excel (.xls, .xlsx, .ods)

Форматы Excel (.xls для старых версий, .xlsx для новых) и OpenDocument Spreadsheet (.ods) очень популярны. Однако R не может работать с ними “из коробки”, как с простыми текстовыми файлами. Для этого требуются специальные пакеты.

2.2.1 Способ 1: Пакет readxl (Рекомендуемый)

Пакет readxl является частью “экосистемы Tidyverse” и является золотым стандартом для чтения файлов Excel. Он прост в использовании и не требует установки Java.

Сначала установите и загрузите пакет:

# Устанавливаем пакет, если делаем это впервые
# install.packages("readxl")

# Загружаем пакет в текущую сессию
library(readxl)

Теперь вы можете читать данные:

# Укажите путь к вашему файлу. 
# Для примера создадим путь к несуществующему файлу.
path_to_xlsx <- file.path("data", "clinical_trial_data.xlsx")

# Читаем первый лист из файла
df_excel <- read_excel(path_to_xlsx)

# Читаем конкретный лист по имени
df_excel_sheet <- read_excel(path_to_xlsx, sheet = "Baseline")

# Читаем конкретный лист по номеру (нумерация с 1)
df_excel_sheet2 <- read_excel(path_to_xlsx, sheet = 2)

2.2.2 Способ 2: Пакет openxlsx

Это отличная альтернатива readxl. Его главное преимущество — он позволяет не только читать, но и записывать (создавать) файлы .xlsx без зависимостей от Java.

# install.packages("openxlsx")
library(openxlsx)

# Чтение файла
df_openxlsx <- read.xlsx(path_to_xlsx, sheet = "Baseline")

2.2.3 Способ 3: Пакет readODS для файлов .ods

Если ваши данные в формате OpenDocument (.ods), который используется в LibreOffice Calc, используйте пакет readODS.

# install.packages("readODS")
library(readODS)

path_to_ods <- file.path("data", "experiment_data.ods")

df_ods <- read_ods(path_to_ods)

2.3 2.3 Работа с текстовыми файлами (.csv, .txt)

Текстовые файлы, такие как CSV (Comma-Separated Values) и обычный текст (.txt), являются одним из самых распространенных и надежных способов обмена данными. R имеет мощные встроенные функции для их чтения.

2.3.1 Загрузка CSV-файлов

Файлы .csv — это, по сути, таблицы, где столбцы разделены запятой. Для их загрузки используется функция read.csv().

# Укажите путь к вашему CSV-файлу
path_to_csv <- file.path("data", "patient_measurements.csv")

# Читаем CSV-файл
df_csv <- read.csv(path_to_csv)

Важное замечание: Функция read.csv() по умолчанию использует запятую (,) в качестве разделителя и точку (.) в качестве десятичного знака. Это стандарт для англоязычных стран. Если ваш файл был создан в русскоязычной версии Excel, скорее всего, разделителем будет точка с запятой (;), а десятичным знаком — запятая (,). В этом случае используйте более гибкую функцию read.csv2().

# Читаем CSV-файл с региональными настройками (разделитель ;, десятичный знак ,)
df_csv_eu <- read.csv2(path_to_csv)

2.3.2 Загрузка TXT-файлов и работа с разделителями

Файлы .txt могут использовать любой символ в качестве разделителя столбцов: табуляцию, точку с запятой, пробел и т.д. Если вы не знаете, какой разделитель используется, вам придется его определить.

Как найти разделитель? Самый простой способ — открыть файл в любом текстовом редакторе (например, Блокнот, Notepad++, VS Code) и визуально определить, какой символ разделяет значения в строках.

Для загрузки таких файлов используется универсальная функция read.table(). Ключевой аргумент здесь — sep.

# Укажите путь к вашему TXT-файлу
path_to_txt <- file.path("data", "lab_results.txt")

# --- Пример 1: Разделитель - точка с запятой ---
df_txt_semicolon <- read.table(path_to_txt, sep = ";", header = TRUE)

# --- Пример 2: Разделитель - табуляция ---
# Табуляция обозначается как "\t"
df_txt_tab <- read.table(path_to_txt, sep = "\t", header = TRUE)

# --- Пример 3: Разделитель - один или несколько пробелов ---
# Для этого используется sep = "" (пустая строка)
df_txt_space <- read.table(path_to_txt, sep = "", header = TRUE)

Совет: Аргумент header = TRUE указывает R, что первая строка в файле содержит названия столбцов. Это почти всегда так в исследовательских данных. Если заголовка нет, используйте header = FALSE.

2.3.3 Современный подход: Пакет readr

Пакет readr (также часть Tidyverse) предлагает более быстрые и удобные аналоги базовых функций R (read.csv, read.table). Его главные преимущества: * Скорость: Значительно быстрее на больших файлах. * Удобство: Не преобразует строки в факторы (factors) по умолчанию, что часто более ожидаемо. * Информативность: При загрузке выводит краткую информацию о типах данных в колонках.

# install.packages("readr")
library(readr)

# Аналог read.csv(), но быстрее и умнее
df_csv_readr <- read_csv(path_to_csv)

# Аналог read.table()
# readr пытается угадать разделитель автоматически, но лучше указать его явно
df_txt_readr <- read_tsv(path_to_txt) # Для файлов с табуляцией
df_txt_delim <- read_delim(path_to_txt, delim = ";") # Для произвольного разделителя

2.3.4 Итог по текстовым файлам

  1. Для .csv с запятой: read.csv() или read_csv().
  2. Для .csv с точкой с запятой: read.csv2() или read_delim(delim = ";").
  3. Для .txt или других файлов: Определите разделитель и используйте read.table(sep = "...") или read_delim(delim = "...").
  4. Всегда используйте header = TRUE, если в вашем файле есть строка с названиями столбцов.

2.4 2.4 Сохранение таблиц в файлы

После того как вы провели обработку и анализ данных, часто возникает необходимость сохранить результат в файл для дальнейшего использования, отчетности или обмена с коллегами. Рассмотрим, как это сделать для самых популярных форматов.

2.4.1 Сохранение в файлы Excel (.xlsx)

Сохранение в формат .xlsx идеально подходит, если вы хотите отправить данные коллегам, которые не работают в R, или если вам нужно сохранить несколько таблиц на разных листах в одном файле.

Для этого мы снова воспользуемся пакетами openxlsx и writexl.

2.4.1.1 Способ 1: Пакет openxlsx

Пакет openxlsx позволяет не только сохранять один фрейм данных, но и создавать сложные Excel-файлы с несколькими листами, форматированием и даже графиками.

# install.packages("openxlsx")
library(openxlsx)

# Создадим пример данных для сохранения
final_data <- data.frame(
  patient_id = 1:5,
  group = c("Control", "Treatment", "Treatment", "Control", "Treatment"),
  blood_pressure = c(120, 115, 118, 122, 110)
)

# Укажем путь для сохранения
output_xlsx_path <- file.path("results", "final_results.xlsx")

# Сохраняем один фрейм данных на лист по умолчанию
write.xlsx(final_data, file = output_xlsx_path)

# Создадим второй набор данных
summary_stats <- data.frame(
  Statistic = c("Mean_BP", "SD_BP"),
  Value = c(mean(final_data$blood_pressure), sd(final_data$blood_pressure))
)

# Сохраняем два фрейма данных на разные листи в ОДНОМ файле
write.xlsx(
  final_data,    # Первый фрейм
  file = output_xlsx_path,
  sheetName = "Raw Data", # Имя первого листа
  append = FALSE  # Перезаписываем файл, если он существует
)

write.xlsx(
  summary_stats, # Второй фрейм
  file = output_xlsx_path,
  sheetName = "Summary", # Имя второго листа
  append = TRUE   # Добавляем лист к существующему файлу
)

2.4.1.2 Способ 2: Пакет writexl

Пакет writexl — это легкая и очень быстрая альтернатива, специально созданная для записи файлов .xlsx. Она особенно удобна для сохранения нескольких листов за один вызов.

# install.packages("writexl")
library(writexl)

# Сохраняем несколько фреймов данных в разные листы за один раз
# Имена элементов списка станут именами листов
list_of_dataframes <- list(
  "Raw Data" = final_data,
  "Summary" = summary_stats
)

write_xlsx(list_of_dataframes, path = output_xlsx_path)

2.4.2 Сохранение в файлы CSV (.csv)

Формат .csv (Comma-Separated Values) — это универсальный стандарт для обмена табличными данными. Он прост, компактен и может быть открыт практически любой программой для работы с таблицами.

2.4.2.1 Способ 1: Базовая функция write.csv()

Это стандартный способ, встроенный в R.

# Укажем путь для сохранения
output_csv_path <- file.path("results", "final_results.csv")

# Сохраняем фрейм данных в CSV
write.csv(final_data, file = output_csv_path, row.names = FALSE)

Важно: Аргумент row.names = FALSE почти всегда следует использовать. По умолчанию R записывает в файл номера строк, что в большинстве случаев является лишней и ненужной информацией.

2.4.2.2 Способ 2: Функция write_csv() из пакета readr

Как и в случае с чтением, пакет readr (часть Tidyverse) предлагает более современную и быструю альтернативу.

# install.packages("readr")
library(readr)

# Сохраняем фрейм данных в CSV
# Функция по умолчанию НЕ записывает имена строк, что удобно
write_csv(final_data, file = output_csv_path)

2.4.3 Сохранение в текстовые файлы (.txt)

Сохранение в .txt полезно, когда вы используете разделитель, отличный от запятой, например, точку с запятой или табуляцию. Для этого используется универсальная функция write.table().

# Укажем путь для сохранения
output_txt_path <- file.path("results", "final_results.txt")

# --- Пример 1: Сохранение с разделителем-точкой с запятой ---
write.table(
  final_data, 
  file = output_txt_path, 
  sep = ";",        # Указываем разделитель
  row.names = FALSE,
  quote = FALSE     # Не заключать текстовые данные в кавычки
)

# --- Пример 2: Сохранение с разделителем-табуляцией ---
# Часто используется для файлов с расширением .tsv (Tab-Separated Values)
output_tsv_path <- file.path("results", "final_results.tsv")

write.table(
  final_data,
  file = output_tsv_path,
  sep = "\t",       # Символ табуляции
  row.names = FALSE,
  quote = FALSE
)

2.5 2.5 Сохранение и загрузка R-объектов (форматы .RData и .rds)

Помимо универсальных форматов, таких как .csv или .xlsx, R имеет собственные форматы для сохранения объектов. Они позволяют сохранить R-объект (например, таблицу, модель, список - ВСЁ, ЧТО УГОДНО!) “как есть”, со всеми его свойствами, такими как типы данных столбцов, уровни факторов и т.д. Это невероятно удобно для сохранения промежуточных результатов анализа. Объекты будут в том же виде воспроизведены на другом компьютере.

2.5.1 Формат .RData (или .Rda)

Формат .RData используется для сохранения нескольких объектов R в один файл. Вы можете выбрать, какие именно объекты из вашей рабочей среды вы хотите сохранить.

Когда использовать: Идеально подходит для сохранения всего состояния вашего рабочего пространства после долгого сеанса анализа, чтобы продолжить работу позже.

2.5.1.1 Сохранение в .RData

Для сохранения используется функция save(). Вы перечисляете имена объектов (в кавычках), которые хотите сохранить.

# Создадим несколько объектов, с которыми мы работали
final_data <- data.frame(
  patient_id = 1:5,
  group = c("Control", "Treatment", "Treatment", "Control", "Treatment"),
  blood_pressure = c(120, 115, 118, 122, 110)
)

summary_stats <- data.frame(
  Statistic = c("Mean_BP", "SD_BP"),
  Value = c(mean(final_data$blood_pressure), sd(final_data$blood_pressure))
)

# Укажем путь для сохранения
output_rdata_path <- file.path("results", "analysis_workspace.RData")

# Сохраняем два объекта (final_data и summary_stats) в один файл
save(final_data, summary_stats, file = output_rdata_path)

# Также можно сохранить ВСЕ объекты из рабочей среды командой save.image()
# Это то, что R предлагает сделать при закрытии сессии
# save.image(file = "results/full_workspace.RData")

2.5.1.2 Загрузка из .RData

Для загрузки используется функция load(). Важно помнить: load() восстанавливает объекты в вашу рабочую среду с их исходными именами.

# Очистим рабочую среду для наглядности
rm(list = ls()) # Удаляем все объекты

# Проверим, что среда пуста
ls()
# character(0)

# Теперь загрузим наш файл
load(output_rdata_path)

# Проверим, какие объекты появились
ls()
# [1] "final_data"    "summary_stats"

# Объекты доступны по их именам
print(head(final_data))

2.5.2 Формат .rds

Формат .rds используется для сохранения только одного объекта R в файл.

Когда использовать: Идеально подходит для сохранения одного конкретного результата, например, большой обработанной таблицы, сложной модели или функции. Это делает ваш код более модульным.

2.5.2.1 Сохранение в .rds

Для сохранения используется функция saveRDS().

# Укажем путь для сохранения
output_rds_path <- file.path("results", "final_data_table.rds")

# Сохраняем только один объект - final_data
saveRDS(final_data, file = output_rds_path)

2.5.2.2 Загрузка из .rds

Для загрузки используется функция readRDS(). Ключевое отличие от load() в том, что readRDS() не загружает объект в рабочую среду, а возвращает его. Это значит, что вы должны присвоить результат переменной (с новым или старым именем).

# Очистим рабочую среду
rm(list = ls())

# Загружаем объект из .rds файла и присваиваем его новой переменной
loaded_data <- readRDS(output_rds_path)

# Проверим, что в среде есть только один объект с новым именем
ls()
# [1] "loaded_data"

print(head(loaded_data))

2.5.3 В чем разница и какой формат выбрать? (.RData vs .rds)

Это главный вопрос, который нужно понять.

Характеристика .RData .rds
Количество объектов Несколько (или все) объектов в одном файле. Только один объект на файл.
Имена объектов Имена “зашиты” в файле. При загрузке load() они восстанавливаются. Имя не сохраняется. Вы сами решаете, как назвать переменную при загрузке (readRDS()).
Использование Сохранение состояния всего рабочего пространства. Сохранение одного конкретного результата (таблицы, модели).
Безопасность Рискованно. load() может перезаписать существующие в вашей среде объекты с такими же именами без предупреждения. Безопасно. Вы присваиваете результат переменной, нет риска случайно перезаписать что-то важное.

Рекомендация для практики:

Для большинства задач в рамках одного проекта формат .rds является предпочтительным. Он более безопасен, гибок и способствует написанию чистого, модульного кода. Вы сохраняете результат одного шага анализа в .rds-файл, а затем в следующем скрипте загружаете его с помощью readRDS() и продолжаете работу. Это делает ваш анализ воспроизводимым и легким для отладки.

3 Глава 3: Типы данных и операции с ними

Понимание типов данных — это ключ к эффективной работе в R. Каждый объект в R, будь то число, текст или целая таблица, имеет определенный тип, который определяет, какие операции можно с ним выполнять. В этой главе мы разберем основные типы данных, научимся их проверять и преобразовывать.

3.1 3.1 Переменные, константы и операция присвоения

В R, как и в других языках программирования, мы храним данные в объектах. Давайте разберемся с основными понятиями.

3.1.1 Переменные

Переменная — это, по сути, именованное “хранилище” для данных. Вы можете присвоить переменной значение (число, текст, таблицу), а затем обращаться к ней по имени, чтобы использовать это значение в дальнейших вычислениях.

Значение, хранящееся в переменной, может изменяться в ходе выполнения программы — отсюда и название “переменная”.

3.1.2 Операция присвоения значения

Чтобы создать переменную и сохранить в ней значение, используется оператор присвоения. В R есть два основных оператора для этого: <- и =.

  • <- (стрелка): Это традиционный и наиболее рекомендуемый способ присвоения в R. Он наглядно показывает, что значение “передается” в переменную, стоящую слева.
  • = (равно): Этот оператор также работает для присвоения, но его использование может быть менее очевидным в некоторых контекстах (например, внутри функций). Для простого присвоения переменной он абсолютно корректен.

Рекомендация: Для ясности и следования общепринятым стандартам в R-сообществе, используйте <- для присвоения значений переменным.

# Создаем переменную 'patient_age' и присваиваем ей числовое значение 45
patient_age <- 45

# Создаем переменную 'patient_group' и присваиваем ей текстовое значение "Treatment"
patient_group <- "Treatment"

# Теперь мы можем использовать эти переменные
print(patient_age)
print(patient_group)

# Можно выполнять операции с переменными
patient_age_in_months <- patient_age * 12
print(patient_age_in_months)

Синтаксический сахар: В RStudio (и многих других редакторах) можно использовать сочетание клавиш Alt + - (на Windows/Linux) или Option + - (на Mac), чтобы автоматически вставить оператор <-.

3.1.3 Константы

Константа — это переменная, значение которой не предполагается изменять после ее создания. В отличие от некоторых других языков программирования, в R нет строгого механизма для создания “защищенных” констант.

Однако, по соглашению, программисты используют имена в ВЕРХНЕМ РЕГИСТРЕ для обозначения констант. Это не техническое ограничение, а скорее стилистическое правило, которое говорит другим программистам (и вам в будущем): “Этот объект не следует изменять”.

# Создаем константу для значения гравитационной постоянной
GRAVITATIONAL_CONSTANT <- 6.674e-11

# Создаем константу для максимальной допустимой дозы препарата
MAX_DOSE_MG <- 500

# Технически мы можем ее изменить, но это плохая практика
# GRAVITATIONAL_CONSTANT <- 0 # Не делайте так!

Использование констант делает код более читаемым и помогает избежать случайных изменений важных значений, таких как пороговые уровни, конверсионные коэффициенты или пути к ключевым файлам.

Конечно! Продолжаем главу. Этот раздел — ключ к пониманию того, как “думает” R.

3.2 3.2 Определение типов и структуры объектов

Когда вы работаете с данными, особенно если они были загружены из файла или созданы другим скриптом, часто возникает вопрос: “Что это за объект передо мной?”. В этом разделе мы разберем основные функции R, которые помогают исследовать объекты: определить их тип, структуру и содержимое.

3.2.1 Основные типы (классы) объектов

В R существует несколько фундаментальных типов (или классов) объектов. Вот самые важные для начинающего:

  • Вектор (vector): Самый базовый объект в R. Упорядоченная последовательность элементов одного и того же типа (например, все числа или все строки).
  • Матрица (matrix): Двумерный массив (таблица) элементов одного и того же типа.
  • Фрейм данных (data.frame): Самый важный тип для анализа данных. Таблица, где столбцы могут быть разных типов (например, один столбец с числами, другой с текстом, третий с датами). Каждый столбец во фрейме данных — это вектор.
  • Список (list): Упорядоченный набор объектов, которые могут быть любого типа и структуры. В списке может храниться число, вектор, целый фрейм данных и даже другой список.
  • Функция (function): Объект, который содержит набор инструкций для выполнения определенной задачи.

3.2.2 Функции для инспектирования объектов

Представьте, что вы — детектив, а объект — это улика. Вам нужны инструменты, чтобы ее изучить. Вот ваш основной набор инструментов.

3.2.2.1 class() — Узнать класс объекта

Функция class() отвечает на вопрос: “Что это за объект?”. Она возвращает высокоуровневый класс объекта.

# Создадим несколько объектов для примера
my_vector <- c(10, 20, 30)
my_matrix <- matrix(1:6, nrow = 2)
my_df <- data.frame(id = 1:3, name = c("A", "B", "C"))
my_list <- list(numbers = my_vector, table = my_df)
my_function <- mean

# Проверяем их классы
class(my_vector)    # "numeric" (или "integer", "character" и т.д.)
class(my_matrix)    # "matrix" "array"
class(my_df)        # "data.frame"
class(my_list)      # "list"
class(my_function)  # "function"

3.2.2.2 typeof() — Узнать низкоуровневый тип

Функция typeof() показывает, как объект хранится в памяти R. Это более “техническая” информация.

typeof(my_vector)  # "double" (для чисел с плавающей точкой)
typeof(my_matrix)  # "integer" (если числа целые) или "double"
typeof(my_df$name) # "character" (для текста)
typeof(my_list)    # "list"

3.2.2.3 str() — Структура объекта (самая полезная!)

Функция str() (от structure) — ваш лучший друг. Она показывает и класс объекта, и его внутреннюю структуру: количество элементов, их типы и первые несколько значений. Это незаменимо для быстрого понимания фреймов данных и списков.

# Посмотрим на структуру нашего фрейма данных
str(my_df)

Вывод покажет, что my_df — это data.frame с 3 наблюдениями (строками) и 2 переменными (столбцами), а также покажет тип каждого столбца и его первые значения.

# Посмотрим на структуру списка
str(my_list)

Вывод покажет, что my_list — это list из 2 элементов с именами numbers и table, и покажет структуру каждого из этих элементов.

3.2.2.4 dim(), length(), nrow(), ncol() — Размеры объекта

Эти функции помогают быстро понять размер объекта.

  • dim(): Возвращает количество строк и столбцов для двумерных объектов (матриц, фреймов данных).
  • length(): Возвращает общее количество элементов (полезно для векторов и списков).
  • nrow() и ncol(): Возвращают количество строк и столбцов отдельно (работают для фреймов данных и матриц).
# Размеры фрейма данных
dim(my_df)     # 3 2 (3 строки, 2 столбца)
nrow(my_df)    # 3
ncol(my_df)    # 2

# Длина вектора
length(my_vector) # 3

3.2.2.5 head() и tail() — Показать начало и конец объекта

Эти функции выводят первые (head) или последние (tail) несколько строк объекта. Очень удобно для быстрого просмотра содержимого таблицы, не распечатывая ее целиком.

# Создадим более длинный фрейм данных
long_df <- data.frame(
  patient_id = 1:100,
  blood_pressure = rnorm(100, mean = 120, sd = 10)
)

# Показываем первые 6 строк (по умолчанию)
head(long_df)

# Показываем последние 3 строки
tail(long_df, n = 3)

3.2.2.6 summary() — Статистическое резюме

Функция summary() предоставляет быстрый статистический обзор объекта. Для числовых столбцов она вычисляет минимум, медиану, среднее, максимум и квартили. Для текстовых — частоту встречаемости значений.

# Статистическое резюме для нашего фрейма данных
summary(long_df)

Это один из самых быстрых способов выполнить первичную разведку данных (Exploratory Data Analysis).

3.2.3 Краткий итог: ваш алгоритм инспектирования

Когда вы сталкиваетесь с новым объектом x, рекомендуемый порядок действий:

  1. class(x): Понять, с чем вы имеете дело в целом.
  2. str(x): Детально изучить его внутреннюю структуру. Это часто уже дает 90% нужной информации.
  3. head(x): Быстро взглянуть на сами данные.
  4. summary(x): Получить первое статистическое представление о данных.
  5. dim(x) или length(x): Уточнить размеры.

3.2.4 3.2.1 В чем разница между вектором и массивом?

Хотя оба термина относятся к структурам данных, в R между ними есть фундаментальное различие. Если говорить кратко: вектор — это одномерная структура, а массив — это многомерная структура.

Давайте разберем это подробно.

3.2.4.1 Вектор (Vector)

  • Размерность: Всегда одномерный. Это просто последовательность элементов.
  • Структура: У него есть только длина (length()), но нет понятий “строки” или “столбцы”.
  • Представление: Можно представить как одну длинную строку данных.
# Создаем числовой вектор
my_vector <- c(1, 2, 3, 4, 5, 6)

# Проверяем его свойства
print(my_vector)
length(my_vector) # Длина вектора - 6
class(my_vector)  # Класс - "numeric" (тип данных, а не структура)
dim(my_vector)    # Размерность - NULL, т.к. у вектора нет измерений

3.2.4.2 Массив (Array)

  • Размерность: Может быть многомерным (двумерным, трехмерным и т.д.).
  • Структура: Имеет атрибут dim (размерность), который определяет его форму. Например, c(2, 3) означает 2 строки и 3 столбца.
  • Матрица: Матрица — это частный случай массива, который имеет ровно два измерения.

Важное уточнение: array() — это функция для создания многомерных массивов. Объекты, которые она создает, имеют класс “array”.

# Создаем тот же набор данных, но в виде матрицы (2D массива) с помощью функции matrix()
my_matrix <- matrix(data = c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3)

# Проверяем его свойства
print(my_matrix)
dim(my_matrix)    # Размерность - 2 3 (2 строки, 3 столбца)
class(my_matrix)  # Класс - "matrix" "array"
length(my_matrix) # Длина - 6 (общее количество элементов)

Теперь создадим трехмерный массив с помощью функции array(). Это наглядно покажет разницу между матрицей (частным случаем) и полноценным массивом.

# Создаем трехмерный массив (2x3x2) с помощью функции array()
# Это как 2 матрицы размером 2x3, сложенные вместе
my_3d_array <- array(data = 1:12, dim = c(2, 3, 2))

# Проверяем его свойства
print(my_3d_array)
dim(my_3d_array)    # Размерность - 2 3 2
class(my_3d_array)  # Класс - "array"
length(my_3d_array) # Длина - 12 (общее количество элементов)

3.2.4.3 Ключевые отличия в одной таблице

Характеристика Вектор (vector) Массив (array) / Матрица (matrix)
Размерность 1D (одномерный) 2D и более (многомерный)
Атрибут dim Отсутствует (NULL) Обязателен (например, c(3, 2) или c(2, 3, 2))
Функция dim() Вернет NULL Вернет вектор размерностей
Аналогия Список или одна строка данных Таблица с рядами и колонками, куб данных
Функция для создания c(), : matrix(), array()

3.2.4.4 Как R это “видит”?

Функция class() четко показывает эту разницу:

class(my_vector)  # "numeric" (тип данных элементов)
class(my_matrix)  # "matrix" "array" (структура данных)
class(my_3d_array) # "array" (структура данных)

Для вектора R сообщает о типе хранящихся в нем данных, а для массива/матрицы — о его структуре.

Всегда помните, что вектор — это фундаментальный одномерный “строительный блок” в R, а массив (включая матрицу) — это более сложная, многомерная структура, созданная с помощью функций matrix() или array() на основе вектора.

Отличная тема! Именование и атрибуты — это то, что делает код в R читаемым и удобным. Добавляем этот подраздел.

3.2.5 3.2.2 Атрибуты и именование элементов

Помимо самих данных, объекты в R могут иметь дополнительную мета-информацию. Эта мета-информация называется атрибутами. Самый распространенный и важный атрибут — это имена элементов.

Атрибуты — это пары “имя-значение”, которые прикреплены к объекту. Они не влияют на сами данные, но предоставляют дополнительную информацию о них. Вы можете посмотреть все атрибуты объекта с помощью функции attributes().

Давайте рассмотрим, как именование работает для разных структур данных.

3.2.5.1 Именование в векторах

Вы можете присвоить имена каждому элементу вектора. Это делает код гораздо понятнее, так как вы можете обращаться к элементу по его имени, а не по числовому индексу.

# Создадим вектор с показателями крови
blood_pressure <- c(120, 80)

# Присвоим имена с помощью функции names()
names(blood_pressure) <- c("systolic", "diastolic")

# Теперь вектор стал гораздо информативнее
print(blood_pressure)

# Мы можем получить имена всех элементов
names(blood_pressure)

# И обратиться к элементу по его имени
blood_pressure["systolic"]

3.2.5.2 Именование в матрицах и массивах

Для двумерных массивов (матриц) мы можем задавать имена для строк и столбцов с помощью функций rownames() и colnames().

# Создадим матрицу с результатами анализов
lab_results <- matrix(
  c(4.5, 140, 98, 5.1, 155, 101),
  nrow = 2,
  byrow = TRUE
)

# Зададим имена столбцов
colnames(lab_results) <- c("Hemoglobin", "Glucose", "HeartRate")

# Зададим имена строк (например, ID пациентов)
rownames(lab_results) <- c("Patient_01", "Patient_02")

# Посмотрим на результат
print(lab_results)

# Теперь можно обращаться к данным по именам
lab_results["Patient_01", "Glucose"]

3.2.5.3 Именование в списках

Списки — это структура, где именование особенно полезно. Имена элементов списка позволяют создавать легко читаемые “контейнеры” для разных объектов.

# Создадим список, содержащий разную информацию о пациенте
patient_info <- list(
  id = "P-001",
  age = 45,
  diagnosis = "Hypertension",
  measurements = blood_pressure # Используем именованный вектор из предыдущего примера
)

# Имена элементов списка
names(patient_info)

# Обращение к элементам списка по имени
patient_info$age
patient_info$diagnosis

3.2.5.4 Именование во фреймах данных (таблицах)

Фрейм данных — это, по сути, список векторов-столбцов. Поэтому именование столбцов работает так же, как и в списках. Имена столбцов — это самый важный вид именования в анализе данных.

# Создадим фрейм данных
clinical_trial <- data.frame(
  patient_id = c("P-001", "P-002", "P-003"),
  treatment_group = c("Drug", "Placebo", "Drug"),
  baseline_bp = c(145, 138, 152)
)

# Имена столбцов
names(clinical_trial)
# или
colnames(clinical_trial)

# Обращение к столбцу по имени
clinical_trial$treatment_group

3.2.6 Оператор $ для доступа к элементам

Оператор доллара ($) — это один из самых часто используемых операторов в R. Его основная задача — быстрый и удобный доступ к элементам списка или столбцам фрейма данных по их имени.

Как он работает: имя_объекта$имя_элемента

# Получаем столбец 'treatment_group' из фрейма данных
clinical_trial$treatment_group

# Получаем элемент 'measurements' из списка
patient_info$measurements

Ключевые особенности оператора $:

  1. Только для имен: Вы можете использовать его только с именами элементов, а не с числовыми индексами.
  2. Упрощение доступа: Это самый короткий и читаемый способ получить доступ к столбцу таблицы.
  3. Частичное совпадение: Оператор $ допускает частичное совпадение имен. Если имя элемента уникально, вы можете написать только его начало.
# R поймет, что мы имеем в виду 'treatment_group'
clinical_trial$treat

Предупреждение: Хотя частичное совпадение удобно, оно может быть источником ошибок, если у вас есть похожие имена (например, treatment_group и treatment_date). Для надежности всегда используйте полные имена.

3.2.7 Итог

  • Атрибуты — это мета-информация об объекте. Самый важный атрибут — names.
  • Именование делает код самодокументируемым и легким для понимания.
  • Используйте names(), rownames(), colnames() для установки и получения имен.
  • Оператор $ — ваш основной инструмент для доступа к именованным элементам списков и столбцам фреймов данных.

3.2.8 3.2.3 Типы данных внутри объектов и их преобразование

Мы разобрались со структурой объектов (вектор, список, фрейм данных). Теперь посмотрим “внутрь” этих объектов на типы данных, которые они содержат. Правильное определение и преобразование типов данных — основа чистого анализа и отсутствия ошибок.

3.2.8.1 Основные типы данных (атомарные векторы)

В R существует несколько базовых типов данных. Большинство столбцов в ваших таблицах будут относиться к одному из них:

Тип данных Описание Пример Как проверить (is.*())
character Текстовые строки. "Treatment", "Patient A" is.character()
numeric (или double) Числа с плавающей точкой (дробные). 12.5, 0.001 is.numeric()
integer Целые числа. 10L, -5L (суффикс L указывает на целое число) is.integer()
logical Логические значения. TRUE, FALSE is.logical()
factor Категориальные данные (факторы). "Control", "Low", "Medium" is.factor()

Факторы (factor) — это специальный тип данных в R для хранения категориальных переменных. Внутри они хранятся как числа (уровни), но отображаются как строки. Это очень эффективно для переменных с ограниченным набором повторяющихся значений (например, группа лечения, пол, стадия заболевания).

3.2.8.2 Проверка типа данных

Чтобы узнать тип данных вектора или столбца во фрейме данных, используйте функции class() или typeof(). Для проверки на конкретный тип удобнее использовать функции is.*().

# Создадим фрейм данных с разными типами
sample_data <- data.frame(
  patient_id = 1:3,
  group = c("Control", "Drug", "Drug"),
  recovered = c(TRUE, FALSE, TRUE),
  stringsAsFactors = FALSE # Важно: не преобразовывать строки в факторы автоматически
)

# Проверяем типы столбцов
class(sample_data$patient_id)  # "integer"
class(sample_data$group)       # "character"
class(sample_data$recovered)   # "logical"

# Используем функции is.*() для проверки
is.numeric(sample_data$patient_id)  # TRUE
is.character(sample_data$group)      # TRUE
is.logical(sample_data$recovered)    # TRUE

3.2.8.3 Преобразование типов (Coercion)

Очень часто данные загружаются не в том формате. Например, числовой столбец может быть прочитан как текстовый (character), если в нем есть опечатка (например, запятая вместо точки). В таких случаях тип нужно явно преобразовать.

Для этого используются функции-преобразователи вида as.*():

  • as.character() -> в текстовый
  • as.numeric() -> в числовой
  • as.integer() -> в целый
  • as.logical() -> в логический
  • as.factor() -> в фактор

Важно: Если преобразование невозможно, R вернет значение NA (Not Available) и выдаст предупреждение.

# --- Пример 1: Текст в число ---
# Допустим, столбец возраста был загружен как текст
age_text <- c("25", "30", "45")

# Преобразуем его в числа
age_numeric <- as.numeric(age_text)
print(age_numeric)

# --- Пример 2: Число в фактор (категорию) ---
# Допустим, у нас есть столбец с кодами групп
group_codes <- c(1, 2, 2, 1)

# Преобразуем в осмысленные факторы
group_factor <- as.factor(group_codes)
print(group_factor)

# Можно сразу задать имена уровням
group_factor_named <- factor(group_codes, levels = c(1, 2), labels = c("Placebo", "Treatment"))
print(group_factor_named)

# --- Пример 3: Проблемное преобразование ---
# В данных есть ошибка ("N/A" вместо числа)
bad_numbers <- c("10", "20", "N/A", "40")

# Попытка преобразования
as.numeric(bad_numbers)
# R выдаст предупреждение: "NAs introduced by coercion"
# Результат будет: 10 20 NA 40

3.2.8.4 Неявное преобразование (Coercion)

Иногда R пытается преобразовать типы автоматически, без явного запроса. Это называется неявным преобразованием. Чаще всего это происходит при объединении векторов разных типов в один.

R использует следующую иерархию: logical -> integer -> numeric -> character.

Это означает, что если вы объедините вектор чисел с вектором текста, все элементы станут текстом.

# Объединяем число и текст
mixed_vector <- c(10, 20, "some text")

# Какого типа получился вектор?
class(mixed_vector) # "character"
print(mixed_vector)  # "10" "20" "some text"

Неявное преобразование — распространенный источник ошибок. Например, если в вашем числовом столбце случайно оказалось одно текстовое значение, весь столбец может превратиться в текстовый, и математические операции перестанут работать. Поэтому всегда проверяйте типы данных после загрузки с помощью str() или summary().

3.2.9 Итог

  1. Проверяйте типы данных ваших столбцов с помощью str() или class().
  2. Используйте is.*() для проверки на конкретный тип.
  3. Используйте as.*() для явного преобразования типов. Будьте готовы к появлению NA, если преобразование не удалось.
  4. Остерегайтесь неявного преобразования при объединении данных, особенно при работе с текстом.

3.2.10 3.2.4 Работа с факторами (категориальными переменными)

Факторы (factor) — это специальный тип данных в R, предназначенный для хранения категориальной информации. Внешне они выглядят как строки (например, "Мужской", "Женский"), но внутри хранятся как целые числа (уровни), что делает работу с ними очень эффективной. Факторы критически важны для статистического моделирования, так как R будет правильно обрабатывать их как категориальные предикторы.

Давайте рассмотрим ключевые аспекты работы с факторами, включая преобразование бинарных данных.

3.2.10.1 1. Референсная категория (Reference Level)

При построении статистических моделей (например, линейной регрессии) R выбирает одну категорию фактора в качестве референсной (или базовой). Все остальные категории сравниваются с ней. По умолчанию в качестве референсной выбирается первая категория по алфавиту.

Это не всегда то, что нам нужно. Например, при сравнении лечения с плацебо логично было бы выбрать “Плацебо” в качестве референсной категории, а не “Лечение” (что может случиться по алфавиту).

Как посмотреть и изменить референсную категорию?

Для этого используется функция relevel().

# Создадим фактор с группами лечения
treatment <- factor(c("Drug", "Placebo", "Drug", "Placebo"))

# Посмотрим на уровни фактора
levels(treatment) # "Drug" "Placebo" (по алфавиту)

# Текущая референсная категория - "Drug"
# Давайте сделаем "Placebo" референсной
treatment_releveled <- relevel(treatment, ref = "Placebo")

# Проверим уровни снова
levels(treatment_releveled) # "Placebo" "Drug"

Теперь в любой модели, использующей treatment_releveled, за основу будет браться группа “Placebo”.

3.2.10.2 2. Порядковые факторы (Ordinal Factors)

Обычные факторы являются номинальными — их категории не имеют внутреннего порядка (например, “Мужской” и “Женский”). Но часто мы имеем дело с порядковыми данными, где важна последовательность категорий (например, “Легкая”, “Средняя”, “Тяжелая” степень заболевания).

Для создания порядкового фактора используется тот же конструктор factor(), но с дополнительным аргументом ordered = TRUE.

# Данные о стадии заболевания
disease_stage <- c("Moderate", "Mild", "Severe", "Mild")

# Создаем порядковый фактор, явно указывая правильный порядок уровней
ordinal_stage <- factor(
  disease_stage,
  levels = c("Mild", "Moderate", "Severe"),
  ordered = TRUE
)

# Посмотрим на результат
print(ordinal_stage)

# Обратите внимание на вывод: уровни теперь упорядочены
# Levels: Mild < Moderate < Severe

Указание ordered = TRUE важно не только для корректного отображения, но и для некоторых статистических тестов, которые учитывают порядок категорий.

3.2.10.3 3. Преобразование бинарных переменных (0 и 1) в факторы и обратно

Очень часто в наборах данных бинарные переменные (например, “да/нет”, “болен/здоров”, “мужчина/женщина”) представлены числами 0 и 1. Для анализа и визуализации их часто полезно преобразовать в осмысленный фактор.

Преобразование числового вектора (0, 1) в фактор

Используйте функцию factor(), указав labels для понятных названий категорий.

# Числовой вектор, где 1 = "Умер", 0 = "Выжил"
mortality_numeric <- c(0, 1, 0, 0, 1)

# Преобразуем в фактор с осмысленными метками
mortality_factor <- factor(
  mortality_numeric,
  levels = c(0, 1),
  labels = c("Выжил", "Умер")
)

print(mortality_factor)
# Levels: Выжил Умер

Обратное преобразование: из фактора в числовой вектор (0, 1)

Это более тонкий момент. Если вы попробуете применить as.numeric() напрямую к фактору, вы получите не исходные 0 и 1, а внутренние числовые коды уровней фактора (обычно 1 и 2).

# НЕПРАВИЛЬНЫЙ СПОСОБ
mortality_factor <- factor(
  mortality_numeric,
  levels = c(0, 1)
)

as.numeric(mortality_factor)

# [1] 1 2 1 1 2  (Это внутренние коды: "Выжил"=1, "Умер"=2)

Чтобы получить исходные числовые значения (0 и 1), нужно сначала преобразовать фактор в текстовый вектор (as.character()), а затем — в числовой (as.numeric()).

# ПРАВИЛЬНЫЙ СПОСОБ: через промежуточное преобразование в текст
mortality_numeric_back <- as.numeric(as.character(mortality_factor))

print(mortality_numeric_back)
# [1] 0 1 0 0 1  (Это именно те значения, которые мы хотели)

Почему это работает? as.character() возвращает вектор текстовых меток ("Выжил", "Умер"). Затем as.numeric() пытается преобразовать этот текстовый вектор в числа. Но "Выжил" и "Умер" не являются числами, и R вернет NA. Чтобы получить 0 и 1, мы должны сначала преобразовать фактор в его уровни (as.numeric(as.factor(x))), а затем, если нужно, сопоставить их с нужными значениями. Однако самый надежный способ — это хранить исходное соответствие или использовать логическую индексацию.

Самый надежный и явный способ получить 0 и 1 из фактора — это сравнение:

# САМЫЙ НАДЕЖНЫЙ СПОСОБ
# Сравниваем фактор с нужным уровнем. Результатом будет логический вектор TRUE/FALSE.
mortality_factor <- factor(
  mortality_numeric,
  levels = c(0, 1),
  labels = c("Выжил", "Умер")
)

# as.numeric() преобразует TRUE в 1, а FALSE в 0.
mortality_binary <- as.numeric(mortality_factor == "Умер")

print(mortality_binary)
# [1] 0 1 0 0 1

3.2.10.4 4. Переопределение категорий (изменение уровней)

Часто возникает необходимость изменить названия категорий, объединить некоторые из них или исправить опечатки.

Способ 1: Изменение имен всех уровней сразу

Если вы хотите переименовать все категории, вы можете просто присвоить новый вектор имен функции levels().

# Вернемся к нашему фактору лечения
levels(treatment) # "Drug" "Placebo"

# Переименуем категории на русский язык
levels(treatment) <- c("Лечение", "Плацебо")

print(treatment)
# Levels: Лечение Плацебо

Способ 2: Переименование или объединение конкретных категорий

Более гибкий способ — создать новый вектор с нужными значениями, используя старые. Это позволяет не только переименовывать, но и объединять категории.

Представим, что в данных есть опечатка и несколько категорий нужно объединить.

# Создадим фактор с "грязными" данными
group <- factor(c("Control", "Treatment", "Treatmnt", "Control", "Placebo"))

# Посмотрим на уникальные уровни
levels(group) # "Control" "Placebo" "Treatmnt" "Treatment"

# Исправим опечатку "Treatmnt" и объединим "Placebo" и "Control" в одну группу "No Treatment"
# Для этого создадим новый вектор, используя логическую индексацию
cleaned_group <- group
levels(cleaned_group) <- c("No Treatment", "No Treatment", "Treatment", "Treatment")

# Посмотрим на результат
print(cleaned_group)
# Levels: No Treatment Treatment

Объяснение: Мы создали новый вектор имен c("No Treatment", "No Treatment", "Treatment", "Treatment"). R сопоставил его со старыми уровнями в алфавитном порядке: - Control стал "No Treatment" - Placebo стал "No Treatment" - Treatmnt стал "Treatment" - Treatment стал "Treatment"

Таким образом, мы исправили опечатку и объединили две категории в одну.

4 Глава 4: Функции

Функции — это фундаментальные строительные блоки в R. По сути, вы уже использовали множество функций: c(), mean(), read.csv(), str() и т.д. Они позволяют упаковывать последовательность команд в единый блок, который можно вызывать по имени, передавая ему разные входные данные. Это делает код модульным, читаемым и переиспользуемым.

В этой главе мы научимся создавать собственные функции.

4.1 4.1 Устройство и создание функций

Функция в R — это объект, который принимает на вход некоторые аргументы (входные данные), выполняет определенный набор инструкций и возвращает результат.

4.1.1 Анатомия функции

Базовый синтаксис для создания функции выглядит так:

my_function <- function(arg1, arg2, ...) {
  # Тело функции: блок кода, который выполняется
  # ...
  
  # Возвращаемое значение (результат работы функции)
  return(result)
}

Разберем по частям:

  1. Имя функции (my_function): Это имя, по которому вы будете вызывать свою функцию.
  2. Присвоение (<-): Мы создаем объект типа “функция” и присваиваем его имени.
  3. Ключевое слово function: Сообщает R, что мы создаем функцию.
  4. Аргументы (arg1, arg2, ...): Это входные данные, которые функция будет использовать. Их может быть любое количество (или ноль).
  5. Тело функции ({ ... }): Это код, который выполняется при вызове функции. Все, что находится внутри фигурных скобок, является телом функции.
  6. Возвращаемое значение (return()): Это результат, который функция “отдает” наружу после завершения работы.

4.1.2 Создание и вызов простой функции

Давайте создадим нашу первую функцию. Она будет принимать два числа и возвращать их сумму.

# Создаем функцию для сложения двух чисел
add_numbers <- function(num1, num2) {
  sum_result <- num1 + num2
  
  # Явно возвращаем результат с помощью return()
  return(sum_result)
}

Теперь вызовем нашу функцию:

# Вызываем функцию с аргументами 5 и 7
result <- add_numbers(5, 7)

print(result)

4.1.3 Неявное возвращение значения

В R есть удобная особенность: если вы не используете return(), функция автоматически вернет результат последнего выполненного выражения.

Перепишем нашу функцию в более идиоматичном для R стиле:

# Более короткий и распространенный вариант
add_numbers_short <- function(num1, num2) {
  num1 + num2 # Последнее выражение, его результат будет возвращен
}

# Вызываем новую функцию
result_short <- add_numbers_short(10, 20)
print(result_short)

4.1.4 Аргументы по умолчанию

Вы можете задать для аргументов значения по умолчанию. Это делает функцию более гибкой: если пользователь не предоставил значение для такого аргумента, будет использовано значение по умолчанию.

# Функция для приветствия с аргументом по умолчанию
greet <- function(name, greeting = "Hello") {
  paste(greeting, name, "!")
}

# Вызываем с двумя аргументами
greet("Alice")

# Вызываем, переопределяя аргумент по умолчанию
greet("Bob", greeting = "Goodbye")

4.1.5 Область видимости переменных (Scope)

Важно понимать, что переменные, созданные внутри функции, существуют только внутри этой функции. Они не видны извне. Это называется локальной областью видимости.

my_multiplier <- function(x) {
  internal_var <- 10 # Эта переменная существует только внутри функции
  result <- x * internal_var
  return(result)
}

# Вызываем функцию
my_multiplier(5)

# Попытка обратиться к internal_var извне вызовет ошибку
# print(internal_var) 
# Error: object 'internal_var' not found

Это защищает вашу рабочую среду от “засорения” временными переменными из функций.

4.2 4.2 Анонимные (лямбда) функции

Анонимная функция (часто называемая лямбда-функцией) — это функция без имени. Она создается “на лету” и обычно используется там, где нужна небольшая функция на один раз, например, в качестве аргумента для другой функции.

Они особенно полезны в пакетах purrr (часть Tidyverse) или в базовых функциях семейства apply для итерации по спискам и векторам.

В R есть два основных способа создавать такие функции “в одну строку”.

4.2.1 Способ 1: Классический синтаксис function()

Это самый явный и универсальный способ. Вы просто пишете ключевое слово function и ее тело.

function(arg1, arg2) {
  # тело функции
  arg1 + arg2
}

Пример использования с lapply():

Представим, что у нас есть список числовых векторов, и мы хотим вычислить среднее значение для каждого из них, вычтя перед этим определенное значение.

# Создаем список данных
data_list <- list(
  group_a = c(10, 12, 15),
  group_b = c(5, 6, 7, 8),
  group_c = c(100, 110)
)

# Используем lapply для применения функции к каждому элементу списка
# Мы хотим вычислить mean(x - 10) для каждого вектора x
# Вместо того чтобы создавать отдельную именованную функцию, мы используем анонимную
result_list <- lapply(data_list, function(x) {
  mean(x - 10)
})

print(result_list)

В этом примере function(x) { mean(x - 10) } — это лямбда-функция. Мы не дали ей имени, а сразу передали в lapply. lapply поочередно подставляет каждый элемент списка в аргумент x.

4.2.2 Способ 2: Короткий синтаксис с обратным слэшем \

Для очень простых анонимных функций, которые состоят из одного вызова другой функции, в базовом R есть удобный сокращенный синтаксис.

Синтаксис: \(аргументы) выражение

Обратный слэш \ — это просто короткий псевдоним для слова function.

Давайте перепишем предыдущий пример, используя этот синтаксис. Наша анонимная функция function(x) mean(x - 10) не подходит, так как она выполняет две операции (вычитание и mean). Но если бы нам нужно было просто применить одну функцию, этот синтаксис был бы идеален.

# Предположим, мы просто хотим применить mean() к каждому элементу списка
result_list_short <- lapply(data_list, \(x) mean(x))

print(result_list_short)

Этот код абсолютно идентичен следующему:

# То же самое, но написано полноценным синтаксисом
result_list_long <- lapply(data_list, function(x) mean(x))

Короткий синтаксис \(x) делает код более лаконичным и читаемым для простых случаев.

4.2.3 Лямбда-функции в Tidyverse (purrr)

В пакете purrr для этого есть свой, еще более популярный синтаксис с помощью “тильды” (~).

# install.packages("purrr")
library(purrr)

# Аналогичная операция с помощью map() из purrr
# ~ создает анонимную функцию, а .x является ее аргументом
result_purrr <- map(data_list, ~ mean(.x - 10))

print(result_purrr)

Синтаксис ~ из purrr очень мощный, так как он позволяет легко выполнять более сложные операции, например, с помощью “магических” точек ..1, ..2 для нескольких аргументов.

4.2.4 Какой способ выбрать?

  • function(...): Используйте для любых анонимных функций, особенно если их тело состоит из нескольких строк или сложной логики. Это самый универсальный и понятный способ.
  • \(...): Идеально подходит для быстрых, однострочных анонимных функций в базовом R. Делает код короче.
  • ~ ...: Стандарт для работы в экосистеме tidyverse. Очень лаконичен и мощен, но является частью purrr, а не базового R.

4.2.5 4.3 Функции семейства apply и их аналоги map

В программировании часто возникает задача применить одну и ту же операцию к каждому элементу сложного объекта (например, к каждой строке или столбцу таблицы, или к каждому элементу списка). Вместо того чтобы писать циклы for, в R предпочтительнее использовать функциональное программирование с помощью функций apply (из базового R) или map (из пакета purrr).

Эти функции принимают на вход: 1. Объект для обработки. 2. Функцию, которую нужно применить. 3. (Опционально) Дополнительные аргументы для этой функции.

4.2.5.1 Функции семейства apply (Базовый R)

Это “классический” способ итерации в R. Разные функции семейства предназначены для разных типов данных.

4.2.5.1.1 apply() — для матриц и массивов

Функция apply() применяется к строкам или столбцам матрицы.

  • MARGIN = 1: применить функцию к строкам.
  • MARGIN = 2: применить функцию к столбцам.
# Создадим матрицу с данными (например, 3 пациента, 4 показателя)
patient_matrix <- matrix(rnorm(12, mean = 100, sd = 15), nrow = 3, ncol = 4)
colnames(patient_matrix) <- c("Test1", "Test2", "Test3", "Test4")
rownames(patient_matrix) <- c("Patient_A", "Patient_B", "Patient_C")

print(patient_matrix)

# Рассчитаем средний балл для КАЖДОГО ПАЦИЕНТА (по строкам)
mean_scores_by_patient <- apply(patient_matrix, MARGIN = 1, mean)
print(mean_scores_by_patient)

# Рассчитаем среднее значение по КАЖДОМУ ТЕСТУ (по столбцам)
mean_by_test <- apply(patient_matrix, MARGIN = 2, mean)
print(mean_by_test)
4.2.5.1.2 lapply() — для списков и векторов

Функция lapply() (list-apply) применяет функцию к каждому элементу списка (или вектора) и возвращает список.

# Используем список данных из предыдущего примера
data_list <- list(
  group_a = c(10, 12, 15),
  group_b = c(5, 6, 7, 8),
  group_c = c(100, 110)
)

# Вычислим среднее для каждой группы
# lapply вернет список
lapply_result <- lapply(data_list, mean)
print(lapply_result)
4.2.5.1.3 sapply() — упрощенный lapply()

Функция sapply() (simplified apply) работает так же, как lapply(), но пытается упростить результат. Если результаты всех итераций имеют одинаковую длину (например, это все одиночные числа), sapply() вернет вектор или матрицу вместо списка.

# sapply "угадает", что результат можно упростить до вектора
sapply_result <- sapply(data_list, mean)
print(sapply_result)
class(sapply_result) # "numeric" (вектор)

Предупреждение: sapply() удобен, но его поведение может быть непредсказуемым, если результаты имеют разную длину или тип. В надежном коде часто предпочитают lapply() или функции из purrr, которые всегда возвращают предсказуемый тип объекта.

4.2.5.2 Функции map() из пакета purrr (Tidyverse)

Пакет purrr предоставляет современный, более последовательный и мощный аналог семейства apply. Все функции map() всегда возвращают один и тот же тип объекта: map() возвращает список, map_dbl() — числовой вектор, map_df() — фрейм данных и т.д. Это делает код более надежным.

4.2.5.2.1 map() и его варианты
  • map(.x, .f): Возвращает список. Прямой аналог lapply().
  • map_dbl(.x, .f): Возвращает числовой вектор (double).
  • map_chr(.x, .f): Возвращает символьный вектор (character).
  • map_df(.x, .f, .id = "col_name"): Возвращает фрейм данных, объединяя результаты по строкам.
# install.packages("purrr")
library(purrr)

# Используем тот же список data_list

# Аналог lapply() - возвращает список
map_result_list <- map(data_list, mean)
print(map_result_list)

# Аналог sapply() - возвращает числовой вектор. Тип результата гарантирован.
map_result_vector <- map_dbl(data_list, mean)
print(map_result_vector)
4.2.5.2.2 Итерация по нескольким объектам одновременно

Одно из главных преимуществ purrr — это удобная работа с несколькими векторами одновременно с помощью map2() и pmap().

map2(.x, .y, .f) принимает два списка/вектора и применяет функцию к парам их элементов.

# Допустим, у нас есть списки значений и порогов
values <- list(a = c(1, 2, 3), b = c(10, 20))
thresholds <- list(a = 2, b = 15)

# Мы хотим для каждого списка посчитать, сколько элементов больше порога
# .x - элемент из `values`, .y - элемент из `thresholds`
map2_dbl(values, thresholds, ~ sum(.x > .y))

pmap(.l, .f) (parallel map) делает то же самое, но для любого количества списков. Она принимает список списков (.l).

# Допустим, мы хотим сгенерировать данные для линейной модели y = a*x + b
params <- list(
  a = c(1, 2, 3),
  x = list(c(1, 2, 3), c(4, 5, 6), c(7, 8, 9)),
  b = c(0, 5, 10)
)

# Функция для расчета
calculate_y <- function(a, x, b) {
  a * x + b
}

# pmap автоматически передает элементы из каждого списка в функцию calculate_y
pmap(params, calculate_y)

4.2.6 Итог: apply vs map

Характеристика apply (Базовый R) map (purrr)
Надежность sapply() может возвращать разные типы (вектор, матрица, список), что непредсказуемо. Каждый вариант (map_dbl, map_chr и т.д.) гарантированно возвращает определенный тип.
Синтаксис Разные функции для разных задач (apply, lapply, mapply). Единообразный синтаксис (map, map2, pmap).
Лямбда-функции function(x) ... или короткий \(x) .... ~ .x + 1 (очень лаконично).
Итерация по >1 объекту mapply() (синтаксис может быть запутанным). map2() и pmap() (очень интуитивно).
Рекомендация Полезно знать для понимания старого кода. Предпочтительный современный подход. Более надежный, читаемый и мощный.

4.2.7 4.4 Пайп-операторы (|> и %>%) и другие инфиксные операторы

Пайп (от англ. pipe — “труба”) — это мощный инструмент, который позволяет делать ваш код более линейным и читаемым. Вместо того чтобы вкладывать функции друг в друга (что приводит к “пирогам” из скобок) или создавать множество промежуточных переменных, вы можете передавать результат одной функции напрямую в качестве первого аргумента для следующей.

4.2.7.1 Оператор |> (Native Pipe, Пайп базового R)

Начиная с версии R 4.1.0, в базовый R был встроен собственный пайп-оператор |>. Он работает очень просто: он берет результат выражения слева и передает его в качестве первого аргумента в функцию справа.

Синтаксис: данные |> функция(другие_аргументы)

Это эквивалентно записи: функция(данные, другие_аргументы).

Давайте сравним “традиционный” подход и подход с пайпом.

# Создадим фрейм данных для примера
df <- data.frame(
  id = 1:5,
  score = c(88, 92, 75, 99, 68),
  group = c("A", "B", "A", "B", "A")
)

# --- Традиционный подход (вложенные функции) ---
# Это сложно читать из-за вложенности и скобок
result_nested <- head(df[order(df$score, decreasing = TRUE), ], 3)

print(result_nested)

# --- Подход с пайпом |> ---
# Код читается как последовательность действий: "взять df, отсортировать, взять первые 3 строки"
result_piped <- df |>
  arrange(desc(score)) |>
  head(3)

Ошибка в коде выше: Функция arrange и desc не являются базовыми функциями R. Они принадлежат пакету dplyr. Для корректной работы этого примера нужно загрузить пакет dplyr. ::: {.cell}

# install.packages("dplyr")
library(dplyr)

result_piped <- df |>
  arrange(desc(score)) |>
  head(3)

print(result_piped)

:::

4.2.7.2 Оператор %>% (Magrittr Pipe, Пайп из Tidyverse)

Оператор %>% был популяризирован пакетом magrittr и является неотъемлемой частью экосистемы Tidyverse. Он работает аналогично |>, но с некоторыми дополнительными возможностями.

Синтаксис: данные %>% функция(другие_аргументы)

Основное отличие от |> заключается в том, что %>% по умолчанию также подставляет данные в качестве первого аргумента, но позволяет явно указать другое место с помощью специального символа . (точка).

# install.packages("magrittr")
library(magrittr)
library(dplyr)

# Пример использования точки (.)
# Мы хотим построить график, где данные - это df, но заголовок должен быть основан на df
df %>%
  {ggplot(., aes(x = group, y = score)) +
    geom_boxplot() +
    ggtitle(paste("Распределение очков для", nrow(.), "пациентов"))}

В этом примере . внутри ggtitle() представляет собой все еще исходный фрейм данных df, что позволяет нам динамически создать заголовок на основе его свойств.

4.2.7.3 %>% vs |>: Что выбрать?

  • |> (Native Pipe):
    • Доступен “из коробки” в современных версиях R.
    • Быстрее.
    • Простой и предсказуемый синтаксис.
    • Рекомендация для новых проектов, если вы не используете специфические возможности %>%.
  • %>% (Magrittr Pipe):
    • Требует загрузки пакета (library(magrittr) или library(dplyr)).
    • Более гибкий благодаря использованию точки ..
    • Имеет другие полезные операторы, такие как %T>% (tee pipe) и %$% (exposition pipe).
    • Необходим для работы со старым кодом Tidyverse.

В большинстве случаев они взаимозаменяемы.

4.2.7.4 Другие инфиксные операторы %(operator)%

Синтаксис %(operator)% не уникален для пайпов. Это общий способ создания инфиксных операторов — операторов, которые размещаются между своими операндами (как +, *, ==). Вы можете создавать свои собственные операторы!

Функция, которая определяет такой оператор, должна называться `%оператор%`.

Давайте создадим свой собственный оператор, например, для возведения в степень, который будет называться %^%.

# Создаем свою функцию-оператор
`%^%` <- function(base, exponent) {
  base ^ exponent
}

# Теперь мы можем использовать его как встроенный оператор
5 %^% 3 # Эквивалентно 5 ^ 3

# Можно использовать с векторами
c(2, 3, 4) %^% c(1, 2, 3)

Это демонстрирует гибкость R: вы можете расширять язык, создавая собственные синтаксические конструкции для удобства и выразительности кода. Пайп-операторы %>% и |> — самые известные примеры таких инфиксных операторов.

4.2.8 4.5 Возврат значений и вывод информации

Когда мы пишем функции или просто выполняем код, нам нужно не только производить вычисления, но и управлять результатами: возвращать значения из функций и выводить полезную информацию для пользователя или для отладки. В R для этого есть несколько инструментов.

4.2.8.1 4.5.1 Функция return(): Явный возврат значения

Как мы уже обсуждали в подразделе 4.1, return() — это стандартный способ явно указать, какой результат должна вернуть функция. Когда R встречает return(), она немедленно прекращает выполнение функции и возвращает указанное значение.

my_function <- function(x) {
  # ... какие-то вычисления ...
  result <- x * 2
  
  # Явно возвращаем результат
  return(result)
  
  # Этот код никогда не выполнится
  print("Это не будет напечатано")
}

return() особенно полезен в функциях с условной логикой, где результат может быть вычислен в разных местах.

# Функция, которая определяет тип числа
check_number <- function(num) {
  if (num > 0) {
    return("Положительное")
  } else if (num < 0) {
    return("Отрицательное")
  } else {
    return("Ноль")
  }
}

print(check_number(10))
print(check_number(-5))

4.2.8.2 4.5.2 Функции вывода: print() и cat()

Для вывода информации на консоль (в панель “Console” в RStudio) существуют две основные функции, которые ведут себя по-разному.

4.2.8.2.2 cat(): Конкатенация и вывод

cat() (от concatenate and print) — это функция для вывода текста. Она преобразует свои аргументы в строки, “склеивает” их (если указан разделитель sep) и выводит на консоль.

Ключевые отличия cat() от print():

  1. Нет кавычек: cat() выводит текст “как есть”.
  2. Нет переноса строки: cat() не добавляет перенос строки в конце, если вы не попросите это явно с помощью \n.
  3. Аргумент sep: Позволяет указать разделитель между аргументами.
  4. Возвращает NULL: cat() ничего не возвращает, ее нельзя использовать в промежуточных вычислениях.
# cat() выводит текст без кавычек
cat(my_text)

# Добавляем перенос строки вручную
cat(my_text, "\n")

# Используем разделитель
name <- "Иван"
age <- 30
cat("Пациент:", name, "Возраст:", age, "лет.\n", sep = " ")

4.2.8.4 4.5.4 Другие полезные функции вывода

  • message(): Выводит информационное сообщение. По умолчанию выводится оранжевым цветом в RStudio. Используется для информирования пользователя о важных, но не критических событиях. Сообщения можно подавлять с помощью suppressMessages().
  • warning(): Выводит предупреждение. По умолчанию выводится желтым цветом. Говорит о потенциальной проблеме, но не останавливает выполнение кода. Предупреждения можно подавлять с помощью suppressWarnings().
  • stop(): Выводит сообщение об ошибке и полностью останавливает выполнение текущей операции. Используется, когда дальнейшее выполнение невозможно или бессмысленно (например, если функция получила неверный тип данных).
# Пример с warning()
log_transform <- function(x) {
  if (any(x < 0)) {
    warning("В данных есть отрицательные значения. Результат для них будет NaN.")
  }
  return(log(x))
}
log_transform(c(1, 5, -2, 10))
Warning in log_transform(c(1, 5, -2, 10)): В данных есть отрицательные
значения. Результат для них будет NaN.
Warning in log(x): созданы NaN
[1] 0.000000 1.609438      NaN 2.302585
# Пример с stop()
divide_numbers <- function(a, b) {
  if (b == 0) {
    stop("Деление на ноль невозможно!")
  }
  return(a / b)
}
# divide_numbers(10, 0) # Эта строка вызовет ошибку и остановит выполнение

5 Глава 5: Циклы

Циклы позволяют выполнять один и тот же блок кода многократно. Это полезно, когда нужно повторить операцию для каждого элемента вектора, списка или для определенного числа итераций.

Хотя в R для многих задач предпочтительнее использовать функциональное программирование (функции apply или map из Главы 4), понимание циклов необходимо для написания сложных алгоритмов и для случаев, когда порядок выполнения операций имеет значение.

В R есть два основных типа циклов: for и while.

5.1 5.1 Цикл for

Цикл for используется, когда вы заранее знаете, сколько раз нужно выполнить итерацию. Он проходит по каждому элементу вектора или последовательности.

5.1.1 Синтаксис

for (variable in sequence) {
  # Код, который выполняется на каждой итерации
}
  • variable: Переменная, которая на каждой итерации будет принимать очередное значение из sequence.
  • sequence: Последовательность, по которой идет итерация (например, вектор c(1, 2, 3) или последовательность 1:10).

5.1.2 Примеры использования

1. Итерация по числовой последовательности

Самый частый случай — повторить что-то N раз.

# Создадим пустой вектор, в который будем складывать результаты
results <- c()

# Цикл от 1 до 5
for (i in 1:5) {
  # На каждой итерации i будет равно 1, затем 2, 3, 4, 5
  results[i] <- i * 10
  print(paste("Итерация:", i, "Результат:", results[i]))
}

print(results)

2. Итерация по элементам вектора

Цикл for может перебирать не только числа, но и любые элементы вектора.

# Вектор с названиями групп
groups <- c("Control", "Treatment A", "Treatment B")

for (group_name in groups) {
  # На каждой итерации group_name будет равно "Control", затем "Treatment A" и т.д.
  print(paste("Анализируем группу:", group_name))
}

3. Итерация по именам столбцов фрейма данных

Это очень распространенная задача. Мы можем получить имена столбцов с помощью names() и перебрать их в цикле.

# Используем наш фрейм данных
df <- data.frame(
  id = 1:5,
  score = c(88, 92, 75, 99, 68),
  group = c("A", "B", "A", "B", "A")
)

# Итерируемся по именам столбцов
for (col_name in names(df)) {
  # Получаем доступ к столбцу с помощью [[ ]]
  column_data <- df[[col_name]]
  
  print(paste("Столбец:", col_name, "Тип данных:", class(column_data)))
}

Важно: Для доступа к столбцу по его имени в виде строки (которое хранится в переменной col_name) используется оператор [[ ]], а не $.

5.2 5.2 Цикл while

Цикл while выполняется до тех пор, пока истинно некоторое логическое условие. Он используется, когда вы не знаете точное число итераций заранее.

5.2.1 Синтаксис

while (condition) {
  # Код, который выполняется, пока condition == TRUE
}

Внимание! Убедитесь, что внутри цикла есть код, который в конечном итоге сделает условие ложным. Иначе вы создадите бесконечный цикл.

5.2.2 Пример использования

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

# Начальное значение
current_value <- 1

# Счетчик итераций, чтобы избежать бесконечного цикла
iteration <- 0

# Цикл будет работать, пока current_value меньше 100
while (current_value < 100) {
  # Увеличиваем значение
  current_value <- current_value * 2
  
  # Увеличиваем счетчик
  iteration <- iteration + 1
  
  print(paste("Итерация:", iteration, "Текущее значение:", current_value))
}

print(paste("Цикл завершен. Финальное значение:", current_value))

5.3 5.3 Управление циклом: break и next

Иногда нужно прервать цикл досрочно или пропустить одну итерацию. Для этого служат операторы break и next.

  • break: Немедленно прерывает выполнение всего цикла.
  • next: Немедленно прекращает текущую итерацию и переходит к следующей.

5.3.1 Пример с break

Найдем первый элемент в векторе, который больше 50, и остановим поиск.

numbers <- c(10, 23, 55, 41, 88, 12)

for (num in numbers) {
  if (num > 50) {
    print(paste("Нашли число больше 50:", num))
    break # Прерываем цикл
  }
}

5.3.2 Пример с next

Напечатаем квадраты всех чисел от 1 до 10, но пропустим числа, которые делятся на 3 без остатка.

for (i in 1:10) {
  if (i %% 3 == 0) { # %% - оператор остатка от деления
    next # Пропускаем остаток итерации и переходим к следующему i
  }
  
  print(paste("Квадрат числа", i, "равен", i^2))
}

5.4 5.4 Когда использовать циклы, а когда apply/map?

  • Используйте apply/map, когда вы хотите применить одну и ту же функцию к каждому элементу списка/вектора/столбцу и результат не зависит от предыдущих итераций. Этот подход обычно быстрее и приводит к более читаемому “функциональному” стилю кода.
  • Используйте циклы, когда:
    • Результат на шаге i зависит от результата на шаге i-1 (например, итеративные процессы, симуляции).
    • Внутри цикла сложная логика с множеством условий (if/else).
    • Вам нужно досрочно прервать выполнение на основе некоторого условия (break).

В анализе данных старайтесь отдавать предпочтение функциям apply/map, но не бойтесь использовать циклы там, где они более уместны.

Конечно! Это очень важный аспект при работе с циклами, особенно при отладке или обработке больших объемов данных. Добавляем этот подраздел в Главу 5.

5.5 5.5 Некоторые функции вывода внутри циклов

Циклы выполняют итеративные процессы, и часто бывает полезно или даже необходимо получать информацию на каждом шаге. Это может быть для отладки (чтобы понять, что происходит) или для информирования пользователя (например, о прогрессе длительной операции). Здесь функции вывода, которые мы рассмотрели в Главе 4, становятся незаменимыми.

5.5.0.1 Отладка и мониторинг с помощью print() и cat()

Самый частый сценарий использования вывода внутри циклов — это отслеживание значений переменных на каждой итерации.

  • print() удобен, чтобы быстро посмотреть на структуру объекта, который изменяется в цикле.
  • cat() идеален для создания форматированных, информативных строк, которые описывают текущий шаг.

Пример: Допустим, мы моделируем процесс накопления капитала с процентами и хотим видеть, как меняется баланс каждый год.

# Начальные параметры
balance <- 1000
interest_rate <- 0.07
years <- 5

cat("Начальный баланс:", balance, "\n")
cat("Годовая ставка:", interest_rate * 100, "%\n\n")

# Цикл для расчета баланса на 5 лет
for (year in 1:years) {
  balance <- balance * (1 + interest_rate)
  
  # Используем cat() для вывода информационного сообщения каждый год
  cat("В конце года", year, "баланс составит:", round(balance, 2), "\n")
}

В этом примере cat() создает понятный и читаемый отчет о процессе, который происходит внутри цикла.

5.5.0.2 Создание простого прогресс-бара

Для длительных циклов, которые могут выполняться секунды или минуты, полезно показывать пользователю, что программа не зависла, а работает.

Пример: Допустим, мы обрабатываем 10 файлов, и каждый обрабатывается некоторое время.

total_files <- 10
processed_files <- list() # Пустой список для результатов

for (i in 1:total_files) {
  # Имитация длительной обработки
  Sys.sleep(0.3) 
  
  # Добавляем "обработанные" данные в список
  processed_files[[i]] <- paste("Результат для файла", i)
  
  # Выводим прогресс
  cat("\rОбработано файлов:", i, "из", total_files) # \r возвращает курсор в начало строки
}

cat("\nОбработка завершена!\n")

Здесь \r (carriage return) используется, чтобы каждая новая строка вывода заменяла предыдущую, создавая эффект обновляющегося счетчика.

5.5.0.3 Предотвращение бесконечного вывода

Будьте осторожны с print() внутри циклов, который может выполняться много тысяч раз. Если вы распечатываете целый фрейм данных на каждой из 10 000 итераций, ваша консоль мгновенно заполнится, и производительность упадет.

Используйте условные операторы, чтобы выводить информацию только на определенных шагах.

# Выводим информацию только каждые 100 шагов
for (i in 1:500) {
  # ... какая-то обработка ...
  
  if (i %% 100 == 0) { # %% - оператор остатка от деления
    cat("Выполнено", i, "итераций из 500\n")
  }
}

5.5.0.4 Использование message() и warning()

  • message() — отличный способ сообщать о прогрессе, не загромождая основной вывод. В RStudio такие сообщения выделяются цветом и их можно легко скрыть, если они не нужны.
  • warning() — полезен, если в ходе итерации обнаружилась потенциальная проблема, но цикл должен продолжить работу.

Пример: Читаем список файлов, и некоторые из них могут отсутствовать.

# Создадим список файлов, где один файл "отсутствует"
file_list <- c("data1.csv", "data_missing.csv", "data2.csv")

for (file_name in file_list) {
  # Имитация проверки наличия файла
  if (file_name == "data_missing.csv") {
    warning(paste("Файл", file_name, "не найден и будет пропущен."))
    next # Переходим к следующей итерации
  }
  
  # Имитация успешной обработки
  message(paste("Файл", file_name, "успешно обработан."))
}

5.5.1 Итог

  • Используйте cat() для создания кастомных информационных сообщений и простых прогресс-баров.
  • Используйте print() для отладки, когда вам нужно быстро посмотреть на структуру объекта на определенном шаге.
  • Используйте message() для стандартных информационных сообщений о прогрессе.
  • Используйте warning() для уведомления о некритических проблемах.
  • Используйте stop() внутри цикла (часто с break), если обнаружена критическая ошибка, требующая немедленной остановки.

6 Глава 6: Элементарные операции с данными

В этой главе мы рассмотрим основные операции для изменения структуры и содержимого данных. Мы будем работать с тремя основными структурами: вектором, матрицей и фреймом данных (таблицей).

6.1 6.1 Манипуляции с векторами, матрицами и таблицами

6.1.1 Операции с векторами

Вектор — это самая простая структура. Большинство операций с ним интуитивно понятны.

# Создадим пример вектора
scores <- c(88, 92, 75, 99, 68)
names(scores) <- c("Anna", "Peter", "Maria", "John", "Sofia")
print(scores)

6.1.1.1 Сортировка

Для сортировки используется функция sort().

# Сортировка по значению (по возрастанию)
sorted_scores <- sort(scores)
print(sorted_scores)

# Сортировка по убыванию
sorted_scores_desc <- sort(scores, decreasing = TRUE)
print(sorted_scores_desc)

6.1.1.2 Поиск по индексу / имени

Доступ к элементам вектора осуществляется с помощью квадратных скобок [].

# Поиск по числовому индексу (нумерация с 1)
scores[2] # Второй элемент (92)

# Поиск по имени (если есть)
scores["Maria"] # Элемент с именем "Maria" (75)

# Поиск по нескольким индексам или именам
scores[c(1, 3, 5)] # Элементы 1, 3 и 5

6.1.1.3 Добавление нового элемента

Для добавления элемента используется функция c() (конкатенация).

# Добавляем новый элемент в конец вектора
scores["David"] <- 95
print(scores)

6.1.1.4 Замена элемента

Чтобы заменить элемент, обратитесь к нему по индексу или имени и присвойте новое значение.

# Заменяем элемент по индексу
scores[1] <- 90

# Заменяем элемент по имени
scores["John"] <- 100

print(scores)

6.1.1.5 Удаление элемента

Для удаления элемента используется знак минус - внутри квадратных скобок.

# Удаляем третий элемент (Maria)
scores <- scores[-3]
print(scores)

# Удаляем несколько элементов
scores <- scores[-c(1, 3)] # Удаляем элементы с индексами 1 и 3 (Anna и Sofia)
print(scores)

6.1.2 Операции с матрицами и таблицами (фреймами данных)

Матрицы и фреймы данных (далее — таблицы) имеют два измерения: строки и столбцы. Для доступа к ним используется синтаксис [строки, столбцы].

# Создадим пример матрицы
mat <- matrix(1:9, nrow = 3, byrow = TRUE)
rownames(mat) <- c("R1", "R2", "R3")
colnames(mat) <- c("C1", "C2", "C3")
print(mat)

# Создадим пример фрейма данных
df <- data.frame(
  id = 1:4,
  name = c("A", "B", "C", "D"),
  score = c(88, 92, 75, 99),
  passed = c(TRUE, TRUE, FALSE, TRUE)
)
print(df)

6.1.2.1 Сортировка

Сортировка таблиц — это частая задача. Лучший способ — использовать функции из пакета dplyr.

# install.packages("dplyr")
library(dplyr)

# Сортировка таблицы df по столбцу 'score' по возрастанию
df_sorted_asc <- df %>% arrange(score)
print(df_sorted_asc)

# Сортировка по нескольким столбцам и по убыванию
df_sorted_desc <- df %>% arrange(desc(passed), score)
print(df_sorted_desc)

6.1.2.2 Поиск по индексу / имени

Синтаксис [строки, столбцы] очень гибок.

# --- Работа со столбцами ---

# Получить один столбец как вектор (самый простой способ)
df$name

# Получить один столбец как фрейм данных
df["name"]

# Получить несколько столбцов
df[, c("id", "score")]

# --- Работа со строками ---

# Получить одну строку
df[2, ]

# Получить несколько строк
df[c(1, 3), ]

# --- Комбинированный доступ ---

# Получить значение из второй строки, третьего столбца
df[2, 3]

# Получить столбцы 'name' и 'score' для строк, где score > 80
df[df$score > 80, c("name", "score")]

6.1.2.3 Добавление и замена строк и колонок

Добавление/замена колонок:

# Добавить новую колонку 'grade' на основе существующей 'score'
df$grade <- ifelse(df$score >= 90, "A", "B")
print(df)

# Заменить существующую колонку
df$score <- df$score * 1.05 # Добавим 5% к каждому баллу
print(df)

Добавление строк:

Для добавления строк используется функция rbind() (row bind). Новая строка должна иметь такое же количество и типы столбцов.

# Создаем новую строку как фрейм данных
new_row <- data.frame(id = 5, name = "E", score = 85, passed = TRUE, grade = "B")

# Добавляем строку в конец таблицы
df <- rbind(df, new_row)
print(df)

6.1.2.4 Удаление строк и колонок

Удаление колонок:

# Удалить колонку 'grade' путем присвоения NULL
df$grade <- NULL
print(df)

# Удалить несколько колонок по номерам
df <- df[, -c(4, 5)] # Удаляем 4-й и 5-й столбцы (passed и grade)
print(df)

Удаление строк:

# Удалить строку по номеру
df <- df[-2, ] # Удаляем вторую строку

# Удалить строки по условию (самый частый случай)
# Удаляем все строки, где score < 90
df <- df[df$score >= 90, ]
print(df)

6.1.3 Сравнение: Базовый R vs. dplyr (Tidyverse)

Хотя базовый R позволяет выполнять все необходимые манипуляции с данными, пакет dplyr из экосистемы Tidyverse предлагает альтернативный синтаксис, который часто оказывается более читаемым, последовательным и менее подверженным ошибкам, особенно в сложных цепочках преобразований.

Давайте сравним оба подхода на конкретных задачах. Мы будем использовать наш фрейм данных df.

# Восстановим исходный фрейм данных для чистоты эксперимента
df <- data.frame(
  id = 1:4,
  name = c("A", "B", "C", "D"),
  score = c(88, 92, 75, 99),
  passed = c(TRUE, TRUE, FALSE, TRUE)
)

6.1.3.1 1. Выбор колонок (Select Columns)

Задача: Выбрать только колонки name и score.

  • Базовый R: Используется числовая индексация или вектор имен внутри [].

    # По номерам колонок
df_base_1 <- df[, c(2, 3)]

# По именам колонок (предпочтительнее)
df_base_2 <- df[, c("name", "score")]

print(df_base_2)

  • dplyr: Используется функция select().

    # install.packages("dplyr")
library(dplyr)

df_dplyr <- df %>% select(name, score)
print(df_dplyr)

Преимущество dplyr: Код более декларативен (“выбрать name и score”), а не процедурен (“взять 2-й и 3-й столбцы”). Не нужно запоминать синтаксис с запятой.

6.1.3.2 2. Фильтрация строк (Filter Rows)

Задача: Оставить только тех студентов, кто сдал тест (passed == TRUE).

  • Базовый R: Используется логическая индексация. Нужно дважды написать имя фрейма данных.

    df_base <- df[df$passed == TRUE, ]
print(df_base)

  • dplyr: Используется функция filter().

    df_dplyr <- df %>% filter(passed == TRUE)
print(df_dplyr)

Преимущество dplyr: Код короче и читабельнее. Не нужно использовать $ и запятые. В filter() можно легко комбинировать условия через запятую (которая работает как логическое “И”): filter(passed == TRUE, score > 90).

6.1.3.3 3. Сортировка (Arrange / Sort)

Задача: Отсортировать студентов по убыванию балла (score).

  • Базовый R: Используется функция order(), которая возвращает индексы для правильного порядка.

    df_base <- df[order(df$score, decreasing = TRUE), ]
print(df_base)

  • dplyr: Используется функция arrange().

    df_dplyr <- df %>% arrange(desc(score))
print(df_dplyr)

Преимущество dplyr: Синтаксис более интуитивен. arrange(desc(score)) прямо говорит “отсортировать по убыванию score”, в то время как order() требует большего количества “шумовых” символов (df$, decreasing = TRUE, , ]).

6.1.3.4 4. Создание/изменение колонок (Create / Modify Columns)

Задача: Создать новую колонку grade (оценка), где “A” за балл >= 90 и “B” в остальных случаях.

  • Базовый R: Прямое присвоение нового столбца.

    df_base <- df
df_base$grade <- ifelse(df_base$score >= 90, "A", "B")
print(df_base)

  • dplyr: Используется функция mutate().

    df_dplyr <- df %>% mutate(grade = ifelse(score >= 90, "A", "B"))
print(df_dplyr)

Преимущество dplyr: mutate() позволяет создавать несколько новых колонок одновременно, при этом новые колонки можно использовать в том же вызове mutate(). Это делает код более компактным.

6.1.3.5 5. Суммирование / Агрегация (Summarize / Aggregate)

Задача: Рассчитать средний балл (score) и количество студентов (id) в каждой группе passed.

  • Базовый R: Используется функция aggregate(). Синтаксис может быть непривычным.

    df_base <- aggregate(
  cbind(score, id) ~ passed, # Формула: зависимые переменные ~ группирующая переменная
  data = df,
  FUN = function(x) c(mean = mean(x), count = length(x))
)
print(df_base)

  • dplyr: Используется комбинация group_by() и summarize(). Это “святой грааль” dplyr.

    df_dplyr <- df %>%
  group_by(passed) %>%
  summarize(
    avg_score = mean(score),
    student_count = n() # n() - специальная функция для подсчета строк в группе
  )
print(df_dplyr)

Преимущество dplyr: Код читается как последовательность мыслей: “взять df, сгруппировать по passed, затем суммировать, рассчитав средний score и количество строк”. Результат получается в аккуратном “плоском” фрейме данных.

6.1.3.6 6. Цепочка операций (The Pipe)

Задача: Выполнить несколько операций подряд: отфильтровать студентов, создать новую колонку и отсортировать.

  • Базовый R: Требует вложенности или создания множества промежуточных переменных.

    # Вариант 1: Вложенность (трудно читать)
df_base_nested <- df[order(
  ifelse(df[df$passed == TRUE, ]$score >= 90, "A", "B"),
  decreasing = TRUE
), ]

# Вариант 2: Промежуточные переменные (длинно)
df_temp_1 <- df[df$passed == TRUE, ]
df_temp_1$grade <- ifelse(df_temp_1$score >= 90, "A", "B")
df_base_chain <- df_temp_1[order(df_temp_1$grade, decreasing = TRUE), ]
print(df_base_chain)

  • dplyr: С помощью пайпа (%>% или |>) код становится линейным и легко читаемым.

    df_dplyr_chain <- df %>%
  filter(passed == TRUE) %>%
  mutate(grade = ifelse(score >= 90, "A", "B")) %>%
  arrange(desc(grade))

print(df_dplyr_chain)

Преимущество dplyr: Пайп позволяет писать код в логической последовательности действий, избегая вложенных скобок и временных переменных. Это делает его невероятно читаемым и легким для отладки.

6.1.4 Итог: Когда что использовать?

Критерий Базовый R dplyr (Tidyverse)
Производительность Очень быстр для простых операций. Очень быстр (использует C++), но может быть медленнее на очень маленьких данных из-за накладных расходов.
Читаемость Может быть громоздким, особенно при сложных операциях. Высокая, особенно при использовании пайпов. Код читается как повествование.
Последовательность Разные функции для похожих задач ([], $, subset, aggregate). Единообразный набор “глаголов” (select, filter, mutate, arrange, summarize).
Надежность Легко допустить ошибку с запятыми в []. Меньше “магии” и синтаксических ловушек.
Рекомендация Полезно знать для понимания основ R, работы с матрицами и написания функций. Стандарт де-факто для анализа данных. Предпочтительный выбор для большинства задач по манипуляции таблицами.

6.2 6.2 Транспонирование и “сводные таблицы” (Pivot)

Помимо базовых манипуляций со строками и столбцами, часто требуется кардинально изменить структуру таблицы. Две ключевые операции для этого — это транспонирование и “пивотинг” (изменение формы данных из длинной в широкую и обратно).

6.2.1 Транспонирование

Транспонирование — это простая операция, которая меняет строки и столбцы местами. Строка i становится столбцом i, а столбец j становится строкой j.

6.2.1.1 Базовый R: t()

В базовом R для этого используется функция t().

# Создадим матрицу
m <- matrix(1:6, nrow = 2, byrow = TRUE)
rownames(m) <- c("R1", "R2")
colnames(m) <- c("C1", "C2", "C3")
print("Исходная матрица:")
print(m)

# Транспонируем матрицу
m_t <- t(m)
print("Транспонированная матрица:")
print(m_t)

Важное замечание для фреймов данных: При транспонировании фрейма данных с помощью t() результатом будет матрица, а не фрейм данных. Все данные в ней будут приведены к одному общему типу (обычно к символьному character), чтобы избежать потери информации.

# Создадим фрейм данных
df <- data.frame(
  id = 1:3,
  name = c("A", "B", "C"),
  score = c(88, 92, 75)
)
print("Исходный фрейм данных:")
print(str(df))

# Транспонируем его
df_t <- t(df)
print("Транспонированный объект (матрица):")
print(str(df_t))

Как видите, все значения стали текстовыми (chr). Это не всегда желаемое поведение.

6.2.2 Pivot: Изменение формы данных

“Пивотинг” — это более сложная и мощная операция, чем простое транспонирование. Она меняет форму данных с “широкой” на “длинную” и наоборот.

  • Широкий формат (Wide): Каждая переменная имеет свой столбец. Удобно для чтения человеком.
  • Длинный формат (Long/Tidy): Каждая строка — это одно наблюдение. Удобно для анализа и визуализации в R.

Для этих операций лучше всего использовать функции из пакетов tidyr и data.table. Мы сосредоточимся на tidyr как на части Tidyverse.

Представим, что у нас есть данные о продажах разных типов продуктов в разные годы.

# Данные в ШИРОКОМ формате
sales_wide <- data.frame(
  year = c(2021, 2022, 2023),
  product_A = c(100, 120, 150),
  product_B = c(80, 90, 110),
  product_C = c(50, 55, 65)
)

print("Широкий формат:")
print(sales_wide)

6.2.2.1 pivot_longer() — из широкого в длинный

Эта операция “собирает” несколько столбцов в два: один с названиями исходных столбцов (ключ), а другой — с их значениями.

  • cols: Столбцы, которые нужно “собрать”.
  • names_to: Название нового столбца для имен (ключей).
  • values_to: Название нового столбца для значений.
# install.packages("tidyr")
library(tidyr)

sales_long <- sales_wide %>%
  pivot_longer(
    cols = c(product_A, product_B, product_C), # Собираем столбцы с продуктами
    names_to = "product",                      # Имена столбцов (product_A) пойдут в столбец 'product'
    values_to = "sales"                        # Значения пойдут в столбец 'sales'
  )

print("Длинный формат:")
print(sales_long)

Теперь данные в “аккуратном” (tidy) формате: одна строка на одну продажу продукта за один год.

6.2.2.2 pivot_wider() — из длинного в широкий

Это обратная операция. Она “разворачивает” два столбца (ключ и значение) в множество столбцов.

  • names_from: Столбец, из которого будут взяты имена для новых столбцов.
  • values_from: Столбец, из которого будут взяты значения для новых столбцов.
# Вернемся к широкому формату из нашего длинного
sales_wide_again <- sales_long %>%
  pivot_wider(
    names_from = product,  # Имена новых столбцов берем из 'product'
    values_from = sales    # Значения для них берем из 'sales'
  )

print("Снова широкий формат:")
print(sales_wide_again)

6.2.3 Когда это полезно?

  • pivot_longer:
    • Визуализация в ggplot2: ggplot2 работает лучше всего с длинными данными. Чтобы построить график продаж всех продуктов на одном графике, вам понадобится длинный формат.
    • Статистическое моделирование: Большинство функций для моделирования ожидают данные в длинном формате.
    • Объединение данных: Легче объединить (join) два набора данных по одному общему ключу, когда они в длинном формате.
  • pivot_wider:
    • Создание отчетов для людей: Широкий формат часто легче читать и воспринимать в виде таблицы.
    • Сводные таблицы: Создание матриц, похожих на те, что можно сделать в Excel или Google Sheets.

6.2.4 Итог

Операция Функция Назначение
Транспонирование t() Простая замена строк и столбцов.
Из широкого в длинный pivot_longer() Сборка нескольких столбцов в два столбца (ключ, значение).
Из длинного в широкий pivot_wider() Разворачивание двух столбцов (ключ, значение) в несколько столбцов.

Освоение pivot_longer() и pivot_wider() — это важный шаг на пути к свободному владению R для анализа данных.

6.3 6.3 Элементарные математические операции

R — это в первую очередь статистическая среда, и ее возможности для выполнения математических вычислений очень широки. Давайте рассмотрим основные операции и то, как они применяются к разным структурам данных.

6.3.1 6.3.1 Операции над отдельными значениями (скалярами)

С одиночными числами (в R они являются векторами длины 1) все работает интуитивно понятно.

Оператор Описание Пример
+ Сложение 5 + 3
- Вычитание 10 - 4
* Умножение 7 * 8
/ Деление 20 / 5
^ или ** Возведение в степень 3 ^ 2 или 3 ** 2
%% Остаток от деления 10 %% 3 (результат 1)
%/% Целочисленное деление 10 %/% 3 (результат 3) 
# Примеры со скалярами
a <- 10
b <- 3

sum_ab <- a + b
diff_ab <- a - b
prod_ab <- a * b
div_ab <- a / b
power_ab <- a ^ b
mod_ab <- a %% b
int_div_ab <- a %/% b

print(paste("Сумма:", sum_ab))
print(paste("Разность:", diff_ab))
print(paste("Произведение:", prod_ab))
print(paste("Деление:", div_ab))
print(paste("Степень:", power_ab))
print(paste("Остаток:", mod_ab))
print(paste("Целочисленное деление:", int_div_ab))

6.3.2 6.3.2 Операции над векторами

Здесь проявляется одна из главных сил R — векторизация. Когда вы применяете операцию к двум векторам, она выполняется поэлементно.

# Создадим два вектора
v1 <- c(10, 20, 30)
v2 <- c(1, 2, 3)

# Поэлементное сложение
v_sum <- v1 + v2 # (10+1, 20+2, 30+3)
print(v_sum)

# Поэлементное вычитание
v_diff <- v1 - v2
print(v_diff)

# Поэлементное умножение
v_prod <- v1 * v2
print(v_prod)

# Поэлементное деление
v_div <- v1 / v2
print(v_div)

6.3.2.1 Правило рециклирования (Recycling)

Что произойдет, если векторы разной длины? R использует правило рециклирования: короткий вектор “растягивается” до длины длинного вектора.

v_long <- c(10, 20, 30, 40, 50, 60)
v_short <- c(1, 2, 3)

# v_short будет "рециклирован" в c(1, 2, 3, 1, 2, 3)
result_recycle <- v_long + v_short
print(result_recycle)

Внимание: Если длина длинного вектора не является кратной длине короткого, R выполнит операцию, но выдаст предупреждение.

v_long2 <- c(10, 20, 30, 40, 50)
v_short2 <- c(1, 2)

# Длина 5 не кратна длине 2. Будет предупреждение.
result_warning <- v_long2 + v_short2
print(result_warning)

6.3.3 6.3.3 Операции над матрицами

С матрицами можно выполнять два типа операций: поэлементные и матричные.

6.3.3.1 Поэлементные операции

Если вы используете стандартные операторы (+, -, *, /), R будет выполнять операции поэлементно, так же как и с векторами. Для этого матрицы должны иметь одинаковый размер.

m1 <- matrix(1:4, nrow = 2, byrow = TRUE)
m2 <- matrix(5:8, nrow = 2, byrow = TRUE)

print("Матрица 1:")
print(m1)
print("Матрица 2:")
print(m2)

# Поэлементное умножение
elementwise_prod <- m1 * m2
print("Поэлементное произведение:")
print(elementwise_prod)

6.3.3.2 Матричные операции (линейная алгебра)

Для выполнения “настоящих” матричных операций из линейной алгебры используются специальные функции.

  • Матричное умножение: %*%
  • Транспонирование: t()
  • Обратная матрица: solve()
  • Определитель: det()
# Используем те же матрицы m1 и m2
# m1:
#      [,1] [,2]
# [1,]    1    2
# [2,]    3    4

# m2:
#      [,1] [,2]
# [1,]    5    6
# [2,]    7    8

# Матричное умножение
matrix_prod <- m1 %*% m2
print("Матричное произведение:")
print(matrix_prod)

6.3.4 6.3.4 Возведение в степень и логарифмирование

Для этих операций в R есть встроенные функции.

6.3.4.1 Возведение в степень

Как мы уже видели, для возведения в степень используется оператор ^ или **.

# Квадрат и куб числа
5^2
2^3

# Поэлементное возведение в степень для векторов
v <- c(1, 2, 3, 4)
v_squared <- v ^ 2
print(v_squared)

6.3.4.2 Логарифмирование

Для логарифмов есть функция log(), которая может принимать основание в качестве аргумента.

# Натуральный логарифм (основание e)
log(10)

# Десятичный логарифм (основание 10)
log10(100)

# Двоичный логарифм (основание 2)
log2(8)

# Логарифм с произвольным основанием (например, 5)
log(25, base = 5)

6.3.4.3 Другие полезные математические функции

  • sqrt(x): Квадратный корень из x.
  • abs(x): Модуль (абсолютное значение) x.
  • round(x, digits): Округление x до digits знаков после запятой.
  • ceiling(x): Округление вверх до ближайшего целого.
  • floor(x): Округление вниз до ближайшего целого.
  • Тригонометрические функции: sin(), cos(), tan() и т.д. (углы в радианах).
# Примеры других функций
sqrt(25)
abs(-15)
round(3.14159, 2)
ceiling(4.2)
floor(4.8)

6.4 6.4 Работа с пропущенными значениями (NA)

В реальных наборах данных почти всегда встречаются пропуски. Это может быть связано с ошибками при вводе данных, отказом респондента отвечать на вопрос или невозможностью провести измерение. В R для обозначения пропущенных значений используется специальное значение NA (Not Available).

Понимание того, как R обрабатывает NA, критически важно для правильного анализа.

6.4.1 Что такое NA?

NA — это не строка и не число, а специальный “логический” объект-заполнитель, который говорит: “здесь нет данных”. Ключевая особенность NA в том, что любая арифметическая или логическая операция с участием NA дает в результате NA.

# Арифметика с NA
10 + NA
100 / NA

# Сравнения с NA
10 > NA
NA == NA # Даже NA не "равно" самому себе!

Это логично: если мы не знаем одно из значений, мы не можем знать и результат операции.

6.4.2 Проверка на наличие NA

Поскольку NA == NA возвращает NA, для проверки на пропущенное значение используются специальные функции.

  • is.na(x): Возвращает TRUE для значений NA и FALSE для всех остальных. Это самый распространенный способ проверки.
  • anyNA(x): Быстро проверяет, есть ли хотя бы одно значение NA в объекте.
# Создадим вектор с пропущенными значениями
data_with_na <- c(10, 20, NA, 40, NA, 60)

# Проверяем каждый элемент
is.na(data_with_na)

# Проверяем, есть ли пропуски в векторе вообще
anyNA(data_with_na)

6.4.3 Типы NA

Хотя мы обычно видим просто NA, на самом деле это “семейство” объектов. Тип NA наследуется от типа вектора, в котором он находится. Это помогает R отслеживать, какие данные были пропущены.

  • NA — для логических и числовых векторов (по умолчанию).
  • NA_character_ — для символьных (текстовых) векторов.
  • NA_integer_ — для целочисленных векторов.
  • NA_complex_ — для комплексных чисел.

В большинстве случаев вам не нужно об этом беспокоиться, но это полезно знать при отладке.

# Создадим вектор с текстовыми данными
text_vector <- c("hello", "world", NA_character_, "R")

# Проверим тип каждого элемента
typeof(text_vector)
# [1] "character" "character" "character" "character"

# is.na() корректно находит специальный текстовый NA
is.na(text_vector)

6.4.4 Другие специальные значения

Кроме NA, в R есть и другие “особые” значения, которые важно не путать.

  • NaN (Not a Number): Возникает при математически неопределенных операциях (например, деление нуля на ноль).

    0 / 0

    NaN считается разновидностью NA, поэтому is.na(0/0) вернет TRUE. Для проверки именно на NaN есть функция is.nan().

  • Inf и -Inf (Infinity): Возникают при делении на ноль.

    1 / 0  # Infinity
-1 / 0 # -Infinity

6.4.5 Обработка пропущенных значений

Пропущенные значения могут мешать вычислениям. Большинство функций в R по умолчанию не могут обработать данные с NA и вернут NA.

# Функция mean() вернет NA, если в векторе есть пропуски
mean(data_with_na)

Чтобы это исправить, многие функции имеют аргумент na.rm = TRUE (NA remove), который говорит функции временно удалить все NA перед выполнением вычислений.

# Считаем среднее, удалив NA
mean(data_with_na, na.rm = TRUE)

Другие функции с аргументом na.rm: sum(), min(), max(), sd() (стандартное отклонение) и многие другие.

6.4.6 Удаление пропусков из данных

Часто возникает необходимость удалить строки, содержащие NA, из таблицы.

6.4.6.1 Базовый R: na.omit()

Функция na.omit() удаляет все строки, которые содержат хотя бы одно значение NA.

# Создадим фрейм данных с пропусками
df_na <- data.frame(
  id = 1:4,
  score = c(100, NA, 80, 95),
  group = c("A", "B", NA, "A")
)

print("Исходный фрейм данных:")
print(df_na)

# Удаляем строки с NA
df_clean <- na.omit(df_na)

print("Фрейм данных без строк с NA:")
print(df_clean)

6.4.6.2 Tidyverse: tidyr::drop_na()

Пакет tidyr предлагает более гибкую функцию drop_na().

# install.packages("tidyr")
library(tidyr)

# Удаляем строки с NA в любых столбцах (аналог na.omit)
df_clean_1 <- df_na %>% drop_na()

# Удаляем строки, где NA находится только в конкретном столбце 'score'
df_clean_2 <- df_na %>% drop_na(score)

print("Удалены строки с NA в столбце 'score':")
print(df_clean_2)

6.4.7 Замена пропусков (Imputation)

Удаление строк — не всегда лучший вариант, так как можно потерять много данных. Альтернатива — замена (импутация) пропусков, например, на среднее значение или медиану.

# Заменим пропуски в столбце 'score' на среднее значение по этому столбцу
mean_score <- mean(df_na$score, na.rm = TRUE)

# Используем ifelse() для замены: если значение NA, ставим среднее, иначе оставляем исходное
df_na$score_imputed <- ifelse(is.na(df_na$score), mean_score, df_na$score)

print("Фрейм данных с замененными пропусками в 'score':")
print(df_na)

6.4.8 Итог

  1. NA — это специальное значение для пропущенных данных.
  2. Любая операция с NA обычно дает NA.
  3. Используйте is.na() для проверки на пропуски.
  4. Используйте na.rm = TRUE в функциях, чтобы проигнорировать NA при вычислениях.
  5. Для удаления строк с NA используйте na.omit() (базовый R) или drop_na() (tidyr).
  6. Для замены пропусков используйте условную логику, например, ifelse().

6.5 6.5 Элементарные преобразования переменных с tidyverse

В предыдущих разделах мы уже познакомились с функцией mutate() из пакета dplyr для создания новых колонок. Здесь мы рассмотрим ее более подробно и в сочетании с case_when(), которая является “швейцарским ножом” для условной логики.

6.5.1 6.5.1 Мощь mutate()

Функция mutate() — это основной инструмент в dplyr для добавления или изменения колонок. Ее главная сила заключается в том, что она позволяет создавать новые переменные на основе существующих, а также использовать только что созданные переменные в том же вызове mutate().

Давайте создадим более сложный пример для демонстрации. Представим, что у нас есть данные о пациентах.

# Создадим фрейм данных с данными пациентов
patients <- data.frame(
  patient_id = 1:5,
  first_name = c("Иван", "Мария", "Петр", "Анна", "Сергей"),
  last_name = c("Иванов", "Петрова", "Сидоров", "Козлова", "Смирнов"),
  age = c(45, 34, 67, 29, 52),
  systolic_bp = c(130, 125, 155, 118, 142),
  diastolic_bp = c(85, 80, 95, 78, 90)
)

print("Исходные данные:")
print(patients)

Теперь используем mutate() для нескольких преобразований:

  1. Создать полное имя (full_name).
  2. Рассчитать пульсовое давление (pulse_pressure).
  3. Создать категорию возраста (age_category).
  4. Рассчитать среднее артериальное давление (map), используя пульсовое давление, созданное на шаге 2.
library(dplyr)

patients_transformed <- patients %>%
  mutate(
    # 1. Создаем полное имя, склеивая строки
    full_name = paste(first_name, last_name),
    
    # 2. Рассчитываем пульсовое давление
    pulse_pressure = systolic_bp - diastolic_bp,
    
    # 3. Создаем категорию возраста с помощью ifelse (для простого случая)
    age_category = ifelse(age >= 65, "Старший", "Младший"),
    
    # 4. Рассчитываем среднее артериальное давление, используя pulse_pressure
    # Формула: (диастолическое + пульсовое/3)
    mean_arterial_pressure = diastolic_bp + pulse_pressure / 3
  )

print("Преобразованные данные:")
print(patients_transformed)

Как видите, мы смогли создать несколько новых колонок в одном вызове, и mean_arterial_pressure была рассчитана с использованием переменной pulse_pressure, которая не существовала в исходном фрейме данных.

6.5.2 6.5.2 Условная логика с case_when()

Функция ifelse() отлично подходит для простых бинарных условий (если/иначе). Но что если у вас несколько условий? Вложенные ifelse(ifelse(...)) быстро становятся нечитаемыми.

Для таких случаев идеально подходит case_when(). Это обобщенная версия ifelse, которая работает как switch или if-elif-else в других языках.

Синтаксис case_when():

case_when(
  условие1 ~ результат1,
  условие2 ~ результат2,
  ...
  TRUE ~ результат_по_умолчанию # Важно! Что делать, если ни одно условие не выполнено
)
  • условие: Логическое выражение, которое возвращает TRUE или FALSE.
  • ~: Специальный оператор, который разделяет условие и результат.
  • результат: Значение, которое будет возвращено, если соответствующее условие истинно.
  • TRUE ~ ...: Строка по умолчанию, которая срабатывает, если все предыдущие условия были ложными. Это хорошая практика для обработки непредвиденных случаев.

6.5.2.1 Пример использования case_when()

Давайте перепишем создание категории возраста из предыдущего примера, сделав ее более детальной.

patients_detailed_age <- patients %>%
  mutate(
    age_category = case_when(
      age < 18 ~ "Ребенок",
      age >= 18 & age < 35 ~ "Взрослый (молодой)",
      age >= 35 & age < 65 ~ "Взрослый (средний возраст)",
      age >= 65 ~ "Пожилой",
      TRUE ~ "Неизвестная категория" # На всякий случай
    )
  )

print("Данные с детальной категорией возраста:")
print(patients_detailed_age)

case_when() делает код гораздо чище и понятнее, чем цепочка из ifelse().

6.5.2.2 Комбинированный пример: mutate() + case_when()

Теперь давайте создадим более сложную переменную — уровень риска для сердечно-сосудистых заболеваний — на основе нескольких условий.

patients_with_risk <- patients %>%
  mutate(
    risk_level = case_when(
      # Высокий риск: возраст > 60 И систолическое > 140
      age > 60 & systolic_bp > 140 ~ "Высокий",
      
      # Средний риск: возраст > 45 ИЛИ систолическое > 130
      age > 45 | systolic_bp > 130 ~ "Средний",
      
      # Низкий риск: во всех остальных случаях
      TRUE ~ "Низкий"
    )
  )

print("Данные с рассчитанным уровнем риска:")
print(patients_with_risk)

6.5.3 Итог: mutate() + case_when() vs. Базовый R

  • Базовый R: Для создания новых колонок используется df$new_col <- .... Для условной логики — вложенные ifelse(). r # Пример на базовом R patients$risk_level <- ifelse(patients$age > 60 & patients$systolic_bp > 140, "Высокий", ifelse(patients$age > 45 | patients$systolic_bp > 130, "Средний", "Низкий"))
  • tidyverse: Используется mutate() для преобразований и case_when() для сложной условной логики. r # Пример на tidyverse (см. код выше) patients %>% mutate(risk_level = case_when(...))

Подход tidyverse выигрывает за счет: * Читаемости: case_when намного нагляднее вложенных ifelse. * Модульности: Легко добавлять, удалять или изменять условия. * Интеграции: mutate() позволяет создавать переменные последовательно в одном блоке.

Отлично! Переходим к самой сути анализа данных. Создаем Главу 7.

7 Глава 7: Элементарный статистический анализ

После того как мы загрузили и очистили данные, наступает время для их изучения. Описательная статистика помогает нам суммировать и понять основные характеристики данных, их распределение и основные тенденции. В этом разделе мы рассмотрим, как рассчитывать базовые статистические показатели в R.

7.1 7.1 Базовые статистики для числовых данных

Для числовых векторов (или столбцов во фрейме данных) мы обычно хотим знать о центральной тенденции (типичное значение), разбросе (насколько данные разнообразны) и позиции (квантили).

7.1.1 7.1.1 Меры центральной тенденции

Эти показатели описывают “типичное” или “центральное” значение в наборе данных.

7.1.1.1 Среднее арифметическое (Mean)

Самый известный показатель центральной тенденции. Это сумма всех значений, деленная на их количество.

Формула для среднего арифметического набора данных \(x_1, x_2, ..., x_n\):

\[ \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i \]

В R для расчета используется функция mean().

# Создадим вектор с данными о росте в см
heights <- c(165, 170, 175, 180, 182, 158, 190)

# Рассчитаем средний рост
mean_height <- mean(heights)
print(paste("Средний рост:", round(mean_height, 2), "см"))

7.1.1.2 Медиана (Median)

Медиана — это значение, которое находится ровно посередине упорядоченного набора данных. Половина значений меньше медианы, а другая половина — больше.

Медиана менее чувствительна к выбросам (экстремально большим или малым значениям), чем среднее.

# Рассчитаем медиану
median_height <- median(heights)
print(paste("Медианный рост:", median_height, "см"))

7.1.2 7.1.2 Меры разброса (дисперсии)

Эти показатели описывают, насколько далеко значения в наборе данных расположены от их центрального значения.

7.1.2.1 Дисперсия (Variance)

Дисперсия измеряет средний квадрат отклонения каждого значения от среднего. Чем больше дисперсия, тем сильнее данные разбросаны.

Формула для выборочной дисперсии (\(s^2\)):

\[ s^2 = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^2 \]

В R для расчета используется функция var().

# Рассчитаем дисперсию
var_height <- var(heights)
print(paste("Дисперсия роста:", round(var_height, 2)))

7.1.2.2 Стандартное отклонение (Standard Deviation)

Стандартное отклонение — это корень квадратный из дисперсии. Оно выражено в тех же единицах, что и исходные данные, что делает его более интерпретируемым, чем дисперсию.

Формула для выборочного стандартного отклонения (\(s\)):

\[ s = \sqrt{s^2} = \sqrt{\frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^2} \]

В R для расчета используется функция sd().

# Рассчитаем стандартное отклонение
sd_height <- sd(heights)
print(paste("Стандартное отклонение роста:", round(sd_height, 2), "см"))

7.1.3 7.1.3 Меры позиции (Квантили)

Квантили делят упорядоченный набор данных на равные части.

7.1.3.1 Минимум и максимум

Самые простые позиционные меры.

# Рассчитаем минимум и максимум
min_height <- min(heights)
max_height <- max(heights)

print(paste("Минимальный рост:", min_height, "см"))
print(paste("Максимальный рост:", max_height, "см"))

7.1.3.2 Процентили и квартили (Percentiles and Quartiles)

Процентиль — это значение, ниже которого падает определенный процент всех наблюдений. * 25-й процентиль — это первый квартиль (\(Q_1\)). * 50-й процентиль — это медиана (\(Q_2\)). * 75-й процентиль — это третий квартиль (\(Q_3\)).

В R для расчета квантилей используется универсальная функция quantile(), где аргумент probs задает нужные процентили в виде дроби от 0 до 1.

# Рассчитаем ключевые квартили (25%, 50%, 75%)
quartiles <- quantile(heights, probs = c(0.25, 0.5, 0.75))
print(quartiles)

# Рассчитаем 90-й процентиль
p90 <- quantile(heights, probs = 0.9)
print(paste("90-й процентиль роста:", p90, "см"))

7.1.4 7.1.4 Функция summary() для быстрого обзора

Вместо того чтобы вызывать каждую функцию по отдельности, можно использовать функцию summary(). Она автоматически рассчитывает и выводит набор самых важных описательных статистик для числового вектора или для каждого числового столбца во фрейме данных.

# Получаем сводку по всем основным статистикам
summary_stats <- summary(heights)
print(summary_stats)

Как видите, summary() предоставляет минимальное значение, первый квартиль, медиану, среднее, третий квартиль и максимальное значение. Это самый быстрый способ получить первое представление о распределении ваших данных.

7.1.5 7.1.5 Оценка нормальности распределения

Многие статистические тесты (например, t-тест, ANOVA) предполагают, что данные следуют нормальному распределению (также известному как распределение Гаусса). Прежде чем применять такие тесты, важно проверить это предположение. Есть два основных подхода: визуальный и формальный.

7.1.5.1 Визуальная оценка

Самый простой и интуитивный способ — посмотреть на данные.

  1. Гистограмма (Histogram): Показывает форму распределения данных. Для нормального распределения гистограмма должна иметь симметричную колоколообразную форму.
  2. График квантиль-квантиль (Q-Q Plot): Сравнивает квантили ваших данных с теоретическими квантилями нормального распределения. Если данные распределены нормально, точки на графике должны лежать примерно вдоль прямой линии.

Для построения этих графиков мы будем использовать пакет ggplot2.

# install.packages("ggplot2")
library(ggplot2)

# Создадим два набора данных: один нормальный, один с отклонением
set.seed(42) # для воспроизводимости
normal_data <- rnorm(100, mean = 50, sd = 10)
skewed_data <- rexp(100, rate = 0.1) # Экспоненциальное распределение (скошенное)

# --- Визуализация для НОРМАЛЬНЫХ данных ---
ggplot(data.frame(value = normal_data), aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 15, fill = "skyblue", color = "black") +
  stat_function(fun = dnorm, args = list(mean = mean(normal_data), sd = sd(normal_data)), color = "red", size = 1) +
  labs(title = "Гистограмма нормальных данных", x = "Значение", y = "Плотность")

# Q-Q plot для нормальных данных
ggplot(data.frame(value = normal_data), aes(sample = value)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "Q-Q Plot для нормальных данных")

# --- Визуализация для СКОШЕННЫХ данных ---
ggplot(data.frame(value = skewed_data), aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 15, fill = "salmon", color = "black") +
  labs(title = "Гистограмма скошенных данных", x = "Значение", y = "Плотность")

# Q-Q plot для скошенных данных
ggplot(data.frame(value = skewed_data), aes(sample = value)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "Q-Q Plot для скошенных данных")

На графиках видно, что для normal_data точки на Q-Q графике лежат на прямой, а для skewed_data они сильно отклоняются от нее, особенно по краям.

7.1.5.2 Формальные тесты на нормальность

Эти тесты проверяют нулевую гипотезу (\(H_0\)) о том, что данные взяты из нормального распределения.

  • Если p-value < 0.05, мы отвергаем \(H_0\) и считаем, что данные не распределены нормально.
  • Если p-value >= 0.05, у нас нет оснований отвергать \(H_0\), и мы можем считать, что данные распределены нормально.

7.1.5.3 Тест Шапиро-Уилка (Shapiro-Wilk Test)

Один из самых популярных и мощных тестов на нормальность, хорошо работает на выборках до 5000 наблюдений. Он особенно чувствителен к отклонениям в хвостах распределения.

# Тест для нормальных данных
shapiro.test(normal_data)
# p-value = 0.5363 > 0.05. Нет оснований отвергать нормальность.

# Тест для скошенных данных
shapiro.test(skewed_data)
# p-value < 2.2e-16 < 0.05. Данные не распределены нормально.

7.1.5.4 Тест Д’Агостино и Пирсона (D’Agostino’s K-squared Test)

Этот тест проверяет нормальность на основе асимметрии (skewness) и эксцесса (kurtosis). Асимметрия мера скошенности распределения, а эксцесс — мера его “островершинности”. Тест Д’Агостино проверяет, совпадают ли эти показатели с теми, что ожидаются у нормального распределения.

Доступен в пакете fBasics.

# install.packages("fBasics")
library(fBasics)

# Тест для нормальных данных
dagoTest(normal_data)
# p-value = 0.6355 > 0.05. Нет оснований отвергать нормальность.

# Тест для скошенных данных
dagoTest(skewed_data)
# p-value < 2.2e-16 < 0.05. Данные не распределены нормально.

7.1.5.5 Тест Жарка-Бера (Jarque–Bera Test)

Подобно тесту Д’Агостино, этот тест также использует асимметрию и эксцесс для проверки нормальности. Он очень популярен в эконометрике. Доступен в пакете tseries.

# install.packages("tseries")
library(tseries)

# Тест для нормальных данных
jarque.bera.test(normal_data)
# p-value = 0.843 > 0.05. Нет оснований отвергать нормальность.

# Тест для скошенных данных
jarque.bera.test(skewed_data)
# p-value < 2.2e-16 < 0.05. Данные не распределены нормально.

7.1.5.6 Тест Андерсона-Дарлинга (Anderson–Darling Test)

Это модификация теста Колмогорова-Смирнова, которая придает больший вес значениям в хвостах распределения. Это делает его особенно чувствительным к отклонениям от нормальности на краях выборки. Доступен в пакете nortest.

# install.packages("nortest")
library(nortest)

# Тест для нормальных данных
ad.test(normal_data)
# p-value = 0.4231 > 0.05. Нет оснований отвергать нормальность.

# Тест для скошенных данных
ad.test(skewed_data)
# p-value < 2.2e-16 < 0.05. Данные не распределены нормально.

Важное замечание: На больших выборках (тысячи наблюдений) тесты на нормальность становятся слишком чувствительными и могут отвергать \(H_0\) из-за незначительных отклонений. В таких случаях визуальная оценка (Q-Q plot) часто является более надежным индикатором.

7.1.6 Какой тест выбрать?

  • shapiro.test(): Отличный выбор по умолчанию для небольших и средних выборок (n < 5000).
  • ad.test(): Используйте, если вас особенно беспокоят отклонения в хвостах распределения.
  • dagoTest() и jarque.bera.test(): Полезны, когда вы хотите понять, какой именно аспект (асимметрия или эксцесс) вызывает отклонение от нормальности.

На практике лучше использовать комбинацию методов: начните с визуализации (Q-Q plot), а затем подтвердите свои подозрения с помощью одного или двух формальных тестов.

7.1.7 7.1.6 Трансформация данных для достижения нормальности

Если ваши данные не распределены нормально, а для анализа это требование критично, можно попробовать применить математическую трансформацию. Цель трансформации — сделать распределение данных более симметричным и близким к нормальному.

7.1.7.1 Логарифмическая трансформация

Наиболее часто используемая трансформация, особенно для правосторонне скошенных (positively skewed) данных (длинный хвост справа). Она эффективна для данных, которые не могут быть отрицательными (например, концентрации, доходы, время).

# Применяем логарифмическую трансформацию к скошенным данным
# Используем log1p(x) = log(x + 1), чтобы избежать проблем с нулями
skewed_data <- rexp(100, rate = 0.5)

log_transformed_data <- log1p(skewed_data)

# Проверим результат визуально
ggplot(data.frame(value = log_transformed_data), aes(sample = value)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "Q-Q Plot после логарифмической трансформации")

# Проведем тест Шапиро-Уилка снова
shapiro.test(log_transformed_data)
# p-value теперь может быть > 0.05, что указывает на успех трансформации.

7.1.7.2 Трансформация Бокса-Кокса (Box-Cox)

Более мощный метод, который автоматически находит наилучшую степень трансформации из семейства степенных преобразований. Трансформация определяется параметром \(\lambda\):

\[ y(\lambda) = \begin{cases} \frac{x^\lambda - 1}{\lambda} & \text{if } \lambda \neq 0 \\ \log(x) & \text{if } \lambda = 0 \end{cases} \]

Для этого используется функция boxcox() из пакета MASS.

# install.packages("MASS")
library(MASS)

# Создадим модель, чтобы найти оптимальную лямбда
# Мы используем линейную модель с одним предиктором (просто константой)
# Это эквивалентно поиску лямбды для самого вектора
bc <- boxcox(lm(skewed_data ~ 1))

# Находим лямбду, которая максимизирует лог-правдоподобие
lambda <- bc$x[which.max(bc$y)]
print(paste("Оптимальная лямбда:", round(lambda, 2)))

# Применяем трансформацию с найденной лямбдой
boxcox_transformed_data <- (skewed_data^lambda - 1) / lambda

# Проверяем результат
ggplot(data.frame(value = boxcox_transformed_data), aes(sample = value)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "Q-Q Plot после трансформации Бокса-Кокса")

shapiro.test(boxcox_transformed_data)

7.1.8 Итог

  1. Всегда начинайте с визуализации (гистограмма, Q-Q plot), чтобы получить интуитивное понимание распределения.
  2. Используйте формальные тесты (например, shapiro.test()) для статистического подтверждения.
  3. Если данные не нормальны, попробуйте логарифмическую трансформацию или более мощную трансформацию Бокса-Кокса.
  4. После трансформации всегда проверяйте нормальность снова.

Конечно! Переходим от числовых данных к категориальным. Это не менее важный аспект анализа.

7.2 7.2 Базовые статистики для категориальных и бинарных данных

В отличие от числовых данных, для категориальных переменных (например, пол, группа лечения, диагноз) мы не можем рассчитать среднее или стандартное отклонение. Вместо этого мы интересуемся частотами и пропорциями каждой категории.

Бинарные данные — это частный случай категориальных данных, имеющих всего две категории (например, “да/нет”, “1/0”, “болен/здоров”).

7.2.1 7.2.1 Таблица частот

Самый первый и основной шаг анализа категориальной переменной — это подсчет того, сколько раз каждая категория встречается в наборе данных. Для этого используется функция table().

# Создадим вектор с категориальными данными
treatment_group <- factor(c("A", "B", "A", "A", "B", "C", "B", "A", "C"))

# Создадим бинарный вектор (1 = событие произошло, 0 = не произошло)
event_occurred <- c(1, 0, 1, 1, 0, 0, 1, 1, 0)

# --- Анализ категориальных данных ---

# Создаем таблицу частот для групп лечения
group_counts <- table(treatment_group)
print(group_counts)

# --- Анализ бинарных данных ---

# Создаем таблицу частот для событий
event_counts <- table(event_occurred)
print(event_counts)

Результатом функции table() является объект класса table, который показывает абсолютное количество наблюдений в каждой категории.

7.2.2 7.2.2 Пропорции и проценты

Часто абсолютные числа менее информативны, чем пропорции (доли) или проценты. Чтобы их рассчитать, нужно использовать функцию prop.table().

Функция prop.table() принимает на вход таблицу частот и возвращает таблицу пропорций (сумма всех элементов равна 1).

# Рассчитаем пропорции для групп лечения
group_proportions <- prop.table(group_counts)
print(group_proportions)

# Чтобы получить проценты, умножим на 100
group_percentages <- group_proportions * 100
print(group_percentages)

# То же самое для бинарных данных
event_proportions <- prop.table(event_counts)
print(event_proportions)

7.2.3 7.2.3 Функция summary() для категориальных данных

Как и в случае с числовыми данными, функция summary() предоставляет быстрый и удобный обзор. Для фактора она выведет таблицу частот.

# summary() для фактора
summary(treatment_group)

# summary() для бинарных данных, представленных как числа
summary(event_occurred)

Как видите, для фактора treatment_group summary() автоматически вызвала table(). Для числового вектора event_occurred она показала стандартные числовые статистики (квартили, медиану), что может быть не всегда то, что нам нужно. Это подчеркивает важность хранения бинарных данных в виде фактора с осмысленными метками.

7.2.4 7.2.4 Визуализация категориальных данных

Лучший способ представить частоты — это визуально.

7.2.4.1 Столбчатая диаграмма (Bar Plot)

Столбчатая диаграмма идеально подходит для отображения частот или пропорций категорий.

# install.packages("ggplot2")
library(ggplot2)

# Создадим фрейм данных для удобства построения графика
plot_data <- data.frame(
  group = names(group_counts),
  count = as.numeric(group_counts)
)

# Строим столбчатую диаграмму
ggplot(plot_data, aes(x = group, y = count, fill = group)) +
  geom_bar(stat = "identity") + # stat="identity" т.к. мы уже задали y=count
  labs(title = "Частота распределения по группам лечения",
       x = "Группа лечения",
       y = "Количество пациентов") +
  theme_minimal() +
  guides(fill = "none") # Убираем лишнюю легенду

7.2.4.2 Круговая диаграмма (Pie Chart)

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

# --- ШАГ 1: Создание исходных данных ---
# Этот код находится в том же фрагменте, что и график,
# поэтому переменные будут доступны для `ggplot()`.

# Создаем вектор с категориальными данными
treatment_group <- factor(c("A", "B", "A", "A", "B", "C", "B", "A", "C"))

# --- ШАГ 2: Расчет статистик ---

# Создаем таблицу частот
group_counts <- table(treatment_group)

# Рассчитываем проценты
group_percentages <- prop.table(group_counts) * 100

# --- ШАГ 3: Подготовка данных для ggplot2 ---

# Создаем фрейм данных, который будет использоваться для построения графика.
# Этот фрейм содержит две колонки: 'group' (для заливки) и 'percentage' (для высоты столбца).
pie_data <- data.frame(
  group = names(group_percentages),
  percentage = as.numeric(group_percentages) # Преобразуем в числовой формат
)

# --- ШАГ 4: Построение круговой диаграммы ---

# Теперь переменная 'pie_data' и ее колонки гарантированно существуют.
library(ggplot2)

ggplot(pie_data, aes(x = "", y = percentage, fill = group)) +
  geom_bar(stat = "identity", width = 1, color = "white") +
  coord_polar("y", start = 0) +
  labs(
    title = "Процентное распределение по группам лечения",
    x = NULL, # Убираем подпись оси X
    y = NULL  # Убираем подпись оси Y
  ) +
  theme_void() + # Убираем фон, сетку и оси
  theme(legend.position = "right") # Размещаем легенду справа

7.2.5 Итог

  1. Начинайте с table(), чтобы получить абсолютные частоты каждой категории.
  2. Используйте prop.table(), чтобы рассчитать пропорции и сравнивать категории относительно друг друга.
  3. Для быстрого обзора применяйте summary().
  4. Для визуализации используйте столбчатые диаграммы (bar plots) для точного сравнения и круговые диаграммы (pie charts) для демонстрации долей.

8 Глава 8: Интервальные оценки

9 Глава 9: Проверка статистических гипотез для одной выборки

10 Глава 10: Описательные статистики таблиц

11 Глава 11: Сравнительный анализ и проверка соответствующих гипотез

Сравнительный анализ — это ядро большинства биомедицинских исследований. Мы хотим сравнить эффект нового препарата с плацебо, оценить различия в показателях крови между группами пациентов или выяснить, изменился ли уровень маркера после лечения. В основе всех этих задач лежит проверка статистических гипотез.

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

Ключевые концепции:

  • Нулевая гипотеза (\(H_0\)): Это утверждение об отсутствии эффекта или различий. Например, “среднее артериальное давление в группе лечения и в группе плацебо не отличается”. Наша цель — найти доказательства, чтобы отвергнуть \(H_0\).
  • Альтернативная гипотеза (\(H_1\) или \(H_A\)): Это утверждение о наличии эффекта или различий. Например, “среднее артериальное давление в группе лечения отличается от такового в группе плацебо”. Это то, что мы хотим доказать.
  • p-value (p-значение): Это вероятность получить наблюдаемые (или еще более выраженные) результаты, если бы нулевая гипотеза была верна.
    • Маленькое p-value (обычно < 0.05) говорит о том, что такие данные были бы очень маловероятны при верной \(H_0\). Это дает нам основание отвергнуть \(H_0\) и считать, что различия статистически значимы.
    • Большое p-value (>= 0.05) говорит о том, что наши данные согласуются с \(H_0\). У нас нет оснований отвергать \(H_0\). Важно: это не означает, что \(H_0\) верна, а лишь то, что у нас не хватило доказательств, чтобы ее отвергнуть.
  • Уровень значимости (\(\alpha\)): Это порог, который мы устанавливаем заранее (чаще всего \(\alpha = 0.05\)). Если p-value < \(\alpha\), мы отвергаем \(H_0\).

Алгоритм выбора статистического теста

Выбор правильного теста зависит от нескольких факторов. Ваша схема прекрасно иллюстрирует этот процесс. Давайте пройдемся по ключевым вопросам, которые нужно себе задать.

Шаг 1: Какие у нас данные?

  1. Количественные (непрерывные) данные: Числовые данные, которые можно измерить (рост, вес, концентрация вещества в крови, артериальное давление).
  2. Категориальные данные: Данные, которые можно разделить на группы.
    • Бинарные: Две категории (пол: муж/жен, исход: выздоровел/не выздоровел).
    • Порядковые (ординальные): Категории с естественным порядком (стадия заболевания: легкая/средняя/тяжелая).
    • Номинальные: Категории без порядка (группа крови: A/B/O/AB).

Шаг 2: Сколько у нас групп для сравнения?

  • Одна выборка: Мы сравниваем показатель выборки с известным или теоретическим значением (например, средний рост пациентов в нашей клинике со средним по стране).
  • Две выборки: Мы сравниваем две группы между собой (группа лечения vs. группа плацебо).
  • Более двух выборок: Мы сравниваем три или более групп (доза препарата 10мг, 20мг, 30мг).

Шаг 3: Независимые или зависимые (парные) выборки?

  • Независимые выборки: Наблюдения в одной группе не связаны с наблюдениями в другой (например, разные пациенты в контрольной и экспериментальной группах).
  • Зависимые (парные) выборки: Наблюдения в одной группе связаны с наблюдениями в другой (например, измерение показателя у одних и тех же пациентов “до” и “после” лечения).

Шаг 4: Каковы параметры распределения (для количественных данных)?

  • Нормальное распределение: Данные в каждой группе распределены нормально. Мы проверяем это с помощью тестов на нормальность (Шапиро-Уилка, etc.) и визуально (Q-Q plot).
  • Не нормальное распределение: Данные не следуют нормальному распределению.

Итоговая схема выбора теста

Основываясь на ответах на эти вопросы, мы можем выбрать подходящий тест, следуя вашей логической схеме:

flowchart TD
    A["Выбор статистического теста"] --> B["Тип данных: Количественные"]
    A --> C["Тип данных: Категориальные"]

    B --> D["Распределение: Нормальное"]
    B --> E["Распределение: Не нормальное"]

    D --> F["Одна выборка"]
    D --> G["Две выборки"]
    D --> H["Более двух выборок"]

    F --> F1["t-тест Стьюдента"]
    G --> G1["t-тест Стьюдента (независ.)"]
    G --> G2["t-тест Стьюдента (завис.)"]
    H --> H1["ANOVA"]
    H --> H2["Дисперсионный анализ с повторными измерениями"]

    E --> I["Одна выборка"]
    E --> J["Две выборки"]
    E --> K["Более двух выборок"]

    I --> I1["Критерий хи-квадрат (χ²)"]
    I --> I2["Критерий Колмогорова-Смирнова"]
    J --> J1["U-критерий Манна-Уитни"]
    J --> J2["Критерий Вилкоксона"]
    K --> K1["H-критерий Краскела-Уоллиса"]
    K --> K2["Критерий Фридмана"]

    C --> L["Одна выборка"]
    C --> M["Две выборки"]
    C --> N["Более двух выборок"]

    L --> L1["Критерий хи-квадрат (χ²)"]
    L --> L2["Биномиальный тест"]
    M --> M1["Критерий хи-квадрат (χ²)"]
    M --> M2["Точный критерий Фишера"]
    N --> N1["Критерий хи-квадрат (χ²)"]

11.0.1 11.1 Логика сравнения выборок

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

В основе этого процесса лежит проверка статистических гипотез. Представим, что мы хотим сравнить среднее артериальное давление в двух группах пациентов: одной давали новый препарат, другой — плацебо. Мы видим, что средние значения немного отличаются. Наша задача — выяснить, является ли это различие статистически значимым.

Процесс проверки гипотез можно свести к следующему алгоритму (представлен на слайде):

  1. Формулировка гипотез. Мы всегда формулируем две гипотезы:

    • Нулевая гипотеза (\(H_0\)): Это гипотеза об отсутствии различий. В нашем примере: “Среднее артериальное давление в группе лечения не отличается от среднего в группе плацебо”. Это консервативное предположение, которое мы пытаемся опровергнуть.
    • Альтернативная гипотеза (\(H_1\)): Это гипотеза о наличии различий. В нашем примере: “Среднее артериальное давление в группе лечения отличается от среднего в группе плацебо”. Это то, что мы хотим доказать.

  2. Выбор статистического теста. Мы не можем просто так сравнить средние значения. Нам нужен статистический инструмент, который поможет оценить значимость различий. Выбор теста — это ключевой шаг, который зависит от трех основных факторов, как показано на схеме:

    • Тип данных: Количественные (например, давление, рост) или категориальные (например, пол, исход лечения).
    • Количество выборок: Одна, две или более групп для сравнения.
    • Параметры распределения: Следуют ли данные нормальному закону распределения или нет. От этого зависит, будем ли мы использовать параметрические тесты (например, t-тест) или их непараметрические аналоги (например, U-критерий Манна-Уитни).

  3. Вычисление p-value. После выбора и проведения теста мы получаем на выходе p-value (p-значение). Это и есть главный результат нашего анализа. p-value — это вероятность получить наблюдаемые (или еще более сильные) различия между группами, если бы на самом деле верна нулевая гипотеза (то есть, если бы реальных различий не было).

  4. Принятие статистического решения. Мы сравниваем полученный p-value с заранее заданным порогом — уровнем значимости (\(\alpha\)). В большинстве биомедицинских исследований \(\alpha = 0.05\).

    • Если p-value < 0.05: Мы считаем, что получить такие различия случайно было бы очень маловероятно. Поэтому мы отвергаем нулевую гипотезу (\(H_0\)) и делаем вывод, что различия между группами статистически значимы.
    • Если p-value ≥ 0.05: Мы считаем, что такие различия вполне могли возникнуть случайно. У нас нет оснований отвергать нулевую гипотезу. Важно: это не значит, что \(H_0\) верна, а лишь то, что у нас не хватило доказательств, чтобы ее опровергнуть.

Таким образом, вся логика сравнения выборок — это структурированный процесс, который позволяет нам объективно оценить, являются ли наблюдаемые в исследовании эффекты реальными, а не случайными.

Отлично! Это очень важный практический аспект, который делает применение t-критерия более осознанным. Давайте интегрируем эту логику в подраздел.

11.0.2 11.2 Сравнение количественных показателей в 2 независимых выборках

11.0.3 11.2.1 Использование t-критерия Стьюдента

t-критерий Стьюдента — это параметрический тест, используемый для сравнения средних значений двух независимых групп. Основная идея теста — определить, является ли наблюдаемая разница между средними значениями двух выборок статистически значимой, или она могла возникнуть случайно.

Процесс сравнения включает вычисление t-статистики, а затем нахождение соответствующего p-value на основе t-распределения.

11.0.3.1 Оценка равенства дисперсий

Выбор между классическим t-критерием Стьюдента и критерием Уэлча зависит от того, можно ли считать дисперсии в двух группах равными. Хотя формальный тест на равенство дисперсий (например, F-тест) существует, на практике часто используют более простое и надежное эмпирическое правило: сравнение отношения дисперсий.

Пусть \(D(x)\) и \(D(y)\) — это выборочные дисперсии двух групп.

  1. Сходные дисперсии: Если отношение дисперсий находится в диапазоне от 0.5 до 2, то есть \(0.5 < \frac{D(x)}{D(y)} < 2\), можно считать, что дисперсии не сильно различаются. В этом случае можно использовать классический t-критерий Стьюдента.
  2. Значительно разные дисперсии: Если отношение дисперсий выходит за пределы этого диапазона (\(\frac{D(x)}{D(y)} \le 0.5\) или \(\frac{D(x)}{D(y)} \ge 2\)), дисперсии считаются существенно различными. В этом случае необходимо использовать поправку Уэлча.

Важное замечание: Критерий Уэлча является более надежным и универсальным, так как он хорошо работает как при равных, так и при неравных дисперсиях. Поэтому он является методом по умолчанию в R и часто является предпочтительным выбором.

11.0.3.2 Формулы для вычислений

В основе t-критерия лежит сравнение разницы средних (\(\bar{x} - \bar{y}\)) с стандартной ошибкой разницы средних.

\[ t = \frac{\bar{x} - \bar{y}}{SE_{\bar{x} - \bar{y}}} \]

1. Предположение о равенстве дисперсий (критерий Стьюдента)

Если мы считаем, что дисперсии двух генеральных совокупностей равны (\(\sigma_x^2 = \sigma_y^2\)), мы используем объединенное (взвешенное) стандартное отклонение (\(s_p\)).

Формула для объединенного стандартного отклонения:

\[ s_p = \sqrt{\frac{(n_x - 1)s_x^2 + (n_y - 1)s_y^2}{n_x + n_y - 2}} \]

Формула для стандартной ошибки разницы средних:

\[ SE_{\bar{x} - \bar{y}} = s_p \sqrt{\frac{1}{n_x} + \frac{1}{n_y}} \]

Степени свободы для этого случая:

\[ df = n_x + n_y - 2 \]

2. Предположение о неравенстве дисперсий (критерий Уэлча)

Если дисперсии существенно различаются, мы используем коррекцию Уэлча.

Формула для стандартной ошибки разницы средних по Уэлчу:

\[ SE_{Welch} = \sqrt{\frac{s_x^2}{n_x} + \frac{s_y^2}{n_y}} \]

Формула для степеней свободы по Уэлчу (более сложная):

\[ df = \frac{\left( \frac{s_x^2}{n_x} + \frac{s_y^2}{n_y} \right)^2}{\frac{(s_x^2/n_x)^2}{n_x-1} + \frac{(s_y^2/n_y)^2}{n_y-1}} \]

3. Расчет 95% Доверительного интервала (ДИ)

95% ДИ для разницы средних показывает диапазон, в котором с 95% уверенностью находится истинная разница между средними в генеральных совокупностях.

Формула для 95% ДИ:

\[ (\bar{x} - \bar{y}) \pm t_{(\alpha/2, df)} \times SE_{\bar{x} - \bar{y}} \]

где: * \((\bar{x} - \bar{y})\) — разница выборочных средних. * \(t_{(\alpha/2, df)}\) — критическое значение t-распределения для двустороннего теста с \(df\) степенями свободы и уровнем значимости \(\alpha=0.05\). * \(SE_{\bar{x} - \bar{y}}\) — стандартная ошибка разницы средних.

11.0.3.3 Расчет в R и вручную

Давайте сравним две выборки x и y с помощью t-критерия.

# Данные из примера
x <- c(4.33, 6.13, 5.40, 7.41, 4.83, 6.51, 8.26, 2.93, 1.15, 5.14)
y <- c(15.79, 9.80, 4.74, 2.84, 1.17, 1.67, -0.82, 12.78, 11.16, 5.86)

# --- Шаг 0: Оценим равенство дисперсий ---
var_x <- var(x)
var_y <- var(y)
variance_ratio <- var_x / var_y

cat("Дисперсия X:", var_x, "\n")
cat("Дисперсия Y:", var_y, "\n")
cat("Отношение дисперсий (D(x)/D(y)):", variance_ratio, "\n")

# Отношение дисперсий 1.41, что находится в диапазоне (0.5, 2).
# Можно предположить, что дисперсии сходны.
# Тем не менее, для надежности мы продолжим с критерием Уэлча (по умолчанию в R).

1. Расчет с помощью R (рекомендуемый способ)

Функция t.test() в R выполняет все вычисления автоматически. По умолчанию она использует критерий Уэлча (var.equal = FALSE).

# Проводим t-тест для двух выборок
t_test_result <- t.test(x, y)

# Выводим результат
print(t_test_result)

Вывод функции t.test() содержит всю необходимую информацию: * t: вычисленное значение t-статистики. * df: количество степеней свободы. * p-value: p-value, соответствующее данной t-статистике и df. * conf.int: 95% доверительный интервал для разницы средних. * mean of x, mean of y: средние значения для каждой выборки.

2. Расчет вручную (для понимания процесса)

Давайте воспроизведем расчеты, чтобы понять, как R получает эти значения. Мы будем использовать формулы для критерия Уэлча, так как он применяется по умолчанию.

# --- Шаг 1: Вычисляем базовые статистики для каждой выборки ---
n_x <- length(x)
n_y <- length(y)
mean_x <- mean(x)
mean_y <- mean(y)

# --- Шаг 2: Вычисляем t-статистику по Уэлчу ---
se_diff <- sqrt(var_x / n_x + var_y / n_y) # SE_{Welch}
t_stat <- (mean_x - mean_y) / se_diff

# --- Шаг 3: Вычисляем степени свободы по Уэлчу ---
numerator <- (var_x / n_x + var_y / n_y)^2
denominator <- (var_x / n_x)^2 / (n_x - 1) + (var_y / n_y)^2 / (n_y - 1)
df_welch <- numerator / denominator

# --- Шаг 4: Вычисляем p-value ---
# Для двустороннего теста умножаем на 2
p_value_manual <- 2 * pt(abs(t_stat), df = df_welch, lower.tail = FALSE)

# --- Шаг 5: Вычисляем 95% Доверительный интервал ---
# Находим критическое значение t для 95% ДИ (квантиль 0.975)
t_critical <- qt(0.975, df = df_welch)

# Рассчитываем разницу средних и границы интервала
mean_diff <- mean_x - mean_y
ci_lower <- mean_diff - t_critical * se_diff
ci_upper <- mean_diff + t_critical * se_diff

# --- Шаг 6: Сравниваем результаты ---
cat("\nРезультаты ручного расчета:\n")
cat("Разница средних:", mean_diff, "\n")
cat("t-статистика:", t_stat, "\n")
cat("Степени свободы:", df_welch, "\n")
cat("p-value:", p_value_manual, "\n")
cat("95% ДИ:", ci_lower, "to", ci_upper, "\n")

Теперь результаты ручного расчета (t-статистика, степени свободы, p-value и 95% ДИ) полностью совпадают с результатами, полученными с помощью функции t.test(). Это подтверждает, что мы правильно понимаем механику работы теста. На практике всегда следует использовать встроенную функцию t.test(), так как она менее подвержена ошибкам и вычисляет все необходимые метрики автоматически.

Конечно! Рассказываем про U-критерий Манна-Уитни — основной непараметрический аналог t-критерия.

11.0.4 11.2.2 Использование U-критерия Манна-Уитни

U-критерий Манна-Уитни (также известный как критерий суммы рангов Уилкоксона) — это непараметрический тест, который используется для сравнения двух независимых выборок. Он является альтернативой t-критерию Стьюдента.

11.0.4.1 Основная идея и когда его использовать

В отличие от t-критерия, который сравнивает средние значения, U-критерий Манна-Уитни сравнивает распределения двух выборок. Его основная гипотеза (\(H_0\)) заключается в том, что обе выборки взяты из одного и того же распределения.

U-критерий используется, когда: 1. Данные не распределены нормально. Это основное условие для его применения. 2. Данные являются порядковыми (ординальными), когда нельзя рассчитать среднее, но можно упорядочить значения (например, оценки “плохо”, “удовлетворительно”, “хорошо”).

11.0.4.2 Логика вычислений

Логика теста основана на рангах, а не на самих значениях.

  1. Объединение и ранжирование: Данные из обеих выборок объединяются в один общий список. Затем всем значениям присваиваются ранги от 1 (самое маленькое значение) до N (самое большое значение, где N — общее число наблюдений). Если есть одинаковые значения (связки), им присваивается средний ранг.
  2. Суммирование рангов: Ранги, принадлежащие каждой из исходных выборок, суммируются отдельно.
  3. Вычисление U-статистики: U-статистика показывает, сколько раз значение из одной выборки предшествует значению из другой выборки в общем ранжированном списке. Фактически, U-статистика для одной выборки (\(U_x\)) вычисляется как:

\[ U_x = R_x - \frac{n_x(n_x + 1)}{2} \]

где: * \(R_x\) — сумма рангов для выборки \(x\). * \(n_x\) — количество наблюдений в выборке \(x\).

U-статистика для второй выборки (\(U_y\)) вычисляется аналогично. В отчетах обычно используется меньшая из двух U-статистик.

  1. Определение p-value: Полученное значение U-статистики сравнивается с критическими значениями из таблиц распределения Манна-Уитни (для малых выборок) или аппроксимируется нормальным распределением (для больших выборок, \(n > 20\)), чтобы найти p-value.

11.0.4.3 Интерпретация результата

  • Маленькое p-value (p < 0.05): Мы отвергаем нулевую гипотезу. Это означает, что распределения двух выборок статистически значимо различаются. Часто это интерпретируют как различие в медианах, хотя формально тест проверяет именно равенство распределений (стохастическое доминирование).
  • Большое p-value (p ≥ 0.05): У нас нет оснований отвергать нулевую гипотезу. Мы не можем утверждать, что распределения выборок различаются.

11.0.4.4 Расчет в R

В R для проведения U-критерия Манна-Уитни используется функция wilcox.test(). Она работает так же, как и t.test().

Давайте используем те же данные, но представим, что они не распределены нормально.

# Данные из примера
x <- c(4.33, 6.13, 5.40, 7.41, 4.83, 6.51, 8.26, 2.93, 1.15, 5.14)
y <- c(15.79, 9.80, 4.74, 2.84, 1.17, 1.67, -0.82, 12.78, 11.16, 5.86)

# Проводим U-критерий Манна-Уитни
wilcox_result <- wilcox.test(x, y)

# Выводим результат
print(wilcox_result)

Вывод функции wilcox.test() содержит: * W: значение U-статистики (в R она обозначается как W). * p-value: p-value, соответствующее этой статистике. * alternative: тип альтернативной гипотезы (по умолчанию two.sided).

11.0.4.5 Сравнение с t-критерием

Характеристика t-критерий Стьюдента U-критерий Манна-Уитни
Тип данных Количественные, интервальные Количественные, порядковые
Предположение Нормальное распределение Никаких предположений о распределении
Что сравнивает Средние значения Ранги (т.е. распределения)
Мощность Выше, если данные нормальны Ниже, если данные нормальны
Мощность Ниже, если данные не нормальны Выше, если данные не нормальны

Итог: U-критерий Манна-Уитни — это мощный и надежный инструмент для сравнения двух групп, когда предположение о нормальности распределения данных не выполняется.

11.0.5 11.2.3 Односторонние и двусторонние тесты

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

11.0.5.1 Двусторонний тест (Two-tailed test)

Это наиболее распространенный и консервативный тип теста. Он используется, когда у нас нет предварительных предположений о том, какая группа покажет более высокий результат.

  • Нулевая гипотеза (\(H_0\)): Средние значения (или распределения) двух групп равны. \(H_0: \mu_x = \mu_y\)
  • Альтернативная гипотеза (\(H_1\)): Средние значения (или распределения) двух групп не равны. Различие может быть в любую сторону. \(H_1: \mu_x \neq \mu_y\)

Когда использовать: В большинстве научных исследований, особенно на начальных этапах. Например, мы сравниваем новый препарат с плацебо и хотим узнать, есть ли какая-либо разница в эффективности, не предполагая заранее, лучше он или хуже.

Как это работает: p-value рассчитывается как вероятность получить наблюдаемую (или еще более экстремальную) разницу в любом направлении. То есть мы суммируем площади под кривой распределения в обоих “хвостах” (слева и справа от наблюдаемого значения).

В R это является поведением по умолчанию.

# Двусторонний t-тест (по умолчанию)
t_test_two_sided <- t.test(x, y)
print(t_test_two_sided)

11.0.5.2 Односторонний тест (One-tailed test)

Этот тест используется, когда у нас есть четкое априорное предположение о направлении эффекта. Мы проверяем, не является ли одна группа специфически больше (или специфически меньше) другой.

  • Нулевая гипотеза (\(H_0\)): Среднее значение группы X не больше (или не меньше) среднего значения группы Y.
    • Для проверки “X > Y”: \(H_0: \mu_x \le \mu_y\)
    • Для проверки “X < Y”: \(H_0: \mu_x \ge \mu_y\)
  • Альтернативная гипотеза (\(H_1\)): Среднее значение группы X больше (или меньше) среднего значения группы Y.
    • Для проверки “X > Y”: \(H_1: \mu_x > \mu_y\)
    • Для проверки “X < Y”: \(H_1: \mu_x < \mu_y\)

Когда использовать: * Когда у вас есть сильные теоретические или предыдущие экспериментальные данные, предсказывающие направление эффекта. * В клинических испытаниях, где вы хотите доказать, что новый препарат не хуже стандартного (non-inferiority trial).

Как это работает: p-value рассчитывается как вероятность получить наблюдаемую (или еще более экстремальную) разницу только в одном направлении. Мы смотрим на площадь под кривой распределения только в одном “хвосте”.

В R это задается с помощью аргумента alternative.

# Односторонний тест: проверяем гипотезу, что среднее X МЕНЬШЕ среднего Y
# (т.е. Y больше X)
t_test_one_sided <- t.test(x, y, alternative = "less")
print(t_test_one_sided)

# Односторонний тест: проверяем гипотезу, что среднее X БОЛЬШЕ среднего Y
# t.test(x, y, alternative = "greater")

11.0.5.3 Сравнение p-value

Обратите внимание на p-value в результатах.

p-value для одностороннего теста ровно в два раза меньше, чем для двустороннего. Это логично, так как мы “ловим” экстремальные значения только в одном хвосте распределения, а не в двух.

11.0.5.4 Применение к U-критерию Манна-Уитни

Те же самые принципы применимы и к непараметрическому U-критерию Манна-Уитни. Функция wilcox.test() также имеет аргумент alternative.

  • alternative = "two.sided" (по умолчанию): \(H_1\): распределения не равны.
  • alternative = "less": \(H_1\): распределение X сдвинуто в сторону меньших значений, чем Y (значения из X систематически меньше, чем из Y).
  • alternative = "greater": \(H_1\): распределение X сдвинуто в сторону больших значений, чем Y.
# Односторонний U-критерий: проверяем, что значения X систематически МЕНЬШЕ значений Y
wilcox_test_one_sided <- wilcox.test(x, y, alternative = "less")
print(wilcox_test_one_sided)

11.0.6 Предостережение и рекомендации

  1. Не “подгоняйте” результат. Решение о использовании одностороннего или двустороннего теста должно быть принято до просмотра данных, на этапе планирования исследования. Выбор одностороннего теста после того, как вы увидели, что p-value для двустороннего чуть больше 0.05, является плохой научной практикой (p-hacking).
  2. Обосновывайте свой выбор. В отчете или статье всегда указывайте, почему вы использовали односторонний тест, и четко формулируйте вашу альтернативную гипотезу.
  3. Когда сомневаетесь, используйте двусторонний тест. Он является более строгим и общепринятым стандартом в большинстве областей.

11.0.7 11.3 Сравнение 2-х независимых выборок с бинарными исходами

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

11.0.8 11.3.1 Сравнительный анализ для пропорций (Z-тест) и вычисление 95% ДИ

Представим, что мы проводим клиническое исследование нового препарата. У нас есть две группы: лечение (treatment) и плацебо (placebo). Мы хотим сравнить долю пациентов, у которых наступило улучшение.

11.0.8.1 Основные метрики

Для каждой группы мы рассчитываем: * \(n\) — общее количество пациентов в группе. * \(k\) — количество пациентов с положительным исходом (улучшением). * \(\hat{p} = \frac{k}{n}\) — выборочная пропорция (доля) успехов. * \(\hat{q} = 1 - \hat{p}\) — доля неудач.

Пример данных: * Группа лечения (treatment): \(n_1 = 100\), \(k_1 = 45\) -> \(\hat{p}_1 = 0.45\) * Группа плацебо (placebo): \(n_2 = 105\), \(k_2 = 30\) -> \(\hat{p}_2 = 0.286\)

11.0.8.2 Статистика Z и проверка гипотезы

Для сравнения двух пропорций используется Z-критерий для разницы пропорций. Он основан на нормальном распределении.

Нулевая гипотеза (\(H_0\)): Пропорции успехов в двух группах равны (\(p_1 = p_2\)). Альтернативная гипотеза (\(H_1\)): Пропорции успехов в двух группах не равны (\(p_1 \neq p_2\)).

1. Расчет объединенной пропорции

В соответствии с нулевой гипотезой, мы предполагаем, что обе группы происходят из одной генеральной совокупности с общей пропорцией успехов. Мы оцениваем эту общую пропорцию (\(\hat{p}_{pooled}\)), объединяя данные из обеих групп:

\[ \hat{p}_{pooled} = \frac{k_1 + k_2}{n_1 + n_2} \]

2. Расчет стандартной ошибки разницы пропорций

Стандартная ошибка разницы пропорций (\(SE_{\hat{p}_1 - \hat{p}_2}\)) вычисляется с использованием этой объединенной пропорции:

\[ SE_{\hat{p}_1 - \hat{p}_2} = \sqrt{ \hat{p}_{pooled} \hat{q}_{pooled} \left( \frac{1}{n_1} + \frac{1}{n_2} \right) } \] где \(\hat{q}_{pooled} = 1 - \hat{p}_{pooled}\).

3. Расчет Z-статистики

Z-статистика показывает, насколько наблюдаемая разница пропорций (\(\hat{p}_1 - \hat{p}_2\)) отклоняется от нуля (разницы, ожидаемой при \(H_0\)), измеренная в единицах стандартной ошибки.

\[ Z = \frac{\hat{p}_1 - \hat{p}_2}{SE_{\hat{p}_1 - \hat{p}_2}} \]

Полученное значение Z сравнивается со стандартным нормальным распределением для нахождения p-value.

11.0.8.3 95% Доверительный интервал для разницы пропорций

Для оценки величины эффекта мы строим 95% доверительный интервал для разницы пропорций (\(p_1 - p_2\)). Здесь есть важный нюанс: для расчета ДИ мы не используем объединенную пропорцию, так как мы оцениваем разницу, а не проверяем гипотезу о равенстве. Вместо этого мы используем отдельные пропорции каждой группы.

Формула для стандартной ошибки разницы пропорций для ДИ (метод Вальда):

\[ SE_{diff} = \sqrt{ \frac{\hat{p}_1 \hat{q}_1}{n_1} + \frac{\hat{p}_2 \hat{q}_2}{n_2} } \]

Формула для 95% ДИ:

\[ (\hat{p}_1 - \hat{p}_2) \pm Z_{(\alpha/2)} \times SE_{diff} \]

где \(Z_{(\alpha/2)}\) — критическое значение стандартного нормального распределения (для 95% ДИ это 1.96).

11.0.8.4 Расчет в R

В R для сравнения пропорций используется функция prop.test().

# --- Данные для анализа ---
# Количество успехов в каждой группе
successes <- c(45, 30)
# Общее количество наблюдений в каждой группе
totals <- c(100, 105)

# --- Расчет с помощью R ---
# prop.test() по умолчанию выполняет тест с коррекцией непрерывности Йейтса.
# Чтобы получить "чистый" Z-тест, отключаем ее: correct = FALSE
prop_test_result <- prop.test(successes, totals, correct = FALSE)

# Выводим результат
print(prop_test_result)

Анализ вывода prop.test(): * X-squared: Значение статистики хи-квадрат. Обратите внимание, что для двух пропорций Z-статистика и \(\chi^2\) связаны: \(\chi^2 = Z^2\). * df: Степени свободы (для двух групп всегда 1). * p-value: p-value теста. Если p < 0.05, мы отвергаем \(H_0\). * estimate of probability: Выборочные пропорции для каждой группы (\(\hat{p}_1\) и \(\hat{p}_2\)). * 95 percent confidence interval: 95% доверительный интервал для разницы пропорций.

11.0.8.5 Расчет вручную (для понимания)

Давайте воспроизведем расчеты для нашего примера.

# --- Шаг 1: Основные метрики ---
n1 <- 100; n2 <- 105
k1 <- 45;  k2 <- 30
p1_hat <- k1 / n1
p2_hat <- k2 / n2

# --- Шаг 2: Расчет Z-статистики ---
p_pooled <- (k1 + k2) / (n1 + n2)
q_pooled <- 1 - p_pooled
se_diff <- sqrt(p_pooled * q_pooled * (1/n1 + 1/n2))
z_stat <- (p1_hat - p2_hat) / se_diff

# --- Шаг 3: Расчет p-value ---
# Для двустороннего теста
p_value_manual <- 2 * pnorm(abs(z_stat), lower.tail = FALSE)

# --- Шаг 4: Расчет 95% ДИ ---
se_ci <- sqrt((p1_hat * (1-p1_hat) / n1) + (p2_hat * (1-p2_hat) / n2))
diff_prop <- p1_hat - p2_hat
ci_lower <- diff_prop - 1.96 * se_ci
ci_upper <- diff_prop + 1.96 * se_ci

# --- Сравнение результатов ---
cat("Результаты ручного расчета:\n")
cat("Разница пропорций:", diff_prop, "\n")
cat("Z-статистика:", z_stat, "\n")
cat("p-value:", p_value_manual, "\n")
cat("95% ДИ для разницы пропорций:", ci_lower, "to", ci_upper, "\n")

Как видите, результаты ручного расчета (p-value и 95% ДИ) практически совпадают с результатами функции prop.test(). Использование prop.test() значительно проще и является стандартной практикой.

11.0.9 11.3.2 Таблицы сопряженности, отношение шансов и относительный риск

Когда мы сравниваем два бинарных исхода (например, “улучшение есть/нет” в группах лечения и плацебо), самый удобный способ представить данные — это таблица сопряженности 2x2.

11.0.9.1 Таблица сопряженности 2x2

Давайте представим наши данные в виде таблицы, где строки — это группы (лечение/плацебо), а столбцы — исходы (улучшение/отсутствие улучшения).

Улучшение (+) Нет улучшения (-) Всего
Лечение a b a + b
Плацебо c d c + d
Всего a + c b + d N

Для нашего примера (45/100 в лечении и 30/105 в плацебо) таблица выглядит так:

Улучшение (+) Нет улучшения (-) Всего
Лечение 45 55 100
Плацебо 30 75 105
Всего 75 130 205

11.0.9.2 Меры эффекта: Отношение шансов и Относительный риск

Чтобы понять, насколько эффективно лечение, мы рассчитываем меры эффекта. Две самые распространенные — это Отношение Шансов (ОШ) и Относительный Риск (ОР).

11.0.9.2.1 1. Отношение Шансов (Odds Ratio, OR)

Шанс (Odds) события — это отношение вероятности того, что событие произойдет, к вероятности того, что оно не произойдет. * Шанс улучшения в группе лечения: \(Odds_1 = \frac{a}{b}\) * Шанс улучшения в группе плацебо: \(Odds_2 = \frac{c}{d}\)

Отношение шансов (ОШ) — это отношение шансов в двух группах.

\[ OR = \frac{Odds_1}{Odds_2} = \frac{a/b}{c/d} = \frac{a \cdot d}{b \cdot c} \]

Интерпретация: * OR = 1: Шансы события в обеих группах одинаковы (эффекта нет). * OR > 1: Шансы события в группе лечения выше, чем в группе плацебо (позитивный эффект). * OR < 1: Шансы события в группе лечения ниже, чем в группе плацебо (негативный эффект).

11.0.9.2.2 2. Относительный Риск (Relative Risk, RR)

Риск (Risk) события — это просто его вероятность (доля). * Риск улучшения в группе лечения: \(Risk_1 = \frac{a}{a+b}\) * Риск улучшения в группе плацебо: \(Risk_2 = \frac{c}{c+d}\)

Относительный риск (ОР) — это отношение рисков в двух группах.

\[ RR = \frac{Risk_1}{Risk_2} = \frac{a/(a+b)}{c/(c+d)} \]

Интерпретация: * RR = 1: Риск события в обеих группах одинаков (эффекта нет). * RR > 1: Риск события в группе лечения выше, чем в группе плацебо (позитивный эффект). * RR < 1: Риск события в группе лечения ниже, чем в группе плацебо (негативный эффект).

11.0.9.3 Расчет стандартных ошибок и 95% ДИ

Поскольку OR и RR рассчитываются на основе выборочных данных, мы должны оценить их статистическую неопределенность с помощью стандартной ошибки (SE) и 95% доверительного интервала (ДИ).

Ключевой момент: Распределение OR и RR асимметрично. Поэтому ДИ для них рассчитывается в логарифмической шкале (log-scale), а затем результат обратно преобразуется.

11.0.9.3.1 1. Для Отношения Шансов (OR)
  • Стандартная ошибка логарифма OR: \[ SE(\ln(OR)) = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \]
  • 95% ДИ для OR: \[ 95\% CI_{OR} = \left( e^{\ln(OR) - 1.96 \times SE(\ln(OR))}, e^{\ln(OR) + 1.96 \times SE(\ln(OR))} \right) \]
11.0.9.3.2 2. Для Относительного Риска (RR)
  • Стандартная ошибка логарифма RR: \[ SE(\ln(RR)) = \sqrt{ \frac{1-a}{a \cdot (a+b)} + \frac{1-c}{c \cdot (c+d)} } \]
  • 95% ДИ для RR: \[ 95\% CI_{RR} = \left( e^{\ln(RR) - 1.96 \times SE(\ln(RR))}, e^{\ln(RR) + 1.96 \times SE(\ln(RR))} \right) \]

11.0.9.4 Оценка статистической значимости (p-value)

После расчета OR или RR и их 95% ДИ нам нужно определить, является ли полученный эффект статистически значимым. Есть три основных способа это сделать.

Способ 1: Через доверительный интервал (самый простой и интуитивный)

Это самый распространенный и рекомендуемый способ. Правило простое: * Если 95% ДИ для OR (или RR) включает 1, то эффект не является статистически значимым (p ≥ 0.05). * Если 95% ДИ для OR (или RR) не включает 1, то эффект является статистически значимым (p < 0.05).

В нашем примере: * 95% ДИ для OR: [1.15, 3.65] (не включает 1) -> значимо. * 95% ДИ для RR: [1.17, 2.14] (не включает 1) -> значимо.

Способ 2: Через Z-тест (аппроксимация нормальным распределением)

Этот метод заключается в проверке, отличается ли логарифм вашей меры эффекта (\(\ln(OR)\) или \(\ln(RR)\)) от нуля (поскольку \(\ln(1)=0\)) на величину, большую, чем ожидается случайностью.

  1. Рассчитайте Z-статистику: \[ Z = \frac{\ln(\text{мера эффекта})}{SE(\ln(\text{меры эффекта}))} \]
  2. Найдите p-value: Полученное значение Z сравнивается со стандартным нормальным распределением.

Способ 3: Через критерий хи-квадрат Пирсона

Это третий, классический способ проверки значимости для таблицы 2x2. Он проверяет нулевую гипотезу об отсутствии ассоциации между строками и столбцами (т.е. о равенстве пропорций).

Важнейший факт: Для таблицы 2x2 все три способа (анализ ДИ, Z-тест и критерий хи-квадрат) математически эквивалентны и дадут практически одинаковое p-value.

Связь между Z-тестом и хи-квадрат прямая: \[ \chi^2 = Z^2 \]

Таким образом, вы можете рассчитать OR и RR, а затем проверить их значимость с помощью простого и интуитивного критерия хи-квадрат.

11.0.9.5 Расчет в R

# --- Данные из таблицы сопряженности ---
a <- 45; b <- 55
c <- 30; d <- 75

# --- Расчет OR и его 95% ДИ ---
OR <- (a * d) / (b * c)
SE_log_OR <- sqrt(1/a + 1/b + 1/c + 1/d)
OR_ci_lower <- exp(log(OR) - 1.96 * SE_log_OR)
OR_ci_upper <- exp(log(OR) + 1.96 * SE_log_OR)

# --- Расчет RR и его 95% ДИ ---
RR <- (a / (a+b)) / (c / (c+d))
SE_log_RR <- sqrt( (1/a) - (1/(a+b)) + (1/c) - (1/(c+d)) )
RR_ci_lower <- exp(log(RR) - 1.96 * SE_log_RR)
RR_ci_upper <- exp(log(RR) + 1.96 * SE_log_RR)

# --- Проверка значимости с помощью хи-квадрат ---
observed_matrix <- matrix(c(a, c, b, d), nrow = 2, byrow = FALSE) # Важен порядок!
chi_sq_test <- chisq.test(observed_matrix, correct = FALSE)
p_value <- chi_sq_test$p.value


# --- Вывод результатов ---
cat("Отношение Шансов (OR):", round(OR, 2), "\n")
cat("95% ДИ для OR:", round(OR_ci_lower, 2), "to", round(OR_ci_upper, 2), "\n\n")

cat("Относительный Риск (RR):", round(RR, 2), "\n")
cat("95% ДИ для RR:", round(RR_ci_lower, 2), "to", round(RR_ci_upper, 2), "\n\n")

cat("p-value (по критерию Хи-квадрат):", p_value, "\n")

Интерпретация результата: * OR = 2.05: Шансы улучшения в группе лечения в 2.05 раза выше, чем в группе плацебо. * RR = 1.58: Риск улучшения в группе лечения в 1.58 раза выше, чем в группе плацебо. * Поскольку 95% ДИ для обеих мер не включает 1, а p-value < 0.05, результат является статистически значимым.

11.0.10 11.3.3 Критерий хи-квадрат Пирсона

Критерий хи-квадрат Пирсона (\(\chi^2\)) — это один из самых распространенных непараметрических методов для анализа категориальных данных. Он используется для проверки, существует ли статистически значимая связь между двумя категориальными переменными.

В нашем контексте (сравнение двух групп с бинарным исходом) он проверяет нулевую гипотезу о том, что пропорции успехов в двух группах равны.

11.0.10.1 Основная идея и ожидаемые частоты

Логика теста заключается в сравнении наблюдаемых (Observed, O) частот в ячейках таблицы сопряженности с ожидаемыми (Expected, E) частотами.

  • Наблюдаемые частоты (O): Это реальные данные из вашего исследования (ячейки a, b, c, d).
  • Ожидаемые частоты (E): Это теоретические частоты, которые мы бы ожидали увидеть, если бы нулевая гипотеза об отсутствии связи была верна.

Ожидаемая частота для каждой ячейки рассчитывается по формуле:

\[ E_{ij} = \frac{(\text{Сумма по строке } i) \times (\text{Сумма по столбцу } j)}{\text{Общий итог (N)}} \]

Для нашей таблицы:

Улучшение (+) Нет улучшения (-) Всего
Лечение O11 = a O12 = b R1 = a+b
Плацебо O21 = c O22 = d R2 = c+d
Всего C1 = a+c C2 = b+d N

Ожидаемые частоты (\(E_{ij}\)) будут:

  • \(E_{11} = \frac{R_1 \cdot C_1}{N}\)
  • \(E_{12} = \frac{R_1 \cdot C_2}{N}\)
  • \(E_{21} = \frac{R_2 \cdot C_1}{N}\)
  • \(E_{22} = \frac{R_2 \cdot C_2}{N}\)

11.0.10.2 Вычисление статистики \(\chi^2\)

Статистика хи-квадрат суммирует квадраты стандартизированных разностей между наблюдаемыми и ожидаемыми частотами по всем ячейкам таблицы.

\[ \chi^2 = \sum_{i=1}^{R} \sum_{j=1}^{C} \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]

Для таблицы 2x2 эта формула развернется в сумму четырех слагаемых: \[ \chi^2 = \frac{(a - E_{11})^2}{E_{11}} + \frac{(b - E_{12})^2}{E_{12}} + \frac{(c - E_{21})^2}{E_{21}} + \frac{(d - E_{22})^2}{E_{22}} \]

Чем больше расхождение между наблюдаемыми и ожидаемыми частотами, тем больше будет значение \(\chi^2\).

11.0.10.3 Степени свободы и p-value

Чтобы определить, является ли полученное значение \(\chi^2\) статистически значимым, мы сравниваем его с критическим значением из распределения хи-квадрат.

Количество степеней свободы (\(df\)) для таблицы сопряженности рассчитывается как:

\[ df = (R - 1) \times (C - 1) \] где R — количество строк, а C — количество столбцов. Для таблицы 2x2: \[ df = (2 - 1) \times (2 - 1) = 1 \]

p-value — это вероятность получить значение \(\chi^2\) равное или большее, чем наше наблюдаемое значение, при условии, что верна нулевая гипотеза.

11.0.10.4 Условия применения

  1. Случайная выборка: Данные должны быть получены случайным образом.
  2. Независимость наблюдений: Каждое наблюдение должно быть независимым от других.
  3. Размер ожидаемых частот: Тест надежно работает, если ожидаемые частоты во всех ячейках таблицы 2x2 равны 5 или более. Если это условие не выполняется, следует использовать точный критерий Фишера.

11.0.10.5 Расчет в R

В R для этого используется функция chisq.test().

# --- Данные из таблицы сопряженности ---
# Создаем матрицу с наблюдаемыми частотами
observed_matrix <- matrix(c(45, 30, 55, 75), nrow = 2, byrow = TRUE)
rownames(observed_matrix) <- c("Лечение", "Плацебо")
colnames(observed_matrix) <- c("Улучшение", "Нет улучшения")

print("Наблюдаемая таблица:")
print(observed_matrix)

# --- Расчет с помощью R ---
# chisq.test() по умолчанию применяет коррекцию непрерывности Йейтса для таблиц 2x2.
# Чтобы получить "чистый" критерий Пирсона, отключаем ее: correct = FALSE
chi_sq_result <- chisq.test(observed_matrix, correct = FALSE)

# Выводим результат
print(chi_sq_result)

Анализ вывода chisq.test(): * X-squared: Рассчитанное значение статистики \(\chi^2\). * df: Степени свободы (для таблицы 2x2 всегда 1). * p-value: p-value теста. Если p < 0.05, мы отвергаем \(H_0\) и делаем вывод о наличии статистически значимой связи между переменными. * Observed: Показывает таблицу с наблюдаемыми частотами. * Expected: Показывает таблицу с ожидаемыми частотами, рассчитанными по формуле.

Как видите, p-value для критерия хи-квадрата (0.0069) практически идентично p-value для Z-теста пропорций (0.0068), что подтверждает их математическую эквивалентность для таблиц 2x2.

11.0.11 Сравнение Z-теста для пропорций и критерия хи-квадрат Пирсона

Хотя Z-тест для пропорций и критерий хи-квадрат Пирсона могут показаться разными (один сравнивает две пропорции, другой — целые таблицы), для таблиц сопряженности размером 2x2 они математически эквивалентны и дадут практически идентичные p-value.

Давайте разберем их сходства и различия.

11.0.11.1 Математическая эквивалентность

Ключевая связь между этими двумя тестами заключается в том, что статистика Z-теста, возведенная в квадрат, равна статистике хи-квадрат.

\[ \chi^2 = Z^2 \]

Это означает, что если вы проведете оба теста для одних и тех же данных, вы получите одно и то же статистическое решение (отвергнуть или не отвергнуть нулевую гипотезу).

Давайте проверим это на нашем примере:

# --- Данные ---
successes <- c(45, 30)
totals <- c(100, 105)

# --- Z-тест для пропорций ---
# prop.test() возвращает X-squared, но мы можем извлечь Z-статистику
# Z = sqrt(X-squared) с правильным знаком
z_stat <- sqrt(prop_test_result$statistic)
if (prop_test_result$estimate[1] < prop_test_result$estimate[2]) {
  z_stat <- -z_stat
}
p_value_z <- prop_test_result$p.value


# --- Критерий хи-квадрат ---
chi_sq_stat <- chi_sq_result$statistic
p_value_chi <- chi_sq_result$p.value

# --- Сравнение ---
cat("Z-статистика:", z_stat, "\n")
cat("p-value (Z-тест):", p_value_z, "\n\n")

cat("Статистика Хи-квадрат:", chi_sq_stat, "\n")
cat("p-value (Хи-квадрат):", p_value_chi, "\n\n")

cat("Z^2:", z_stat^2, "\n")
cat("Хи-квадрат:", chi_sq_stat, "\n")

Как видите, Z^2 и Хи-квадрат совпадают (с точностью до округления), и p-value также идентичны.

11.0.11.2 Различия в подходе и интерпретации

Несмотря на математическую эквивалентность, концептуально эти тесты подходят к проблеме с разных сторон.

Аспект Z-тест для пропорций Критерий хи-квадрат Пирсона
Основной фокус Сравнение двух пропорций (\(\hat{p}_1\) vs \(\hat{p}_2\)). Анализ ассоциации/связи между двумя категориальными переменными.
Нулевая гипотеза (\(H_0\)) \(H_0: p_1 = p_2\) (пропорции равны). \(H_0\): Переменные независимы (нет связи между группой и исходом).
Статистика Z-статистика, которая следует стандартному нормальному распределению. Статистика \(\chi^2\), которая следует распределению хи-квадрат с 1 степенью свободы.
Основной результат Разница пропорций и ее 95% ДИ. Статистика \(\chi^2\) и p-value.
Сфера применения Прямое сравнение двух групп (например, лечение vs. плацебо). Более универсален. Может применяться для таблиц любого размера (2x2, 2x3, 3x3 и т.д.), а не только для двух пропорций.

11.0.11.3 Когда какой тест предпочесть?

  1. Для сравнения двух пропорций (таблица 2x2):
    • Оба теста дадут одинаковый p-value.
    • Z-тест часто предпочтительнее, потому что он напрямую предоставляет разницу пропорций и ее 95% доверительный интервал. Эта мера эффекта (Risk Difference) часто более интерпретируема, чем просто наличие/отсутствие связи.
    • Критерий хи-квадрат просто говорит вам “есть ли связь”, но не характеризует ее величину так же удобно.
  2. Для таблиц сопряженности большего размера (например, 2x3 или 3x3):
    • Здесь Z-тест неприменим, так как он сравнивает только две пропорции.
    • Критерий хи-квадрат является единственным выбором для проверки наличия ассоциации между переменными.

11.0.12 Итог

  • Для таблиц 2x2: Z-тест и критерий хи-квадрат — это два способа сделать одно и то же. Z-тест часто более информативен, так как дает оценку разницы пропорций.
  • Для больших таблиц: Критерий хи-квадрат является стандартным инструментом для анализа связи между категориальными переменными.

11.0.13 11.3.4 Точный тест Фишера

Точный тест Фишера (Fisher’s Exact Test) — это непараметрический статистический тест, который используется для анализа таблиц сопряженности, особенно когда размеры выборок малы. В отличие от критерия хи-квадрата, который является аппроксимацией, тест Фишера вычисляет точное p-value.

11.0.13.1 Когда использовать тест Фишера?

Основное правило: используйте точный тест Фишера, если критерий хи-квадрат неприменим.

Критерий хи-квадрат является надежной аппроксимацией, только если ожидаемые частоты во всех ячейках таблицы 2x2 равны 5 или более.

Если это условие нарушается (например, одна или несколько ожидаемых частот меньше 5), p-value, полученное с помощью хи-квадрата, может быть неточным. В таких случаях точный тест Фишера является предпочтительным методом.

11.0.13.2 Основная идея теста

Тест Фишера отвечает на вопрос: “Какова вероятность получить наблюдаемую (или еще более экстремальную) таблицу при условии, что итоговые суммы по строкам и столбцам фиксированы, а нулевая гипотеза (об отсутствии связи) верна?”

Логика: 1. Фиксируются итоги: Тест предполагает, что общее количество наблюдений в каждой группе (суммы по строкам) и общее количество каждого исхода (суммы по столбцам) являются фиксированными. 2. Перебираются все возможные таблицы: Рассматриваются все возможные комбинации чисел в ячейках (a, b, c, d), которые удовлетворяют этим фиксированным итогам. 3. Рассчитывается вероятность каждой таблицы: Для каждой возможной таблицы рассчитывается ее вероятность с помощью гипергеометрического распределения. 4. Суммируются вероятности: Точечное p-value — это сумма вероятностей наблюдаемой таблицы и всех “еще более экстремальных” таблиц (т.е. тех, где связь между переменными выражена еще сильнее).

“Экстремальность” определяется по величине отклонения от того, что ожидалось бы при отсутствии связи.

11.0.13.3 Расчет в R

В R для проведения точного теста Фишера используется функция fisher.test(). Она работает так же, как и chisq.test().

Пример 1: Маленькая выборка, где тест Фишера необходим

Представим, что у нас гораздо меньше данных.

Улучшение (+) Нет улучшения (-) Всего
Лечение 8 2 10
Плацебо 3 7 10
Всего 11 9 20

Здесь ожидаемые частоты будут меньше 5, поэтому хи-квадрат не подходит.

# --- Данные для маленькой выборки ---
small_matrix <- matrix(c(8, 3, 2, 7), nrow = 2, byrow = FALSE)
rownames(small_matrix) <- c("Лечение", "Плацебо")
colnames(small_matrix) <- c("Улучшение", "Нет улучшения")

print("Таблица для маленькой выборки:")
print(small_matrix)

# --- Расчет с помощью точного теста Фишера ---
fisher_result_small <- fisher.test(small_matrix)

# Выводим результат
print(fisher_result_small)

Анализ вывода fisher.test(): * p-value: Точное p-value. В нашем примере оно равно 0.06979. * odds ratio: Тест Фишера также рассчитывает отношение шансов (OR) и его 95% доверительный интервал.

Пример 2: Сравнение с хи-квадрат на наших исходных данных

Давайте применим тест Фишера к нашим основным данным (где хи-квадрат был применим), чтобы сравнить результаты.

# --- Данные из нашего основного примера ---
observed_matrix <- matrix(c(45, 30, 55, 75), nrow = 2, byrow = FALSE)
rownames(observed_matrix) <- c("Лечение", "Плацебо")
colnames(observed_matrix) <- c("Улучшение", "Нет улучшения")

# --- Расчет с помощью точного теста Фишера ---
fisher_result_main <- fisher.test(observed_matrix)

# Выводим результат
print(fisher_result_main)

Сравните p-value с результатом chisq.test: * p-value (Хи-квадрат): 0.0069 * p-value (Фишер): 0.0062

Как видите, при больших выборках результаты очень близки. Но при малых выборках они могут сильно отличаться.

11.0.14 Итог: Хи-квадрат vs. Тест Фишера

Характеристика Критерий хи-квадрат Точный тест Фишера
Тип Аппроксимация Точный расчет
Условие применения Ожидаемые частоты во всех ячейках ≥ 5 Нет ограничений на размер выборки
Вычислительная сложность Очень низкая Высокая (особенно для больших таблиц)
Когда использовать Стандартный выбор для больших выборок Обязателен для малых выборок или когда есть ячейки с ожидаемой частотой < 5

Практическая рекомендация: Если сомневаетесь, используйте точный тест Фишера. Современные компьютеры легко справляются с вычислениями, и он является более точным и универсальным методом.

11.0.15 11.4 Сравнение количественных показателей в нескольких независимых выборках

Когда нам нужно сравнить средние значения более чем двух независимых групп, использование множественных t-тестов некорректно. Каждое сравнение увеличивает вероятность совершить ошибку первого рода (ложноположительный результат). Для решения этой проблемы используется дисперсионный анализ (ANOVA).

Отличная идея! Такой наглядный график сделает концепцию разложения дисперсии гораздо более понятной. Давайте создадим этот трехчастный график.

11.0.16 11.4.1 Однофакторный дисперсионный анализ (One-Way ANOVA)

Однофакторный ANOVA — это статистический метод, который сравнивает средние значения трех или более групп. Фактически, он проверяет, являются ли различия между средними значимыми в целом.

11.0.16.1 Исследовательские гипотезы в ANOVA

  • Нулевая гипотеза (\(H_0\)): Средние значения во всех группах равны. \(H_0: \mu_1 = \mu_2 = ... = \mu_k\) (где \(k\) — количество групп).
  • Альтернативная гипотеза (\(H_1\)): Хотя бы одно среднее значение отличается от остальных. \(H_1\): Не все \(\mu\) равны.

Важно понимать, что если ANOVA показывает значимый результат (мы отвергаем \(H_0\)), это не говорит нам, какие именно группы различаются между собой. Для этого требуются пост-хок тесты (например, Тьюки).

11.0.16.2 Расчет сумм квадратов (Sum of Squares)

Чтобы количественно оценить дисперсии, ANOVA использует суммы квадратов отклонений (Sum of Squares, SS).

  • Общая сумма квадратов (Total Sum of Squares, \(SS_T\)): Общая изменчивость всех данных относительно общего среднего (\(\bar{X}_{grand}\)). \[ SS_T = \sum_{i=1}^{k} \sum_{j=1}^{n_i} (X_{ij} - \bar{X}_{grand})^2 \]
    • \(k\) - количество групп.
    • \(n_i\) - количество наблюдений в \(i\)-й группе.
    • \(X_{ij}\) - \(j\)-е наблюдение в \(i\)-й группе.
    • \(\bar{X}_{grand}\) - общее среднее по всем данным.
  • Межгрупповая сумма квадратов (Between-Group Sum of Squares, \(SS_B\)): Изменчивость, обусловленная различиями между средними групп (\(\bar{X}_i\)) и общим средним. \[ SS_B = \sum_{i=1}^{k} n_i (\bar{X}_i - \bar{X}_{grand})^2 \]
  • Внутригрупповая сумма квадратов (Within-Group Sum of Squares, \(SS_W\)): Сумма изменчивостей внутри каждой группы. \[ SS_W = \sum_{i=1}^{k} \sum_{j=1}^{n_i} (X_{ij} - \bar{X}_i)^2 \]

Эти три суммы квадратов связаны простым соотношением: \[ SS_T = SS_B + SS_W \]

11.0.16.3 Концепция разложения дисперсии

Основная идея ANOVA заключается в том, чтобы разложить общую изменчивость (дисперсию) в данных на две компоненты:

  1. Внутригрупповая дисперсия (Within-Group Variance):
    • Это изменчивость внутри каждой группы. Она отражает случайную вариацию данных, не связанную с фактором группировки.
    • Малая внутригрупповая дисперсия означает, что данные внутри каждой группы сгруппированы плотно вокруг своего среднего.
  2. Межгрупповая дисперсия (Between-Group Variance):
    • Это изменчивость между средними значениями разных групп. Она отражает эффект от фактора группировки (например, от разных видов лечения).
    • Большая межгрупповая дисперсия означает, что средние значения групп сильно различаются.

11.0.16.4 F-статистика и ее распределение

Чтобы проверить гипотезу, мы сравниваем межгрупповую и внутригрупповую дисперсии. Для этого мы превращаем суммы квадратов в средние квадраты (Mean Squares, MS), разделив их на соответствующее количество степеней свободы (\(df\)).

  • Степени свободы:
    • \(df_B = k - 1\) (между группами)
    • \(df_W = N - k\) (внутри групп, где N - общее число наблюдений)
  • Средние квадраты:
    • \(MS_B = \frac{SS_B}{df_B}\)
    • \(MS_W = \frac{SS_W}{df_W}\)

F-статистика — это отношение межгрупповой дисперсии к внутригрупповой:

\[ F = \frac{MS_B}{MS_W} = \frac{\text{Вариативность между группами}}{\text{Вариативность внутри групп}} \]

Как соотносятся дисперсии и F-статистика:

  • Если верна нулевая гипотеза (\(H_0\)): Все средние равны. Межгрупповая дисперсия (\(MS_B\)) оценивает ту же самую случайную изменчивость, что и внутригрупповая (\(MS_W\)). Их отношение (\(F\)) будет близко к 1.
  • Если верна альтернативная гипотеза (\(H_1\)): Хотя бы одно среднее отличается. Межгрупповая дисперсия (\(MS_B\)) будет больше внутригрупповой (\(MS_W\)), и их отношение (\(F\)) будет значительно больше 1.

Распределение F-статистики:

Если нулевая гипотеза верна, F-статистика следует F-распределению с двумя параметрами: степенями свободы для числителя (\(df_1 = k-1\)) и знаменателя (\(df_2 = N-k\)).

Мы вычисляем наблюдаемое значение F и сравниваем его с этим распределением, чтобы найти p-value — вероятность получить такое (или большее) значение F случайно.

11.0.16.5 Расчет в R

В R для проведения однофакторного ANOVA используется функция aov().

# Создадим данные для 3-х групп
set.seed(42)
group1 <- rnorm(30, mean = 10, sd = 2)
group2 <- rnorm(30, mean = 15, sd = 2)
group3 <- rnorm(30, mean = 12, sd = 2)

# Создаем фрейм данных в "длинном" формате, который требуется для aov()
anova_data <- data.frame(
  value = c(group1, group2, group3),
  group = factor(rep(c("Группа 1", "Группа 2", "Группа 3"), each = 30))
)

# Проводим однофакторный ANOVA
anova_result <- aov(value ~ group, data = anova_data)

# Выводим результат в виде ANOVA-таблицы
summary(anova_result)

Анализ ANOVA-таблицы: * Df: Степени свободы (df_B и df_W). * Sum Sq: Суммы квадратов (SS_B и SS_W). * Mean Sq: Средние квадраты (MS_B и MS_W). * F value: Рассчитанное значение F-статистики. * Pr(>F): p-value. Если p < 0.05, мы отвергаем \(H_0\) и делаем вывод, что средние хотя бы двух групп статистически значимо различаются.

11.0.17 11.4.1.1 Допущения для ANOVA и способы их проверки

Однофакторный дисперсионный анализ (ANOVA) является параметрическим методом, а это значит, что его корректность зависит от выполнения нескольких ключевых допущений относительно данных. Если эти допущения серьезно нарушены, результаты ANOVA (особенно p-value) могут быть неточными.

Вот три основных допущения и способы их проверки.

11.0.17.1 1. Независимость наблюдений (Independence of Observations)

Это самое важное допущение. Оно означает, что значение одной переменной не должно быть связано со значением другой переменной. Другими словами, данные в одной группе не должны быть связаны с данными в другой группе, и внутри группы наблюдения также должны быть независимы.

  • Пример нарушения: Измерение артериального давления у одних и тех же пациентов несколько раз без учета этого фактора; измерение показателей у студентов одной группы, которые общаются друг с другом и могут влиять на ответы друг друга.
  • Как проверить: Это допущение проверяется в первую очередь через понимание дизайна исследования. Нет простого статистического теста для проверки независимости. Вы должны быть уверены, что сбор данных проводился корректно (например, использовалась случайная выборка, не было “заражения” данных между группами).
  • Что делать, если нарушено: Если допущение нарушено, результаты ANOVA недействительны. Нужно пересмотреть дизайн исследования или использовать статистические модели, которые учитывают зависимость (например, дисперсионный анализ с повторными измерениями или смешанные модели).

11.0.17.2 2. Нормальность распределения остатков (Normality of Residuals)

ANOVA не требует, чтобы данные в каждой группе были распределены нормально. Требование заключается в том, чтобы остатки (ошибки модели), то есть отклонения каждого наблюдения от среднего своей группы, были распределены нормально.

  • Как проверить:
    1. Визуальный анализ (основной способ): Постройте график квантиль-квантиль (Q-Q plot) для остатков модели. Если остатки распределены нормально, точки на графике должны лежать примерно вдоль прямой линии.
    2. Формальные тесты: Можно использовать тесты на нормальность, например, тест Шапиро-Уилка (shapiro.test()), примененный к остаткам.

    Важное замечание: На больших выборках (например, > 30-50 наблюдений в группе) тесты на нормальность становятся слишком чувствительными и могут отвергать гипотезу о нормальности из-за незначительных отклонений. В таких случаях визуальный анализ (Q-Q plot) является более надежным индикатором.

  • Что делать, если нарушено:
    • Преобразование данных: Попробуйте применить математическую трансформацию к зависимой переменной (например, логарифмическую log() или преобразование Бокса-Кокса boxcox()), чтобы сделать распределение остатков более близким к нормальному.
    • Использовать непараметрический аналог: Если преобразование не помогает, используйте H-критерий Краскела-Уоллиса. Он не требует предположения о нормальности и является прямым аналогом ANOVA для непараметрических данных.

11.0.17.3 3. Равенство дисперсий (Homogeneity of Variances)

Это допущение, также известное как гомоскедастичность, означает, что дисперсии (изменчивость) внутри каждой из сравниваемых групп должны быть примерно равны.

  • Как проверить:
    1. Визуальный анализ: Постройте график “ящик с усами” (boxplot) для каждой группы. Визуально сравните длины “ящиков” и “усов”. Если они примерно одинаковы, допущение, вероятно, выполняется. Сильные различия в длине указывают на неравенство дисперсий.
    2. Формальные тесты:
      • Тест Левена (Levene’s Test): Является наиболее популярным и надежным тестом. Он менее чувствителен к отклонениям от нормальности, чем тест Барлетта. В R реализован в функции leveneTest() из пакета car.
      • Тест Барлетта (Bartlett’s Test): Более чувствителен к отклонениям от нормальности. Используйте его, если вы уверены, что распределение остатков нормально.
  • Что делать, если нарушено:
    • Использовать ANOVA с поправкой Уэлча: R по умолчанию использует поправку Уэлча в функции oneway.test(), которая не требует равенства дисперсий. Это самый простой и распространенный способ.
    • Преобразование данных: Иногда трансформация (например, логарифмическая) может также помочь выровнять дисперсии.
    • Использовать непараметрический аналог: H-критерий Краскела-Уоллиса также не требует предположения о равенстве дисперсий.

11.0.17.4 Пример с гетероскедастичными (неравными дисперсиями) данными

Давайте создадим данные, где средние могут быть одинаковы, но дисперсии (разброс) сильно различаются.

# install.packages("car")
# install.packages("ggplot2")
library(car)
library(ggplot2)

# Создаем данные с разными дисперсиями
set.seed(42)
hetero_data <- data.frame(
  value = c(rnorm(30, mean = 10, sd = 2),  # Группа 1: sd = 2
             rnorm(30, mean = 10, sd = 10), # Группа 2: sd = 10 (большой разброс)
             rnorm(30, mean = 10, sd = 2)), # Группа 3: sd = 2
  group = factor(rep(c("Группа 1", "Группа 2", "Группа 3"), each = 30))
)

# --- 1. Визуальная проверка равенства дисперсий (Boxplot) ---
ggplot(hetero_data, aes(x = group, y = value, fill = group)) +
  geom_boxplot() +
  labs(title = "График 'Ящик с усами' для данных с неравными дисперсиями",
       subtitle = "Видно, что 'ящик' для Группы 2 гораздо выше",
       x = "Группа", y = "Значение") +
  theme_minimal() +
  guides(fill = "none")

# --- 2. Формальная проверка равенства дисперсий (Тест Левена) ---
levene_test_hetero <- leveneTest(value ~ group, data = hetero_data)
print(levene_test_hetero)

Анализ результатов:

  • На графике: Ящик с усами для “Группы 2” значительно выше и шире, чем для “Группы 1” и “Группы 3”. Это явный визуальный признак гетероскедастичности (неравенства дисперсий).
  • Тест Левена: p-value = 0.0006355. Поскольку p-value < 0.05, мы отвергаем нулевую гипотезу о равенстве дисперсий. Тест формально подтверждает то, что мы увидели на графике.

Что делать в этом случае?

Поскольку допущение о равенстве дисперсий нарушено, стандартный ANOVA (функция aov()) может дать неверный p-value. Мы должны использовать более надежный метод.

# --- 3. Используем ANOVA с поправкой Уэлча ---
# Функция oneway.test() по умолчанию использует поправку Уэлча
oneway_test_hetero <- oneway.test(value ~ group, data = hetero_data)

print(oneway_test_hetero)

Результат oneway.test() показывает p-value (p-value = 0.4241), которое корректно учитывает неравенство дисперсий. В данном случае, p-value > 0.05, и мы делаем вывод, что, несмотря на разный разброс, средние значения в группах статистически значимо не различаются.

11.0.18 Итог

  1. Независимость: Проверяется через дизайн исследования.
  2. Нормальность остатков: Проверяется визуально (Q-Q plot) и с помощью теста Шапиро-Уилка.
  3. Равенство дисперсий: Проверяется визуально (boxplot) и с помощью теста Левена.

Если допущения нарушены, не паникуйте. Существуют надежные альтернативы: ANOVA с поправкой Уэлча или непараметрический H-критерий Краскела-Уоллиса.