Около data.table

Тази бележка ще бъде интересна за тези, които използват библиотеката за обработка на таблични данни за R — data.table, и вероятно ще се зарадват да видят гъвкавостта на нейното приложение в различни примери.

Вдъхновени от добър пример колеги, и надявайки се, че вече сте прочели неговата статия, предлагам да се задълбочим в оптимизацията на кода и производителността на основата на data.table.

Въведение: откъде идва data.table?

Най-добре е да започнете запознаването с библиотеката отдалеч, а именно, от структурите данни, от които може да бъде получен обектът data.table (по-нататък, ДТ).

Масив

Код

## arrays ---------

arrmatr <- array(1:20, c(4,5))

class(arrmatr)

typeof(arrmatr)

is.array(arrmatr)

is.matrix(arrmatr)

Една от тези структури е масив (?base::array). Както в други езици, масивите тук са многомерни. Въпреки това, интересно е, че например двуизмерен масив започва да наследява свойства от класа матрица (?base::matrix), а едноизмерен масив, което също е важно, не наследява от вектор (?base::vector).

При това трябва да се разбира, че типът данни, съдържащи се в някакъв обект, следва да се проверява с функцията base::typeof, която връща вътрешно описание на типа съгласно R Internals — общият протокол на езика, свързан с първородния C.

Още една команда за определяне на класа на обекта, base::class, връща във случаи на вектори векторен тип (той се различава по име от вътрешния, но позволява също да се разбере типа на данните).

Списък

От двуизмерния масив, който е и матрица, може да се премине към списък (?base::list).

Код

## lists ------------------

mylist <- as.list(arrmatr)

is.vector(mylist)

is.list(mylist)

При това стават няколко неща наведнъж:

• Второто измерение на матрицата се свива, тоест получаваме едновременно и списък, и вектор.
• Списъкът, така, наследява от тези класове. Трябва да се има предвид, че на елемента на списъка ще отговаря едно (скалярно) значение от клетката на матрицата-обект.

Благодарение на факта, че списъкът е също и вектор, към него могат да се прилагат определени функции за вектори.

Датафрейм

От списък, матрица или вектор може да се премине към датафрейм (?base::data.frame).

Код

## data.frames ------------

df <- as.data.frame(arrmatr)
df2 <- as.data.frame(mylist)

is.list(df)

df$V6 <- df$V1 + df$V2

Интересното в него е: датафреймът наследява от списък! Колоните на датафрейма са клетки на списъка. Това ще бъде важно в бъдеще, когато използваме функции, приложими за списъци.

data.table

Може да се получи ДТ (?data.table::data.table) от датафрейм, списък, вектор или матрица. Например, ето как (in place).

Код

## data.tables -----------------------
library(data.table)

data.table::setDT(df)

is.list(df)

is.data.frame(df)

is.data.table(df)

Полезно е, че, както и датафреймът, ДТ наследява свойства на списъка.

ДТ и памет

В отличие от всички останали обекти в R base, DТ се предават по референция. Ако е необходимо да се направи копие в нова област на паметта, е необходима функция data.table::copy или е необходимо да се направи селекция от стария обект.

Код

df2 <- df

df[V1 == 1, V2 := 999]

data.table::fsetdiff(df, df2)

df2 <- data.table::copy(df)

df[V1 == 2, V2 := 999]

data.table::fsetdiff(df, df2)

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

Някои примери за използване на свойствата на data.table

Като списък…

Итерирането по редовете на датафрейма или DТ не е най-добрата идея, тъй като кодът на цикъла на езика R е значително по-бавен C, но преминаването в цикъл по колоните, които обикновено са значително по-малко, е напълно допустимо. Когато преминаваме през колоните, помним, че всяка колона е елемент от списък, съдържащ, като правило, вектор. А операциите с вектори са добре векторизирани в основните функции на езика. Освен това можем да използваме операторите за избор, типични за списъци и вектори: `[[`, `$`.

Код

## operations on data.tables ------------

#using list properties

df$'V1'[1]

df[['V1']]

df[[1]][1]

sapply(df, class)

sapply(df, function(x) sum(is.na(x)))

Векторизация

Ако е необходимо да преминем през редовете на голям DТ, най-доброто решение ще бъде написването на функция с векторизация. Но ако това не е възможно, трябва да помним, че цикълът вътре DТ е все пак по-бърз от цикъла в R, тъй като се изпълнява на C.

Нека пробваме на по-голям пример с 100К реда. Ще извлечем първата буква от думите, входящи в колоната-вектор w.

Актуализирано

Код

library(magrittr)
library(microbenchmark)

## По-голям пример ----

rown <- 100000

dt %
	.[, d := 1 + b + c + rnorm(nrow(.))]

# векторизация

microbenchmark({
	dt[
		, first_l := unlist(strsplit(w, split = ' ', fixed = T))[1]
		, by = 1:nrow(dt)
	   ]
})

# второ

first_l_f %
		do.call(rbind, .) %>%
		`[`(,1)
}

dt[, first_l := NULL]

microbenchmark({
	dt[
		, first_l := .(first_l_f(w))
		]
})

# трето

first_l_f2 %
		unlist %>%
		matrix(nrow = 3) %>%
		`[`(1,)
}

dt[, first_l := NULL]

microbenchmark({
	dt[
		, first_l := .(first_l_f2(w))
		]
})

Първо изпълнение с итериране по редовете:

Единица: милисекунди
expr min
{ dt[, `:=`(first_l, unlist(strsplit(w, split = " ", fixed = T))[1]), by = 1:nrow(dt)] } 439.6217
lq mean median uq max neval
451.9998 460.1593 456.2505 460.9147 621.4042 100

Втори цикъл, в който векторизацията се извършва чрез преобразуване на списък в матрица и взимане на елементи от сечението с индекс 1 (последното е именно векторизация). Ще се поправя: векторизация на ниво функция strsplit, който може да приема вектор на вход. Оказва се, че процедурата за преобразуване на списък в матрица е много по-тежка от самата векторизация, но и в този случай е значително по-бърза от невекторизирания вариант.

Единица: милисекунди
expr min lq mean median uq max neval
{ dt[, `:=`(first_l, .(first_l_f(w)))] } 93.07916 112.1381 161.9267 149.6863 185.9893 442.5199 100

Ускорение по медиана в 3 пъти.

Трети цикъл, където схемата за преобразуване в матрица е променена.

Единица: милисекунди
expr min lq mean median uq max neval
{ dt[, `:=`(first_l, .(first_l_f2(w)))] } 32.60481 34.13679 40.4544 35.57115 42.11975 222.972 100

Ускорение по медиана в 13 пъти.

С това трябва да се експериментира, колкото повече — толкова по-добре.

Още един пример с векторизация, където също има текст, но той е близък до реални условия: различна дължина на думите, различно количество думи. Необходимо е да се извлекат първите 3 думи. Ето как:

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

Код

# fourth

rown <- 100000

words <-
	sapply(
		seq_len(rown)
		, function(x){
			nwords <- rbinom(1, 10, 0.5)
			paste(
				sapply(
					seq_len(nwords)
					, function(x){
						paste(sample(letters, rbinom(1, 10, 0.5), replace = T), collapse = '')
					}
				)
				, collapse = ' '
			)
		}
	)

dt <- 
	data.table(
		w = words
		, a = sample(letters, rown, replace = T)
		, b = runif(rown, -3, 3)
		, c = runif(rown, -3, 3)
		, e = rnorm(rown)
	) %>%
	.[, d := 1 + b + c + rnorm(nrow(.))]

first_l_f3 <- function(sd, n)
{
	l <- strsplit(sd, split = ' ', fixed = T)
	
	maxl <- max(lengths(l))
	
	sapply(l, "length<-", maxl) %>%
		`[`(n,) %>%
		as.character
}

microbenchmark({
	dt[
		, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
		]
})

dt[
	, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
	]

Единица: милисекунди
expr min lq mean median

{ dt[, `:=`((paste0(«w_», 1:3)), strsplit(w, split = " ", fixed = T))] } 851.7623 916.071 1054.5 1035.199
uq max neval
1178.738 1356.816 100

Скриптът работи със средна скорост от 1 секунда. Не е зле.

Свързани с една верига…

С обектите на DT може да се работи, използвайки свързване. Това изглежда като прикрепване на синтаксис на скобки отдясно, всъщност, захар.

Код

# chaining

res1 <- dt[a == 'a'][sample(.N, 100)]

res2 <- dt[, .N, a][, N]

res3 <- dt[, coefficients(lm(e ~ d))[1], a][, .(letter = a, coef = V1)]

Тече през тръбите…

Подобни операции могат да се извършат чрез piping, изглежда подобно, но функционално е по-богато, тъй като може да се използват всякакви методи, а не само DT. Ще изведем коефициентите на логистичната регресия за нашите синтетични данни с редица филтри за DT.

Код

# piping

samplpe_b <- dt[a %in% head(letters), sample(b, 1)]

res4 <- 
	dt %>%
	.[a %in% head(letters)] %>%
	.[, 
	  {
	  	dt0 <- .SD[1:100]
	  	
	  	quants <- 
	  		dt0[, c] %>%
	  		quantile(seq(0.1, 1, 0.1), na.rm = T)
	  	
	  	.(q = quants)
	  }
	  , .(cond = b > samplpe_b)
	  ] %>%
	glm(
		cond ~ q -1
		, family = binomial(link = "logit")
		, data = .
	) %>%
	summary %>%
	.[[12]]

Статистика, машинно обучение и други в DT

Може да се използват лямбда-функции, но понякога е по-добре да ги създадете отделно, да напишете цялото проучване на данни, и напред — те работят в DT. Примерът е обогатен с всички гореспоменати функции, плюс няколко полезни неща от арсенала на DT (като обращение към самия DT в DT по референция, понякога въведени не последователно, но така че да е удобно).

Код

# function

rm(lm_preds)

lm_preds <- function(
	sd, by, n
)
{
	
	if(
		n < 100 | 
		!by[['a']] %in% head(letters, 4)
	   )
	{
		
		res <-
			list(
				low = NA
				, mean = NA
				, high = NA
				, coefs = NA
			)
		
	} else {

		lmm <- 
			lm(
				d ~ c + b
				, data = sd
			)
		
		preds <- 
			stats::predict.lm(
				lmm
				, sd
				, interval = "prediction"
				)
		
		res <-
			list(
				low = preds[, 2]
				, mean = preds[, 1]
				, high = preds[, 3]
				, coefs = coefficients(lmm)
			)
	}

	res
	
}

res5 <- 
	dt %>%
	.[e < 0] %>%
	.[.[, .I[b > 0]]] %>%
	.[, `:=` (
		low = as.numeric(lm_preds(.SD, .BY, .N)[[1]])
		, mean = as.numeric(lm_preds(.SD, .BY, .N)[[2]])
		, high = as.numeric(lm_preds(.SD, .BY, .N)[[3]])
		, coef_c = as.numeric(lm_preds(.SD, .BY, .N)[[4]][1])
		, coef_b = as.numeric(lm_preds(.SD, .BY, .N)[[4]][2])
		, coef_int = as.numeric(lm_preds(.SD, .BY, .N)[[4]][3])
	)
	, a
	] %>%
	.[!is.na(mean), -'e', with = F]


# plot

plo <- 
	res5 %>%
	ggplot +
	facet_wrap(~ a) +
	geom_ribbon(
		aes(
			x = c * coef_c + b * coef_b + coef_int
			, ymin = low
			, ymax = high
			, fill = a
		)
		, size = 0.1
		, alpha = 0.1
	) +
	geom_point(
		aes(
			x = c * coef_c + b * coef_b + coef_int
			, y = mean
			, color = a
		)
		, size = 1
	) +
	geom_point(
		aes(
			x = c * coef_c + b * coef_b + coef_int
			, y = d
		)
		, size = 1
		, color = 'black'
	) +
	theme_minimal()

print(plo)

Заключение

Надявам се, че успях да създам цялостна, но разбира се, не пълна, картина на такъв обект като data.table, започвайки от неговите свойства, свързани с наследяването от класовете R и завършвайки с неговите уникални функции и околна среда от елементи на tidyverse. Надявам се, че това ще ви помогне по-добре да изучите и прилагате тази библиотека за работа и развлечение.

Спасибо!

Пълен код

Код

## load libs ----------------

library(data.table)
library(ggplot2)
library(magrittr)
library(microbenchmark)


## arrays ---------

arrmatr <- array(1:20, c(4,5))

class(arrmatr)

typeof(arrmatr)

is.array(arrmatr)

is.matrix(arrmatr)


## lists ------------------

mylist <- as.list(arrmatr)

is.vector(mylist)

is.list(mylist)


## data.frames ------------

df <- as.data.frame(arrmatr)

is.list(df)

df$V6 <- df$V1 + df$V2


## data.tables -----------------------

data.table::setDT(df)

is.list(df)

is.data.frame(df)

is.data.table(df)

df2 <- df

df[V1 == 1, V2 := 999]

data.table::fsetdiff(df, df2)

df2 <- data.table::copy(df)

df[V1 == 2, V2 := 999]

data.table::fsetdiff(df, df2)


## operations on data.tables ------------

#using list properties

df$'V1'[1]

df[['V1']]

df[[1]][1]

sapply(df, class)

sapply(df, function(x) sum(is.na(x)))


## Bigger example ----

rown <- 100000

dt <- 
	data.table(
		w = sapply(seq_len(rown), function(x) paste(sample(letters, 3, replace = T), collapse = ' '))
		, a = sample(letters, rown, replace = T)
		, b = runif(rown, -3, 3)
		, c = runif(rown, -3, 3)
		, e = rnorm(rown)
	) %>%
	.[, d := 1 + b + c + rnorm(nrow(.))]

# vectorization

# zero - for loop

microbenchmark({
	for(i in 1:nrow(dt))
		{
		dt[
			i
			, first_l := unlist(strsplit(w, split = ' ', fixed = T))[1]
		]
	}
})

# first

microbenchmark({
	dt[
		, first_l := unlist(strsplit(w, split = ' ', fixed = T))[1]
		, by = 1:nrow(dt)
	   ]
})

# second

first_l_f <- function(sd)
{
	strsplit(sd, split = ' ', fixed = T) %>%
		do.call(rbind, .) %>%
		`[`(,1)
}

dt[, first_l := NULL]

microbenchmark({
	dt[
		, first_l := .(first_l_f(w))
		]
})

# third

first_l_f2 <- function(sd)
{
	strsplit(sd, split = ' ', fixed = T) %>%
		unlist %>%
		matrix(nrow = 3) %>%
		`[`(1,)
}

dt[, first_l := NULL]

microbenchmark({
	dt[
		, first_l := .(first_l_f2(w))
		]
})

# fourth

rown <- 100000

words <-
	sapply(
		seq_len(rown)
		, function(x){
			nwords <- rbinom(1, 10, 0.5)
			paste(
				sapply(
					seq_len(nwords)
					, function(x){
						paste(sample(letters, rbinom(1, 10, 0.5), replace = T), collapse = '')
					}
				)
				, collapse = ' '
			)
		}
	)

dt <- 
	data.table(
		w = words
		, a = sample(letters, rown, replace = T)
		, b = runif(rown, -3, 3)
		, c = runif(rown, -3, 3)
		, e = rnorm(rown)
	) %>%
	.[, d := 1 + b + c + rnorm(nrow(.))]

first_l_f3 <- function(sd, n)
{
	l <- strsplit(sd, split = ' ', fixed = T)
	
	maxl <- max(lengths(l))
	
	sapply(l, "length<-", maxl) %>%
		`[`(n,) %>%
		as.character
}

microbenchmark({
	dt[
		, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
		]
})

dt[
	, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
	]


# chaining

res1 <- dt[a == 'a'][sample(.N, 100)]

res2 <- dt[, .N, a][, N]

res3 <- dt[, coefficients(lm(e ~ d))[1], a][, .(letter = a, coef = V1)]

# piping

samplpe_b <- dt[a %in% head(letters), sample(b, 1)]

res4 <- 
	dt %>%
	.[a %in% head(letters)] %>%
	.[, 
	  {
	  	dt0 <- .SD[1:100]
	  	
	  	quants <- 
	  		dt0[, c] %>%
	  		quantile(seq(0.1, 1, 0.1), na.rm = T)
	  	
	  	.(q = quants)
	  }
	  , .(cond = b > samplpe_b)
	  ] %>%
	glm(
		cond ~ q -1
		, family = binomial(link = "logit")
		, data = .
	) %>%
	summary %>%
	.[[12]]


# function

rm(lm_preds)

lm_preds <- function(
	sd, by, n
)
{
	
	if(
		n < 100 | 
		!by[['a']] %in% head(letters, 4)
	   )
	{
		
		res <-
			list(
				low = NA
				, mean = NA
				, high = NA
				, coefs = NA
			)
		
	} else {

		lmm <- 
			lm(
				d ~ c + b
				, data = sd
			)
		
		preds <- 
			stats::predict.lm(
				lmm
				, sd
				, interval = "prediction"
				)
		
		res <-
			list(
				low = preds[, 2]
				, mean = preds[, 1]
				, high = preds[, 3]
				, coefs = coefficients(lmm)
			)
	}

	res
	
}

res5 <- 
	dt %>%
	.[e < 0] %>%
	.[.[, .I[b > 0]]] %>%
	.[, `:=` (
		low = as.numeric(lm_preds(.SD, .BY, .N)[[1]])
		, mean = as.numeric(lm_preds(.SD, .BY, .N)[[2]])
		, high = as.numeric(lm_preds(.SD, .BY, .N)[[3]])
		, coef_c = as.numeric(lm_preds(.SD, .BY, .N)[[4]][1])
		, coef_b = as.numeric(lm_preds(.SD, .BY, .N)[[4]][2])
		, coef_int = as.numeric(lm_preds(.SD, .BY, .N)[[4]][3])
	)
	, a
	] %>%
	.[!is.na(mean), -'e', with = F]


# plot

plo <- 
	res5 %>%
	ggplot +
	facet_wrap(~ a) +
	geom_ribbon(
		aes(
			x = c * coef_c + b * coef_b + coef_int
			, ymin = low
			, ymax = high
			, fill = a
		)
		, size = 0.1
		, alpha = 0.1
	) +
	geom_point(
		aes(
			x = c * coef_c + b * coef_b + coef_int
			, y = mean
			, color = a
		)
		, size = 1
	) +
	geom_point(
		aes(
			x = c * coef_c + b * coef_b + coef_int
			, y = d
		)
		, size = 1
		, color = 'black'
	) +
	theme_minimal()

print(plo)

Източник: habr.com