Здесь приведён мой код на языке программирования R для проведения статистических расчётов по всем примерам из книги «Занимательная статистика. Манга», опубликованной издательством Додэка-XXI. Документ был подготовлен при помощи языка разметки R Markdown, который позволяет представлять в красивом виде код на языке R вместе с результатами его выполнения. Все числа округлены так же, как и в книге.

Я хочу поблагодарить Джеймса Чёрча (James Church) за разрешение использовать свой код в качестве отправной точки для моего варианта. Я практически полностью переработал его код на свой лад, но кое-какие элементы ещё остались. Его код к американскому изданию книги, содержащий только базовые функции R без использования дополнительных пакетов, опубликован на сайте издательства No Starch Press: www.nostarch.com/download/MangaGuideToStatistics.R.

В моём варианте кода я использовал три довольно популярных в последнее время пакета, не входящих в базовую поставку R:

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

library(magrittr)
library(dplyr)
library(ggplot2)

Глава 1: Разберёмся с типами данных

В этой главе содержится только теория без каких-либо статистических расчётов.

Глава 2: Знакомимся с количественными данными

Стр. 33: Таблицу с ценами на лапшу рамэн можно представить в виде простого вектора (так в R называются одномерные массивы):

ramen.prices <- c(700, 850, 600, 650, 980, 750, 500, 890, 880, 700,
                  890, 720, 680, 650, 790, 670, 680, 900, 880, 720,
                  850, 700, 780, 850, 750, 780, 590, 650, 580, 750,
                  800, 550, 750, 700, 600, 800, 800, 880, 790, 790,
                  780, 600, 670, 680, 650, 890, 930, 650, 777, 700)

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

Названия столбцов в R всё-таки лучше задавать латинскими буквами. Для наглядности, перед выводом на экран я буду менять названия на русские и округлять числовые значения. Для этого будет служить функция Output:

# функция для смены названий столбцов и округления численных значений перед выводом на экран
Output <- function(df, newnames, digits = 1){
  nums <- sapply(df, is.numeric)
  df[nums] <- round(df[nums], digits)
  names(df) <- newnames[names(df)]
  return(df)
}

# функция для создания таблицы аналогичной приведённым на стр. 38 и 56
DistributionTable <- function(hist.data){
  interval.lower <- hist.data$breaks[-length(hist.data$breaks)]
  interval.upper <- hist.data$breaks[-1]
  hist.table <- data.frame(
    Interval = paste(interval.lower, "-", interval.upper),
    Midpoint = hist.data$mids,
    Frequency = hist.data$counts,
    Rel.Frequency = hist.data$counts / sum(hist.data$counts)
  )
  return(hist.table)
}

ramen.hist <- hist(ramen.prices, breaks = 5, right = FALSE, plot = FALSE)

# используем вектор соответсвий названий столбцов для возможности лёгкого переименования
hist.table.rus <- c(Interval = "Интервал", Midpoint = "Середина интервала", Frequency = "Частота", Rel.Frequency = "Относительная частота")

DistributionTable(ramen.hist) %>% Output(hist.table.rus, 2)
Интервал Середина интервала Частота Относительная частота
500 - 600 550 4 0.08
600 - 700 650 13 0.26
700 - 800 750 18 0.36
800 - 900 850 12 0.24
900 - 1000 950 3 0.06

Стр. 39: Гистограмма цен на рамэн построена при помощи функции ggplot из пакета ggplot2. График был построен таким образом, чтобы он был похожим на книжный вариант.

ggplot(NULL, aes(x = ramen.prices)) + theme_bw() +
  geom_histogram(binwidth = 100, colour = "black", fill = "gray") +
  ggtitle("Гистограмма цен на рамэн 50-ти лучших ресторанов") +
  xlab("Цена на рамэн") + ylab("Частота (количество ресторанов") +
  coord_cartesian(xlim = c(490, 1010), ylim = c(0, 20)) +
  scale_x_continuous(breaks = seq(550, 950, 100))

Стр. 41: Введём таблицу (фрейм) данных, содержащую результаты игры в боулинг для трёх команд:

bowling <- data.frame(
  team.A = c(86, 73, 124, 111, 90, 38),
  team.B = c(84, 71, 103, 85, 90, 89),
  team.C = c(229, 77, 59, 95, 70, 88)
)

teams.rus <- c(team.A = "Команда А", team.B = "Команда Б", team.C = "Команда В")

bowling %>% Output(teams.rus)
Команда А Команда Б Команда В
86 84 229
73 71 77
124 103 59
111 85 95
90 90 70
38 89 88

Стр. 42: Средние значения можно рассчитать при помощи встроенной функции mean. А функция summarise_each из пакета dplyr позволяет применить какую-либо функцию для каждого столбца таблицы. Поэтому можно легко получить средние результаты сразу для всех трёх команд:

bowling %>% summarise_each(funs(mean)) %>% Output(teams.rus)
Команда А Команда Б Команда В
87 87 103

Стр. 43: В R нет встроенной функции для расчёта cредней геометрической, но её очень легко можно написать самим:

GeometricMean <- function(values) {
  gmean <- prod(values) ^ (1 / length(values))
  return(gmean)
}

bowling %>% summarise_each(funs(GeometricMean)) %>% Output(teams.rus, 2)
Команда А Команда Б Команда В
81.61 86.48 92.06

Стр. 43: Средние гармонические результатов игры в боулинг для трёх команд:

HarmonicMean <- function(values) {
  hmean <- 1 / (sum(1 / values) / length(values))
  return(hmean)
}

bowling %>% summarise_each(funs(HarmonicMean)) %>% Output(teams.rus, 2)
Команда А Команда Б Команда В
75.16 85.95 85.13

Стр. 45: Медианные значения результатов трёх команд (для расчёта медианы существует встроенная функция median):

bowling %>% summarise_each(funs(median)) %>% Output(teams.rus)
Команда А Команда Б Команда В
88 87 82.5

Стр. 50: Стандартное отклонение (\(\sigma\)). Существует встроенная функция sd для расчёта несмещённой оценки стандарного отклонения генеральной совокупности по данным выборки: \(\sigma = \sqrt{\frac{1}{N-1} \sum_{i=1}^N (x_i-\overline{x})^2}\) (см. стр. 52). Но опять же можно легко написать функцию для расчётов по формуле для генеральной совокупности \(\sigma = \sqrt{\frac{1}{N} \sum_{i=1}^N (x_i-\overline{x})^2}\):

StdDev <- function(values) {
    m <- mean(values)
    stddev <- sqrt(sum((values - m) ^ 2) / length(values))
    return(stddev)
}

bowling %>% summarise_each(funs(StdDev)) %>% Output(teams.rus, 1)
Команда А Команда Б Команда В
27.5 9.5 57.5

Стр. 55: Количество интервалов группировки по правилу Стерджеса (используя встроенную функцию nclass.Sturges) и величина интервала:

cat("Применим правило Стерджеса на данные о ценах на рамэн в 50-ти лучших ресторанах:\n")
nclass.Sturges(ramen.prices) %>% paste("Количество интервалов:", ., "\n") %>% cat()

BreakRange <- function(values) {
    return(ceiling((max(values) - min(values)) / nclass.Sturges(values)))
}

BreakRange(ramen.prices) %>% paste("Величина интервала:", ., "\n") %>% cat()
Применим правило Стерджеса на данные о ценах на рамэн в 50-ти лучших ресторанах:
Количество интервалов: 7 
Величина интервала: 69

Стр. 56: Таблица распределения цены на рамэн в 50-ти лучших ресторанах (группировка по правилу Стерджеса):

# вектор с границами интервалов группировки
ramen.breaks <- seq(min(ramen.prices), by = BreakRange(ramen.prices),
                    length.out = nclass.Sturges(ramen.prices) + 1)
ramen.breaks
[1] 500 569 638 707 776 845 914 983
ramen.histS <- hist(ramen.prices, breaks = ramen.breaks, right = FALSE, plot = FALSE)
DistributionTable(ramen.histS) %>% Output(hist.table.rus, 2)
Интервал Середина интервала Частота Относительная частота
500 - 569 534.5 2 0.04
569 - 638 603.5 5 0.10
638 - 707 672.5 15 0.30
707 - 776 741.5 6 0.12
776 - 845 810.5 10 0.20
845 - 914 879.5 10 0.20
914 - 983 948.5 2 0.04

Стр. 57: Упражнение. Данные по результатам забега на 100 метров:

race.times <- c(16.3, 22.4, 18.5, 18.7, 20.1);
mean(race.times) %>% paste("Среднее значение:", ., "\n") %>% cat()
median(race.times) %>% paste("Медиана:", ., "\n") %>% cat()
StdDev(race.times) %>% round(2) %>% paste("Стандартное отклонение:", ., "\n") %>% cat()
Среднее значение: 19.2 
Медиана: 18.7 
Стандартное отклонение: 2.01

Глава 3: Знакомимся с качественными данными

Стр. 61: Результаты школьного анкетирования по вопросу «Нравится ли вам новая форма?»:

school.uniform.survey <- factor(
  c("Нравится", "Так себе", "Нравится", "Так себе", "Не нравится",
    "Нравится", "Нравится", "Нравится", "Нравится", "Нравится",
    "Нравится", "Нравится", "Так себе", "Нравится", "Нравится",
    "Так себе", "Нравится", "Нравится", "Нравится", "Нравится",
    "Нравится", "Нравится", "Не нравится", "Так себе", "Нравится",
    "Нравится", "Не нравится", "Нравится", "Нравится", "Нравится",
    "Так себе", "Так себе", "Нравится", "Не нравится", "Нравится",
    "Нравится", "Нравится", "Нравится", "Так себе", "Нравится"))

Стр. 62: Частотную таблицу можно получить при помощи функции table. Итоговые суммы можно добавить, используя функцию addmargins:

school.uniform.table <- table(school.uniform.survey)

school.uniform.table %>% addmargins(FUN=list(`Итого` = sum))
Не нравится Нравится Так себе Итого
4 28 8 40

Для получения процентной таблицы применим функцию prop.table:

school.uniform.ptable <- prop.table(school.uniform.table) * 100

school.uniform.ptable %>% addmargins(FUN=list(`Итого` = sum))
Не нравится Нравится Так себе Итого
10 70 20 100

Пользуясь этими функциями, можно сделать таблицу, аналогичную таблице на стр. 62:

CrossTab <- function(fact){
  tab <- addmargins(table(fact), FUN=list(`Итого` = sum))
  ptab <- addmargins(prop.table(table(fact)), FUN=list(`Итого` = sum))
  ctab <- data.frame(Response = unlist(unname(dimnames(tab))),
                     Frequency = as.vector(tab),
                     Procent = as.vector(ptab) * 100)
  return(ctab)
}

ctab.rus <- c(Response = "Вариант ответа", Frequency = "Количество", Procent = "Процент")

CrossTab(school.uniform.survey) %>% Output(ctab.rus)
Вариант ответа Количество Процент
Не нравится 4 10
Нравится 28 70
Так себе 8 20
Итого 40 100

Графическое представление качественной переменной (гистограмма):

ggplot(NULL, aes(x = school.uniform.survey)) + theme_bw() +
  geom_bar(colour = "black", fill = "gray") +
  ggtitle("Результаты анкетирования о новой школьной форме") +
  xlab("Ответ") + ylab("Частота") +
  coord_cartesian(ylim = c(0, 30))

Стр. 64: Упражнение. Какая из двух политических партий победит на следующих выборах?

party.victory <- factor(c("Победит Б", "Победит Б", "Победит Б", "Не знаю", "Победит А", "Победит Б", "Победит А", "Не знаю", "Победит Б", "Победит Б"))

CrossTab(party.victory) %>% Output(ctab.rus)
Вариант ответа Количество Процент
Не знаю 2 20
Победит А 2 20
Победит Б 6 60
Итого 10 100

Глава 4: Нормированное отклонение и рейтинг успеваемости

Стр. 67-69: Результаты тестов в баллах по четырём школьным предметам:

RUS.LETTERS <- c("А", "Б", "В", "Г", "Д", "Е", "Ж", "З", "И", "К", "Л",
                 "М", "Н", "О", "П", "Р", "С", "Т", "У", "Ф", "Х", "Ц",
                 "Ч", "Ш", "Щ", "Ы", "Э", "Ю", "Я")

test.scores <- data.frame(
  Student = c("Руи", "Юми", RUS.LETTERS[1:16]),
  English = c(90, 81, 73, 97, 85, 60, 74, 64, 72, 67, 87, 78, 85, 96, 77 ,100 ,92 ,86),
  Literature = c(71, 90, 79, 70, 67, 66, 60, 83, 57, 85, 93, 89, 78, 74, 65, 78, 53, 80),
  History = c(73, 61, 14, 41, 49, 87, 69, 65, 36, 7, 53, 100, 57, 45, 56, 34, 37, 70),
  Biology = c(59, 73, 47, 38, 63, 56, 15, 53, 80, 50, 41, 62, 44, 26, 91, 35, 53, 68)
)

# вектор для переименования названий столбцов на русский перед выводом на экран
test.rus <- c(Student = "Ученик", English = "Английский язык", Literature = "Литература",
                     History = "История", Biology = "Биология")

test.scores %>% Output(test.rus)
Ученик Английский язык Литература История Биология
Руи 90 71 73 59
Юми 81 90 61 73
А 73 79 14 47
Б 97 70 41 38
В 85 67 49 63
Г 60 66 87 56
Д 74 60 69 15
Е 64 83 65 53
Ж 72 57 36 80
З 67 85 7 50
И 87 93 53 41
К 78 89 100 62
Л 85 78 57 44
М 96 74 45 26
Н 77 65 56 91
О 100 78 34 35
П 92 53 37 53
Р 86 80 70 68

Стр. 68: Средний балл по каждому предмету можно рассчитать при помощи функции summarise_each, отбрасывая при этом столбец с именами учеников (-Student):

test.scores %>% summarise_each(funs(mean), -Student) %>% Output(test.rus)
Английский язык Литература История Биология
81.3 74.3 53 53

Стр. 72: Для нормирования вектора значений применяется функция scale, но по умолчанию она использует выборочное стандарное отклонение (со знаменателем \(n-1\)), поэтому в данном случае придётся поменять параметр scale:

Normalize <- function(values) {
  norm.v <- scale(values, scale = StdDev(values))
  return(norm.v)
}

# Normalize each column except Student
test.scores %>% mutate_each(funs(Normalize), -Student) %>% Output(test.rus, 2)
Ученик Английский язык Литература История Биология
Руи 0.77 -0.30 0.88 0.33
Юми -0.03 1.39 0.35 1.09
А -0.74 0.41 -1.71 -0.33
Б 1.39 -0.39 -0.53 -0.82
В 0.33 -0.65 -0.18 0.55
Г -1.90 -0.74 1.49 0.16
Д -0.65 -1.27 0.70 -2.08
Е -1.54 0.77 0.53 0.00
Ж -0.83 -1.54 -0.75 1.48
З -1.27 0.95 -2.02 -0.16
И 0.50 1.66 0.00 -0.66
К -0.30 1.30 2.07 0.49
Л 0.33 0.33 0.18 -0.49
М 1.30 -0.03 -0.35 -1.48
Н -0.39 -0.83 0.13 2.08
О 1.66 0.33 -0.84 -0.98
П 0.95 -1.90 -0.70 0.00
Р 0.41 0.50 0.75 0.82

Стр. 74: Рейтинг успеваемости:

DeviationScore <- function(values){
  return(Normalize(values) * 10 + 50)
}

test.scores %>% mutate_each(funs(DeviationScore), -Student) %>% Output(test.rus)
Ученик Английский язык Литература История Биология
Руи 57.7 47.0 58.8 53.3
Юми 49.7 63.9 53.5 60.9
А 42.6 54.1 32.9 46.7
Б 63.9 46.1 44.7 41.8
В 53.3 43.5 48.2 55.5
Г 31.0 42.6 64.9 51.6
Д 43.5 37.3 57.0 29.2
Е 34.6 57.7 55.3 50.0
Ж 41.7 34.6 42.5 64.8
З 37.3 59.5 29.8 48.4
И 55.0 66.6 50.0 43.4
К 47.0 63.0 70.7 54.9
Л 53.3 53.3 51.8 45.1
М 63.0 49.7 46.5 35.2
Н 46.1 41.7 51.3 70.8
О 66.6 53.3 41.6 40.2
П 59.5 31.0 43.0 50.0
Р 54.1 55.0 57.5 58.2

Стр. 78: Упражнение. Среднее значение и стандартное отклонение, рассчитанные по нормированным отклонениям результатов забега на 100 метров:

mean(Normalize(race.times)) %>% round(4) %>% paste("Среднее значение:", ., "\n") %>% cat()
Среднее значение: 0 
StdDev(Normalize(race.times)) %>% paste("Стандартное отклонение:", ., "\n") %>% cat()
Стандартное отклонение: 1

Глава 5: Вычислим вероятность

Стр. 93: Найдём вероятность того, что случайная величина, распределённая по стандартному нормальному закону попадёт в интервал от \(0\) до \(1.96\). В таблице стандартного нормального распределения из книги затабулирована площадь под кривой от \(0\) до \(x\), а функция pnorm в R по умолчанию использует интервал от \(-\infty\) до \(x\). Чтобы получить правильный ответ можно просто вычесть pnorm(0) из pnorm(1.96).

(pnorm(1.96) - pnorm(0)) %>% round(4)
[1] 0.475

Стр. 97: Необходимо найти площадь под стандартной нормальной кривой от \(1.8\) to \(\infty\). Чтобы использовать интервал от \(x\) до \(\infty\) нужно использовать параметр lower.tail = FALSE.

pnorm(1.8, lower.tail = FALSE) %>% round(4)
[1] 0.0359

Стр. 104: Найдём \(\chi^2\)-значение для соответсвующей вероятности \(0.05\) с \(1\) степенью свободы. Функция qchisq по умолчанию возвращает такое значение \(q\), для которого \(P(X < q)=p\). А таблица, приведённая в книге подразумевает \(P(X \geq q)=p\). Поэтому опять необходимо использовать параметр lower.tail = FALSE:

qchisq(0.05, df = 1, lower.tail = FALSE) %>% round(4)
[1] 3.8415

Стр. 108: Упражнение. Найдём площадь под стандартной нормальной кривой в интервале от \(-\infty\) до \(-0.29\):

pnorm(-0.29) %>% round(4)
[1] 0.3859

Стр. 108: Упражнение. Найдём величину \(\chi^2\) для вероятности \(P=0.05\) и \(2\) степенями свободы:

qchisq(0.05, df = 2, lower.tail = FALSE) %>% round(4)
[1] 5.9915

Глава 6: Что может связывать две переменные

Стр. 116: Результаты анкетирования о расходах на косметику и одежду из журнала P.girls:

street.survey <- data.frame(
  Respondent = c(RUS.LETTERS[1:10]),
  Makeup = c(3000,  5000, 12000,  2000,  7000, 15000,  5000,  6000,  8000, 10000),
  Clothes = c(7000,  8000, 25000,  5000, 12000, 30000, 10000, 15000, 20000, 18000)
)

street.rus <- c(Respondent = "Респондент", Makeup = "Расходы на косметику, йены", Clothes = "Расходы на одежду, йены")

street.survey %>% Output(street.rus)
Респондент Расходы на косметику, йены Расходы на одежду, йены
А 3000 7000
Б 5000 8000
В 12000 25000
Г 2000 5000
Д 7000 12000
Е 15000 30000
Ж 5000 10000
З 6000 15000
И 8000 20000
К 10000 18000

Стр. 116: Точечная диаграмма результатов анкетирования (при помощи пакета ggplot2):

ggplot(street.survey, aes(Makeup, Clothes)) + geom_point(size=3) + theme_bw() +
  ggtitle("Точечная диаграмма расходов на косметику и одежду (в месяц)") +
  xlab("Расходы на косметику (¥)") + ylab("Расходы на одежду (¥)") +
  coord_cartesian(xlim = c(0, 30000), ylim = c(0, 32000))

Стр. 118: Коэффициент линейной корреляции между расходами на косметику и расходами на одежду можно вычислить, используя функцию cor:

street.survey %$% cor(Makeup, Clothes) %>% round(4)
[1] 0.968

Стр. 121: Результаты анкетирования «Возраст и любимый бренд одежды»:

brand.survey <- data.frame(
  Respondent = c(RUS.LETTERS[1:15]),
  Age = c(27, 33, 16, 29, 32, 23, 25, 28, 22, 18, 26, 26, 15, 29, 26),
  Brand = as.factor(c("Benetton", "Zara", "O'STIN", "O'STIN", "Zara",
                      "Benetton", "Zara", "Benetton", "O'STIN", "O'STIN",
                      "Zara", "Benetton", "O'STIN", "Zara", "O'STIN"))
)

brand.rus <- c(Respondent = "Респондент", Age = "Возраст", Brand = "Бренд")

brand.survey %>% Output(brand.rus)
Респондент Возраст Бренд
А 27 Benetton
Б 33 Zara
В 16 O’STIN
Г 29 O’STIN
Д 32 Zara
Е 23 Benetton
Ж 25 Zara
З 28 Benetton
И 22 O’STIN
К 18 O’STIN
Л 26 Zara
М 26 Benetton
Н 15 O’STIN
О 29 Zara
П 26 O’STIN

Стр. 122: Точечная диаграмма по результатам опроса:

ggplot(brand.survey, aes(Brand, Age)) + geom_point(size=3) +
  theme_bw() + ggtitle("Точечная диаграмма «Любимый бренд одежды» и «Возраст»") +
  xlab("Предпочитаемый бренд одежды") + ylab("Возраст")

Стр. 123-124: Вычислим корреляционное отношение между возрастом девушек и предпочитаемым брендом одежды. В R нет встроенной функции для расчёта корреляционного отношения, так что выполним все вычисления по формуле. Здесь хорошо проявляются преимущества пакетов dplyr и magrittr, при помощи которых для вычисления внутригрупповой и межгрупповой дисперсий потребуется буквально по одной строчке кода:

# Внутригрупповая дисперсия
intraclass.var <- brand.survey %>% group_by(Brand) %>% 
                  summarise(S = sum((Age - mean(Age))^2)) %>%
                  summarise(sum(S)) %>% as.numeric()

# Средний возраст респондентов
age.mean <- brand.survey %>% summarise(mean(Age))

# Межгрупповая дисперсия
interclass.var <- brand.survey %>% group_by(Brand) %>%
                  summarise(S = n()*(mean(Age) - age.mean)^2) %>%
                  summarise(sum(S)) %>% as.numeric()

# Корреляционное отношение
cor.ratio <- interclass.var / (intraclass.var + interclass.var)

Осталось только вывести результаты вычислений:

intraclass.var %>% round(4) %>% paste("Внутригрупповая дисперсия:", ., "\n") %>% cat()
interclass.var %>% round(4) %>% paste("Межгрупповая дисперсия:", ., "\n") %>% cat()
cor.ratio %>% round(4) %>% paste("Корреляционное отношение:", ., "\n") %>% cat()
Внутригрупповая дисперсия: 224 
Межгрупповая дисперсия: 180 
Корреляционное отношение: 0.4455

Стр. 128: Таблица взаимной сопряжённости, полученная по результатам опроса 300 школьников о предпочитаемом способе признания в любви:

asked.out.table <- matrix(c(34, 61, 53, 38, 40, 74), byrow=T, nrow=2)
dimnames(asked.out.table) = list(Sex=c("Женский", "Мужской"),
  DesiredWay = c("По телефону", "По SMS", "При встрече"))

asked.out.table %>% addmargins(FUN=list(`Итого` = sum), quiet = TRUE)
По телефону По SMS При встрече Итого
Женский 34 61 53 148
Мужской 38 40 74 152
Итого 72 101 127 300

Стр. 132: Критерий согласия Пирсона (\(\chi^2_0\)), применяемый для расчёта коэффициента корреляции Крамера, можно вычислить при помощи функции chisq.test:

asked.out.chisq <- chisq.test(asked.out.table)$statistic
asked.out.chisq %>% round(4)
X-squared 
   8.0091

Стр. 133: Осталось лишь определить величину самого коэффициента корреляции Крамера:

Cramers <- function(ctable) {
    chi2 <- as.numeric(chisq.test(ctable)$statistic)
    s <- sum(ctable)
    d <- min(dim(ctable))
    return(sqrt(chi2 / (s * (d - 1))))
}

Cramers(asked.out.table) %>% round(4)
[1] 0.1634

Стр. 138: Упражнение. Вычислим коэффициент корреляции Крамера, чтобы установить взаимосвязь между «обычно заказываемой кухней» и «предпочитаемым напитком»:

drink.table <- matrix(c(43, 33, 51, 53, 29, 41), byrow=T, nrow=3)
dimnames(drink.table) = list(FoodType=c("Японская", "Европейская", "Китайская"),
  Preference = c("Кофе", "Чай"))

drink.table %>% addmargins(FUN=list(`Итого` = sum), quiet = TRUE)
Кофе Чай Итого
Японская 43 33 76
Европейская 51 53 104
Китайская 29 41 70
Итого 123 127 250

Стр. 141: Упражнение. Рассчитаем значение критерия согласия Пирсона \(\chi^2_0\):

chisq.test(drink.table)$statistic %>% round(4)
X-squared 
   3.3483

Стр. 141: Упражнение. Определим коэффициент корреляции Крамера:

Cramers(drink.table) %>% round(4)
[1] 0.1157

Глава 7: А что это за проверка гипотезы о независимости?

Попробуем провести эксперимент, описанный на стр. 152-153 и касающийся критерия согласия Пирсона (\(\chi^2_0\)) и распределения \(\chi^2\):

Стр. 153: Предполагаемая таблица взаимной сопряжённости о предпочитаемом способе признания в любви для всех школьников Японии:

asked.out.table.all <- as.table(matrix(c(400, 1600, 2000, 600, 2400, 3000),
                                   byrow=T, nrow=2))
dimnames(asked.out.table.all) = list(Sex=c("Женский", "Мужской"),
  DesiredWay = c("По телефону", "По SMS", "При встрече"))

asked.out.table.all %>% addmargins(FUN=list(`Итого` = sum), quiet = TRUE)
По телефону По SMS При встрече Итого
Женский 400 1600 2000 4000
Мужской 600 2400 3000 6000
Итого 1000 4000 5000 10000

Сначала сделаем таблицу с частотами всех возможных комбинаций заначений двух признаков:

asked.out.freq <- as.data.frame(asked.out.table.all)
asked.out.rus <- c(Sex = "Пол респондента", DesiredWay = "Способ признания в любви",
                   Freq = "Частота")

asked.out.freq %>% Output(asked.out.rus)
Пол респондента Способ признания в любви Частота
Женский По телефону 400
Мужской По телефону 600
Женский По SMS 1600
Мужской По SMS 2400
Женский При встрече 2000
Мужской При встрече 3000

Теперь, чтобы получить таблицу по всем 10000 школьникам (аналогичную таблице 7.1. на стр. 153) необходимо «повторить» все строки с соответсвующей частотой:

FreqToCases <- function(row) {
  n <- row$Freq
  df <- data.frame(
    Sex = rep(row$Sex, n),
    DesiredWay = rep(row$DesiredWay, n)
  )
  return(df)
}

# Вызываем функцию FreqToCases для размножения каждой строки
# таблицы частот asked.out.freq и "склеиваем" результаты
# в одну большую таблицу на 10000 строк
asked.out.data <- asked.out.freq %>% rowwise() %>% do(FreqToCases(.))

Теперь нужно получить 20000 различных случайных выборок по 300 школьников и рассчитать для каждой из них величину критерия согласия Пирсона \(\chi^2_0\):

chisq.data <- data.frame(Sample = 1:20000)

SampleChiSq <- function(n) {
  # случайная выборка из n школьников
  asked.out.sample <- sample_n(asked.out.data, n)
  # таблица взаимной сопряжённости для этой выборки
  asked.out.sample.table <- table(asked.out.sample)
  # возвращаем величину критерия согласия Пирсона
  return(chisq.test(asked.out.sample.table)$statistic)
}

# Векторизируем функцию, чтобы не нужно было вызывать функцию отдельно для каждой из 20000 выборок
SampleChiSq %<>% Vectorize()

# Вычисляем 20000 критериев согласия Пирсона по различным выборкам из 300 школьников
chisq.data$ChiSq <- SampleChiSq(rep(300, 20000))

# Выведем первые 10 строк полученной таблицы
chisq.rus <- c(Sample = "Выборка", ChiSq = "Критерий согласия Пирсона")
chisq.data %>% head(10) %>% Output(chisq.rus, 5)
Выборка Критерий согласия Пирсона
1 1.62330
2 2.89592
3 1.13590
4 1.06070
5 2.53846
6 0.39047
7 0.70536
8 0.60186
9 0.11794
10 0.58609

Стр. 154: Гистограмма распределения полученных 20000 значений \(\chi^2_0\) по 20000 различным случайным выборкам (величина интервала равна 0.5 для большей наглядности):

# название графика с математическами аннотациями
title <- expression(paste("Гистограмма распределения значений ", chi[0]^2, " по данным 20000 выборок"))

# рисуем график при помощи пакета ggplot2
chi.plot1 <- ggplot(chisq.data, aes(x = ChiSq)) + theme_bw() +
  geom_histogram(aes(y = ..density..), binwidth = 0.5, colour = "black", fill = "gray") +
  xlab(expression(paste(chi[0]^2, "-значение"))) + ylab("Плотность распределения") +
  ggtitle(title) + coord_cartesian(xlim = c(0, 15), ylim = c(0, 0.52))

chi.plot1

Стр. 158: Критическое значение \(\chi^2_{крит}\) для уровня значимости \(0.05\) и двумя степенями свободы:

chi.crit <- qchisq(0.05, df = 2, lower.tail = FALSE)

chi.crit
[1] 5.991465

Стр. 159: Добавим график плотности вероятности \(\chi^2\)-распределения к полученной гистограмме и покажем критическую область:

# название графика с математическами аннотациями
title <- expression(atop(paste("Гистограмма распределения значений ", chi[0]^2, " по данным 20000 выборок"),
         paste("и график ", chi^2, "-распределения с 2 степенями свободы")))

# добавляем элементы к графику chi.plot1
chi.plot2 <- chi.plot1 + ggtitle(title) +
  stat_function(fun = dchisq, args = list(df=2), colour = "blue") +
  geom_vline(aes(xintercept = chi.crit), colour = "red", linetype = "dashed") +
  annotate("text", label = "{chi[крит]^2}(alpha == 0.05)", parse = TRUE, x = chi.crit + 2,
           y = 0.07, color = "red", size = 6)

chi.plot2

Гистограмма распределения практически совпадает с графиком \(\chi^2\)-распределения, так как коэффициент корреляции Крамера равен нулю:

Cramers(asked.out.table.all) %>% cat()
0

Посмотрим, для какого числа выборок значение \(\chi^2_0\) оказалось больше или равно \(\chi^2_{крит}\):

chisq.data %>% filter(ChiSq >= chi.crit) %>% nrow()
[1] 897

Это составляет 4.5% от всех 20000 случайных выборок, т.е. около 5%.

Стр. 177: P-значение:

chisq.test(asked.out.table)$p.value %>% round(4)
[1] 0.0182

Или можно получить сразу полную сводку по критерию согласия Пирсона:

chisq.test(asked.out.table)

    Pearson's Chi-squared test

data:  asked.out.table
X-squared = 8.0091, df = 2, p-value = 0.01823

Стр. 188: Упражнение. Проверить гипотезу о независимости по данным анкетирования об обычно заказываемой кухне и предпочитаемым напитком:

chisq.test(drink.table)

    Pearson's Chi-squared test

data:  drink.table
X-squared = 3.3483, df = 2, p-value = 0.1875
qchisq(0.01, df = 2, lower.tail = FALSE) %>% round(4) %>%
  paste("Критическое χ²-значение:", ., "\n") %>% cat()
Критическое χ²-значение: 9.2103