Документ взят из кэша поисковой машины. Адрес оригинального документа : http://herba.msu.ru/shipunov/software/r/r-intro-ru.txt Дата изменения: Wed Dec 5 15:08:40 2007 Дата индексирования: Tue Oct 2 02:49:17 2012 Кодировка: Windows-1251 Поисковые слова: eruption

Введение в R

Это представление R ( "GNU S"), языка и среды для статистических
вычислений и графики. R близка к популярной системе S, которая
была разработана в Bell Laboratories Джоном Чамберсом и др.. Она
предусматривает широкий спектр статистических и графических методов
(линейное и нелинейное моделирование, статистические тесты, анализ
временных рядов, классификация, кластеризация, ...).

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

Текущая версия этого документа - 2.6.0 (2007-10-03).

ISBN 3-900051-12-7

Предисловие

Это введение в R возникло из первоначального сборника заметок описывающих
среды S и S-Plus написанных Венаблес Билл и Дэвид М. Смит (Инсайтфул
корпорация). Мы внесли ряд небольших изменений, чтобы отразить различия
между R и S кодом, и расширить некоторые материалы.

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

Замечания и исправления всегда приветствуются. Пожалуйста адресуйте
электронную почту на R-core@R-project.org.

Примечания для Читателя

Большинству R новичков следует начать с вводной сессии в Приложении
А. Это должно познакомить со стилем R сессий и, что более важно дать
некоторое впечатление о том, что действительно происходит.

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

1.1 R среда

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

* эффективная обработка и хранение данных

* набор операторов для обработки массивов, в частности матриц,

* цельная, непротиворечивая, комплексная коллекция утилит для анализа
данных,

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

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

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

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

1.2 Сопутствующее программное обеспечение и дополнительная документация

R можно рассматривать как реализацию языка S, который был разработан
в Bell Laboratories Риком Бекером, Джоном Чамбером и Алланом Вилксом,
и который собственно лежит в основе S-Plus систем.

Эволюция S языке характеризуется в четырех книгах Джона Чамбера с
соавторами. Для R, основой является "Новый S Язык: Программная среда
для анализа данных и визуализации" Ричарда А. Беккера, Джон М. Чамбера и
Аллана Р. Вилкса. Новые возможности реализации S в 1991 году освещаются
в "Статистические модели в S" под редакцией Джон М. Чамбера и Тревора
Дж. Хастие. Формальные методы и классы из пакета методов основаны на
описаниях в "Кодирование с данными" Джона М. Чамбера. См. раздел Ссылки
для точных обращений к источникам.

Сейчас есть целый ряд книг, которые описывают, как использовать R для
анализа данных и статистических вычислений, также документация S / S-Plus,
как правило, может быть использована с R, если помнить о различиях между
реализациями S. См. "Какая существует документация по R"? .

1.3 R и статистика

В нашем описании R среды не упоминается о статистике, но многие люди
используют R в качестве статистической системы. Мы предпочитаем думать
о нем как о окружении, в котором были реализованы многие классические и
современные статистические методы. Некоторые из них встроены в базовое
окружение R, но многие из них поставляется как пакеты. Есть около 25
пакетов поставляемых с R (так называемые "стандарт" и "рекомендовано"
наборы пакетов), много больше можно получить через CRAN семейство интернет
сайтов (через http://CRAN.R-project.org) и из других источников. Более
подробную информацию о пакетах рассмотрим позднее (см. Пакеты).

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

Существует важное различие между философией в S (и, следовательно, R)
и других основных статистических системах. В S статистический анализ,
как правило, выполняется в виде серии шагов, с сохранением промежуточных
результатов в объектах. Так когда SAS и SPSS выведут обильные результаты
процедуры регрессионного или дискриминантного анализа, R выведет минимум
результатов и сохранит вывод в подогнанный объект для последующего
использования в дальнейшем в вызываемых R функциях.

1.4 R и графический десктоп

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

Большинство пользователей сталкиваются время от времени с необходимостью
взаимодействовать непосредственно с операционной системы на своем
компьютере. В этом руководстве, мы в основном обсуждаем взаимодействие
с операционной системой UNIX машин. Если вы запускаете R под Windows
или MacOS нужно будет сделать несколько небольших корректировок.

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

1.5 Используем R интерактивно

Когда Вы используете R программу, она выводит приглашение в момент
ожидания ввода команд. По умолчанию приглашение это ">", что в UNIX
может совпадать с приглашением оболочки, и поэтому может показаться,
что ничего не произошло. Однако, как мы увидим, можно, при желании,
легко выбрать другое R приглашение. В дальнейшем мы предполагаем, что
приглашение оболочки UNIX является "$".

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

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

$ mkdir work $ cd work


2. Запустить R программу с помощью команды

$ R


3. С этого момента можно вводить R команды (см. ниже). 4. Выходом из R
программы служит команда

> q()


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

Последующие R сессии требуют меньше действий.

1. Сделать work рабочим каталогом и запустить программу, как делали ранее:

$ cd work $ R


2. Использовать R программу, которую завершить при помощи q() команды
в конце сессии.

Для использования под Windows R процедура в основе своей та же.
Создайте папку как рабочий каталог, и установите ее в поле Start Вашего
R ярлыка. Затем запустите R двойным нажатием на иконку.

1.6 Первый сеанс

Читателям, желающим испытать R на компьютере, прежде чем приступить
настоятельно рекомендуется проработать вводную сессию приведенную в
сессии примера А.

1.7 Получение помощи по функциям и возможностям

R имеет встроенную команду help по аналогии с командой man в UNIX.
Чтобы получить дополнительную информацию о какой-либо конкретной функции,
например solve, нужно ввести

> help(solve)

Альтернативный вариант

>?solve

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

> help("[[")

Любая форма кавычек может быть использована, чтобы экранировать другие
кавычки, как в строке "It's important". Мы будем придерживаться по
возможности использования двойных кавычек.

На большинстве установленных R систем справка доступна в формате HTML,
достаточно выполнить

> help.start()

будет запущен веб-браузер, что позволит просматривать страницы помощи в
режиме гипертекста. В UNIX, последующие help запросы будут направляться в
HTML вариант системы помощи. Ссылка `Поисковая система и ключевые слова"
доступная на странице загружаемой после выполнения help.start () является
особенно полезной, поскольку содержит высокоуровневое оглавление, для
поиска доступных функций. Это может быть очень удобно, чтобы быстро
получить понимание вопроса, и понять широкий спектр возможностей
предлагаемых R.

help.search команда позволяет искать подсказку различными способами:
выполните ?help.search для подробностей и конкретных примеров.

Примеры по теме страницы помощи можно, как правило, запустить коммандой

> example(тема)

Windows версии R имеют другие дополнительные системы помощи:
воспользуйтесь

>?help

для дальнейших подробностей.

1.8 Команды R, учет регистра и т.п.

Технически R является языком выражений с очень простым синтаксисом.
Он учитывает регистр как и большинство других программ UNIX, так
"А" и "а" различные символы и будут обозначать различные переменные.
Набор символов, которые могут быть использованы как имена в R зависит
от операционной системы и страны, в которых будет запущен R (технически
говоря, от используемой локали). Обычно разрешены все алфавитно-цифровые
символы сноска1 (а в некоторых странах включая и акцентированные буквы)
вместе с "." и "_", с ограничением, что имя должно начинаться с "." или
буквы, а если оно начинается с "." второй символ не должен быть цифрой.

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

Команды разделяются либо точкой с запятой (";"), или переводом
строки. Простые команды могут быть сгруппированы в единое составное
выражение фигурными скобками ( "{" и "}"). Комментарии могут быть
практически сноска2 где угодно, начинаются с символа решетки ("#"),
при этом все до конца строки является комментарием.

Если команда не завершена в конце строки, R выдаст особое приглашение,
по умолчанию

+

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

Командные строки набираемые в консоли являются ограниченными сноска3 в
размере до 1024 байт (не символов!).

1.9 Повтор и коррекция предыдущих команд

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

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

Кроме того, текстовый редактор Emacs предоставляет более полные механизмы
поддержки (посредством ESS, Emacs Speaks Statistics) для работы в
интерактивном режиме с R. См. R и Emacs.

1.10 Выполнение команд записанных в файл или перенаправление вывода в файл

Если команды сноска4 хранятся во внешнем файле, например commands.R
в рабочем каталоге work, они могут быть выполнены в любой момент в R
сессии с помощью команды

> source("commands.R")

Для Windows Source можно найти в меню "Файл". Функция sink,

> sink("record.lis")

будет направлять весь последующий вывод консоли во внешний файл,
record.lis. Команда

> sink()

восстанавливает вывод в консоли еще раз.

1.11 Сохранение данных и удаление объектов

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

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

> objects()

(также как ls()) может быть использована для отображения названий (в
основном) объектов, которые в настоящее время хранятся в R. Набор объектов
который в настоящее время хранится называется рабочее пространство. Чтобы
удалить объекты доступна функция rm:

> rm(x, y, z, ink, junk, temp, foo, bar)

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

Если R будет запущена позже из этого каталога рабочее пространство
будет перезагружено из этого файла. Одновременно загрузится связанная
история команд.

Мы рекомендуем Вам использовать отдельную рабочую директорию для каждого
из анализов проведенных в R. Очень распространено использовать для
объектов созданых в ходе анализа имена х и у. Имена, подобные этим имеют
смысл в контексте конкретного анализа, но может быть довольно трудно
определить, что они означают, когда несколько анализов были проведены
в одном и том же каталоге.

2.1 Вектора и присваивания

R оперирует именованными структурами данных. Простейшая такая структура
это численный вектор, который представляет собой совокупность, состоящую
из упорядоченного набора чисел. Чтобы создать вектор с именем х, например
состоящий из пяти чисел, а именно 10.4, 5.6, 3.1, 6.4 и 21.7, используют
такую команду R

> х <- с(10.4, 5.6, 3.1, 6.4, 21.7)

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

Одиночное число входящее в выражение трактуется как вектор единичной
длины.

Учтите, что оператор присваивания ("<-"), который состоит из
двух символов "<" ( "меньше") и "-" ( "минус") выполняется строго
однонаправлено и "указывает" на объект получающий значение выражения.
В большинстве случаев в качестве альтернативы может быть использован
оператор "=". Присваивание можно также сделать с помощью функции
assign(). Эквивалентный способ присвоения, приведенному выше, выглядит:

> assign("x", c(10.4, 5.6, 3.1, 6.4, 21.7))

Обычно оператор, <-, можно рассматривать как синтаксическое сокращение
приведенного.

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

> с(10.4, 5.6, 3.1, 6.4, 21.7) -> х

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

> 1/х

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

Дальнейшее присваивание

> y <- c(x, 0, x)

создаст вектор y с 11 элементами, состоящий из двух копий х с нулем
между ними.

2.2 Действия над векторами

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

> v <- 2*х + у + 1

создает новый вектор v длинной 11 составленный путем сложения, элемент
с элементом, 2*х повторенного 2,2 раза, y повторенного только раз,
и 1 повторенной 11 раз.

Элементарные арифметические операторы это +, -, *, / и ^ для возведения
в степень. Дополнительно присутсвуют все простые арифметические
функции. log, exp, sin, cos, tan, sqrt, и т.д., все имеют обычное
значение. max и min выделяют соответственно наибольший и наименьший
из элементов вектора. range это функция, значение которой вектор из
двух элеменов, а именно c(min(x), max(x)). length(x) - количество
элементов в х, sum(x) дает сумму значений элементов в х, а prod(x) их
произведение. Две статистические функции это mean(x), которая вычисляет
выборочное среднее, что соответсвует sum(x)/length(x), и var(х),
которая дает

sum((x-mean(x))^2)/(length(x)-1)

или выборочную дисперсию. Если аргумент var() - это n-на-р матрица,
значение - р-на-р матрица выборочной ковариации, полученная путем
интерпретации строк как p независимых векторов-выборок.

sort(x) возвращает вектор того же размера как х с элементами расположеными
в порядке возрастания; однако есть и другие более гибкие возможности
сортировки (см. order() или sort.list(), которые для сортировки
производят перестановки). Заметим, что max и min выделяют наибольшей и
наименьшей значение в своих аргументах, даже если они получат несколько
векторов. Параллельные максимальные и минимальные функции pmax и pmin
возвращают вектор (длиной, равной длине длиннейшего аргумента), который
содержит в каждом своем элементе наибольшее (наименьшее) значение на
этой позиции во всех входных векторах. В большинстве случаев пользователю
не важно, являются ли "числа" в числовых векторах целыми, реальными или
даже комплексными. Внутренние расчеты осуществляются в реальных числах
двойной точности, или комплексных числах двойной точности, если входные
данные являются комплексными.

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

sqrt(-17)

даст NaN и предупреждение, но

sqrt(-17+0i)

будет обработано, как комплексное число.

2,3 Генерация регулярных последовательностей

R имеет ряд возможностей для получения общеупотребительных
последовательностей чисел. Например 1:30 - вектор с(1, 2, ..., 29,
30). Оператор двоеточие имеет высший приоритет в выражении, так, например
2*1:15 это вектор с(2, 4, ..., 28, 30). Присвойте n <- 10 и сравните
последовательности 1:n-1 и 1:(n-1).

Выражение 30:1 может быть использовано для создания обратной
последовательности.

Функция seq() представляет более общий механизм для генерации
последовательностей. Она имеет пять аргументов, лишь часть из которых
может быть определена одновременно. Первые два аргумента, если они заданы,
определяют начало и конец последовательности, и если указаны только два
аргумента возвращаемый результат такой же, как у оператора двоеточие. Так
seq(2,10) - это тот же вектор, как и 2:10.

Параметры seq(), а также для многих других функций R, также можно
задавать по имени, тогда порядок, в котором они появляются, не имеет
значения. Первые два параметра могут быть обозначены from=значение и
to=значение; таким образом, seq(1,30), seq(from=1, to=30) и seq(to=30,
from=1), тоже самое что и 1:30. Следующие два параметра seq(), можно
присвоить как by=значение и length=значение, что определяет размер шага
и длину последовательности соответственно. Если ни один из них не задан,
применяется заданное по by=1.

Например

> seq(-5, 5, by=.2) -> s3

генерирует в s3 вектор c(-5.0, -4.8, -4.6, ..., 4.6, 4.8, 5.0).
Точно так же

> s4 <- seq(length=51, from=-5, by=.2)

порождает тот же вектор в s4.

Пятый параметр можно задать как along=вектор, когда он задан другие
параметры недопустимы, и он создает последовательность 1, 2, ...,
length(вектор), или пустую последовательность, если вектор пуст (как
может случится).

Дополнительно можно использовать функцию rep(), которая может повторять
объект различными сложными способами. Простейшей формой является

> s5 <- rep(x, times=5)

, что поместит в s5 пять копий х край-в-край. Другой полезный вариант

> s6 <- rep(x, each=5)

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

2.4 Логические векторы

Также как числовыми векторами, R может манипулировать логическими
значениями. Элементы логического вектора могут иметь значения TRUE,
FALSE и NA (для "пропущено", см. ниже). Первые два часто сокращают
до T и F, соответственно. Заметьте, однако, что T и F всего лишь
переменные, которые установлены в TRUE и FALSE по умолчанию, но
не зарезервированные идентификаторы и, следовательно, могут быть
переназначены пользователем. Поэтому вы всегда должны использовать TRUE
и FALSE. Логические векторы порождаются сравнениями. Например

> temp <- x > 13

задает temp как вектор той же длины, как х со значениями FALSE
соответствующими элементам х, когда условие не выполняется, и TRUE,
когда имеет место.

Логические операторы <, <=, >, >=, == для точного равенства и !=
для неэквивалентности. Кроме того, если c1 и c2 являются логическими
выражениями, тогда c1 & c2 является их умножением ("и"), c1 | c2 их
сложением ("или") и !c1 является отрицанием c1. Логические векторы могут
быть использованы в обычных арифметических операциях, и в этом случае они
приводятся к числовым векторам, FALSE становится 0, а TRUE становится
1. Однако есть ситуации, когда логические векторы и их приведения к
численному виду отличаются, для примера см. следующий подраздел.

2.5 Пропушенные значения

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

Функция is.na(х) дает логический вектор того же размера, как x со
значением TRUE там и только там, где соответствующий элемент х является
NA.

> z <- c(1:3,NA); ind <- is.na(z)

Учтите, что логическое выражения х == NA немного отличается от is.na(х),
поскольку NA на самом деле не значение, а признак отсутствующей
величины. Так х == NA это вектор той же длины что и х, все значения
которого будут NA, так как само логическое выражение является неполным и,
следовательно, неразрешимым.

Заметим, что существует второй тип "пропущенных" значений, которые
возникают в численных вычислениях, так называемое значение не-число,
NaN. Примерами являются

> 0/0

или

> Inf - Inf

, которые оба дают NaN поскольку результат не может быть определен
по правилам.

В целом, is.na(хх) дает TRUE как для NA и значения NaN. Что бы отличить
их, is.nan(хх) дает TRUE только для NaNs. Отсутствующие значения иногда
выводятся как , тогда когда как символьные вектора будут выводится
без квотирования.

2.6 Символьные вектора

Символьные величины и символьные вектора часто используются в R, например,
как метки графиков. При необходимости они описываются последовательностю
символов, ограниченных двойными кавычками, например, "x-values", "New
iteration results".

строки символов вводятся с использованием либо двойные (") или одиночные
(') кавычки, но печатаются с использованием двойных кавычек (или иногда
без кавычек). При этом используются управляющие последовательности в
стиле C, с использованием \ как esc-символ, так \\ вводится и выводится
как \\, а внутри двойных кавычек "вводится как \". Другие полезные
являются управляющие последовательности \n, новая строка, \t, табуляция
и \b, забой.

Символьные вектора можно объединять в вектор с помощью функции с();
случаи такого использования будут встречатся часто. Функция paste()
принимает произвольное количество аргументов и объединяет их по очереди в
одну строку символов. Все числа включенные в аргументы очевидным образом
преобразуются в символьные строки, так же, как если бы они были выведены
на печать. Аргументы по умолчанию разделены в результате одним пробелом,
но это может быть изменено именованым параметром, sep="строка", который
заменяет разделитель на "строка", возможно нулевой длины.

Например

> labs <- paste(c("X","Y"), 1:10, sep="")

преобразует labs в символьный вектор

c("X1", "Y2", "X3", "Y4", "X5", "Y6", "X7", "Y8", "X9", "Y10")

Обратите внимание что имеет место повторное использование коротких
списков; так с("X", "Y") повторяется 5 раз что бы совпасть с
последовательностью 1:10. сноска8

2.7 Индексные вектора; выбор и изменение подмножества набора данных

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

Такой индексный вектор может иметь любой из четырех различных типов.

1. Логический вектор. В этом случае индексный вектор должен быть
той же длины, как и вектор из которого которые должны быть выбраны
элементы. Значения соответствущие TRUE в индексном векторе выбираются,
а соответствующие FALSE пропускаются. Например

> y <- x[!is.na(x)]

создает (или пересоздает) объект y в котором будут содержаться, "не
отсутствующие" значения х, с сохранением порядка. Заметим, что если х
имеет отсутствующие значения, y будет короче, чем х. Итак

> (x+1)[(!is.na(x)) & x>0] -> z

создает объект z и помещает в него значения вектора х+1, для которых
соответствующие значения в х было одновременно не отсутсвующим и
положительным. 2. Вектор положительных целых величин. В этом случае
значения в индексном векторе должны попадать в набор {1, 2, ...,
length(х)}. В результат выбираются и объединяются, в указанном порядке,
соответствующие элементы вектора. Индексный вектор может быть любой длины,
а результат имеет такую же длину как индексный вектор. Например х[6]
является шестым элементом х и

> х[1:10]

выбирает первые 10 элементов х (предполагая, что length(x) не менее
10). Итак

> c("x","y")[rep(c(1,2,2,1), times=4)]

(безусловно маловероятно чтобы так пришлось делать) выводит символьный
вектор длиной 16, содержащий "х", "у", "у", "х" повторенные четыре
раза. 3. Вектор отрицательных целых величин. Такой индексный вектор
определяет значения которые должны быть исключены, а не включены.
Таким образом

> y <- x[-(1:5)]

передает в у все, кроме первых пяти элементов х. 4. Вектор символьных
строк. Эта возможность применяется только тогда, когда объект имеет
такое свойство как имена для обращения к его компонентам. В этом случае
под-вектор вектора имен может быть использован так же, как положительные
целые метки в пункте 2 выше.

> fruit <- c(5, 10, 1, 20) > names(fruit) <- c("orange", "banana",
"apple", "peach") > lunch <- fruit[c("apple","orange")]

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

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

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

Например

> x[ is.na(x)] <- 0

заменяет все отсутствующие значения в х нулями, а

> y[y < 0] <- -y[y < 0]

действует так же, как

> y <- abs(y)

2.8 Другие типы объектов

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

* матрицы или в целом массивы как многомерное обобщение векторов.
Фактически они являются векторами, которые могут быть проиндексированы
двумя или более индексами, и должны выводится в специальном формате.
См. Массивы и матрицы.

* факторы представляют компактный способ работы с качественными
данными. См. Факторы.

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

* таблицы данных это матрице-подобные структуры, в которой столбцы
могут быть различных типов. Думайте о фреймах данных как о "матрицах
данных" с одной строкой на каждое наблюдение, но включающих (возможно)
как числовые переменные, так и качественные. Многие эксперименты лучше
всего записываются во фреймах данных: воздействие качественное, а отклик
колличественный. См. Таблицы данных.

* функции сами также объекты R, которые можно сохранить в окружении памяти
рабочего проекта. Это обеспечивает простой и удобный способ расширения
R. См. Написание собственных функций.


3.1 Обязательные атрибуты: режим и продолжительность

R оперирует записями с технической точки известными как объекты.
Примерами могут быть векторы с численными (реальными) или комплексными
величинами, векторы с логическими значениями и векторы строк символов.
Они известны как структуры "атомы", поскольку их компоненты все
имеют одинаковый тип, или вид, а именно численным (numeric) сноска 9,
комплексным (complex), логическим (logical), символьным (character)
и прямого доступа (raw).

Все значения вектора должны иметь один и тот же тип. Таким образом,
любой вектор должен быть однозначно либо логическим, либо числовым,
либо комплексным, либо символьным или прямого доступа. (Единственным
допустимым исключением из этого правила является специальное "значение"
обозначаемое как NA для отсутствующих значений, хотя в реальности
различают несколько видов NA). Заметим, что вектор может быть пустым,
и все еще иметь тип. Например пустой вектор символьных строк обозначается
как character(0) а пустой численный вектор, как numeric(0).

R также работает с объектами, называемыми списки, которые имеют тип
list. Они являются упорядоченными последовательностями объектов, каждый
из которых может иметь любой тип. Списки известны как "рекурсивные", а не
"атомные" структуры, поскольку их компоненты могут сами быть списками
в своем собственном формате.

Другие рекурсивные структуры имеют типы функция (function) и выражение
(expression). Функции являются объектами, которые являются составной
частью системы R, наряду с аналогичными написанными самим пользователем
функциями, что мы обсудим детально позже. Выражения, как объекты
составляют самую сложную часть R, которая не будет обсуждаться в этом
руководстве, кроме как косвенно, когда мы будем обсуждать формулы,
используемые при моделировании в R.

Под типом объекта, мы имеем в виду базовый тип его основных
составляющих. Это особый случай "свойств" объекта. Другим свойством
каждого объекта является его длина. Функции mode(object) и length(object)
могут использоваться для определения типа и длины любых определенных
структур сноска10.

Другие свойства объекта обычно получают посредством attributes(object),
см. Получение и настройки свойств. Из-за этого, тип и длинну также
называют "внутренними характеристиками" объекта. Например, если z -
это комплексный вектор длины 100, тогда при операции mode(z) выводится
символьная строка "complex", а при length(z) выводится 100.

R заботится об изменении типа практически в любом месте, где это имеет
смысл, (и мало, где - нет). Например имея

> z <- 0:9

мы можем ввести

> digits <- as.character(z)

после чего digits содержит символьный вектор с("0", "1", "2", ...,
"9"). Далее приведение, или смена типа, реконструирует числовой вектора
снова:

> d <- as.integer (digits)

Теперь d и z являются идентичными. сноска11 Существует большой набор
функций вида as.нечто(), как для приведения от одного типа к другому
так для придания объекту ряда новых свойств которыми он пока не
обладает. Читателю придется проконсультироваться с несколькими файлами
руководства, чтобы освоится с этим.

3,2 Изменение длины объекта

"Пустой" объект все еще будет иметь тип. Например

> e <- numeric()

создает в е пустую структуру вектора численного типа. Аналогично
character() - это пустой символьный вектор, и так далее. Когда создан
объект любого размера, новые компоненты могут добавляться просто, придавая
значению его индекса значение за пределом текущего диапазона. Таким
образом

> е[3] <- 17

создает е как вектор длиной 3, (первые две компоненты, на данный момент,
обе NA). Это применимо к абсолютно любой структуре, при условии что тип
дополнительного компонента(ов) подходит типу объекта впереди.

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

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

> alpha <- alpha[2 * 1:5]

делает его объектом длиной 5 состоящим только из бывших компонент с четным
индексом. (Старые значения индекса при этом конечно не сохраняются.) Также
можно сохранить только первые три значения,

> length(alpha) <- 3

или вектор может быть расширен тем же способом ( отсутствующими
значениями).

3.3 Определение и установка свойств

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

Должна быть проявлена некоторая осторожность при назначении или удалении
атрибутов, поскольку они являются неотъемлемой частью используемой в R
системы объектов.

Когда такое используют в левой части присваивания, оно может быть
использовано, чтобы связать новый атрибут с объектом или изменить
существующий. Например

> attr(z, "dim") <- c(10,10)

позволяет R рассматривать z, как если бы это была матрица 10-на-10.

3.4 Класс объекта

Все объекты в R имеют класс, определяемый при помощи функции class.
Для простых векторов это только тип, например, "numeric", "logical",
"character" или "list", но возможны и другие значения "matrix", "array",
"factor" и " data.frame".

Специальный атрибут известный как класс объекта используется для
обеспечения объектно-ориентированного стиля сноска12 программирования в
R. Например, если объект имеет класс "data.frame", он будет напечатан
определенным способом, plot() функция будет отображать его графически
определенным способом, и другие так называемые общие функции, таких как
summary() будут реагировать на него как параметр, способом применимым
к данному классу.

Чтобы удалить временно эффект класса, используйте функцию unclass().
Например, если winter класса "data.frame", тогда

> winter

будет выведена в виде таблицы данных, которая выглядит подобно матрице,
тогда как

> unclass(winter)

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

4 Упорядоченные и неупорядоченные факторы

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

4.1 Специальный пример

Предположим, к примеру, мы имеем выборку 30 налоговых деклараций из всех
штатов и территорий Австралии сноска13, и их индивидуальное происхождение
определяет символьный вектор абревиатуры штата, как

> state <- c("tas", "sa", "qld", "nsw", "nsw", "nt", "wa", "wa", "qld",
"vic", "nsw", "vic", "qld", "qld", "sa", "tas", "sa", "nt", "wa", "vic",
"qld", "nsw", "nsw", "wa", "sa", "act", "nsw", "vic", "vic", "act")

Заметьте, что в случае символьного вектора, "sorted" означает
отсортировано в алфавитном порядке.

Создаются факторы аналогичным образом с помощью функции factor():

> statef <- factor(state)

Функция print() обрабатывает факторы несколько по другому, чем другие
объекты:

> statef

[1] tas sa qld nsw nsw nt wa wa qld vic nsw vic qld qld sa

[16] tas sa nt wa vic qld nsw nsw wa sa act nsw vic vic act

Levels: act nsw nt qld sa tas vic wa

Чтобы выяснить уровни фактора может быть использована функция levels().

> levels(statef)

[1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa"

4.2 функция tapply() и неровные массивы

Чтобы продолжить предыдущий пример, предположим, что мы имеем доходы
от каждого налогоплательщика в отдельном векторе (в подходящих крупных
денежных единицах)

> incomes <- c(60, 49, 40, 61, 64, 60, 59, 54, 62, 69, 70, 42, 56, 61,
61, 61, 58, 51, 48, 65, 49, 49, 41, 48, 52, 46, 59, 46, 58, 43)

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

> incmeans <- tapply(incomes, statef, mean)

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

act nsw nt qld sa tas vic wa

44.500 57.333 55.500 53.600 55.000 60.500 56.000 52.250

Функция tapply() используется для применения функции, в этом случае
mean(), к каждой группе компонентов первого аргумента, в этом случае
incomes, определенных уровнями второго компонента, в этом случае statef
сноска14, как если бы они были отдельными векторными структурами. В
результате получится структура такой же длины, содержащая результаты,
как длина свойства levels фактора. За подробностями читателю следует
обратится к документации.

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

> stderr <- function(x) sqrt(var(x)/length(x))

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

> incster <- tapply(incomes, statef, stderr)

и вычисленные значения равны

> incster

act nsw nt qld sa tas vic wa

1.5 4.3102 4.5 4.1061 2.7386 0.5 5.244 2.6575

В качестве упражнения позаботьтесь о определении стандартных 95% границ
для средних доходов штатов. Для этого можно использовать tapply() еще раз
с функцией length(), что бы найти размер выборки, и qt() функцией что бы
найти значение точек процентов соответствующего t-распределения. (Можно
также исследовать возможности R в t-тестах.)

Функция tapply() может также использоваться для обработки более сложного
индексирования вектора на несколько категорий. Например, мы могли бы
захотеть разделить налогоплательщиков одновременно по штату и полу. Однако
в данном простом случае (только одного фактора), то что произойдет, можно
представить следующим образом. Значения в векторе собираются в группы,
соответствующие различным позициям в факторе. Затем функция применяется
для каждой из этих групп в отдельности. Результат это вектор значений
функции, маркированный согласно значению levels фактора.

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

4.3 Упорядоченные факторы

Уровни факторов хранятся в алфавитном порядке, или в порядке, в каком
они были определены в факторе, если они были определены принудительно.

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

5.1 Массивы

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

Вектор размерности это вектор из неотрицательных целых. Если его
длина k, тогда у массива k измерений, например матрица это 2х-мерный
массив. Размерности индексируются от единицы до значений, заданных в
векторе размерностей.

Вектор может быть использован в R как массив, только если он в качестве
своего свойства dim имеет вектор размерностей. Предположим, например,
z - вектор из 1500 элементов. Присвоение

> dim(z) <- c(3,5,100)

назначает ему свойство dim, что позволяет вектору использоваться как
массив 3 на 5 на 100.

Другие функции, такие, как matrix() и array(), служат для более просто
и естественно выглядящего присвоения, как мы увидим в функции array().

Значения вектора данных задают значения массива в таком же порядке,
как они происходят в FORTRAN, что означает "столбцы более старшие",
первый индекс сдвигается быстрее всех, а последний медленнее всех.

Например, если вектор размерностей массива, скажем a, это с(3,4,2),
то имеется 3 * 4 * 2 = 24 записи в a, и вектор данных содержит их в
порядке a[1,1,1], a[2,1,1], ..., a[2,4,2], a[3,4,2].

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

5.2 Индексация массивов. Подмножества массивов

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

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

Продолжая предыдущий пример, a[2,,] - это 4 * 2 массив с с вектором
размерности c(4,2) и вектором данных, содержащим значения

c(a[2,1,1], a[2,2,1], a[2,3,1], a[2,4,1], a[2,1,2], a[2,2,2], a[2,3,2],
a[2,4,2])

в указанном порядке. a[,,] соответствует всему массиву, что то же самое,
что опустить индексы полностью и использовать a без них.

Для любого массива, скажем Z, вектор размерностей может быть адресован
напрямую как dim(Z) (по обе стороны операции присвоения).

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

5.3 Матрицы индексов

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

Следующий матричный пример позволит разъяснить процесс. В случае массива
с двойным индексом, индексная матрица может быть задана, состоящей из
двух колонок и стольких строк, сколько необходимо. Записями в индексной
матрице являются строки и столбцы индексов для дважды проиндексированного
массива. Предположим, например, мы имеем 4 на 5 массив X, и мы хотели
бы сделать следующее:

* Извлечь элементы X[1,3], X[2,2] и X[3,1], как векторную структуру, и

* Замените эти записи в массиве X на нулями.

В этом случае нам понадобится 3 на 2 индексный массив, как показано в
следующем примере.

> x <- array(1:20, dim=c(4,5)) # Создает массив 4 на 5.

> x

[,1] [,2] [,3] [,4] [,5]

[1,] 1 5 9 13 17

[2,] 2 6 10 14 18

[3,] 3 7 11 15 19

[4,] 4 8 12 16 20

> i <- array(c(1:3,3:1), dim=c(3,2))

> i # i это индексный массив 3 на 2.

[,1] [,2]

[1,] 1 3

[2,] 2 2

[3,] 3 1

> x[i] # Извлечь указанные элементы

[1] 9 6 3

> x[i] <- 0 # Заменить указанные элементы нулями.

> x

[,1] [,2] [,3] [,4] [,5]

[1,] 1 5 0 13 17

[2,] 2 0 10 14 18

[3,] 0 7 11 15 19

[4,] 4 8 12 16 20

Отрицательные индексы не разрешены в индекс матрицах. NA и нулевые
значения допускаются: строки в индекс матрице, содержащие ноль
игнорируются, а строки, содержащие NA порождают в результате NA.

Как менее тривиальный пример, предположим, что мы хотим создать
(нередуцированную) матрицу дизайна для блочного дизайна заданного блоками
факторов (b уровней) и вариаций (v уровней). Далее предположим, есть n
участков в эксперименте. Мы могли бы действовать следующим образом:

> Xb <- matrix(0, n, b)

> Xv <- matrix(0, n, v)

> ib <- cbind(1:n, blocks)

> iv <- cbind(1:n, varieties)

> Xb[ib] <- 1

> Xv[iv] <- 1

> X <- cbind(Xb, Xv)

Для построения матрицы инцидентности, например N, мы могли бы использовать

> N <- crossprod(Xb, Xv)

Однако простым прямым способом получения этой матрицы будет использовании
table():

> N <- table(blocks, varieties)

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

5.4 Функция array()

Как и присвоением векторной структуре атрибута dim, массивы могут быть
так же изготовлены из вектора с помощью функции array, которая имеет форму

> Z <- array(data_vector, dim_vector)

Например, если вектор h содержит 24 или меньше чисел, тогда команда

> Z <- array(h, dim=c(3,4,2))

будет использовать h, что бы создать 3 на 4 на 2 массив в Z. Если размер
h - точно равен 24, результат тот же, как у

> dim(Z) <- c(3,4,2)

Однако если h короче, чем 24, ее значения рециркулируют с начала вновь что
бы дополнить его до размера 24 (см. Правило рециркуляции). Как крайний,
но распространенный случай

> Z <- array(0, c(3,4,2))

делает Z массивом из одних нулей.

На данный момент dim(Z) обозначает вектор размерностей c(3,4,2), и Z[1:24]
соответствует вектору данных, как это было в h, и Z[] с пустыми индексом
или Z без индекса означает весь массив целиком как массив.

Массивы могут использоваться в арифметических выражениях, результатом
является массив сформированный операциями элемент-за-элементом на
векторе данных. Параметры размерности операндов, как правило, требуются
одинаковые, и наследуются вектором размерности результата. Таким образом,
если A, B и C все подобные массивы, а затем

> D <- 2*A*B + C + 1

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

5.4.1 Смешанная арифметика векторов и массивов. Правило рециркуляции

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

* Выражение разбирается слева направо.

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

* Как только встречаются короткий вектор с массивом, массив должен иметь
тот же значение размерности или произойдет ошибка.

* Любой вектор операнд более длинный, чем операнд матрица или массив
генерирует ошибку.

* Если имеется структура массив и не было ошибки или приведения к вектору,
результатом будет структура массив с значением dim общим с его массивами
операндами.

5.5 Внешнее произведение двух массивов

Важной операцией на массивах является внешнее произведение. Если a и
b два числовых массива, их внешнее произведение это массив, чей вектор
размерностей получают конкатенацией двух исходных векторов размерности
(порядок важен), а вектор данных получают получая все возможные
произведения элементов вектора данных a с элементами вектора данных
b. Внешнее произведение выполняет специальный оператор %о%:

> ab <- a %o% b

Альтернатива

> ab <- outer(a, b, "*")

Функцию умножения может быть заменена произвольной функцией двух
переменных. Например, если бы мы хотели вычислить функцию f(x; y) =
cos(y)/(1 + x^2) на регулярной сетке значений с х и у - координатами,
определенными R векторами х и у соответственно, мы могли бы действовать
следующим образом:

> f <- function(x, y) cos(y)/(1 + x^2)

> z <- outer(x, y, f)

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

Пример: детерминант 2 на 2 одноразрядных матриц

Как искусственный, но емкий пример, рассмотрим детерминант 2 на 2 матрицы
[a, b; c, d], где каждая запись является неотрицательным целым числом
в диапазоне 0, 1, ..., 9, то есть цифрой .

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

Аккуратный способ, это использовать функцию outer() дважды:

> d <- outer(0:9, 0:9)

> fr <- table(outer(d, d, "-"))

> plot(as.numeric(names(fr)), fr, type="h", xlab="Determinant",
ylab="Frequency")

Отметьте приведение значения names таблицы частот к numeric, с тем чтобы
восстановить диапазон значений детерминанта. "Очевидный" путь решения
этой проблемы через циклы, которые будут обсуждаться в "Циклы и условное
выполнение", настолько неэффективен, что просто нецелесообразен.

Кроме того, просто удивительно, что почти одна на 20 этих матриц является
сингулярной.

5.6. Обобщенное транспонирование массивов

Функция aperm(a, perm) может быть использована для перестановок
массива a. Аргумент perm должен быть перестановкой из целых (1, ...,
k), где k - число индексов в a. Вывод функции - массив того же размера,
но с учетом того, что perm[j] задала старым размерностям новые j-ные
измерения. Самый простой способ - думать об этой операции, как о обобщение
для транспонирования матриц. Действительно, если A это матрица, (то есть
массив с двойным индексом), то B полученный путем

> B <- aperm(A, c(2,1))

- это просто транспонированная А. Для этого особого случая есть простая
функция t(), поэтому мы могли бы использовать B <- t().


5.7 Возможности обработки матриц

Как отмечалось выше, таблица - это просто массив с двумя индексами.
Однако он является таким важным особым случаем, что он достоин отдельного
обсуждения. R содержит много операторов и функций, которые доступны
только для матриц. Например, как указано выше, t(X) - это функция
транспонирования матрицы. Функции nrow(A) и ncol(A) выдают число строк
и столбцов в матрице соответственно.

* Умножение

* Линейные уравнения и инверсии

* Собственные значения и собственные вектора

* Сингулярное разложения и детерминант

* Подгонка методом наименьших квадратов и QR разложение

5.7.1 Умножение матриц

Оператор %*% используется для умножения матриц. n на 1 или 1 на n матрицы,
конечно, можно использовать в качестве n-вектора, если контекст будет
подходящий. И наоборот, векторы, которые участвуют в выражении матричного
умножения автоматически полагаются векторами строк или столбцов, как
подойдет для умножения, если это возможно (хотя это не всегда однозначно
возможно, как мы увидим позже).

Если, например, А и В квадратные матрицы одинакового размера, то

> A * B

является матрицей произведений элемент на элемент, а

> A %*% B

является произведением матриц. Если х - вектор, то

> х %*% A %*% х

является квадратичной формой. сноска15

Функция crossprod() формирует "векторное произведение", значит что
crossprod(X, у) тоже самое, что t(X) %*% y, но как операция более
эффективна. Если второй аргумент crossprod() пропущен, он принимается
таким же, как первый.

Смысл diag() зависит от ее аргументов. diag(v), где v - вектор,
дает диагональную матрицу с элементами вектора, в качестве значений
диагонали. С другой стороны diag(М), где M - это матрица, дает вектор
значений главной диагонали М. Это те же умолчания, которые используются
для diag() в Matlab. Кроме того, что несколько неожиданно, если k это
одиночное числовое значение, тогда diag(k) - это k на k единичная матрица!

5.7.2 Линейные уравнения и обращение

Решение линейных уравнений является обращенным произведением матриц.
Когда после

> b <- A %*% x

даны только A и b, вектор х является решением этой системы линейных
уравнений. В R,

> solve(A,b)

решает систему, возвращая х (с некоторой потерей точности). Заметим,
что в линейной алгебре, формально x = A^{-1} %*% b, где A^{-1} обозначает
обращение A, что можно вычислить с помощью

solve(A)

, однако редко требуется. Численно, это неэффективные и потенциально
нестабильно вычислять x <- solve(A) %*% b вместо solve(A,b).

Квадратичная форма x %*% A^{-1} %*% x, которая используется в многомерных
расчетах, следует вычислять как то так сноска16 : x %*% solve(A,x),
а не вычисляя обращенную A.

5.7.3 Собственные значения и собственные вектора

Функция eigen(Sm) вычисляет собственное значение и собственный вектор
у симметричной матрицы Sm. Результатом этой функции является список из
двух компонентов называющихся values и vectors. Присвоение

> ev <- eigen(Sm)

будет присваивать этот список ev. Затем ev$val - вектор на собственных
значений Sm и ev$vec содержит матрицу соответствующих собственных
векторов. Если бы мы нуждались только в собственных значениях мы могли
бы использовать присвоение:

> evals <- eigen(Sm)$values

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

> eigen(Sm)

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

> evals <- eigen(Sm, only.values = TRUE)$values

5.7.4 Метод сингулярного разложения и определитель

Функция svd(M) принимает произвольную матрицу как аргумент, M,
и подсчитывает сингулярное разложение М. Оно состоит из матрицы
ортонормальных столбцов U с теми же столбцами, как M, вторая матрица
ортонормальных столбцов V, чьи столбцы - это строки М и диагональная
матрица положительных записей D, таких, что M = U %*% D %*% t(V).
D фактически возвращается как вектор диагональных элементов. В результате
svd(М) фактически является списком из трех компонентов которые называют d,
u и v, соответственно.

Если M является квадратной, то, не трудно понять, что

> absdetM <- prod(svd(M)$d)

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

> absdet <- function(M) prod(svd(M)$d)

после чего мы можем использовать absdet(), как любую другую функцию R. Как
дальнейший тривиальный, но потенциально полезный пример, Вы возможно
решились бы написать функцию, скажем tr(), для расчета остатков квадратных
матриц. [Примечание: Вам не нужно явно использовать цикл. Взгляните еще
раз на функцию diag().]

R имеет встроенные функции для расчета определителя, включая sign,
и determinant, дающие знак и модуль (опционально в логарифмической шкале).

5.7.5 Подгонка методом наименьших квадратов и QR разложение

Функция lsfit() возвращает список содержащий результаты подгонки методом
наименьших квадратов. Присвоение, такое как

> ans <- lsfit(X, y)

дает результаты подгонки методом наименьших квадратов, где y - вектор
наблюдений, а Х матрица плана. См. руководство для получения более
подробной информации, а также про выполнение, помимо прочего, функции
ls.diag() для диагностики регрессии. Заметим, что общее среднее значение
автоматически включается, и не требует непосредственной вставки как
столбец X. Кроме того, знайте, что Вы почти всегда будете предпочитать
использовать lm(.) (См. линейные модели) вместо lsfit() в регрессионных
моделях.

Другой тесно связаной функцией является qr() и ее сателиты. Рассмотрим
следующие присвоения

> Xplus <- qr(X)

> b <- qr.coef(Xplus, y)

> fit <- qr.fitted(Xplus, y)

> res <- qr.resid(Xplus, y)

Это вычислит ортогональную проекцию y в диапазоне X в fit, проекцию на
ортогональное дополнение в res, и вектор коэффициентов для проецирования
в b, то есть в сущности b является результатом оператора Matlab '\'.

Не предполагается, что X имеет полный диапазон столбцов. Избыточность
будет разыскиваться, и как только найдена, устраняться.

Эта альтернатива является старым, низкоуровневым способом выполнить расчет
методом наименьших квадратов. Хотя это еще и используется в некоторых
случаях, это теперь повсеместно заменяется возможностями аппарата
статистических моделей, как будет показано в "Статистические модели в R".

5.8 Объединение частей матриц, cbind() и rbind()

Как мы уже увидели неформально, матрицы могут строиться из других
векторов и матриц с помощью функции cbind() и rbind(). Можно сказать
что cbind() формирует матрицы соединением вместе матриц по горизонтали,
или по-столбцам, а rbind() вертикально, или по-строчно.

В присвоении

> X <- cbind(arg_1, arg_2, arg_3, ...)

аргументы cbind () должены быть либо векторами любой длины, или матриц
с одинаковым размером колонки, то есть одинаковым количеством строк. В
результате матрица из объединения аргументов arg_1, arg_2, ... образующих
столбцы.

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

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

Предположим, X1 и X2 имеют одинаковое число строк. Чтобы совместить их
по-столбцам в матрице X, а также с первой колонке из единиц мы можем
использовать

> X <- cbind(1, X1, X2)

Результат rbind() или cbind() всегда имеет статус матрицы. Поэтому
cbind(х) и rbind(х) возможно простейший путь чтобы вектор х рассматривался
как столбец или строка матрицы соответственно.

5.9 Применение функции объединения с(), с массивами

Следует отметить, что в то время cbind() и rbind() являются функциями
объединения которые учитывают параметр dim, базовая функция c() не
учитывает, а наоборот снимает со всех численных объектов параметры dim
и dimnames. Это иногда по своему полезно.

Официальный способ привести массив обратно к простому векторному объекту
является использование as.vector()

> vec <- as.vector(X)

Однако аналогичный результат может быть достигнут путем использования с()
с только одним аргументом, просто из-за этого побочного эффекта:

> vec <- с(X)

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

5.10 Частотные таблицы из факторов

Напоминаем, что фактор определяет разделение на группы. Аналогично
пара факторов определяет двумерную перекрестную классификацию, и
так далее. Функция table() дает частотную таблицу, которая может быть
вычислена от равных по длине факторов. Если заданы k аргументов-факторов,
в результате получаем k-мерный массив частот.

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

> statefr <- table(statef)

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

> statefr <- tapply(statef, statef, length)

Кроме того предположим, что incomef является фактором предоставляющим
соответствующее определение "вид доходов" для каждого элемента вектора
данных, например используя функцию cut():

> factor(cut(incomes, breaks = 35+10*(0:7))) -> incomef

Тогда для расчета двумерной таблицы частот:

> table(incomef,statef)

statef

incomef act nsw nt qld sa tas vic wa

(35,45] 1 1 0 1 0 0 1 0

(45,55] 1 1 1 1 2 0 1 3

(55,65] 0 3 1 3 2 2 2 1

(65,75] 0 1 0 0 0 0 1 0

Расширение до многомерных частотных таблиц тривиально.

6.1 Списки

Список R - это объект состоящий из упорядоченного набора объектов
называемых компонентами списка.

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

> Lst <- list(name="Fred", wife="Mary", no.children=3,
child.ages=c(4,7,9))

Компоненты всегда пронумерованы и на них всегда можно сослаться по
номеру. Таким образом, если Lst - имя списка с четырьмя компонентами,
на них можно индивидуально сослаться Lst[[1]], Lst[[2]], Lst[[3]] и
Lst[[4]]. Если, далее, Lst [[4]] - это вектор индексный массив тогда
Lst[[4]][1] является его первой записью.

Если Lst - это список, то функция length(Lst) дает число имеющихся
(максимально) компонентов.

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

> имя$имя_компонента

что даст тот же эффект.

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

Так в простом примере, приведенном выше:

Lst$name то же, что Lst[[1]] и это строка "Fred",

Lst$wife то же, что Lst[[2]] и это строка "Mary",

Lst $ child.ages [1] такая же, как Lst [[4]] [1], и это число 4.

Кроме того, можно также использовать имена компонентов списка в двойных
квадратных скобках, т.е. Lst[["название"]] то же, что Lst$название. Это
особенно полезно, когда имя извлекаемого компонента хранится в другой
переменной, к примеру

> x <- "name"; Lst[[x]]

Очень важно различать Lst[[1]] от Lst[1]. `[[...]]' является оператором
используемым для выбора одиночного элемента, а `[...]' это основной
индекс-оператор. Таким образом, первый случай - это первый объект в списке
Lst, и если это именованый спиок имя не включается. Второй случай дает
подсписок из списка Lst состоящий только из первой записи. Если это был
именованный список, имена передаются в подсписок.

Имена компонентов могут быть сокращенны до минимального числа букв
необходимых чтобы определить их уникально. Так Lst$coefficients может
быть минимально определяется как Lst$coe и Lst$covariance как Lst$cov.

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


6,2 Составление и изменение списков

Новые списки могут создаваться из существующих объектов функций
list(). Назначение вида

> Lst <- list(name_1=object_1, ..., name_m=object_m)

устанавливает список Lst из m компонентов используя object_1, ...,
object_m для компонентов и предоставления им имена, определенные
параметром names, (которые могут быть выбраны произвольно). Если эти имена
не упоминаются, компоненты только нумеруются. Компоненты, использованные
для формирования списка, копируются при создании нового списка и исходные
объекты не затрагиваются.

Списки, как и любой индексированный объект, может быть продлены указанием
дополнительных компонентов. Например

> Lst[5] <- list(matrix=Mat)

6.2.1 Объединение списков

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

> list.ABC <- c(list.A, list.B, list.C)

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


6.3 Таблицы данных

Таблицы данных это список класса "data.frame". На списки, которые могут
быть внесены в таблицы данных есть ограничения, а именно

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

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

* Численные и логические вектора, факторы включаются как есть, символьные
вектора приводят к факторам, уровни которых это уникальные значения,
встречающиеся в векторе.

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

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


6.3.1 Создание таблиц данных

Объекты, удовлетворяющие ограничениям накладываемым на столбцы
(компоненты) в таблицах данных, могут использоваться для создания таблицы
данных, с помощью функции data.frame:

> accountants <- data.frame(home=statef, loot=incomes, shot=incomef)

Список компоненты которого соответствуют ограничениям таблиц
данных может быть приведен к таблице данных, при помощи функции
as.data.frame(). Простейший способ построить таблицу данных с нуля
заключается в использовании функции read.table() для чтения всей таблицы
данных из внешнего файла. Это обсуждается далее в "Чтение данных из
файлов".


6.3.2 attach() и detach()

Нотация с использованием $, такая как accountants$statef, для перечисления
компонентов, не так уж удобна. Полезной возможностью будет каким-то
образом сделать компоненты списка или таблицы данных временно видными,
как переменные под их собственными именами, без необходимости каждый
раз указывать прямо название списка.

Функция attach() принимает "базу данных", такую как список или таблицу
данных, как ее аргумент. Таким образом, предположим, lentils является
данных рамы с тремя переменными lentils$u, lentils$v, lentils$w. Подключая

> attach(lentils)

помещает таблицу данных в путь поиска на позицию 2, и представляя
отсутствие каких-либо показателей u, v или w в позиции 1; у, v и
w делаются доступны в качестве переменных таблицы данных по особым
правилам. С этого момента такие присвоения, как

> u <- v+w

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

> lentils$u <- v+w

Однако новое значение компонента u не видно до тех пор, пока таблица
данных не будет отсоединена и присоединена снова.

Чтобы отсоединить таблицу данных, используют функцию

> detach()

Точнее, это выражение отсоединяет от пути поиска записи находящиеся в
настоящее время в позиции 2. Таким образом, в данном контексте переменные
u, v и w будут больше не видны, за исключением нотации списка, такой
как lentils$u и т.п.. Записи на позициях больше, чем 2 в пути поиска
можно отсоединить путем указания их номера в detach, но гораздо
безопаснее всегда использовать название, например, detach(lentils)
или detach("lentils")

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


6.3.3 Работа с таблицами данных

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

* собрать все переменные для каждой четко определеной и отдельной проблемы
в таблицу данных под достаточно информативным именем;

* при работе с проблемой присоединять соответствующие таблицу данных
на позицию 2 и использовать рабочий каталог на уровне 1 для оперативных
значений и временных переменных;

* перед выходом из проблемы, добавить любые переменные, которые нужно
сохранить для будущего использования в таблицу данных, используя $
вид присвоения, а затем detach();

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

Таким образом очень просто работать со многими проблемами в общем
каталоге, во всех из которых есть переменных с названием х, у и z,
к примеру.


6.3.4 Прикрепление произвольных списков

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

> attach(any.old.list)

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


6.3.5 Управление путем поиска

Функция search показывает текущий путь поиска, и это очень полезный
способ отслеживать, какие таблицы данных и списки (а также пакеты)
были присоединены и отсоединены. Первоначально она дает

> search()

[1] ".GlobalEnv" "Autoloads" "package:base"

, где. GlobalEnv является рабочее пространство. сноска17

После присоединения lentils мы имеем

> search()

[1] ".GlobalEnv" "lentils" "Autoloads" "package:base"

> ls(2)

[1] "u" "v" "w"

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

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

> detach("lentils")

> search()

[1] ".GlobalEnv" "Autoloads" "package:base"

7 Чтение данных из файлов

Большие объекты данных, как правило, будут заполнятся значениями из
внешних файлов, а не вводится с клавиатуры в ходе сессии R. Возможности
ввода R просты, и требования к их использованию достаточно строгие
и даже в определенной стапени негибкие. Существует четкая позиция
разработчиков R, что вы можете изменять ваши исходные данные с помощью
других инструментов, таких, как файловые редакторы или Perl1 сноска8
чтобы привести их к требованиям R. В большинстве случаев это очень просто.

Если переменные содержатся в основном в таблицах данных, как мы
настоятельно рекомендуем, целая таблица данных может быть считана напрямую
функцией read.table(). Существует также более примитивная функция ввода -
scan(), к которую можно обращаться напрямую.

Для получения дополнительной информации об импорте данных в R, а также
о экспорте данных, см. руководство "Импорт/экспорт данных в R".


7.1 Функция read.table()

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

* В первая строка файла должна содержать имя для каждой переменной в
таблице данных.

* Каждая дополнительная строка в файле имеет в качестве первого компонента
метку строки и значения для каждой из переменных.

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


Вид входного файла с именами и метками строк:


Price Floor Area Rooms Age Cent.heat

01 52.00 111.0 830 5 6.2 no

02 54.75 128.0 710 5 7.5 no

03 57.50 101.0 1000 5 4.2 no

04 57.50 131.0 690 6 8.8 no

05 59.75 93.0 900 5 1.9 yes

...

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

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

> HousePrice <- read.table ("houses.data")

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

Входной файл форме без меток строк:

Price Floor Area Rooms Age Cent.heat

52.00 111.0 830 5 6.2 no

54.75 128.0 710 5 7.5 no

57.50 101.0 1000 5 4.2 no

57.50 131.0 690 6 8.8 no

59.75 93.0 900 5 1.9 yes

...

Таблица данных может быть считана так

> HousePrice <- read.table("houses.data", header=TRUE)

, где параметр header=TRUE указывает, что первая строка - строка названий,
и, как следует косвенно из вида файла, нет прямо заданных меток строк.

7.2 Функция scan()

Предположим, что вектора данных имеют одинаковую длину и считываются
параллельно. Далее предположим, что есть три вектора, первый символьного
типа и оставшиеся два численного типа, а файл это input.dat. Первым шагом
является использование scan() для считывания трех векторов в виде списка
следующим образом

> inp <- scan("input.dat", list("",0,0))

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

> label <- inp[[1]]; x <- inp[[2]]; y <- inp[[3]]

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

> inp <- scan("input.dat", list(id="", x=0, y=0))

Если Вы хотите получить доступ к переменным по отдельности, они могут
быть реорганизованы в переменные текущей таблицы данных:

> label <- inp$id; x <- inp$x; y <- inp$y

или список может быть присоединен в позицию 2 пути поиска
(см. "Присоединение произвольных списков").

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

> X <- matrix(scan("light.dat", 0), ncol=5, byrow=TRUE)

Есть и более развитые возможности ввода, они подробно описаны в
руководствах.

7.3 Доступ к встроенным наборам данных

Около 100 наборов данных поставляются с R (в наборах данных пакетов),
доступны в расширениях R (в том числе в рекомендованных расширениях
поставлямых с R). Чтобы просмотреть список доступных в настоящее время
наборов данных используйте

data()

С R версии 2.0.0 всех наборы данных поставляемые с R доступны
непосредственно по имени. Однако многие пакеты все еще используют старый
стиль, в котором data также использовалась для загрузки данных в R,
например

data(infert)

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

7.3.1 Загрузка данных из других пакетов R

Для доступа к данным из конкретного пакета, используйте опцию package,
например

data(package="rpart")

data(Puromycin, package="datasets")

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

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

7.4 Редактирование данных

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

> xnew <- edit(xold)

позволит вам изменить ваш набор данных xold, и по завершении измененый
объект присваивается xnew. Если вы хотите изменить сам набор данных xold,
самым простым способом является использование fix(xold), что эквивалентно
xold <- edit(xold).

Используйте

> xnew <- edit(data.frame())

чтобы вводить данные через табличный интерфейс.

8 Распределение вероятности

* R как набор статистических таблиц

* Изучение распределения выборки данных

* Тесты для одной и двух выборок

8.1 R как набор статистических таблиц

Один удобный способ использовать R - наличие полного набора статистических
таблиц. Предусмотрены функции для вычисления кумулятивных функций
распределения P(X <= х), функции плотности вероятности и квантильные
функции (задано q, найдено наименьший х такой, что P(X <= х) > q),
и генерирование выборки по заданному распределению.

Распределение R-наименование дополнительные аргументы

бета beta shape1, shape2, ncp

биноминальное binom size, prob

Каучи-Лоренса cauchy location, scale

хи-квадрат chisq df, ncp

экспоненциальное exp rate

F f df1, df2, ncp

гамма-распределение gamma shape, scale

геометрическое geom prob

гипергеометрическое hyper m, n, k

лог-нормальное lnorm meanlog, sdlog

логистическое logis location, scale

negative binomial nbinom size, prob

нормальное norm mean, sd

Пуассона pois lambda

t Стьюдента t df, ncp

равномерное unif min, max

Вейбулла weibull shape, scale

Уилкса wilcox m, n

Задать префикс 'd' названию для получения плотности, 'р' для кумулятивной
плотности распредения, 'q' для квантильной функции и 'r' для генерации
(случайная выборка). Первый аргумент х для dxxx, q для pxxx, р для qxxx
и n для rxxx (за исключением rhyper и rwilcox, для которых задается nn).
На настоящий момент не во всех случаях доступен параметр нецентральности
ncp: см. интерактивную справочную информацию.

Все pxxx и qxxx функции имеют логические аргументы lower.tail и log.p и
только dxxx из них имеет параметр log. Это позволяет, например, получение
кумулятивных (или "интегральной") функции риска, H(t) = - log(1 - F(t)),
с помощью

- pxxx(t, ..., lower.tail = FALSE, log.p = TRUE)

или более точно непосредственно log-likelihoods (с помощью dxxx(...,
log = TRUE)).

Кроме того, существуют функции ptukey и qtukey для распределения
преобразованного по Стьюденту разброса выборок из нормального
распределения.

Вот некоторые примеры

> ## 2-стороннее р-значения для t распределения

> 2*pt(-2.43, df = 13)

> ## верхняя 1% точка для распределения F(2, 7)

> qf(0.01, 2, 7, lower.tail = FALSE)


8.2 Изучение распределения набора данных

Получив (одномерный) набор данных, мы можем изучить его распределение
многими способами. Самый простой состоит в изучении чисел. Два
немного разных анализа производятся summary и fivenum с одной стороны,
и отображением чисел в консоли с помощью stem с другой (от названия
графика "stem and leaf").

> attach(faithful)

> summary(eruptions)

Min. 1st Qu. Median Mean 3rd Qu. Max.

1.600 2.163 4.000 3.488 4.454 5.100

> fivenum(eruptions)

[1] 1.6000 2.1585 4.0000 4.4585 5.1000

> stem(eruptions)

Десятичная точка - 1 разряд(ы) с левой стороны от |

16 | 070355555588

18 | 000022233333335577777777888822335777888

20 | 00002223378800035778

22 | 0002335578023578

24 | 00228

26 | 23

28 | 080

30 | 7

32 | 2337

34 | 250077

36 | 0000823577

38 | 2333335582225577

40 | 0000003357788888002233555577778

42 | 03335555778800233333555577778

44 | 02222335557780000000023333357778888

46 | 0000233357700000023578

48 | 00000022335800333

50 | 0370

График "стебли-и-листья" подобен гистограмме, в R есть функция hist для
постоения гистограммы.

> hist(eruptions)

## сделать столбцы меньше, вывести плотность распределения

> hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)

> lines(density(eruptions, bw=0.1))

> rug(eruptions) # show the actual data points

Больше элегантный график плотности распределения может быть сделан
с помощью density, и мы добавили в этом примере линию вычисленную
density. Окно bw было выбрано методом проб и ошибок поскольку окно по
умолчанию дает слишком сильное сглаживание (оно обычно делается для
"интересных" плотности). (Доступны и лучше автоматизированные методы
выбора окна, и в нашем примере bw="SJ" дает хороший результат.)

Мы можем вывести эмпирическую кумулятивную функцию распределения с
помощью функции ecdf.

> plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)

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

> long <- eruptions[eruptions > 3]

> plot(ecdf(long), do.points=FALSE, verticals=TRUE)

> x <- seq(3, 5.4, 0.01)

> lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)

Quantile-quantile (Q-Q) график может помочь нам изучить это более
тщательно.

par(pty="s") # организует квадратную область графического вывода

qqnorm(long); qqline(long)

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

х <- rt (250, df = 5)

qqnorm (х); qqline (х)

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

qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q plot for t dsn")

qqline(x)

Наконец, мы могли бы захотеть более формального критерия принятия
нормальности (или отклонения). R предоставляет тест Шапиро-Уилкса

> shapiro.test(long)

Shapiro-Wilk normality test

data: long

W = 0.9793, p-value = 0.01052

и тест Колмогорова - Смирнова

> ks.test (long, "pnorm", mean = mean(long), sd = sqrt(var(long)))

One-sample Kolmogorov-Smirnov test

data: long

D = 0.0661, p-value = 0.4284

alternative hypothesis: two.sided

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

8.3 Одно- и двух- выборочные тесты

До сих пор мы сравнивали единственную выборку с нормальным
распределением. Гораздо более общая операция это сравнить характеристики
двух выборок. Отметим, что все "классические" критерии в R, в том числе
используемые ниже, содержаться в расширении stats, которое обычно уже
загружено.

Рассмотрим следующие наборы данных об скрытой теплоте таяния льда (калл/г)
по Rice (1995, стр.490)

Method A: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97

80.05 80.03 80.02 80.00 80.02

Method B: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97

Boxplot обеспечивают простое графическое сравнение двух выборок.

A <- scan()

79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97

80.05 80.03 80.02 80.00 80.02

B <- scan()

80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97

boxplot(A, B)

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

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

> t.test(A, B)

Welch Two Sample t-test

data: A and B

t = 3.2499, df = 12.027, p-value = 0.00694

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.01385526 0.07018320

sample estimates:

mean of x mean of y

80.02077 79.97875

который свидетельствует о значительной разнице, принимая нормальность
распределения. По умолчанию функция R не предполагает равенства дисперсий
в обоих выборках (в отличие от аналогичной t.test функции в S-Plus ). Мы
можем использовать F-тест для проверки равенства дисперсий, при условии,
что обе выборки взяты из нормальной генеральной совокупности.

> var.test(A, B)

F test to compare two variances

data: A and B

F = 0.5837, num df = 12, denom df = 7, p-value = 0.3938

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

0.1251097 2.1052687

sample estimates:

ratio of variances

0.5837405

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

> t.test(A, B, var.equal=TRUE)

Two Sample t-test

data: A and B

t = 3.4722, df = 19, p-value = 0.002551

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

0.01669058 0.06734788

sample estimates:

mean of x mean of y

80.02077 79.97875

Все эти тесты требуют нормальность обоих выборок. Двух выборочный тест -
Уилкокса (он же Манна-Уитни) предполагает для нулевой гипотезы только
единое непрерывное распределение.

> wilcox.test(A, B)

Wilcoxon rank sum test with continuity correction

data: A and B

W = 89, p-value = 0.007497

alternative hypothesis: true location shift is not equal to 0

Warning message:

Cannot compute exact p-value with ties in: wilcox.test(A, B)

Обратите внимание на предупреждение: есть несколько пучностей в каждой
выборке, которые сильно свидетельствуют, что эти данные взяты из
дискретного распределения (возможно это из-за округления).

Есть несколько способов визуально сравнить два образца. Мы уже видели
парный boxplot. Следующим

> plot(ecdf(A), do.points=FALSE, verticals=TRUE, xlim=range(A, B))

> plot(ecdf(B), do.points=FALSE , verticals=TRUE, add=TRUE)

будут выведены две эмпирических кумулятивных функции распределения, а
qqplot построит Q-Q график по двум выборкам. Тест Колмогорова-Смирнова
основан на определении максимального вертикального расстояния между
двумя ecdf, предполагая распределение единым и непрерывным:

> ks.test(A, B)

Two-sample Kolmogorov-Smirnov test

data: A and B

D = 0.5962, p-value = 0.05919

alternative hypothesis: two-sided

Warning message:

cannot compute correct p-values with ties in: ks.test(A, B)

9 Группировка, циклы и условное исполнение

* Группировка выражений

* Управляющие операторы

9.1 Группировка выражений

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

Команды могут быть сгруппированы скобками, {expr_1; ...; expr_m},
и в этом случае значением этой группы является результат последнего
выражения вычисленного в группе. Такая группа сама является выражением,
например она может быть включена в скобки, и использоваться как часть
еще большего выражения, и так далее.


9.2 Управляющие операторы

* Условное исполнение

* Повторяющееся выполнение


9.2.1 Условное исполнение: оператор if

В языке имеются условные конструкции в форме

> if (expr_1) expr_2 else expr_3

где expr_1 должно вычислятся в одиночное логическое значение и результат
всего выражения является тогда очевиден.

Операторы "сокращения" && и || часто используются как часть условия if
операторе. Тогда как & и | применяются по-элементно к векторам, && и ||
применяется к векторам единичной длины, и только вычисляют свой второй
аргумент, если это необходимо.

Существует векторизированный вариант if/else конструкции, это функция
ifelse. Она имеет форму ifelse(condition, a, b), и возвращает вектор
длины равной длинне его наибольшему аргументу, с элементами a[i], если
condition[i] истинно, в противном случае b[i].

9.2.2 Повторяющееся выполнение: для циклов, repeat и while

Итак для цикла имеется конструкция форма которой

> for (name in expr_1) expr_2

где name - переменная цикла. expr_1 - векторное выражение, (часто
последовательность, подобная 1:20), часто expr_2 это сгруппированное
выражение, а составляющие его выражения написаны используя name. expr_2
вычисляются покругу следом за тем, как name движется по значениям в
векторе результате expr_1.

В качестве примера предположим, ind - вектор типа indicators, и мы
хотим подготовить отдельные графики y на х в каждой из групп. Одна
из возможностей заключается в использовании coplot(), сноска19,
который выведет таблицу графиков, соответствующих каждому из уровней
фактора. Другой способ сделать это, просто одновременно выведя графики
в одно окно, выглядит следующим образом:

> xc <- split(x, ind)

> yc <- split(y, ind)

> for (i in 1:length(yc)) {

plot(xc[[i]], yc[[i]]);

abline(lsfit(xc[[i]], yc[[i]]))

}

(Примечание: функция split(), которая выводит список векторов, полученных
путем разделения большего вектора в соответствии с разбиением заданным
фактором. Это полезная функция, в основном используемая в связке с
boxplots. Подробнее см. на странице помощи.)

Предупреждение: цикл for() используются в R кода значительно реже,
чем в компилируемых языках. Код, который оперирует объектом "как целым"
одновременно и более прозрачный и более быстрый в R.

К другим возможностям циклов относятся конструкции

> repeat expr

и

> while (condition) expr

Оператор break может быть использован для прекращения любого
цикла, возможно досрочно. Это единственный способ прекратить цикл
repeat. Оператор next может быть использован для прекращения одного
конкретного цикла, и перехода к "следующему". Управляющая конструкция
чаще всего используется в соединении с функциями, которые обсуждаются в
"Написание ваших собственных функций", тогда и появятся другие примеры.

10 Написание вашей собственной функции

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

Следует подчеркнуть, что большинство функций поставляющихся как часть R
системы, например mean(), var(), postscript() и т.п., сами написаны на R
и, следовательно, не существенно отличаются от написанных пользователем
функций.

А функция определяется присвоением вида

> name <- function(arg_1, arg_2, ...) expression

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

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

10.1 Простые примеры

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

Эта функция определяется следующим образом:

> twosam <- function(y1, y2) {

n1 <- length(y1); n2 <- length(y2)

yb1 <- mean(y1); yb2 <- mean(y2)

s1 <- var(y1); s2 <- var(y2)

s <- ((n1-1)*s1 + (n2-1)*s2)/(n1+n2-2)

tst <- (yb1 - yb2)/sqrt(s*(1/n1 + 1/n2))

tst

}

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

> tstat <- twosam(data$male, data$female); tstat

В качестве второго примера рассмотрим функции для прямой эмуляции команды
Matlab "обратный слешь", которая возвращает коэффициенты ортогональной
проекции вектора y на колоночное пространство матрицы X. (Это делают
для оценки методом наименьших квадратов коэффициентов регрессии.) Это
обычно делается функцией qr(), но ее иногда немного сложно использовать
напрямую, так что имеет смысл иметь простую функцию, как указано ниже,
и использовать ее без опаски.

Таким образом, если n на 1 вектор у и n на р матрицу X, тогда X\y
определяется как (X'X )^{-}X'y, где (X'X)^{-} является обобщенная обратной
величиной X'X.

> bslash <- function(X, y) {

X <- qr(X)

qr.coef(X, y)

}

После того как этот объект был создан он может быть использован в таких
командах, как

> regcoeff <- bslash(Xmat, yvar)

и т.п.

Классическая R функция lsfit() выполнит это задание так же хорошо,
и детальнее ссылка20. Выше вместо этого, несколько не интуитивным
образом используются функции qr() и qr.coef(), что бы сделать эту часть
расчета. Отсюда, вероятно, если она будет использоваться часто имеется
определенная ценность в том, что бы иметь изолированно только эту часть
вычислений, в виде простой в использовании функции. Если это так, мы,
возможно, для еще более удобного использования пожелает сделать ее
матричным бинарным оператором.

10.2 Определение новых бинарных операторов

Если бы мы задали функции bslash() другое имя, а именно одну из форм

%anything%

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

> "%!%" <- function(X, y) { ... }

(Заметьте что используются кавычки.) Эта функция может быть использована
в виде X %!% y. (Символ "обратный слешь" сам по себе не удобный выбор,
так как он представляет особые проблеммы в таком контексте.)

Оператор умножения матриц, %*%, и оператор внешнего произведения %о%
это другие примеры определенных таким образом бинарных операторов.

10.3 Именованные аргументы и аргументы по умолчанию

Как впервые отмечается в "Генерация регулярных последовательностей",
если аргументы вызываемых функций, задаются в форме "name=object", то
они могут быть указаны в любом порядке. Кроме того, последовательность
аргументов может начаться не поименованной, позиционной формой, а
поименованные аргументы можно указать после позиционных аргументов.

Таким образом, если есть функция fun1 определена как

> fun1 <- function(data, data.frame, graph, limit) {

[function body omitted]

}

тогда функция может вызываться несколькими способами, к примеру

> ans <- fun1(d, df, TRUE, 20)

> ans <- fun1(d, df, graph=TRUE, limit=20)

> ans <- fun1(data=d, limit=20, graph=TRUE, data.frame=df)

все являются эквивалентными.

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

> fun1 <- function(data, data.frame, graph=TRUE, limit=20) { ... }

ее можно было бы вызвать как

> ans <- fun1(d, df)

что эквивалентно трем примерам выше, или как

> ans <- fun1(d, df, limit=10)

что изменяет одно из значений по умолчанию.

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

10.4 Аргумент `...'

Другая частая потребность разрешить одной функции передать набор
аргументов в другую. Например, многие графические функции используют
функцию par() и функции, такие как plot() позволяют пользователю
передавать графические параметры в par() для контроля графического
вывода. (См. "Функция par()", для более подробной информации по функции
par().) Это можно сделать задавая дополнительный аргумент, буквально
`...', у этой функции, который может затем быть передан. Макет такого
примера приводится ниже.

fun1 <- function(data, data.frame, graph=TRUE, limit=20, ...) {

[пропущены команды]

if (graph)

par(pch="*", ...)

[еще пропущено]

}

10.5 Присвоения в функции

Заметим, что любое обычное присвоение сделаное в функции будет локальным
и временным и теряется после выхода из функции. Таким образом, присвоение
X <- qr(X) не влияет на значение аргумента в вызывающей программе.

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

Если требуются глобальные и постоянные назначения в рамках функций,
то могут быть использованы "superassignment" оператор, <<- или функция
assign(). Подробнее см. документацию. S-Plus пользователи должны понимать,
что <<- имеет отличающуюся семантику в R. Эти вопросы обсуждаются дальше в
"Состояние".

10.6 Более сложные примеры

10.6.1 Влияющие факторы в блочной схеме эксперимента

Как более полный, однако похожий на "закат солнца вручную", пример
функции, рассмотрим нахождение влияющих факторов для блочной схеме
эксперимента. (Некоторые аспекты этой проблемы уже обсуждались в
"Индексные матрицы".)

Блочная схема эксперимента определяется двумя факторами, к примеру blocks
(b уровней) и varieties (v уровней). Если R и K являются соответственно
v на v и b на b матрицами откликов и размеров блоков, и N является b на v
матрицей инцидентности, то влияющие факторы определяются как собственные
значения из матрицы E = I_v - R^{-1/2}N'K^{-1}NR^{-1/2} = I_v - A'A,
where A = K^{-1/2}NR^{-1/2}. Один из способов, чтобы написать такую
функцию, приводится ниже.

> bdeff <- function(blocks, varieties) {

blocks <- as.factor(blocks) # minor безопасно перемещены

b <- length(levels(blocks))

varieties <- as.factor (varieties) # minor безопасно перемещены

v <- length(levels(varieties))

K <- as.vector(table(blocks)) # удаление атрибута dim

R <- as.vector(table(varieties)) # удаление атрибута dim

N <- table(blocks, varieties)

A <- 1/sqrt(K) * N * rep(1/sqrt(R), rep(b, v))

sv <- svd(A)

list(eff=1 - sv$d^2, blockcv=sv$u, varietycv=sv$v)

}

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

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


10.6.2 Сбрасывание всех имен при печати массива

Для целей печати большой матрицы или массива, часто полезно печатать
их форме плотного блока без всяких имен или индексов. Удаление атрибута
dimnames не дает такого эффекта, напротив массиву должен быть назначен
атрибут dimnames, состоящий из пустых строк. Например, чтобы распечатать
матрицу X

> temp <- X

> dimnames(temp) <- list(rep("", nrow(X)), rep("", ncol(X)))

> temp; rm(temp)

Это может быть гораздо более удобно сделать с помощью функции,
no.dimnames(), показаной ниже, как "обертка", чтобы достичь того же
результата. Она также показывает, что некоторые эффективные и полезные
функции пользователь могут быть достаточно короткими.

no.dimnames <- function(a) {

## Удалить все имена размерности из массива для компактной печати.

d <- list()

l <- 0

for(i in dim(a)) {

d[[l <- l + 1]] <- rep("", i)

}

dimnames(a) <- d

a

}

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

> no.dimnames(X)

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

10.6.3 Рекурсивное численное интегрирование

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

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

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

area <- function(f, a, b, eps = 1.0e-06, lim = 10) {

fun1 <- function(f, a, b, fa, fb, a0, eps, lim, fun) {

## функция `fun1" видна только внутри 'area'

d <- (a + b)/2

h <- (b - a)/4

fd <- f(d)

a1 <- h * (fa + fd)

a2 <- h * (fd + fb)

if(abs(a0 - a1 - a2) < eps || lim == 0)

return(a1 + a2)

else {

return(fun(f, a, d, fa, fd, a1, eps, lim - 1, fun) +

fun(f, d, b, fd, fb, a2, eps, lim - 1, fun))

}

}

fa <- f(a)

fb <- f(b)

a0 <- ((fa + fb) * (b - a))/2

fun1(f, a, b, fa, fb, a0, eps, lim, fun1)

}

10.7 Область действия

Обсуждение в этом разделе несколько более техническое, чем в других
частях этого документа. Однако это детали одного из основных различий
между S-Plus и R.

Символы, которые встречаются в теле функции могут быть разделены на
три класса; формальные параметры, локальные переменные и свободные
переменные. Формальные параметры функции это те, которые входят в список
аргументов функции. Их значение определяется в процессе связывания
фактических аргументов функции с формальными аргументами. Локальные это
переменные, значения которых определяются вычислением выражений в теле
функции. Переменные, которые не являются формальными параметрами или
локальными переменными называются свободными переменными. Свободные
переменные становятся локальными переменными, если им они
присваиваются. Рассмотрим следующие определение функции.

f <- function(x) {

y <- 2*x

print(x)

print(y)

print(z)

}

В этой функции, х - формальный параметр, y является локальной переменной,
а z является свободной переменной.

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

cube <- function(n) {

sq <- function() n*n

n*sq()

}

Переменная n в функции sq - это не аргумент для этой функции. Поэтому
она выступает свободной переменной и для определения значения, которое
должно быть ей назначено, должны быть использованы правила области
действия. Согласно статической области действия (S-Plus) значение
определяется глобальной переменной с названием n. В области действия (R)
она является параметром функции cube, поскольку это является активным
связыванием для переменной n в момент определения функции sq. Разница
между вычислением в R и вычислением в S-Plus заключается в том, что S-Plus
ищет глобальную переменную с названием n а R сначала ищет переменную
называемую n в окружении, созданном после определения cube.

## сначала вычислим в S

S> cube(2)

Error in sq(): Object "n" not found

Dumped

S> n <- 3

S> cube(2)

[1] 18

## таже функция вычисляется в R

R> cube(2)

[1] 8

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

В специальный оператора присваивания, <<-, используется для изменения
значения, связанного с total. Этот оператор смотрит назад в окружающие
среды в поисках такого окружения, которое содержит символ total, и когда
он находит такое окружение он заменяет значение, в этом окружении, со
значением с правой стороны. Если глобальное или верхнего уровня окружение
достигается и без нахождения символа total, тогда такая переменная
создается и значение присваивается ей. Для большинства пользователей
<<- создает глобальную переменную и присваивает значение справа на
лево ссылка21. Только тогда, когда <<- была использована в функции,
что была возвращена как значение другой функции происходят описываемое
здесь специальное поведение.

open.account <- function(total) {

list(

deposit = function(amount) {

if(amount <= 0)

stop("Deposits must be positive!\n")

total <<- total + amount

cat(amount, "deposited. Your balance is", total, "\n\n")

},

withdraw = function(amount) {

if(amount > total)

stop("You don't have that much money!\n")

total <<- total - amount

cat(amount, "withdrawn. Your balance is", total, "\n\n")

},

balance = function() {

cat("Your balance is", total, "\n\n")

}

)

}

ross <- open.account(100)

robert <- open.account(200)

ross$withdraw(30)

ross$balance()

robert$balance()

ross$deposit(50)

ross$balance()

ross$withdraw(500)

10.8 Настройка окружения

Пользователи могут настраивать свои окружения несколькими
различными способами. Существует системный файл инициализации и
каждый рабочий каталог может иметь свои собственные специальные файлы
инициализации. Наконец могут быть использованы специальные функции .First
и .Last.

Место расположения системного файл инициализации берется из значения
переменной среды R_PROFILE. Если эта переменная не установлена
используется файл Rprofile.site в подкаталоге etc домашней директории
R. Этот файл должен содержать команды, которые Вы хотите исполнять каждый
раз когда запускаете R в вашей системе. Второй, индивидуальный, профайл
с именем .Rprofile сноска22 может быть размещен в любом каталоге. Если R
запускается в этом каталоге, то файл будет считан. Этот файл предоставляет
индивидуальным пользователям контроль над их рабочим окружением, и
дает в различных рабочих каталогах возможность отличающихся процедур
запуска. Если файла .Rprofile нет в каталоге запуска, тогда R ищет файл
.Rprofile в пользовательском домашнем каталоге, и использует его (если
он существует).

Любая функция с именем .First() в любом из этих двух профайлов или в
образе .RData имеет особый статус. Они автоматически выполняться в начале
R сессии и могут использоваться для инициализации среды. Например,
определение в примере ниже изменяет приглашение системы на $ и
устанавливает различные другие полезные вещи, которые затем могут
использоваться как само собой разумеющееся в остальной части сессии.

Таким образом, последовательность, в которой файлы исполняются это
Rprofile.site, .Rprofile, .RData а затем .First(). Определения в
последующих файлах маскируют определения в предыдущих файлах.

> .First <- function() {

options(prompt="$ ", continue="+\t") # $ как приглашение

options(digits=5, length=999) # формат чисел и область печати

x11() # графический драйвер

par(pch = "+") # рисовать символом

source(file.path(Sys.getenv("HOME"), "R", "mystuff.R"))

# мои собственные функции

library(MASS) # присоединить пакет

}

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

> .Last <- function() {

graphics.off() # небольшая мера безопасности.

cat(paste(date(),"\nAdios\n")) # Время обеда наступило?

)

10.9 Классы, общие функции и объектно-ориентированное программирование

Класс объекта определяет, как он будет обработан так называемыми общими
функциями. Посмотрев с другой стороны, общая функция выполняет действия
и обработку своих аргументов в зависимости от конкретного класса самого
аргумента. Если у аргумента отсутствует атрибут класса, или он имеет
класс не удовлетворяющий формату общей функции, всегда предусмотрено
действие по умолчанию.

Проясним это примером. Класс mechanism предоставляет пользователю
возможность проектирования и подготовки общих функций для специальных
целей. Среди многих других общих функций имеется plot() для вывода
объектов графически, summary() для подведения итогов анализов различных
типов, и anova() для сравнения статистических моделей.

Количество общих функций, которые могут относиться к классу в
специфическим способом может быть очень большим. Например, функции,
которые могут быть применены к некоторым распространенным объектам класса
"data.frame " включают

[ [[<- any as.matrix

[<- mean plot summary

Текущий полный список можно получить используя функцию methods():

> methods(class="data.frame ")

Число классов которые общая функция может обрабатывать также может быть
очень большим. Например, функция plot() метод по умолчанию и варианты
для объектов классов "data.frame", "density", "factor", и многое
другое. Полный список можно получить с помощью функции methods():

> methods(plot)

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

> coef

function (object, ...)

UseMethod("coef")

Присутствие UseMethod указывает что это общая функция. Чтобы узнать,
какие методы имеются мы можем использовать methods()

> methods(coef)

[1] coef.aov* coef.Arima* coef.default* coef.listof*

[5] coef.nls* coef.summary.nls*

Невидимые функции обозначены звездочками

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

> getAnywhere("coef.aov")

A single object matching ' coef.aov' was found

It was found in the following places

registered S3 method for coef from namespace stats

namespace:stats

with value

function (object, ...)

{

z <- object$coef

z[!is.na(z)]

}

> getS3method("coef", "aov")

function (object, ...)

{

z <- object$coef

z[!is.na(z)]

}

Читатель может обратится к "Определение языка R" для более полного
обсуждения этого механизма.

11 Статистические модели в R

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

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

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

11.1 Определение статистических моделей; формулы

Шаблон для статистической модели - это модель линейной регрессии с
независимыми, гомоскедастическими ошибками

y_i = sum_{j=0}^p beta_j x_{ij} + e_i, i = 1, ..., n,

где e_i являются NID(0, sigma^2). В матричной нотации это будет записано,
как

y = X beta + e

, где y является вектором отклика, X - это модель-матрица или матрица
плана имеющая столбцы x_0, x_1, ..., x_p, определяющих переменных. Очень
часто x_0 это столбец из единиц определяющий постоянное слагаемое.

Примеры

Прежде чем дать формальное определение, несколько примеров помогут
составить общее представление.

Предположим, у, х, x0, x1, x2, ... являются численные переменные, X -
это матрица и A, B, C, ... являются факторами. Ниже следующие формулы
задают статистические модели, справа даны описания моделей.

у ~ х

у ~ 1 + х

Обе вводят одну и ту же простую модель линейной регрессии, y по х. Первая
имеет неявно заданный свободный член, а вторая прямо определенный.

у ~ 0 + х

и ~ -1 + х

у ~ х - 1

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

log(y) ~ x1 + x2

Множественная регрессия трансформированной переменной, log(у), по x1 и x2
(со свободным членом).

y ~ poly(x,2)

y ~ 1 + x + I(x^2)

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

y ~ X + poly(x,2)

Множественная регрессия и моделью-матрицей состоящей как из матрицы X,
а также в терминах полинома по х степени 2.

y ~ А

Модель однофакторного дисперсионного анализа y, с классами определяемыми
A.

у ~ А + х

Модель однофакторного дисперсионного анализа y, с классами определяемыми
A, и с ковариациями х

y ~ A*B

y ~ A + B + A:B

y ~ B %in% A

y ~ A/B

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

y ~ (A + B + C)^2

y ~ A*B*C - A:B:C

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

y ~ A * x

y ~ A/x

y ~ A/(1 + x) - 1

Изолированные модели простой линейной регрессии y по х внутри уровней
заданных в А, различными метками. В последнем виде производит четко
столько вычислений различных отсекаемых отрезков и коэффициентов наклона,
сколько имеется уровней A.

y ~ A*B + Error(C)

Эксперимент с двумя факторами воздействия, А и В, и стратифицированной
ошибкой определяемой фактором C. Например, разделить отображение
эксперимента, на участки (и следовательно подрисунки), определяемые
фактором C.

Оператор "~" используется для определения формулы модели в R. Форма,
для простой линейной модели, будет выгдядеть

response ~ op_1 term_1 op_2 term_2 op_3 term_3 .....

где

response

это вектор или матрица (или выражение вычисленное в вектор или матрицу)
задающие переменную(е) отклика.

op_i

является оператором, либо + или -, подразумевая включения или исключения
значений в модель, (первый не обязательный).

term_i

это либо

* вектор либо матричное выражение, или 1,

* фактор, или

* формула в составе факторов, векторов или матриц связаных операторами
формулы.

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

Операторы формула аналогичны по действию операторам в нотации Уилкинсона
и Роджерса, используемых такими программами как Glim и Genstat. Одно из
неизбежных изменений заключается в том, что оператор "." заменяется ":"
поскольку точка является легальным символом имени в R.

Описание нотации приводится ниже (на основе Chambers & Hastie, 1992,
с. 29):

Y ~ M

Y моделируется как М.

M_1 + M_2

Включить M_1 и M_2.

M_1 - M_2

Включить M_1 исключив M_2.

M_1 : M_2

Тензорное произведение M_1 и M_2. Если оба условия являются факторами,
то фактор "подклассов".

M_1 %in% M_2

Как M_1:M_2, но с другим синтаксисом.

M_1 * M_2

M_1 + M_2 + M_1:M_2.

M_1 / M_2

M_1 + M_2 %in% M_1.

M^n

Все условия из М вместе с "взаимодействиями" вплоть до порядка n

I(М)

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

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

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

11.1.1 Контрасты

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

Как насчет k-уровневого фактора А? Ответ отличается для неупорядоченного
и упорядоченного фактора. Для неупорядоченного факторов, k-1 столбцов
генерируются для показатели второго, ..., k-го уровней фактора. (Таким
образом, применяемая параметризация создает на каждом уровне такой
же контраст откликов, что и на первом.) Для упорядоченных факторов,
k-1 столбцы являются ортогональными полиномами по основанию 1, ..., k,
исключая параметры константы.

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

options(contrasts = c("contr.treatment", "contr.poly"))

Основная причина, по которой это упомянуто, это то что R и S имеют
различные значения по умолчанию для неупорядоченного факторов, S
использует контрасты Хелмерта. Так что, если вам нужно сравнить свои
результаты с учебником или статьей, где используется S-Plus, нужно
будет установить

options(contrasts = c("contr.helmert", "contr.poly"))

Это разница преднамеренная, так как считается, что threatment контрасты
(по умолчанию в R) новичкам легче интерпретировать.

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

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

11.2 Линейные модели

Основная функция для установки простых множественных моделей lm(),
а усовершенствованный вариант вызова выглядит следующим образом:

> fitted.model <- lm(formula, data = data.frame)

Например

> fm2 <- lm(y ~ x1 + x2, )

будет соответствовать множественной регрессионной модели y по x1 и x2
(с параметром свободного члена).

Важный (но технически необязательный) параметр data = production
указывает, что все параметры, необходимые для построения модели следует
брать в первую очередь из произведения таблицы данных. Это случай
независимости, подключено произведение таблицы данных в путь поиска,
или нет.

11.3 Общие функции для извлечения информации о модели

Результатом lm() является подогнанный объект model; технически это список
результатов класса "lm". Информация о подогнанной модели может быть
затем отображена, извлечена, визуализирована и т.п., с использованием
общих функций, ориентированных на объекты класса "lm". К ним относятся

add1 deviance formula predict step

alias drop1 kappa print summary

anova effects labels proj vcov

coef family plot residuals

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

anova(object_1, object_2)

Сравнить субмодель с внешней моделью и вывести таблицу дисперсионного
анализа.

coef(object)

Вывести коэффициенты (матрицу) регрессии.

Полный синтаксис: coefficients(object).

deviance(object)

Сумма квадратов остатка, взвешенных в случае необходимости.

formula(object)

Вывести формулу модели.

plot(object)

Выводит четыре графика, показывающих остатки, подогнанные значения и
некоторую диагностическую информацию.

predict(object, newdata= data.frame)

Таблицы данных должны содержать переменные с теми же метками, как и
в оригинале. Результат - вектор или матрица предсказанных значений в
соответствии с подставленными значениями из data.frame.

print(object)

Печать краткого содержания объекта. Наиболее часто используются неявно.

residuals(object)

Извлечение остатков (матрицы остатков), взвешенное если необходимо.

Короткий синтаксис: resid(object).

step(object)

Выбрать подходящую модель, добавив или сбросив параметр с сохранением
иерархии. Возвращается модель с наименьшим значением AIC (Akaike's An
Information Criterion), обнаруженная в пошаговом поиске.

summary(object)

Распечатать всеобъемлющее резюме результатов регрессионного анализа.

vcov(object)

Возвращает вариационно-кавариационную матрицу основных параметров объекта
подогнанной модели.

11.4 Дисперсионный анализ и сравнение моделей

функции подгонки модели aov(formula, data=data.frame) действует на
простейшем уровне очень похоже на функцию lm(), и большинство общих
функций, перечисленных в таблице "Общие функции" применимы для извлечения
информации из модели.

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

response ~ mean.formula + Error(strata.formula)

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

Например, со всеми переменными факторами, формула модели, такая, как
в выражении:

> fm <- aov(yield ~ v + n*p*k + Error(farms/blocks), data=farm.data)

, как правило, используется для описания эксперимента со средним модели
v + n*p*k и тремя слоями ошибок, а именно "между farms", "внутри farms,
между bloks" и "в рамках blocks".

11.4.1 Таблица ANOVA

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

Для многослойного эксперимента начинают процедуру проекцией отклика на
слои ошибки, по порядку, и затем подгоняют модельное среднее для каждой
проекции. Более подробную информацию см. Chambers & Hastie (1992).

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

> anova(fitted.model.1, fitted.model.2, ...)

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

11.5 Обновление подогнанной модели

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

> new.model <- update(old.model, new.formula)

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

> fm05 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = production)

> fm6 <- update(fm05, . ~ . + x6)

> smf6 <- update(fm6, sqrt(.) ~ .)

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

Примечание: в особенности если задано data=аргумент в первоначальном
вызове функции подгонки модели, то эта информация передается с помощью
объекта подогнанной модели в update() и ее союзникам.

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

> fmfull <- lm(y ~ . , data = production)

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

Другие функции для изучения пошаговых последовательностей моделей add1(),
drop1() и step(). Названия эти дают хорошее понятие о их назначении,
но более подробную информацию см. в системе интерактивной помощи.

11.6 Обобщенная линейная модель

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

* Существует отклик, y, и воздействующие переменные x_1, x_2, ...,
значения которых влияют на распределение отклика.

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

eta = beta_1 x_1 + beta_2 x_2 + ... + beta_p x_p,

следовательно x_i, не влияет на распределение y, только если beta_i
равна нулю.

* Распределение y имеет форму

f_Y(y; mu, phi) = exp((A/phi) * (y lambda(mu) - gamma(lambda(mu))) +
tau(y, phi))

где phi - это параметр масштаба (возможно известный), и является
константой для всех наблюдений, А представляет собой предварительный
вес, предположительно известный, но возможно варьирующие с разными
наблюдениями, а $\mu$ - это среднее y. Поэтому предполагается, что
распределение y определяется его средней и, возможно, параметром масштаба.

* Средняя, mu, является гладкой обратимой функцией от линейного
предиктора:

mu = m(eta), eta = m^{-1}(mu) = ell(mu)

и эта обратная функция, ell(), называется связывающей функцией.

Эти предположения достаточно свободны, чтобы охватить широкий класс
моделей полезных в статистической практике, но одновременно достаточно
жесткие чтобы разработать по крайней мере приблизительно единую
методологию оценки и прогноза. Читателя можно направить к любой из работ
по данной тематике за более подробной информации, например, McCullagh &
Nelder (1989) или Dobson (1990).

11.6.1 Семейства

Класс обобщенных линейных моделей обрабатывается инструментарием,
представленным R, в него входят gaussian, binomial, poisson, inverse
gaussian и gamma response распределения, а также квази-правдоподобная
модель, где распределение отклика прямо не задается. В последнем случае
дисперсионная функция должна быть определена как функция среднего,
но в других случаях эта функция вводится как распределение отклика.

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

Имя семейства Функция связи

binomial logit, probit, log, cloglog

gaussian identity, log, inverse

Gamma identity, inverse, log

inverse.gaussian 1/mu^2, identity, inverse, log

poisson identity, log, sqrt

quasi logit, probit, cloglog, identity, inverse, log, 1/mu^2, sqrt

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

11.6.2 Функция glm()

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

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

> fitted.model <- glm(formula, family=family.generator, data=data.frame)

Единственным новым элементом является family.generator, который является
инструментом, с помощью которого описано семейство. Это имя функции,
создающей список функций и выражений, которые в совокупности определяют
и контролируют модель и процесс оценки. Хотя это может показаться на
первый взгляд довольно сложным, его использование очень простое.

Имена стандартных, имеющихся генераторов семейств приводятся в разделе
"Название семейства" в таблице в разделе "Семьи". Если есть возможность
выбора связи, в скобках в качестве параметра, может быть также добавлено
к семейству название связи. В случае квази-семейства, таким способом
может быть также указана дисперсионная функция.

Несколько примеров прояснит этот процесс. В gaussian семействе

Вызов такой, как

> fm <- glm(y ~ x1 + x2, family = gaussian, data = sales)

получит такой же результат, как

> fm <- lm(y ~ x1+x2, data=sales)

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

binomial семейство

Рассмотрим небольшой искусственный пример, Silvey (1970).

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

Возраст: 20 35 45 55 70

Количество испытаний: 50 50 50 50 50

Количество слепых: 6 17 26 37 44

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

Если y это число слепых в возрасте х и n число проверенных, обе модели
имеют вид y ~ B(n, F(beta_0 + beta_1 x)), где для пробит случая F(z)
= Phi(z) является стандартной нормальной функцией распределения, а в
логит случае (по умолчанию), F(z) = e^z/(1+e^z). В обоих случаях LD50
является LD50 = - beta_0/beta_1, то есть той точки, в которой аргумент
функции распределения равен нулю.

Первый шаг - это ввод данных, в виде таблицы данных

> kalythos <- data.frame(x = c(20,35,45,55,70), n = rep(50,5), y =
c(6,17,26,37,44))

Чтобы подогнать binomial модель с использованием glm() есть три варианта
ввода отклика:

* Если отклик это вектор, то предполагается что в нем двоичные данные
и он должен содержать только 0/1.

* Если отклик это двух-колоночная матрица предполагается, что первая
колонка содержит число успешных исходов, а второй содержит число неудач.

* Если отклик является фактором, его первый уровень - воспринимается как
"проигрыш" (0) и все другие уровни как "успех" (1).

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

> kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)

Чтобы подогнать модели мы используем

> fmp <- glm(Ymat ~ x, family = binomial(link=probit), data = kalythos)

> fml <- glm(Ymat ~ x, family = binomial, data = kalythos)

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

> summary(fmp)

> summary(fml)

Обе модели соответствуют хорошо (даже слишком). Чтобы найти LD50 оценки
мы можем использовать простую функцию:

> ld50 <- function(b) -b[1]/b[2]

> ldp <- ld50(coef(fmp)); ldl <- ld50(coef(fml)); c(ldp, ldl)

Фактическая оценка по этим данным дает 43,663 лет и 43,601 лет
соответсвенно.

Пуассон модели

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

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

> fmod <- glm(y ~ A + B + x, family = poisson(link=sqrt), data =
worm.counts)

Квази-максимумправдоподобные модели

Для всех семейств дисперсия отклика будет зависеть от средней и будет
иметь параметр масштаба в качестве множителя. Форма зависимость дисперсии
от средней является характеристикой распределения отклика, например для
распределения Пуассона Var(у) = mu.

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

Например, подгонку нелинейной регрессии y = theta_1 z_1 / (z_2 - theta_2)
+ e, которая может быть записана альтернативно как y = 1 / (beta_1 x_1
+ beta_2 x_2) + e, где x_1 = z_2/z_1, x_2 = -1/z_1, beta_1 = 1/theta_1
и beta_2 = theta_2/theta_1. Предположим что подходящая таблица данных
будет создана, мы могли бы подогнать эту нелинейную регрессию, как

> nlfit <- glm(y ~ x1 + x2 - 1, family = quasi(link=inverse,
variance=constant), data = biochem)

Читатель может обратится к руководству и встроенной документации помощи
за дополнительной информацией, если это необходимо.

11.7 Нелинейные модели, методы наименьших квадратов и максимума
правдоподобия

Некоторые формы нелинейных моделей могут быть подогнаны в Обобщенной
Линейной Моделью (glm()). Но в большинстве случаев мы производим подгонку
нелинейной кривой одним из методов нелинейной оптимизации. Процедуры R в
нелинейной оптимизации это optim(), nlm() и (начиная с R 2.2.0) nlminb(),
которые обеспечивают (и превосходят) функциональность S-Plus-вых - ms()
и nlminb(). Мы подгоняем значения параметров, которые позволяют свести
к минимуму некоторые индекс отклонения-от-цели, и делаем это, пробуя
пошагово различные параметры. В отличие например от линейной регрессии,
нет никаких гарантий, что процедура сойдется к удовлетворяющей оценке. Все
методы требуют первоначальные предположения о значениях параметров
подгонки, и сходимость может зависеть критически от качества выбора
исходных значений параметров.

11.7.1 Метод наименьших квадратов

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

Вот пример из Bates & Watts (1988), стр. 51. Данные:

> x <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.56, 0.56,
1.10, 1.10)

> y <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200)

Моделью для подгонки является:

> fn <- function(p) sum((y - (p[1] * x)/(p[2] + x))^2)

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

> plot(x, y)

> xfit <- seq(.02, 1.1, .05)

> yfit <- 200 * xfit/(0.1 + xfit)

> lines(spline(xfit, yfit))

Мы могли бы выбрать удачнее, но начальные значения 200 и 0,1 выглядят
подходящими. Теперь сделем подгонку:

> out <- nlm(fn, p = c(200, 0.1), hessian = TRUE)

После подгонки, out$minimum это SSE, а out$estimate является оценкой
параметров методом наименьших квадратов. Чтобы получить оценку стандартных
ошибок (SE) по результату подгонки, мы делаем:

> sqrt(diag(2*out$minimum/(length(y) - 2) * solve(out$hessian)))

2 в строке выше, это количество параметров. А 95% доверительный интервал
будет оценен параметром +/- 1,96 SE. Мы можем наложить подгонку методом
наименьших квадратов на следующем графике:

> plot(x, y)

> xfit <- seq(.02, 1.1, .05)

> yfit <- 212.68384222 * xfit/(0.06412146 + xfit)

> lines(spline(xfit, yfit))

Стандартный пакет stats предоставляет гораздо более широкие возможности
для подгонки нелинейных моделей методом наименьших квадратов. Эта модель
которую мы только что подогнали является моделью Michaelis-Menten,
поэтому мы можем использовать

> df <- data.frame(x=x, y=y)

> fit <- nls(y ~ SSmicmen(x, Vm, K), df)

> fit

Nonlinear regression model

model: y ~ SSmicmen(x, Vm, K)

data: df

Vm K

212.68370711 0.06412123

residual sum-of-squares: 1195.449

> summary(fit)

Formula: y ~ SSmicmen(x, Vm, K)

Parameters:

Estimate Std. Error t value Pr(>|t|)

Vm 2.127e+02 6.947e+00 30.615 3.24e-11

K 6.412e-02 8.281e-03 7.743 1.57e-05

Residual standard error: 10.93 on 10 degrees of freedom

Correlation of Parameter Estimates:

Vm

K 0.7651

11.7.2 Метод максимального правдоподобия

Максимальное правдоподобие - это метод подгонки нелинейной модели
который применяется даже если ошибки не являются нормальными. Этот
метод находит значения параметров, которые максимизируют логарифмическую
функцию правдоподобия, или что эквивалентно минимизируют отрицательную
логарифмическую функцию правдоподобия. Вот пример из Dobson (1990),
стр. 108-111. Этот пример подгоняет логистическую модель доза-реакция
данных, в которые четко можно подогнать с помощью glm(). Данные:

> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)

> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)

> n <- c(59, 60, 62, 56, 63, 59, 62, 60)

Отрицательная логарифмическая функция правдоподобия для минимизации это:

> fn <- function(p) sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))
+ log(choose(n, y)) ))

Мы выбираем разумное стартовое значение и делаем подгонку:

> out <- nlm(fn, p = c(-50,20), hessian = TRUE)

После подгонки, out$minimum это отрицательная логарифмическая функция
правдоподобия, а out$estimate является оценкой параметров методом
максимального правдоподобия. Чтобы получить оценку SE по произведенной
подгонке, мы делаем:

> sqrt(diag(solve(out$hessian)))

А 95% доверительный интервал будет оценкой параметра +/- 1,96 SE.

11.8 Некоторые нестандартные модели

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

* Смешанные модели. Рекомендованный nlme пакет предоставляет функции
lme() и nlme() для линейных и нелинейных моделей со смешанными эффектами,
которые является линейными и нелинейными регрессиями, в которых некоторые
коэффициенты соответствуют случайным эффектам. Эти функции интенсивно
используют формулы для описания моделей.

* Локально подогнанная регрессия. Функция loess() подгоняет
непараметрическую регрессию с помощью локальной взвешенной
регрессии. Такие регрессии являются полезными для проявления тенденции
неочевидных данных или для редукции данных для обзора большого набора
данных.

Функция loess входит в стандартный пакет stats, включая код для прогноза
следящей регрессии.

* Робастная регрессия. Есть несколько функций для подгонки регрессионных
моделей устойчивым образом к воздействиям экстремальных выбросов
в данных. Функция lqs в рекомендованном пакете MASS обеспечивает
современные алгоритмы для подгонки с высокой устойчивостью. Имеются в
пакетах менее устойчивые, но статистически более эффективные методы,
например, функция rlm в пакете MASS.

* Аддитивные модели. Этот метод направлен на строительство регрессионной
функции из гладких аддитивных функций независимых переменных, как правило,
по одной для каждой независимой переменной. Функции avas и ace в пакете
acepack и функции bruto и mars в пакете mda предусматривает несколько
примеров этой методики в созданных пользователями пакетах R. Расширением
является Обобщенная аддитивная модель, которая реализована в созданных
пользователями пакетах gam и mgcv.

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

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

Древовидные модели доступны в R через сделанные пользователями пакеты
rpart и tree.

12 Графические функции

Графические средства являются важным и экстремально приспособляемой
частью среды R. Можно использовать эти возможности для показа широкого
спектра статистических графиков, а также что бы создать совершенно новые
виды графиков.

Графические возможности могут быть использованы в интерактивном и пакетном
режимах, но в большинстве случаев, интерактивное использование более
продуктивно. Интерактивное использование также просто в использовании,
поскольку в момент запуска R инициирует драйвер графического устройства,
который открывает специальное окно для отображения интерактивной
графики. Хотя это делается автоматически, полезно знать, что в UNIX
используется команда X11() и команда windows() в Windows.

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

Команды рисования делятся на три основные группы:

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

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

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

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

Данное руководство описывает только, то что известно как "базовая"
графика. Отдельная подсистема графики из пакета grid сосуществует
с базовой - это более мощная, но более сложная в использовании
система. Существует рекомендуемый пакет lattice которая построена на
grid и обеспечивает способы для получения составных графиков наподобие
графиков из Trellis системы в S.

12.1 Высокоуровневые графические команды

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

12.1.1 Функция plot()

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

plot(x, y)

plot(xy)

Если х и у являются векторами, plot(х, у) выдает корреллограмму y по
х. Тот же эффект может быть произведен заданием одного аргумента (вторая
форма), либо как список, содержащий два элемента х и у, либо матрицей
из двух колонок.

plot(x)

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

plot(f)

plot(f, y)

f является фактором, y является вектором чисел. Первый вариант создает
столбчатую диаграмму по f; второй вариант производит коробчатую диаграмму
y для каждого уровня f.

plot(df)

plot(~ expr)

plot(y ~ expr)

df это таблица данных, y - любой объект, expr это список имен объектов,
разделенные "+" (например, a + b + c). В первых двух вариантах выводится
графики распределения по переменным в таблице данных (первый вариант)
или по числу имен объектов (второй вариант). Третий вариант рисует y по
каждому объекту, указанному в expr.

12.1.2 Отображение многомерных данных

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

> pairs(X)

производит матрицу попарных корреллограмм на переменных заданных
в столбцах Х, то есть каждый столбец Х отображается по всем другим
столбцам X и результирующие n(n-1) графиков расположены в виде матрицы
с постоянным масштабом шкал графиков по строках и столбцам этой матрицы.

Когда рассматриваются три или четыре переменных coplot может быть более
полезным. Если a и b вектора чисел и с вектор чисел или фактор-объект
(все одинаковой длинны), тогда команда

> coplot(a ~ b | c)

производит ряд корреллограмм a по b для заданных значений c. Если c
представляет собой фактор, это просто означает, что a выводится по b на
каждом уровне c. Когда c численный вектор, оно делится на ряд непрерывных
интервалов и для каждого интервала a отображается по b для значений c из
данного интервала. Количество и положение интервалов можно контролировать
с помощью given.values=, аргумента coplot() - функция co.intervals ()
полезна для выбора интервалов. Вы также можете использовать два переменных
в команде, например

> coplot(a ~ b | c + d)

, который производит корреллограмму a на b для каждого совместный
интервала группировки по c и d.

Функции coplot() и pairs() обе принимают аргументом panel=, который может
быть использован для настройки типа графика, который появится в каждой
из панелей. По умолчанию применяется points() для вывода корреллограммы,
но путем указания некоторых других низкоуровневых графических функций
двух векторов х и у как значения panel= можно выводить любой нужный
тип графиков. Примером panel функции полезной для coplots является
panel.smooth().

12.1.3 Графический вывод

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

qqnorm(x)

qqline(x)

qqplot(x, y)

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

hist(x)

hist(x, nclass=n)

hist(x, breaks=b, ...)

Выводит гистограмму вектора чисел х. Как правило выбирается разумное
количество интервалов группировки, но можно их задать прямо аргументом
nclass=. Кроме того, точки разбиения можно указать точно аргументом
breaks=. Если задан аргумент probability=TRUE, столбцы представляют
вместо сумм относительные частоты. dotchart(x, ...)

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

image(x, y, z, ...)

contour(x, y, z, ...)

persp(x, y, z, ...)

Графики трех переменных. График image рисует сетку прямоугольников,
используя различные цвета для отображения величины z, график contour
рисует горизонтали представляющие уровни z, и график persp рисует 3D
поверхность.

12.1.4 Параметры высокоуровневых графических функций

Есть ряд параметров, которые могут быть переданы высокоуровневым
графическим функциям, такие как:

add=TRUE

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

axes=FALSE

Подавляет вывод осей - полезно чтобы добавить собственные осей с помощью
функции axis(). По умолчанию, axes=TRUE, что означает отображение осей.

log="x"

log="y"

log="xy"

Приводит х, у или обе оси к логарифмическому масштабу. Это сработает
для многих, но не всех видов графиков.

type=

type= аргумент контролирует тип выводимого графика, а именно:


type="p"

Выводит отдельные точек (принято по умолчанию)

type="l"

Рисует линии

type="b"

Выводит точки, соединенные линиями (вместе)

type="o"

Рисует точки перекрытые линиями

type="h"

Рисует вертикальные линии от точек до нулевого уровня оси (график
high-density)

type="s"

type="S"

Графическая ступенчатая функция. В первом варианте, верх вертикали
определяет точу, а во втором - низ.

type="n"

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

xlab = строка

ylab = строка

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

main=строка

Название рисунка большим шрифтом, расположенное в верхней части рисунка.

sub=строка

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

12.2 Низкоуровневые графические команды

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

Некоторые из наиболее полезных низкоуровневых графических функций:

points(x, y)

lines(x, y)

Добавляет точки или связывающие линии к текущему графику. для plot()
аргументы типа type= также могут быть переданы в эти функции (и по
умолчанию "р" для points() и "l" lines().)

text(x, y, labels, ...)

Добавить текст в рисунок в точке заданных х, у. Обычно метки -
целочисленные или символьные вектора, в этом случае labels[i] - выводится
в точке (x[i], y[i]). Значение по умолчанию - 1:length(x).

Примечание: Эта функция часто используется в сочетании

> plot(x, y, type="n"); text(x, y, names)

Графический параметр type="n" подавляет точки, но создает оси, а функции
text() поставляет специальные символы, которые заданы именами в символьном
векторе для каждой из точек.

abline(a, b)

abline(h=y)

abline(v=x)

abline(lm.obj)

Добавляет в текущий рисунок линию наклона b и отсекающую отрезок a. h=y
может быть использовано для задания y-координаты для высоты горизонтальной
линии проходящей через рисунок, а v=x аналогично для х-координаты задает
вертикальную линию. Также lm.obj может быть списком с компонентами
коэффициентами с length 2 (например результат функции подгонки модели),
которые используются в качестве отсекающего отрезка и наклона, именно
в таком порядке.

polygon(x, y, ...)

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

legend(x, y, legend, ...)

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

legend( , fill=v)

Цвета для заполнения прямоугольников

legend( , col=v)

Цвета, которыми рисуются точки или линии

legend( , lty=v)

Стиль линии

legend( , lwd=v)

Ширина линии

legend( , pch=v)

Выводимые символы (символьный вектор)

title(main, sub)

Большим шрифтом добавляет заголовок main в начало текущего рисунка и
(дополнительно) мелким шрифтом подзаголовок sub внизу.

axis(side, ...)

Добавить оси к текущему рисунку со стороны заданной первым аргументом (от
1 до 4, считая по часовой стрелке снизу.) Другие аргументы контролируют
расположение оси внутри или рядом с рисунком, и позиции отметок и
меток. Полезно для добавления в дальнейшем пользовательских осей вызывать
plot() с аргументом axes=FALSE.

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

Если х и у аргументы обязательны, то достаточно задать один аргумент из
списка с элементами по имени x и y. Аналогично матрица с двумя столбцами
также подходит для ввода. Этим путем функции, такие как locator()
(см. ниже) могут использоваться для интерактивного определения позиции
на рисунке.

12.2.1 Математические формулы

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

> text(x, y, expression(paste(bgroup("(", atop(n, x), ")"), p^x,
q^{n-x})))

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

> help(plotmath)

> example(plotmath)

> demo(plotmath)

12.2.2 Векторные шрифты Hershey

Можно задать векторные шрифты Hershey для рендеринга текста, когда
используются текстовые и контурные функций. Есть три причины для
использования шрифтов Hershey:

* Hershey шрифты могут дать результаты лучше, особенно на экране
компьютера, для ротированного и/или мелкого текста.

* Hershey шрифты предоставляют некоторые символы, которых нет в
стандартных шрифтах. В частности, имеются знаки зодиака, картографические
и астрономические символы.

* Hershey шрифты предоставляют кириллические и японские (Кандзи и Кана)
символы.

Более подробная информация, включая таблицы символов Hershey можно
получить в R с помощью команд:

> help(Hershey)

> demo(Hershey)

> help(Japanese)

> demo(Japanese)

12.3 Интерактивная графика

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

locator(n, type)

Ожидает выбор пользователя локации на текущем рисунке, с использованием
левой кнопки мыши. Это продолжается до тех пор, пока не будет выбрано
n (по умолчанию 512) точек, или нажата другая кнопка мыши. Аргумент
type применим для вывода отмеченных точек и действует так же,
как на высокоуровневые графические команды; по умолчанию вывод
отсутствует. locator() возвращает координаты отдельных точек как список
с двумя компонентами x и y.

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

> text(locator(1), "Outlier", adj=0)

(locator() будет проигнорировано, если текущее устройство, например
postscript не поддерживает интерактивного взаимодействия.)

identify(x, y, labels)

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

Иногда мы хотим идентифицировать избранные точки на графике, а не их
координаты. Например, мы можем пожелать что бы пользователь выбрал
некоторые интересующие наблюдения на графическом дисплее, а затем
обработать эти наблюдения разными способами. Задав ряд (х, у) координат
двумя численными векторами х и у, мы могли бы использовать функцию
identify() следующим образом:

> plot(x, y)

> identify(x, y)

Функция identify() сама не рисует, но просто позволяет пользователю
переместить курсор мыши и нажать левую кнопку мыши вблизи точки. Если
возле указателя мыши есть точка, она будет помечена своим индексом
(то есть, своей позицией в х/y векторах) выведенным поблизости. Можно
также использовать некоторые информационные строки (например, название
случая), в качестве выделения, используя labels аргумент для identify(),
или вообще отключить выделение аргументом plot = FALSE. Когда этот процесс
прекращается (см. выше), identify() возвращает индексы выбранных точек;
вы можете использовать эти индексы для извлечения отдельных точек из
первоначальных векторов х и y.

12.4 Использование графических параметров

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

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

12.4.1 Постоянные изменения параметров: функция par()

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

par()

Без параметров, возвращает список всех графических параметров и их
значения для текущего устройства. par(c("col", "lty"))

Если аргумент символьный вектор, возвращает только названные графические
параметры (тоже как список.) par(col=4, lty=2)

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

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

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

> oldpar <- par(col=4, lty=2)

... команды вывода ...

> par(oldpar)

Чтобы сохранить и восстановить все настраиваемые [сноска23] графические
параметры используем

> oldpar <- par( no.readonly=TRUE)

... команды вывода ...

> par(oldpar)