Факторы

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

Факторы

В статистике существует раздение на качественые и количественные переменные
* Качественные факторы определяются факторами
* Фактор - это гибрид целочиесленного (integer) и строкового (character) типа вектора

Чтобы лучше понять структура факторов, приведу пример того, как они создаются:

levels <- sort(unique(x))
f <- match(x, levels)
levels(f) <- as.character(levels)
class(f) <- "factor"

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

x <- factor(sample(1:100, 30, replace = TRUE))
x

 [1] 20 33 21 74 71 99 53 51 24 62 32 99 25 20 23 48 89 87 65 63 61 7  19 68 66 3  41 19 86 7 
Levels: 3 7 19 20 21 23 24 25 32 33 41 48 51 53 61 62 63 65 66 68 71 74 86 87 89 99

Объединить уровни факторы в один можно при помощи конструкции вида:

levels(x)[1:2] <- levels(x)[1]
x

 [1] 20 33 21 74 71 99 53 51 24 62 32 99 25 20 23 48 89 87 65 63 61 3  19 68 66 3  41 19 86 3 
Levels: 3 19 20 21 23 24 25 32 33 41 48 51 53 61 62 63 65 66 68 71 74 86 87 89 99

RNGkind(sample.kind = "Rounding")

non-uniform 'Rounding' sampler used

set.seed(1337)
f <- factor(sample(LETTERS, 30, replace = TRUE))

f

 [1] O O B L J I Y H G D Z Z V F Z A Z Y I G W S B D O Z Y X M Z
Levels: A B D F G H I J L M O S V W X Y Z

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

length(f)

[1] 30

Фактор - гибрид. Проверим

is.vector(f)

[1] FALSE

is.numeric(f)

[1] FALSE

is.character(f)

[1] FALSE

По формальным признакам ничего гибридного.
Но! проверим, как работает приведение типов:

as.numeric(f)

 [1] 11 11  2  9  8  7 16  6  5  3 17 17 13  4 17  1 17 16  7  5 14 12  2  3 11 17 16 15 10 17

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

as.character(f)

 [1] "O" "O" "B" "L" "J" "I" "Y" "H" "G" "D" "Z" "Z" "V" "F" "Z" "A" "Z" "Y" "I" "G" "W" "S" "B" "D"
[25] "O" "Z" "Y" "X" "M" "Z"

f[17]

[1] Z
Levels: A B D F G H I J L M O S V W X Y Z

Видим, что у каждого целого есть свое некое подставление - метки буквы латинского алфавита в данном случае

# Посмотреть не все элементы фактора ,а уникальные градации фактора
levels(f)

 [1] "A" "B" "D" "F" "G" "H" "I" "J" "L" "M" "O" "S" "V" "W" "X" "Y" "Z"

# то же, что и length для вектора
nlevels(f)

[1] 17

Уровни фактора

Фактор как наследник вектора имеет свои правила индексирования. НО! для обычного фактора определены только обозначения “==”, “!=”

Поэтому к элементам фактора можно обращаться, как к элементам вектора, и не только обращаться, но и производить замену

# исключаем уровень А
# Заменяем А на Z
f[f == "A"] <- "Z"
f

 [1] O O B L J I Y H G D Z Z V F Z Z Z Y I G W S B D O Z Y X M Z
Levels: A B D F G H I J L M O S V W X Y Z

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

# заключение в скобки здесь позволяет сразу вывести результат работы 
(f <- droplevels(f)) 

 [1] O O B L J I Y H G D Z Z V F Z Z Z Y I G W S B D O Z Y X M Z
Levels: B D F G H I J L M O S V W X Y Z

Преобразование уровней фактора

levels(f) <- tolower(levels(f))
#levels(f) <- letters[LETTERS %in% levels(f)]
f

 [1] o o b l j i y h g d z z v f z z z y i g w s b d o z y x m z
Levels: b d f g h i j l m o s v w x y z

Изменились уровни, хотя соответствие осталось таким же

levels(f)[1] <- "bbb"
f

 [1] o   o   bbb l   j   i   y   h   g  
[10] d   z   z   v   f   z   z   z   y  
[19] i   g   w   s   bbb d   o   z   y  
[28] x   m   z  
16 Levels: bbb d f g h i j l m o s ... z

При этом не нужно применять функцию drop levels

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

• Иногда категории обладают отношением порядка - такая конструкция называется упорядоченный фактор, т.е. категория качественной переменной упорядочена (например, стадия заболевания)
• Упорядоченный фактор - как создать?
  

1. функция ordered()
   
2. ordered = TRUE для factor

temp <- c("freezing cold", "cold", "comfortable", "hot", "burning hot")
# в ответах переменной неявным образом прослеживается отношение порядка
ft <- ordered(sample(temp, 14, replace = TRUE), temp)
ft

 [1] comfortable   hot           burning hot   cold          burning hot   burning hot  
 [7] freezing cold freezing cold burning hot   comfortable   cold          comfortable  
[13] cold          hot          
Levels: freezing cold < cold < comfortable < hot < burning hot

Появились дополнительные значки порядке в levels. Теперь можно воспользоваться операторами не только == или != (не равно), но и >/</>= и т.п.
Например, хочу отобрать только те ответы, в которых хотя бы горячо, т.е. “hot”

ft[ft >= "hot"]

[1] hot         burning hot burning hot burning hot burning hot hot        
Levels: freezing cold < cold < comfortable < hot < burning hot

Получили ответы, которые расположены справа от категории hot

Преобразование количественной переменной в качетсвенную

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

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

Например, вызовем cut на 10 нормально распределенных числах

a

 [1] -2.3438702  0.5902631  0.3191343  0.7016064 -1.3774571  0.3126459  1.1208476  0.4718309
 [9]  1.1276853  0.6777231

Указываем вектор интервалов разбиений от -5 до 5 (аргумент breaks)

cut(a, -5:5)

 [1] (-3,-2] (0,1]   (0,1]   (0,1]   (-2,-1] (0,1]   (1,2]   (0,1]   (1,2]   (0,1]  
Levels: (-5,-4] (-4,-3] (-3,-2] (-2,-1] (-1,0] (0,1] (1,2] (2,3] (3,4] (4,5]

Если в качестве аргумента breaks укзать единственное число (скаляр), то данные будут разбиты на интервалы равной длины так, чтобы крайние значения не входили в смежные интервалы (а именно границы смещаются на 0,1% интервала)

cut(a, 3)

 [1] (-2.35,-1.19]  (-0.0295,1.13] (-0.0295,1.13] (-0.0295,1.13] (-2.35,-1.19]  (-0.0295,1.13]
 [7] (-0.0295,1.13] (-0.0295,1.13] (-0.0295,1.13] (-0.0295,1.13]
Levels: (-2.35,-1.19] (-1.19,-0.0295] (-0.0295,1.13]

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

table(cut(rnorm(1000), -5:5))

(-5,-4] (-4,-3] (-3,-2] (-2,-1]  (-1,0]   (0,1]   (1,2]   (2,3]   (3,4]   (4,5] 
      0       0      21     135     340     339     140      21       4       0 

Так мы можем получить элементарную осмысленную статистику - фактически мы получаем нечто похожее на гистограмму, которая в данном случае имеет вид нормального распределния (но мы и вызывали функцию rnorm, так что это ожидаемо)

Важное отступление - функция options

У сессии R есть набор активных настроек, отвечающих за подсчет и вывод результатов вычислений
*?optinons

?options

** digits - количество знаков при печати чисел
** error - повдеение при ошибке
** width - длина строки при печати векторов и матриц
Полный список достаточно большой, но есть ряд достаточно важных

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

Отменить такое поведение можно вызовом
options(stringsAsFactors = FALSE)

Функция tapply

Мы уже не раз поднимали фукнции семейства apply, есть такая и для работы с факторами.
Дело в том, что факторы сами по себе несут мало смысла, они несут смысл как элементу дата фрейма. Одна из наиболее распространенных задач в анализе данных - подсчет некой статистики по группам, которые определяют факторы. Как раз для этих задач удобно пользоваться функцией tapply.
Разберем на примере массива warpbreaks (стародревний массив с данными о поломках прядильных установок в зависимости от типа волоса, с которым работал, и силы натяжения)

str(warpbreaks)

'data.frame':   54 obs. of  3 variables:
 $ breaks : num  26 30 54 25 70 52 51 26 67 18 ...
 $ wool   : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ tension: Factor w/ 3 levels "L","M","H": 1 1 1 1 1 1 1 1 1 2 ...

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

tapply(warpbreaks$breaks, warpbreaks$wool, max)

 A  B 
70 44 

Глоссарий

?facotr
?levels, ?nlevels, ?droplevels
?ordered
?cut, ?table
?options
?tapply

Задача

Используйте показанную мной связку из двух функций, чтобы превратить количественную переменную mag (сила землетрясения в баллах по шкале Рихтера) дата фрейма quakes в качественную. Интервалы должны быть длиной в полбалла, начиная с минимального, при этом левый конец интервала включается. Теперь отсортируйте общее количество случаев, попавших в каждую категорию, в порядке убывания.
1. Посмотрю структуру дата-фрейма

str(quakes)

'data.frame':   1000 obs. of  5 variables:
 $ lat     : num  -20.4 -20.6 -26 -18 -20.4 ...
 $ long    : num  182 181 184 182 182 ...
 $ depth   : int  562 650 42 626 649 195 82 194 211 622 ...
 $ mag     : num  4.8 4.2 5.4 4.1 4 4 4.8 4.4 4.7 4.3 ...
 $ stations: int  41 15 43 19 11 12 43 15 35 19 ...

2. Внесу дата фрейм в переменную, чтобы работать с ним

q <- quakes
summary(q$mag)

   Min. 1st Qu.  Median    Mean 3rd Qu. 
   4.00    4.30    4.60    4.62    4.90 
   Max. 
   6.40 

3. Сгенерирую последовательность от min(q\(mag)* до *max(q\)mag) с шагом 0.5

?ceiling
#ceiling, чтобы включить верхнюю границу
a <- seq(min(q$mag), ceiling(max(q$mag)), by = 0.5)
a

[1] 4.0 4.5 5.0 5.5 6.0 6.5 7.0

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

sort(table(cut(q$mag, br = a, right = FALSE)))

[6.5,7) [6,6.5) [5.5,6) [5,5.5) [4,4.5) [4.5,5) 
      0       5      33     160     377     425 

?cut ?table

Подсказки:

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

?cut

Вернемся к массиву avianHabitat или Дата Фрейм - часть 2

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

str(avian)

'data.frame':   1070 obs. of  17 variables:
 $ Site    : chr  "BunkerHill27" "BunkerHill27" "BunkerHill27" "BunkerHill27" ...
 $ Observer: chr  "RA" "RA" "RA" "RA" ...
 $ Subpoint: int  1 2 3 4 5 6 7 8 9 10 ...
 $ VOR     : num  6 4.5 2 2.5 4 2 5.5 4 3.5 3.5 ...
 $ PDB     : int  3 2 4 3 4 3 3 2 2 2 ...
 $ DBHt    : num  5.2 3.1 5.5 6.2 5.4 4 5.2 4.4 5.7 4.8 ...
 $ PW      : int  0 3 1 0 0 0 2 1 1 0 ...
 $ WHt     : num  0 4.7 5.8 0 0 0 6.3 4.1 5.7 0 ...
 $ PE      : int  4 3 3 3 3 3 2 2 2 1 ...
 $ EHt     : num  2.9 4.1 3.9 4 3.5 4.1 2.6 4.3 5.2 1.7 ...
 $ PA      : int  0 0 0 0 0 0 0 0 0 0 ...
 $ AHt     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ PH      : int  4 3 3 4 4 2 4 5 4 5 ...
 $ HHt     : num  3 3.5 7.5 5 3.7 3.5 5.8 8.2 6.9 5.7 ...
 $ PL      : int  0 2 0 0 0 0 0 0 0 0 ...
 $ LHt     : num  0 1 0 0 0 0 0 0 0 0 ...
 $ PB      : int  0 0 0 0 0 0 0 0 0 0 ...

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

getwd()

[1] "C:/Users/Tony/Desktop/R/Mine"

Ок, я в нужной директории. Из-под Линукса, конечно, работать с диеркториями было бы удобнее. Но так или иначе к файлам из рабочей директории можно обращаться, не указывая полный путь, а просто по имени, т.е.

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

source()

Error in source() : argument "file" is missing, with no default

Удостоверюсь, что массив записан в дата фрейма avian существует

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

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

list.files()

 [1] "avianHabitat.csv"        
 [2] "Data frame.nb.html"      
 [3] "Data frame.Rmd"          
 [4] "Data Intro.nb.html"      
 [5] "Data Intro.Rmd"          
 [6] "Factors.nb.html"         
 [7] "Factors.Rmd"             
 [8] "Intro.nb.html"           
 [9] "Intro.Rmd"               
[10] "Mtrx-and-Lsts---2.html"  
[11] "Mtrx and Lsts - 2.Rmd"   
[12] "Mtrx and Lsts.nb.html"   
[13] "Mtrx and Lsts.Rmd"       
[14] "rsconnect"               
[15] "Strings.nb.html"         
[16] "Strings.Rmd"             
[17] "Vector.nb.html"          
[18] "Vector.Rmd"              
[19] "Матрицы и списки.nb.html"

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

list.files(pattern = ".csv$")

[1] "avianHabitat.csv"

Можно сделать и более хитро, хотя в данном случае большого смысла в этом нет, например:

list.files(pattern = ".*\\.csv$")

[1] "avianHabitat.csv"

# где "." - любой символ, т.е. неважно с какого символа начинается слова
# "*" - он же, т.е. любой символ, но повторенный любое количество раз
# ".*" - значит мне не важно начало
# "\\." - это обычная точка, но чтобы она считывалась как wildcard приходится экранировать двумя бэкслешами
#"\\.csv$" - заканчивается на .csv

Но допустим, этот найденный нам нужный файл по-прежнему очень большой. Открыть его ничем не получается. Что делать?
В такой ситуации можно воспользоваться низкоуровневой функцией readLines, указав количество строк к прочтению

readLines("avianHabitat.csv", 5)

[1] "Site,Observer,Subpoint,VOR,PDB,DBHt,PW,WHt,PE,EHt,PA,AHt,PH,HHt,PL,LHt,PB"
[2] "BunkerHill27,RA,1,6,3,5.2,0,0,4,2.9,0,0,4,3,0,0,0"                        
[3] "BunkerHill27,RA,2,4.5,2,3.1,3,4.7,3,4.1,0,0,3,3.5,2,1,0"                  
[4] "BunkerHill27,RA,3,2,4,5.5,1,5.8,3,3.9,0,0,3,7.5,0,0,0"                    
[5] "BunkerHill27,RA,4,2.5,3,6.2,0,0,3,4,0,0,4,5,0,0,0"                        

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

Результат точно такой же

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

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

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

Теперь вижу, что все правильно

Вообще самая частая проблема, когда наоборот по умолчанию при считывании дата фрейма переменные строкового типа character считаются факторами (factor). Это происходит, т.к. по умолчнаию в сессии R установлена такая опция:
#options(stringsAsFactors = TRUE)
Чтобы этого не происходило перед считыванием дата фрейма лучше менять ее на FALSE, а потом уже вручную строки преобразовывать в факторы (которых обычно заведомо меньше)

Проверка диапазона значений

Возвращаемся к вопросу проверки диапазона значений. Раньше мы это делали по сути вручную

check_percent_range <- function(x) {
  any(x < 0 | x > 100)
}

check_percent_range(avian$PW)

[1] FALSE

У нас уже есть функция, написанная ранее, но что делать, когда переменная не одна, а огромное множество, и все их перепроверять вручную довольно утомительно. Автоматически все это можно сделать с помощью функции семейства apply.
Для начала пересчитаем % total coverage с помощью пакета stringr

library(stringr)

package 㤼㸱stringr㤼㸲 was built under R version 4.0.2

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

#coverage_variables <- names(avian)[-(1:4)][c(T,F)]
#avian$total_coverage <- rowSums(avian[ , coverage_variables])
coverage_variables <- names(avian)[str_detect(names(avian), "^P")]
#где в str_detect мы закинули все имена переменных, начинающихся с буквы P (регулярное выражение "^" - значит начало строки)
coverage_variables

[1] "PDB" "PW"  "PE"  "PA"  "PH"  "PL"  "PB" 

Теперь провернем то, ради чего все это затевалось

sapply(coverage_variables, function(name) check_percent_range(avian[[name]]))

  PDB    PW    PE    PA    PH    PL    PB 
FALSE FALSE FALSE FALSE FALSE FALSE FALSE 

#где function(name) - указание о том, что будем применять анонимную функцию, после которого указываем саму функцию
##[[name]], т.к. функция на вход принимает вектор, поэтому и извлекаем вектор из дата фрейма, а не строку

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

Посмотрим на уникальные названия наблюдений под переменной Site

unique(avian$Site)

  [1] "BunkerHill27"       "BunkerHill26"       "BunkerHill25"      
  [4] "BunkerHill12"       "BunkerHill13"       "BunkerHill14"      
  [7] "BunkerHill15"       "BunkerHill16"       "BunkerHill17"      
 [10] "BunkerHill24"       "BunkerHill23"       "BunkerHill22"      
 [13] "BunkerHill21"       "BunkerHill11"       "CreteCreek18"      
 [16] "CreteCreek16"       "CreteCreek14"       "CreteCreek12"      
 [19] "CreteCreek24"       "CreteCreek22"       "CreteCreek120"     
 [22] "CreteCreek118"      "CreteCreek116"      "CreteCreek114"     
 [25] "CreteCreek112"      "CreteCreek110"      "CreteCreek21"      
 [28] "CreteCreek23"       "CreteCreek29"       "CreteCreek28"      
 [31] "CreteCreek25"       "CreteCreek26"       "CreteCreek27"      
 [34] "HortonCreek16"      "HortonCreek15"      "HortonCreek14"     
 [37] "HortonCreek13"      "HortonCreek12"      "HortonCreek11"     
 [40] "HortonCreek17"      "HortonCreek18"      "HortonCreek19"     
 [43] "HortonCreek110"     "HortonCreek111"     "HortonCreek113"    
 [46] "HortonCreek114"     "HortonCreek115"     "HortonCreek116"    
 [49] "HortonCreek126"     "HortonCreek125"     "HortonCreek124"    
 [52] "HortonCreek123"     "HortonCreek122"     "HortonCreek121"    
 [55] "HortonCreek120"     "HortonCreek119"     "HortonCreek118"    
 [58] "HortonCreek117"     "HortonCreek112"     "LivingstonCreek112"
 [61] "LivingstonCreek113" "LivingstonCreek114" "LivingstonCreek115"
 [64] "LivingstonCreek116" "LivingstonCreek211" "LivingstonCreek210"
 [67] "LivingstonCreek29"  "LivingstonCreek13"  "LivingstonCreek14" 
 [70] "LivingstonCreek15"  "LivingstonCreek16"  "LivingstonCreek17" 
 [73] "LivingstonCreek28"  "LivingstonCreek27"  "LivingstonCreek26" 
 [76] "LivingstonCreek25"  "LivingstonCreek18"  "LivingstonCreek12" 
 [79] "LivingstonCreek11"  "LivingstonCreek19"  "LivingstonCreek110"
 [82] "LivingstonCreek111" "LivingstonCreek24"  "LivingstonCreek23" 
 [85] "LivingstonCreek22"  "LivingstonCreek21"  "McAdamCreek18"     
 [88] "McAdamCreek17"      "McAdamCreek16"      "McAdamCreek15"     
 [91] "McAdamCreek14"      "McAdamCreek13"      "McAdamCreek12"     
 [94] "McAdamCreek11"      "McAdamCreek27"      "McAdamCreek26"     
 [97] "McAdamCreek25"      "McAdamCreek24"      "McAdamCreek23"     
[100] "McAdamCreek22"      "McAdamCreek21"     

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

#1) где первый аргумент - это переменная
#2) далее регулярное выражение, обозначающее все цифры [:digit:]
#а "+" - wildcard, обозначающает, что символ повторен какое-то количество раз
#т.е. "[:digit:]+" - это любая цифра, повторенная какое-то произвольное количество раз
#3) третий аргумент - на что заменяем (на "", т.е. на пустую строку)
str(avian$site_name)

 chr [1:1070] "BunkerHill" "BunkerHill" "BunkerHill" "BunkerHill" "BunkerHill" ...

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

avian$site_name <- factor(str_replace(avian$site, "[:digit:]", ""))
str(avian$site_name)

 Factor w/ 5 levels "BunkerHill","CreteCreek",..: 1 1 1 1 1 1 1 1 1 1 ...

summary(avian$site_name)

     BunkerHill      CreteCreek 
            140             190 
    HortonCreek LivingstonCreek 
            260             290 
    McAdamCreek 
            190 

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

tapply(avian$DBHt, avian$site_name, mean)

     BunkerHill      CreteCreek 
      2.2578571       0.6168421 
    HortonCreek LivingstonCreek 
      0.3753846       0.4648276 
    McAdamCreek 
      0.9042105 

Это уже некий продвинутый анализ данных.
Можем так же оценить, например, регион в котором покрытие растениями (total_coverage) - наименьший

tapply(avian$total_coverage, avian$site_name, summary)

$BunkerHill
   Min. 1st Qu.  Median    Mean 
  4.000   8.000   9.000   9.236 
3rd Qu.    Max. 
 10.000  13.000 

$CreteCreek
   Min. 1st Qu.  Median    Mean 
   6.00    9.00   10.00   10.02 
3rd Qu.    Max. 
  11.00   14.00 

$HortonCreek
   Min. 1st Qu.  Median    Mean 
   3.00    9.00   10.00   10.14 
3rd Qu.    Max. 
  11.00   16.00 

$LivingstonCreek
   Min. 1st Qu.  Median    Mean 
  2.000   8.000  10.000   9.479 
3rd Qu.    Max. 
 11.000  17.000 

$McAdamCreek
   Min. 1st Qu.  Median    Mean 
  2.000   7.000   9.000   8.626 
3rd Qu.    Max. 
 10.000  15.000 

Теперь можем, например, определить, кто их ученых нашел, какие самые высокие растения из представителей своего вида. Ученых у нас три в переменной Observer, которая записана как фактор (JT, RA, RR). А вот видов растений больше: DB (карликовая берёза), W (ива), E (вереск), A (ольха), H (травяные растения), L (лишайники).
1. Найду переменные с высотой видов

names(avian)[str_detect(names(avian), ".Ht$")]

[1] "DBHt" "WHt"  "EHt"  "AHt" 
[5] "HHt"  "LHt" 

2. Запишу их в переменную

Ht_variables

[1] "DBHt" "WHt"  "EHt"  "AHt"  "HHt"  "LHt" 

3. С помощью sapply найду общую статистику

sapply(Ht_variables, function(name) summary(avian[[name]]))

              DBHt       WHt       EHt        AHt      HHt      LHt
Min.     0.0000000  0.000000 0.0000000  0.0000000 0.000000 0.000000
1st Qu.  0.0000000  0.000000 0.3000000  0.0000000 1.400000 0.000000
Median   0.0000000  0.400000 0.8000000  0.0000000 2.300000 0.200000
Mean     0.7827103  1.026822 0.9658879  0.1859813 2.321028 0.268514
3rd Qu.  1.2000000  1.100000 1.4000000  0.0000000 3.100000 0.500000
Max.    10.0000000 24.500000 5.3000000 31.5000000 8.200000 1.300000

Это можно сделать еще так, получив то же самое. Смысл в том, что функции saplly передаю не просто вектор строк, а имена колонок.

sapply(avian[,Ht_variables], summary)

              DBHt       WHt       EHt        AHt      HHt      LHt
Min.     0.0000000  0.000000 0.0000000  0.0000000 0.000000 0.000000
1st Qu.  0.0000000  0.000000 0.3000000  0.0000000 1.400000 0.000000
Median   0.0000000  0.400000 0.8000000  0.0000000 2.300000 0.200000
Mean     0.7827103  1.026822 0.9658879  0.1859813 2.321028 0.268514
3rd Qu.  1.2000000  1.100000 1.4000000  0.0000000 3.100000 0.500000
Max.    10.0000000 24.500000 5.3000000 31.5000000 8.200000 1.300000

Ладно, пока не понимаю, как получить выкладку через одну команду сразу по всем растениям. Сделаем по порядку:

tapply(avian$DBHt, avian$Observer, summary)

$JT
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.700   1.303   1.950   9.900 

$RA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.0000  0.7671  1.2000 10.0000 

$RR
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.0000  0.5829  1.1000  5.0000 

Выигрывает Ra

tapply(avian$WHt, avian$Observer, summary)

$JT
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.000   1.207   0.800  24.500 

$RA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.500   1.013   1.100  18.500 

$RR
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.500   0.972   1.100  22.000 

Выигрывает JT

tapply(avian$EHt, avian$Observer, summary)

$JT
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.400   0.700   1.221   1.700   5.300 

$RA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.4000  0.8000  0.9564  1.4000  4.9000 

$RR
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.2250  0.8000  0.8708  1.3000  4.2000 

tapply(avian$AHt, avian$Observer, summary)

$JT
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.0000  0.8967  0.0000 31.5000 

$RA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.000   0.112   0.000  19.200 

$RR
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000000 0.000000 0.000000 0.000578 0.000000 0.200000 

tapply(avian$HHt, avian$Observer, summary)

$JT
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.400   2.200   2.669   4.050   8.200 

$RA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.400   2.300   2.244   3.000   7.500 

$RR
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.500   2.300   2.298   3.000   7.300 

tapply(avian$LHt, avian$Observer, summary)

$JT
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.200   0.204   0.400   0.800 

$RA
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.2000  0.2532  0.4000  1.3000 

$RR
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   0.300   0.322   0.500   1.100 

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

max_height <- sapply(avian[,Ht_variables], max)
data.frame(max_height, avian$Observer[max_height])

или так

aggregate(list(avian[,Ht_variables]), by = list(avian$Observer), max)

используется функция aggregate - аналог tapply только для нескольких столбцов

https://stackoverflow.com/questions/7029800/how-to-run-tapply-on-multiple-columns-of-data-frame-using-r

Общий глоссарий для этого урока:

• ?paste (?paste0), ?strsplit

• Regular expressions

• ?grep (?grepl, ?gsub)

• library(stringr): ?str_extract, ?str_replace

• ?tolower, ?toupper

• ?getwd, ?setwd, ?list.files, ?list.dirs

• ?formatC, ?cat

• ?factor

• ?levels, ?nlevels, ?droplevels

• ?ordered

• ?cut, ?table

• ?options

• ?tapply