А. Б. Шипунов, Е. М. Балдин, П. А. Волкова, А. И. Коробейников, С. А. Назарова

• Название: А. Б. Шипунов, Е. М. Балдин, П. А. Волкова, А. И. Коробейников, С. А. Назарова
• Анкор: rbook
• Дата: 29.09.2022
• Размер: 3.04 Mb.
• Имя файла: rbook.pdf
• Тип: Документы; #705644
• страница: 14 из 19

1. Какой тип у данных, каким способом (способами) они представ- лены?
2. Однородны ли данные — в каких единицах измерены показатели?

3. Параметрические или непараметрические у нас данные — можно ли предположить нормальное распределение данных?
4. Нужна ли «чистка» данных — есть ли пропущенные данные, вы- бросы, опечатки?
Если в таблице есть, кроме цифр, буквы, то нужно тщательно про- следить, нет ли опечаток. Иногда опечатки бывают почти невидимыми,
скажем, русская буква «с» и английская «c» по виду неразличимы, а вот программа обработки данных посчитает их разными. К счастью,
базовые команды read.table() и summary() позволяют «выловить» по- давляющее большинство подобных проблем. Кстати, если при попытке загрузки данных в R возникает ошибка, не спешите винить программу:
б
0
ольшая часть таких ошибок — это результат неправильного оформле- ния и/или ввода данных.
8.2. Окончательная обработка данных
Для того чтобы помочь читателю разобраться в многочисленных способах анализа, мы составили следующую таблицу (табл. 8.1). Она

Отчет
193
устроена очень просто: надо ответить всего лишь на три вопроса, и спо- соб обработки найдется в соответствующей ячейке. Главный вопрос —
это тип данных (см. соответствующую главу). Если данные количест- венные, то подойдет верхняя половина таблицы, если качественные или порядковые — нижняя половина. Затем надо понять, можно ли посчи- тать ваши данные распределенными нормально, то есть параметриче- скими, или по крайней мере без опаски закрыть глаза на некоторое уклонение от нормальности. И наконец, нужно понять, зависят ли раз- ные колонки данных друг от друга, то есть ст
0
оит или не ст
0
оит приме- нять парные методы сравнения (см. главу о двумерных данных).
Найдя в таблице нужный метод, загляните в указатель и перейдите на страницу с более подробным описанием функции. А можно просто ввести команду help().
8.3. Отчет
Любое полноценное исследование оканчивается отчетом — презен- тацией, статьей или публикацией в Интернете. В R есть множество ав- томатизированных средств подготовки отчетов.
Таблицы, созданные в R, можно сохранить в форматах L
A
TEX или
HTML при помощи пакета xtable. Естественно, хочется пойти дальше и сохранять в каком-нибудь из этих форматов вообще всю R-сессию.
Для HTML такое возможно, если использовать пакет R2HTML:
> library(R2HTML)
> dir.create("example")
> HTMLStart("example")
> 2+2
> plot(1:20)
> HTMLplot()
> HTMLStop()
В рабочей директории будет создана поддиректория example, и туда будут записаны HTML-файлы, содержащие полный отчет о текущей сессии, в том числе и созданный график.
Можно пойти и еще дальше. Что, если создать файл, который будет содержать код R, перемешанный с текстовыми комментариями, и потом
«скормить» этот файл R, так чтобы фрагменты кода заменились на результат их исполнения? Идея эта называется «literate programming»
(грамотное программирование) и принадлежит Дональду Кнуту, созда- телю TEX. В случае R такая система используется для автоматической генерации отчетов — особенности, которая делает R поистине незаме-

194
Статистическая разведка
Данные
Од- на груп- па
Две группы:
различия
Две группы:
связи
Три и более групп:
связи
Три и более групп:
общая картина
Количественные
Параметрические
Независимые t.test()
oneway.
test(),
paiwise.
t.test(),
anova(),
lm()
lda(),
mano- va()
Зависимые t.test(...,
paired =
TRUE)
cor.test(...,
metho d="p e")
–
–
Непараметрические
Независимые wilcox.
test()
kruskal.
test()
pca(),
tree(),
cor(),
hclust(),
isoMDS(),
cmd- scale()
Зависимые summary()
wilcox.test
(...,
paired =
TRUE)
cor.test(...,
metho d="sp")
–
–
Номинальные или шк альные
Непараметрические
Независимые chisq.
test(),
prop.
test(),
binom.
test()
glm(...,
"bino- mial")
–
cor(),
dist(),
hclust(),
isoMDS(),
corresp()
Зависимые mcne- mar.
test()
–
–
–
Таблица 8.1. Варианты статистического анализа и соответствующие функции R

Отчет
195
нимым. Для создания подобного отчета надо вначале набрать простой
L
A
TEX-файл и назвать его, например, test-Sweave.Rnw:
\documentclass[a4paper,12pt]{article}
\usepackage[T2A]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[english,russian]{babel}
\usepackage[noae]{Sweave}
\begin{document} % Тут начинается отчет
\textsf{R} как калькулятор:
<>=
1 + 1 1 + pi sin(pi/2)
@
Картинка:
<>=
plot(1:20)
@
\end{document}
Затем этот файл необходимо обработать в R:
> Sweave("test-Sweave.Rnw")
Writing to file test-Sweave.tex
Processing code chunks ...
1 : echo print term verbatim
2 : echo term verbatim eps pdf
You can now run LaTeX on ’test-Sweave.tex’
При этом создается готовый L
A
TEX-файл test-Sweave.tex. И нако- нец, при помощи L
A
TEX и dvips или pdfL
A
TEX можно получить резуль- тирующий PDF-файл, который изображен на рис. 61.
Такой отчет можно расширять, шлифовать, изменять исходные дан- ные, и при этом усилия по оформлению сводятся к минимуму.

196
Статистическая разведка
R как калькулятор:
> 1 + 1
[1] 2
> 1 + pi
[1] 4.141593
> sin(pi/2)
[1] 1
Картинка:
> plot(1:20)
5 10 15 20 5
10 15 20
Index
1:20 1
Рис. 61. Пример отчета, полученного с помощью команды Sweave()

Отчет
197
Есть и другие системы генерации отчетов, например пакет brew, ко- торый позволяет создавать автоматические отчеты в текстовой форме
(разумеется, без графиков), и пакет odfWeave, который может работать с ODF (формат OpenOffice.org).
И все-таки полностью полагаться на автоматику никогда не ст
0
оит.
Обязательно проверьте, например, полученные графики — достаточен ли размер графических файлов, нет ли проблем с цветами. И не за- будьте запомнить в файл историю ваших команд. Это очень поможет,
если придется опять вернуться к обработке этих данных. А если вы захотите сослаться в своем отчете на тот замечательный инструмент,
который помог вам обработать данные, не забудьте вставить в список литературы то, что выводит строка
> citation()
Удачного вам анализа данных!

Приложение А
Пример работы в R
Это приложение написано для тех, кто не хочет читать длинные тексты, а хочет как можно скорее научиться обрабатывать данные в
R
. Для этого надо «просто» последовательно выполнить приведенные ниже команды. Желательно не копировать их откуда-либо, а набрать с клавиатуры: таким способом запомнить и, стало быть, научиться будет гораздо проще. Хорошо, если про каждую выполненную команду вы прочтете соответствующий раздел справки (напоминаем, это делается командой help(команда)). Мы не стали показывать здесь то, что вы- водит R, и не приводим получающиеся графики: все это вам нужно будет получить и проверить самостоятельно
Все команды будут относиться к файлу данных о воображаемых жуках, состоящему из четырех столбцов (признаков): пол жука (POL:
самки = 0 и самцы = 1), цвет жука: (CVET: красный=1, синий=2, зеле- ный=3), вес жука в граммах (VES) и длина жука в миллиметрах (ROST):
POL CVET
VES
ROST
0 1
10.68 9.43 1
1 10.02 10.66 0
2 10.18 10.41 1
1 8.01 9
0 3
10.23 8.98 1
3 9.7 9.71 1
2 9.73 9.09 0
3 11.22 9.23 1
1 9.19 8.97 1
2 11.45 10.34
Начинаем...
Создаем на жестком диске рабочую директорию, создаем в ней ди- ректорию data; копируем в последнюю файл данных в текстовом фор- мате с расширением *.txt и разделителем-табуляцией (делается из фай- ла Excel или другой подобной прграммы командой меню Save as... /
Сохранить как...
). Пусть ваш файл называется zhuki.txt.

Пример работы в R
199
Внимание!
Проверьте, чтобы десятичный разделитель был «.», то есть точкой, а не запятой: половина — это 0.5, а не 0,5! Если десятич- ный разделитель — запятая, то можно заменить его на точку в любом текстовом редакторе, например в Блокноте/Notepad.
Открываем программу R. Указываем директорию, где находится ваш файл данных, при помощи меню: Файл -> Изменить папку -> выби- раем директорию
, или печатаем команду setwd(...), в аргументе кото- рой должен стоять полный путь к вашей директории (для указания пути надо использовать прямые слэши — «/».
Для проверки местоположения печатаем команду
> dir("data")
...и нажимаем Enter (эту клавишу нужно нажимать каждый раз после ввода команды).
Эта команда должна вывести, среди прочего, имя вашего файла
(zhuki.txt).
Читаем файл данных (создаем в памяти программы объект под на- званием data, который представляет собой копию вашего файла дан- ных). В строке ввода набираем:
> data <- read.table("data/zhuki.txt", h=TRUE)
Внимание!
Будьте осторожны со скобками, знаками препинания и регистром букв! Если десятичный разделитель — запятая и вы почему- либо не смогли заменить его на точку, то есть способ загрузить и такие данные:
> data.2 <- read.table("data/zhuki_zap.txt", h=TRUE, dec=",")
Кстати, похожую таблицу данных можно создать при помощи R сле- дующими командами:
> data.3 <- data.frame(POL=sample(0:1, 100, replace=TRUE),
+ CVET=sample(1:3, 100, replace=TRUE),
+ VES=sample(20:120, 100, replace=TRUE),
+ ROST=sample(50:150, 100, replace=TRUE))
Первая команда создает данные (у нас указана повторность 100)
методом случайной выборки, а вторая записывает их в файл с задан- ным названием. Цифры и графики будут немного другие (выборка-то случайная!), но на наши упражнения это повлиять не должно.
Посмотрим на ваш файл данных (если файл большой, то можно использовать специальную команду head(data)):

200
Пример работы в R
> data
Внимание!
Внутри R вносить изменения в данные не очень удобно.
Разумно вносить их в текстовый файл данных (открыв его, например,
в Excel), а потом заново читать его в R.
Посмотрим на структуру файла данных. Сколько объектов (obs. =
observations), сколько признаков (variables), как названы признаки и в каком порядке они следуют в таблице:
> str(data)
Отметьте для себя, что POL и CVET загружены как числа, в то время как на самом деле это номинальные данные.
Создадим в памяти еще один объект с данными, куда отберем дан- ные только для самок (POL = 0):
> data.f <- data[data$POL == 0,]
А теперь — отдельный объект с данными для крупных (больше
10 мм) самцов:
> data.m.big <- data[data$POL == 1 & data$ROST > 10,]
Кстати, эту команду проще не вводить заново, а получить путем редактирования предыдущей команды (обычное дело в R). Для вызова предыдущей команды воспользуйтесь «стрелкой вверх» на клавиатуре.
Использованные знаки «==» и «&» — это логические выражения «та- ких, что» и «и». Именно они служат критериями отбора. Кроме того,
для отбора данных обязательно нужны квадратные скобки, а если дан- ные табличные (как у нас), то внутри квадратных скобок обязательна запятая, разделяющая выражения для строк и столбцов.
Добавим еще один признак к нашему файлу: удельный вес жука
(отношение веса жука к его длине) — VES.R:
> data$VES.R <- data$VES/data$ROST
Проверьте, что новый признак появился, при помощи уже исполь- зованной нами команды str(data) (стрелка вверх!).
Новый признак был добавлен только к копии вашего файла данных,
который находится в памяти программы. Чтобы сохранить измененный файл данных под именем zhuki_new.txt в вашей директории на жест- ком диске компьютера, нужно написать:
> write.table(data, "data/zhuki_new.txt", quote=FALSE)

Пример работы в R
201
Охарактеризуем выборку...
Посмотрим сначала на основные характеристики каждого признака
(не имеет смысла для категориальных данных):
> summary(data)
Конечно, команду summary(), как и многие прочие, можно приме- нять как к всему файлу данных, так и к любому отдельному признаку:
> summary(data$VES)
А можно вычислять эти характеристики по отдельности. Минимум и максимум:
> min(data$VES)
> max(data$VES)
... медиана:
> median(data$VES)
... среднее арифметическое только для веса и для каждого из при- знаков:
> mean(data$VES)
и
> colMeans(data)
соответственно.
К сожалению, предыдущие команды не работают, если есть пропу- щенные значения. Для того чтобы посчитать среднее для каждого из признаков, избавившись от пропущенных значений, надо ввести
> mean(data, na.rm=TRUE)
Кстати, строки с пропущенными значениями можно удалить из таб- лицы данных так:
> data.o <- na.omit(data)
Иногда бывает нужно вычислить сумму всех значений признака:
> sum(data$VES)

202
Пример работы в R
... или сумму всех значений одной строки (попробуем на примере второй):
> sum(data[2,])
... или сумму значений всех признаков для каждой строки:
> apply(data, 1, sum)
Для номинальных признаков имеет смысл посмотреть, сколько раз встречается в выборке каждое значение признака (заодно узнаем, какие значения признак принимает):
> table(data$POL)
А теперь выразим частоту встречаемости значений признака не в числе объектов, а в процентах, приняв за 100% общее число объектов:
> 100*table(data$POL)/length(data$POL)
И еще округлим значения процентов до целых чисел:
> round(100*table(data$POL)/length(data$POL), 0)
Одна из основных характеристик разброса данных вокруг средне- го значения (наряду с абсолютным и межквартильным разбросом) —
стандартное отклонение (standard deviation). Вычислим его:
> sd(data$VES)
Вычислим напоследок и безразмерный коэффициент вариации (CV):
> 100*sd(data$VES)/mean(data$VES)
Можно вычислять характеристики любого признака отдельно для самцов и для самок. Попробуем на примере среднего арифметического для веса:
> tapply(data$VES, data$POL, mean)
Посмотрим, сколько жуков разного цвета среди самцов и самок:
> table(data$CVET, data$POL)

Пример работы в R
203
(Строки — разные цвета, столбцы — самцы и самки.)
А теперь то же самое, но не в штуках, а в процентах от общего числа жуков:
> 100*table(data$CVET, data$POL)/sum(data$CVET, data$POL)
И наконец, вычислим средние значения веса жуков отдельно для всех комбинаций цвета и пола (для красных самцов, красных самок,
зеленых самцов, зеленых самок...):
> tapply(data$VES, list(data$POL, data$CVET), mean)
Теперь порисуем диаграммы...
Проверим сначала, как распределены данные, нет ли выбросов. Для этого построим гистограммы для каждого признака:
> hist(data$VES, breaks=20)
Детализацию гистограммы можно изменять, варьируя число интер- валов (breaks).
Можно задать точную ширину интервалов значений признака на гистограмме (зададим ширину в 20 единиц, а значения признака пусть изменяются от 0 до 100):
> hist(data$VES, breaks=c(seq(0,100,20)))
Более точно проверить нормальность распределения признака мож- но при помощи двух команд:
> qqnorm(data$VES); qqline(data$VES)
Обе команды «делают» один график (поэтому мы написали их в одну строчку через точку с запятой). Чем больше распределение точек на этом графике отклоняется от прямой линии, тем дальше распреде- ление данных от нормального. Кстати, для того чтобы открыть новое графическое окно (новый график будет нарисован рядом со старым, а не вместо него), можно использовать команду dev.new().
Построим диаграмму рассеяния, на которой объекты будут обозна- чены кружочками. Рост будет по оси абсцисс (горизонтальная ось),
вес — по оси ординат (вертикальная ось):
> plot(data$ROST, data$VES, type="p")
Можно изменять размер кружочков (параметр cex). Сравните:

204
Пример работы в R
> plot(data$ROST, data$VES, type="p", cex=0.5)
и
> plot(data$ROST, data$VES, type="p", cex=2)
Можно изменить вид значка, который обозначает объект (см. номе- ра значков на рис. 62)). Похожую таблицу можно вызвать при помощи команды example(points) и нажать несколько раз Enter, пока «демон- страция» не окончится.
0 1
2 3
4 5
6 7
8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
*
*
+
+
Рис. 62. Номера значков, используемых в стандартных графиках R
Вот, например, обозначим объекты значком 2-го типа (пустым тре- угольником):
> plot(data$ROST, data$VES, type="p", pch=2)
Можно вместо значков обозначить объекты на диаграмме кодом по- ла (0/1):

Пример работы в R
205
> plot(data$ROST, data$VES, type="n")
> text(data$ROST, data$VES, labels=data$POL)
Обе команды здесь действуют на один и тот же график. Первая печатает пустое поле, вторая добавляет туда значки.
Еще можно сделать так, чтобы разные цифры-обозначения имели и разные цвета (мы добавили «+1», иначе бы значки для самок печата- лись нулевым, прозрачным цветом):
> plot(data$ROST, data$VES, type="n")
> text(data$ROST, data$VES, labels=data$POL, col=data$POL+1)
А можно самцов и самок обозначить разными значками:
> plot(data$ROST, data$VES, type="n")
> points(data$ROST, data$VES, pch=data$POL)
Другой, более сложный вариант — с использованием встроенных в
R
шрифтов Hershey (в порядке исключения приводим его на рис. 63):
> plot(data$ROST, data$VES, type="n", xlab="Рост", ylab="Вес")
> text(data$ROST, data$VES,
+ labels=ifelse(data$POL, "\\MA", "\\VE"),
+ vfont=c("serif","plain"), cex=1.5)
... и еще разными цветами:
> plot(data$ROST, data$VES, type="n")
> text(data$ROST, data$VES, pch=data$POL, col=data$POL+1)
Наконец, в случае с разноцветными значками нужно добавить услов- ные обозначения (легенду):
> legend(50, 100, c("male", "female"), pch=c(0,1), col=c(1,2))
Числами указано положение (координаты на графике) верхнего ле- вого угла легенды: первая цифра (50) — абсцисса, вторая цифра (100) —
ордината.
Сохраняем график при помощи меню: Файл -> Сохранить как...
-> PNG -> graph.png или печатаем две команды:
> dev.copy(png, filename="graph.png")
> dev.off()

206
Пример работы в R
9.0 9.5 10.0 10.5 8.
0 8.
5 9.
0 9.
5 10
.0 10
.5 11
.0 11
.5
Рост
В
ес
Рис. 63. Распределение самцов и самок жуков по росту и весу (для символов использованы шрифты Hershey)
Не забудьте напечатать dev.off()!
Чтобы получить график без посторонних надписей, следует дать команду plot(..., main="", xlab="", ylab=""), где многоточие обо- значает все остальные возможные аргументы.
Сохранять график можно, и не выводя его на экран, сразу в графи- ческий файл, опять-таки с помощью двух дополнительных команд:
> png("data.png")
# [... тут печатаем команды для построения графика ...]
> dev.off()
Рисуем коррелограмму:
> plot(data[order(data$ROST), c("ROST", "VES")], type="o")
Здесь мы отсортировали жуков по росту, потому что иначе вместо графика получилось бы множество пересекающихся линий.
А теперь нарисуем две линии на одном графике:

Пример работы в R
207
> plot(data[order(data$CVET), c("CVET", "VES")], type="o",
+ ylim=c(5, 15))
> lines(data[order(data$CVET), c("CVET", "ROST")], lty=3)
Аргумент ylim задает длину оси ординат (от 0 до 50). Аргумент lty задает тип линии («3» — это пунктир).
Рисуем «ящик-с-усами», или, по-другому, боксплот (он покажет вы- бросы, минимум–максимум, квартильный разброс и медиану):
> boxplot(data$ROST)
... а теперь — для самцов и самок по отдельности:
> boxplot(data$ROST

factor(data$POL))
Статистические тесты
Достоверность различий для параметрических данных (тест Стью- дента), для зависимых переменных:
> t.test(data$VES, data$ROST, paired=TRUE)
... и для независимых переменных:
> t.test(data$VES, data$ROST, paired=FALSE)
... если нужно сравнить значения одного признака для двух групп:
> t.test(data$VES data$POL)
Если p-value < 0.05, то различие между выборками достоверно.
В R по умолчанию не требуется проверять, одинаков ли разброс данных относительно среднего.
Достоверность различий для непараметрических данных (тест Вил- коксона):
> wilcox.test(data$VES, data$ROST, paired=TRUE)
Достоверность различий между тремя и более выборками парамет- рических данных (вариант однофакторного дисперсионного анализа):
> oneway.test(data$VES data$CVET)
Посмотрим, какие именно пары выборок достоверно различаются:
> pairwise.t.test(data$VES, data$CVET, p.adj="bonferroni")

208
Пример работы в R
А теперь проверим достоверность различий между несколькими вы- борками непараметрических данных:
> kruskal.test(data$VES data$CVET)
Достоверность соответствия для категориальных данных (тест Пир- сона, хи-квадрат):
> chisq.test(data$CVET, data$POL)
Достоверность различия пропорций (тест пропорций):
> prop.test(c(sum(data$POL)), c(length(data$POL)), 0.5)