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

Для того, чтобы поработать с результатами опроса про спасателей из мультфильма «Чип и Дейл спешат на помощь», нам нужно сначала загрузить их в R. Результаты опроса хранятся в файле chip-n-dale.csv. Расширение .csv расшифровывается как comma separated values, то есть значения, разделенные запятыми. Внешне он немного отличается от более привычного нам формата Excel, но смысл не меняется. Например, такую запись данных

name,age,gender
Anna,25,female
Nick,27,male

компьютер сможет считать как привычную табличку:

##   name age gender
## 1 Anna  25 female
## 2 Nick  27   male

С чтением файлов Excel иногда возникают проблемы, поэтому мы решили преобразовать результаты опроса в формат CSV. Загрузим файл по ссылке с помощью функции read.csv() и сохраним табличку в переменную dat, чтобы потом к ней можно было обратиться:

dat <- read.csv("https://allatambov.github.io/psms/chip-n-dale.csv")

Файлы можно загружать по названию или по ссылке, причем название или ссылка должны быть указаны в кавычках (неважно, одинарных или двойных). Если файл сохранен на компьютере, еще можно воспользоваться функцией file.choose(), после запуска которой откроется окно для выбора файла:

# для примера, 
# если запускали код выше, эту строчку запускать не нужно

dat <- read.csv(file.choose())

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

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

Переменные:

Чтобы избежать проблем при чтении файлов, содержащих текст на кириллице, имена героев были записаны на английском: Chip — Чип, Dale — Дейл, Gadget — Гайка, Monty — Рокки (Рокфор), Zipper — Вжик. Также есть вариант Other, если в вопросе был выбран ответ «Другое».

Описание данных

Посмотрим на «техническое» описание данных в dat, которое выдает нам функция str():

# str — от structure
str(dat)
## 'data.frame':    67 obs. of  11 variables:
##  $ X       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ group   : Factor w/ 6 levels "191","192","193",..: 2 2 2 4 2 2 4 3 1 4 ...
##  $ female  : int  1 1 1 0 1 1 0 0 0 1 ...
##  $ spec    : Factor w/ 2 levels "PA","PG": 1 2 2 2 1 1 2 2 2 2 ...
##  $ fav     : Factor w/ 6 levels "Chip","Dale",..: 1 3 3 2 1 3 1 2 5 3 ...
##  $ assoc   : Factor w/ 6 levels "Chip","Dale",..: 1 3 1 2 1 3 2 3 1 4 ...
##  $ cheer   : int  60 90 75 0 40 50 100 50 75 95 ...
##  $ grump   : int  100 100 90 101 90 60 50 100 60 54 ...
##  $ cheeze  : int  20 50 100 77 80 40 40 0 85 100 ...
##  $ tech    : int  60 60 95 89 70 90 50 40 70 75 ...
##  $ reaction: num  80 90 87 78.7 20 60 80 60 80 89 ...

R сообщает нам, что в dat есть 67 наблюдений (строк) и 11 переменных (столбцов). Среди них есть:

Теперь посмотрим на описательные статистики по каждому столбцу в таблице:

summary(dat)
##        X            group        female       spec        fav        assoc   
##  Min.   : 1.0   191    :16   Min.   :0.0000   PA:24   Chip  :15   Chip  :22  
##  1st Qu.:17.5   192    :17   1st Qu.:0.0000   PG:43   Dale  :10   Dale  :12  
##  Median :34.0   193    :10   Median :1.0000           Gadget:29   Gadget:16  
##  Mean   :34.0   194    :13   Mean   :0.5522           Monty : 6   Monty : 7  
##  3rd Qu.:50.5   195    : 9   3rd Qu.:1.0000           Other : 4   Other : 6  
##  Max.   :67.0   teacher: 2   Max.   :1.0000           Zipper: 3   Zipper: 4  
##                                                                              
##      cheer            grump            cheeze            tech       
##  Min.   :  0.00   Min.   :  5.00   Min.   :  0.00   Min.   :  0.00  
##  1st Qu.: 42.50   1st Qu.: 50.00   1st Qu.: 50.00   1st Qu.: 25.00  
##  Median : 70.00   Median : 70.00   Median : 77.00   Median : 60.00  
##  Mean   : 61.67   Mean   : 68.47   Mean   : 70.59   Mean   : 51.21  
##  3rd Qu.: 80.00   3rd Qu.: 90.00   3rd Qu.:100.00   3rd Qu.: 72.25  
##  Max.   :100.00   Max.   :101.00   Max.   :100.00   Max.   :100.00  
##  NA's   :1        NA's   :1        NA's   :1        NA's   :1       
##     reaction     
##  Min.   :  0.00  
##  1st Qu.: 50.00  
##  Median : 70.00  
##  Mean   : 60.54  
##  3rd Qu.: 78.53  
##  Max.   :100.00  
##  NA's   :1

Уже знакомая нам функция summary() для нечисловых показателей подсчитывает, с какой частотой встречаются те или иные значения, а для числовых выдает набор основных описательных статистик: минимум (Min.), максимум (Max.), среднее (Mean), медиану (Median), нижний и верхний квартили (1st Qu. и 3rd Qu. соответственно).

Что мы здесь видим интересного? Во-первых, мы можем заметить, что в детстве самым любимым героем была Гайечка (29 человек), в то время как сейчас большинство студентов ассоциируют себя с Чипом (22 человека). Во-вторых, бросается в глаза высокий уровень ворчливости. Если посмотреть на медиану grump, можно сделать вывод, что половина студентов оценивает свою ворчливость более, чем на 70. А если посмотреть еще и на верхний квартиль, то и вовсе выяснится, что 25% студентов оценивают свою ворчливость на 90% от Чипа и более. Самоизоляция делает свое дело.

Интересная история также и с любовью с сыру. У переменной cheeze верхний квартиль совпадает с максимальным значением 100! Это означает, что 25% студентов оценивает свою любовь к сыру на 100% от Рокфора! Но при этом в выборке также есть респонденты, которые сыр не любят совсем — минимальное значение по переменной cheeze равно 0.

Если посмотрим на уровень веселости, то обнаружим, что в среднем студенты оценивают свою веселость на 61.67% от Дейла. Медиана тоже несильно отличается от среднего, она равна 70. Следовательно, половина студентов оценивают свою веселость не более, чем на 70% от Дейла.

Нельзя не заметить, что в описании переменных, начиная с cheer, встречается загадочный NA's. Это количество пропущенных значений (NA обычно расшифровывается как Not applicable). У нас в данных есть один человек, который не отвечал на эти вопросы (они не были обязательными). Чтобы пропуски нам не мешали, давайте их уберем! Удалим с помощью функции na.omit() и запишем изменения в датафрейм dat:

dat <- na.omit(dat)

Обратите внимание: na.omit() удаляет строчки целиком, то есть убирает строчки, где есть хотя бы один пропуск, хотя бы одна пустая ячейка.

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

Давайте для описания степени веселости студентов построим гистограмму. Для этого нам нужно выбрать столбец cheer из датафрейма dat. Чтобы выбрать конкретный столбец из таблицы по его названию, нам понадобится оператор $. Например, такой код выведет на экран все значения столбца cheer:

dat$cheer
##  [1]  60  90  75   0  40  50 100  50  75  95  50   0  50  75  50  76  70  97  20
## [20]  80  75  95  60  70  70   0  80  75  80  90  40  40  65  80  70  75 100  70
## [39]  90  40  25  40  95  20  50 100  70  80  50  70  40  75  40  60  50  75  60
## [58]  50 100  30  22  90  50  20  40 100

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

hist(dat$cheer, 
     col = "red", 
     border = "yellow",
     main = "Уровень веселости",
     xlab = "Веселость (в % от Дейла)",
     ylab = "Частота")

Пояснения по коду:

Распределение показателя веселости нельзя назвать симметричным, при этом видно, что в выборке больше всего значений от 75 до 80 процентов от Дейла.

Теперь построим гистограмму для степени ворчливости. И будет эта гистограмма в цветах Чипа:

hist(dat$grump, 
     col = "burlywood1", 
     border = "firebrick",
     main = "Уровень ворчливости",
     xlab = "Ворчливость (в % от Чипа)",
     ylab = "Частота")

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

Чтобы сохранить график себе на компьютер, в окне Plot нужно выбрать ExportSave as image, в предложенном окне выбрать папку и вписать название файла.

Теперь посмотрим на любовь к сыру и построим для cheeze ящик с усами.

boxplot(dat$cheeze, 
        col = "yellow", 
        border = "sienna2")

У нас получился довольно интересный ящик с усами. Ящик, у которого верхнего «уса» нет. Это нормально, так как в нашем случае верхний квартиль и максимум совпали (мы уже обращали на это внимание). Можем сообщить, что нетипичных значений в выборке нет, при этом любовь к сыру у студентов очень велика.

Ящик с усами для степени любви к технике (в процентах от Гайечки) выглядит менее экзотично:

boxplot(dat$tech, 
        col = "plum1")

Нетипичных наблюдений нет. Также по ящику видно, что распределение не является симметричным — медиана находится не в середине ящика, а ближе к верхнему квартилю.

Доверительные интервалы для доли и для среднего

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

Выберем из таблицы только те строчки, где в столбце assoc указано значение "Chip". Для этого воспользуемся функцией subset():

chip <- subset(dat, dat$assoc == "Chip")

В функции subset() на первом месте указывается название таблицы, на втором — условие, которое необходимо проверить. В данном случае условие dat$assoc == "Chip". Здесь мы проверяем равенство каждого элемента в столбце assoc значению "Chip" (один знак = используется для присваивания, например, x = 4, а два знака = используются для проверки условия).

Строчки, которые соответствуют студентам, выбравших Чипа, сохранены в новую таблицу chip. Посчитаем число строк в ней, а также число строк в исходной табличке с помощью nrow():

Nchip <- nrow(chip)
N <- nrow(dat)

Если мы поделим Nchip на N, мы получим ту самую долю, для которой мы собирались построить доверительный интервал.

Nchip / N
## [1] 0.3333333

Теперь перейдем к доверительным интервалам. Для того, чтобы R строил доверительные интервалы, нам понадобится установить библиотеку DescTools:

install.packages("DescTools") 

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

Если библиотека не хочет устанавливаться, выдает сообщения об ошибке и пишет что-то про отсутствие доступа, попробуйте запустить RStudio от имени администратора. Для этого нужно закрыть RStudio, найти его значок в меню или на рабочем столе, кликнуть правой клавишей и выбрать Запуск от имени администратора. Когда RStudio в таком режиме запустится, прогнать строчку выше с install.packages(), должно сработать.

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

library(DescTools)

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

Для построения доверительного интервала для доли нам понадобится функция BinomCI(). Она принимает на вход число успехов и общее число испытаний, а также уровень доверия. В нашем случае число успехов — это Nchip, число тех, кто выбрал Чипа, общее число испытаний — это N. Уровень доверия примем равным 90% или 0.9.

Построим 90%-ный доверительный интервал:

ci90 <- BinomCI(Nchip, N, conf.level = 0.90)
ci90
##            est    lwr.ci    upr.ci
## [1,] 0.3333333 0.2461205 0.4336724

Что есть что? Значение est — это значение доли тех, кто ассоциирует себя с Чипом. Та самая оценка доли (отсюда est — от estimate), полученная по выборке, которую мы на семинарах обозначали \(\hat{p}\). Значение lwr.ci — это нижняя граница доверительного интервала (lwr — от lower), а upr.ci — это верхняя граница доверительного интервала (upr — от upper).

Интерпретация: с 90%-ной уверенностью мы можем утверждать, что доля студентов, которые ассоциируют себя с Чипом, лежит в интервале от 0.25 до 0.43. Если мы будем проводить аналогичное исследование на выборках одного и того же размера много раз, независимо друг от друга, 90% доверительных интервалов будут включать истинное значение доли тех, кто ассоциирует себя с Чипом.

При желании можно умножить результат на 100 и получить вместо долей проценты:

ci90 * 100
##           est   lwr.ci   upr.ci
## [1,] 33.33333 24.61205 43.36724

Теперь построим 95%-ный доверительный интервал:

ci95 <- BinomCI(Nchip, N, conf.level = 0.95)
ci95
##            est    lwr.ci    upr.ci
## [1,] 0.3333333 0.2315643 0.4534366

Интерпретация: с 95%-ной уверенностью мы можем утверждать, что доля студентов, которые ассоциируют себя с Чипом, лежит в интервале от 0.23 до 0.45. Если мы будем проводить аналогичное исследование на выборках одного и того же размера много раз, независимо друг от друга, 95% доверительных интервалов будут включать истинное значение доли любителей рано вставать.

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

Построим 95%-ный доверительный интервал для среднего значения скорости реакции (считаем, что измеряется во «вжиках»):

# для 95% уровня можно conf.level не указывать
# R его выставляет по умолчанию

mci95 <- MeanCI(dat$reaction)
mci95
##     mean   lwr.ci   upr.ci 
## 60.53636 53.97463 67.09810

С 95%-ной уверенностью можно утверждать, что среднее значение скорости реакции лежит в интервале от 53.98 до 67.10.

Давайте построим 99%-ный интервал для среднего значения скорости реакции:

mci99 <- MeanCI(dat$reaction, conf.level = 0.99)
mci99
##     mean   lwr.ci   upr.ci 
## 60.53636 51.81776 69.25496

А теперь сравним их длину 95%-ного и 99%-ного доверительного интервала. Чтобы найти длину доверительного интервала, нам нужно из верхней границы вычесть нижнюю. В mci95 верхняя граница — это третий элемент, а нижняя граница — это второй элемент. Выберем их, выполним вычитание и проделаем то же для 99%-ного интервала:

mci95[3] - mci95[2]
##   upr.ci 
## 13.12347
mci99[3] - mci99[2]
##  upr.ci 
## 17.4372

Как и можно было ожидать, 99%-ный доверительный интервал длинее 95%-ного доверительного интервала.