df = data.frame(a = sample(c('a', 'b', 'c', 'd'), 1000, replace = TRUE),
b = sample(c('x', 'y', 'zizuzi'), 1000, replace = TRUE),
c = rnorm(1000))
Загружаем необходимые библиотеки
library(dplyr) # преобразование данных library(ggplot2) # визуализация library(scales) # на диаграмме значения преобразовывает в %Далее собственно преобразование
library(rvest)
team = c('anaheim-ducks' , 'arizona-coyotes', 'boston-bruins', 'buffalo-sabres',
'calgary-flames', 'carolina-hurricanes', 'chicago-blackhawks', 'colorado-avalanche',
'columbus-blue-jackets', 'dallas-stars', 'detroit-red-wings', 'edmonton-oilers',
'florida-panthers', 'los-angeles-kings', 'minnesota-wild', 'montreal-canadiens',
'nashville-predators', 'new-jersey-devils', 'new-york-islanders', 'new-york-rangers',
'ottawa-senators', 'philadelphia-flyers', 'pittsburgh-penguins',
'san-jose-sharks', 'st-louis-blues', 'tampa-bay-lightning', 'toronto-maple-leafs',
'vancouver-canucks', 'vegas-golden-knights', 'washington-capitals',
'winnipeg-jets')
y = NULL;
for(i in team) {
k=html_table(html_nodes(read_html(paste0('https://www.spotrac.com/nhl/rankings/cash/',i,'/')), "table")[[1]], header = TRUE)
k$team = i
y = rbind(y, k)
}
y
obj = strsplit(as.character(y$Player), "\n\t\t" )
y$Player = sapply(obj, '[', 1)
y$amplua = trimws(lapply(obj, '[', 2))
y$cash = as.numeric(gsub("[^0-9]", "", y$cash))
Кстати, если вы знаете более простой способ скачивания данных - поделитесь!
Спасибо!df = data.frame(id = c(1,1,1,2,2,2),
begin=c('01.01.2018', '01.02.2018','01.03.2018', '01.04.2018','01.05.2018', '01.06.2018' ),
end=c('01.01.2018', '01.03.2018','01.04.2018', '01.04.2018','01.06.2018', '01.07.2018'),
category = c('Прием', 'Повышение', 'Увольнение','Прием', 'Повышение', 'Увольнение' ))
id begin end category
1 1 01.01.2018 01.01.2018 Прием
2 1 01.02.2018 01.03.2018 Повышение
3 1 01.03.2018 01.04.2018 Увольнение
4 2 01.04.2018 01.04.2018 Прием
5 2 01.05.2018 01.06.2018 Повышение
6 2 01.06.2018 01.07.2018 Увольнение
Toy example. Здесь id - табельный работника, и вы видите, что первые три строки - это все про одного работника. Begin - начало как-то периода. End - конец периода. category, соответственно, описание периода. Так, например, работник один был принят 1 января 2018, уволился 1 апреля 2018.
df$number = c(1:nrow(df))Создаем новую переменную, которая содержит номера строк. И
df$end1 = ifelse(df$id == df$id[df$number+1] & df$category == 'Прием', as.character(df$begin[df$number+1]), as.character(df$end)) df id begin end category number end1 1 1 01.01.2018 01.01.2018 Прием 1 01.02.2018 2 1 01.02.2018 01.03.2018 Повышение 2 01.03.2018 3 1 01.03.2018 01.04.2018 Увольнение 3 01.04.2018 4 2 01.04.2018 01.04.2018 Прием 4 01.05.2018 5 2 01.05.2018 01.06.2018 Повышение 5 01.06.2018 6 2 01.06.2018 01.07.2018 Увольнение 6 01.07.2018Можно было сразу результат записывать в переменную df$end, но я создал переменную end1, чтобы можно было сравнить end и end1.
hh
|
Superjob
|
Отраслевые
|
другое
|
1
|
1
|
0
|
0
|
1
|
0
|
0
|
0
|
1
|
0
|
1
|
1
|
1
|
1
|
0
|
0
|
1
|
1
|
0
|
0
|
df$hh = ifelse(grepl('hh',df$На.каких.работных.сайтах.размещали.резюме), 1, 0)
HH = grep("hh",df$На.каких.работных.сайтах.размещали.резюме,value=F)
HH = length(HH)/nrow(df)
Superjob = grep("Superjob",df$На.каких.работных.сайтах.размещали.резюме,value=F)
Superjob = length(Superjob)/nrow(df)
Первой строкой я считаю количество ячеек, содержащих hh, второй строкой я считаю % (делю на количество строк датасета). Если вам не нужны проценты, вы убираете знаменатель и все. Теперь у нас есть два объекта HH и Superjob, нам нужно получить диаграмму (на самом деле у вас будет больше объектов, по числу вариантов ответа, но я просто экономлю пространство, поскольку это просто повторение. И диаграмму я получаю так: объединяю объекты в один датасет js = data.frame(value = c(HH, Superjob, Rabota, job, Росработа, avito, otrasl))
col = c('HH', 'Superjob', 'Rabota', 'job', 'Росработа', 'avito', 'otrasl')
df = cbind(col, js)
И теперь визуализируем
df = df[order(-(df$value)),] # задаю в порядке убывания df library(ggplot2) ggplot(df,aes(x = reorder(df$col, df$value), y = df$value)) + geom_bar(stat="identity", fill="coral1", colour="darkgreen") + geom_text(aes(label=paste(round(df$value*100, 1),"%",sep="") ), hjust = 0.5, size=8) + theme(text = element_text(size=30), plot.title = element_text(hjust = .5)) + coord_flip()+ labs(title = "% респондентов", y = " ", x='сети')+ theme(axis.text = element_text(size = 19))Обратите внимание, я перед загрузкой диаграммы сначала расположил варианты ответа в порядке убывания.
q$hh = ifelse(grepl("hh",q$На.каких.работных.сайтах.размещали.резюме), '1', '0')
hh = q[q$hh == '1', ]
С помощью первой строки я повторяю, что делал выше - создаю новую переменную, а второй строкой фильтром отбираю только те строки полного датасета, где содержится HH. Далее я с помощью пакет dplyr отбираю зарплаты, создаю из них новый датасет hh, и в нем новую переменную js -она - эта переменная - состоит только из значений hh, т.е. мы выбрали зарплаты только тех респондентов, которые размещали свои резюме на HH.RU
library("plyr"); library("dplyr")
hh = select(hh, zp = Размер.зарплаты.на.момент.увольнения)
str(hh)
hh$js = 'hh'
Аналогично я проделываю это с другими сайтами: для Superjob создаю датасет Superjob и так далее
q$Superjob = ifelse(grepl("Superjob",q$На.каких.работных.сайтах.размещали.резюме), '1', '0')
Superjob = q[q$Superjob == '1', ]
Superjob = select(Superjob, zp = Размер.зарплаты.на.момент.увольнения)
str(Superjob)
Superjob$js = 'Superjob'
Теперь нам надо все это объединить в единый датасет. rbind() нам в помощь, поскольку названия переменных во всех датасетах у нас идентичные.
js = rbind(hh, Superjob, Rabota, otrasl ) js$zp = as.numeric(as.character(js$zp)) js$js = as.factor(js$js)Далее мы можем посмотреть значимость различий
kruskal.test(js$zp ~ js$js)Выяснив, что для наших данных таких различий нет. И построим боксплот
ggplot(js, aes(x=js, y=zp/1000, fill=js)) + geom_boxplot(color='blue') +
theme_grey(base_size = 30)+ theme(legend.position="none") +
scale_y_continuous(breaks=seq(0, 250, 25)) + coord_flip(ylim=c(0, 250)) +
ylab('зарплата')
Обратите внимание, что в формуле я использую zp/1000 - т.е. делю зарплату на 1000, чтобы деления не рябили нулями, были не в рублях, а тысячах рублей.job = grep("job",df$На.каких.работных.сайтах.размещали.резюме,value=F)
job = length(job)/nrow(df)
То мы соберем в эту переменную помимо job еще и весь Superjob, в котором также есть job. Я преодолеваю эту засаду таким образом.
job = q[!grepl("Super",q$На.каких.работных.сайтах.размещали.резюме) & grepl("job",q$На.каких.работных.сайтах.размещали.резюме), ]
job = nrow(job)/nrow(q)
Буквально эта формула означает: мы собираем все, что НЕ содержит Super и содержит job.
library(HRanalytics)
data("survey")
Сабсетинг я не буду показывать, если вам, как спецу интересен этот пост, значит вы уже умеете делать это. Я беру полные строки по полу (обратите внимание, там пол указан в разных регистрах), выбираю семейный статус женат и гражданский брак, и полные строки по году рождения. Всего 541 строка. Репрезентативно? Проверим наличие корреляции между возрастом и инстроверсиейcor.test(q1$Год.Рождения, q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.)
Pearson's product-moment correlation
data: q1$Год.Рождения and q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.
t = 4.0552, df = 539, p-value = 5.749e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.08905551 0.25269745
sample estimates:
cor
0.1720633
Действительно, слабая (< 0, 3) связь есть. Чем моложе женщина, тем экстравертированней.
summary(glm(q1$Семейное.положение.на.момент.трудоустройства ~ q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.+ q1$year, family= 'binomial')) q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ. -2.0880 0.5396 -3.870 0.000109 *** q1$year -2.8793 0.6573 -4.380 1.18e-05 ***Показываю уравнение лог регрессии и значимость предикторов. Сами переменные шкалированы в границах 0:1. Значимы обе переменные. Тем не менее, применим формулу sobel(), где зависимая переменная будет as.integer(q1$Семейное.положение.на.момент.трудоустройства)
library(multilevel) sobel(pred=q1$Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.,med=q1$year,out=q1$Семейное.положение.на.момент.трудоустройства)У нас q1$Семейное.положение.на.момент.трудоустройства) - это output или зависимая переменная
$`Mod1: Y~X`
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.967276 0.05170370 38.049047 8.300863e-155
pred -0.432824 0.09323694 -4.642194 4.334812e-06
$`Mod2: Y~X+M`
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.2387412 0.07973573 28.077014 1.725524e-107
pred -0.3620601 0.09306240 -3.890509 1.125716e-04
med -0.4763393 0.10778717 -4.419258 1.198319e-05
$`Mod3: M~X`
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5698981 0.02031517 28.052834 1.948845e-107
pred 0.1485578 0.03663422 4.055165 5.749226e-05
$Indirect.Effect
[1] -0.07076391
Первое, на что обращаем внимание
число уволенных работников за отчетный период / среднесписочная численность работников за отчетный период.И среднесписочная считается так:
Среднесписочная численность работников "за отчетный месяц исчисляется путем суммирования численности работников списочного состава за каждый календарный день отчетного месяца, т.е. с 1 по 30 или 31 число (для февраля - по 28 или 29 число), включая праздничные (нерабочие) и выходные дни, и деления полученной суммы на число календарных дней отчетного месяца.
Численность работников списочного состава за выходной или праздничный (нерабочий) день принимается равной списочной численности работников за предшествующий рабочий день. При наличии двух или более выходных или праздничных (нерабочих) дней подряд численность работников списочного состава за каждый из этих дней принимается равной численности работников списочного состава за рабочий день, предшествовавший выходным или праздничным (нерабочим) дням." Консультант плюс
df$begin = as.Date(df$Дата.поступления, "%d.%m.%Y") df$end = as.Date(df$дата.увольнения, "%d.%m.%Y") df$end[is.na(df$end)] = "2018-05-30 UTC"Последняя строка здесь - я заполнил пустые ячейки дат увольнения (пустые, потому что эти работники продолжают работать) датой выгрузки отчета. И формула
n = NULL;
m = NULL;
for(i in seq(as.Date('2018-01-01'), as.Date('2018-01-31'), by = 'day')) {
x = length(which(i == df$end))
m = rbind(m, x)
y = length(which(df$begin < i))-length(which(i > df$end))
n = rbind(n, y)
}
print(sum(m)/mean(n))
Здесь:seq(as.Date('2018-01-01'), as.Date('2018-01-31'), by = 'day')
это последовательность дней месяцаturnover = function(data, date_begin, date_end) {
n = NULL;
m = NULL;
for(i in seq(as.Date(date_begin), as.Date(date_end), by = 'day')) {
x = length(which(i == df$end))
m = rbind(m, x)
y = length(which(df$begin < i))-length(which(i > df$end))
n = rbind(n, y) }
print(sum(m)/mean(n))
}
И мы можем получить результат, где df - наши данные, первая дата - начала периода, вторая дата - конец периода.turnover(df, '2018-01-01', '2018-01-31')
[1] 0.01064467
turn_plot = function(data, date_begin, date_end) {
n = NULL;
m = NULL;
for(i in seq(as.Date(date_begin), as.Date(date_end), by = 'day')) {
x = length(which(i == df$end))
m = rbind(m, x)
y = length(which(df$begin < i))-length(which(i > df$end))
n = rbind(n, y) }
a =format(seq(as.Date(date_begin), as.Date(date_end), by = 'day'), '%m.%Y')
gh = data.frame(tapply(m/n,a, sum ))
gh$date = row.names(gh)
colnames(gh) = c('num', 'date')
gh$num = as.numeric(gh$num)
library(ggplot2)
library(scales)
ggplot(data=gh, aes(x=date, y=num)) + geom_bar(stat="identity", color="blue", fill="white") +
geom_text(aes(label=scales::percent(num, accuracy = 0.001)), vjust=-0.1, size=6.5) +
theme_grey(base_size = 30)
}
И я задаю даты за 4 месяца
turn_plot(df, '2018-01-01', '2018-04-30') # только имейте ввиду, что дата '2018-04-31' выдаст ошибкуИ мы получаем
library(mlr) library(dplyr) library(data.table) library(mltools) library(ROCR)давайте выберем необходимые переменные
data = dplyr::select(survey, feedback=Как.часто.Ваш.руководитель.давал..дает..Вам.обратную.связь,
iq=Ш1...ОБЩИЙ.БАЛЛ., intro = Ш6...ИНТРОВЕРСИЯ...ЭКСТРАВЕРСИЯ.,
impls = Ш8...ИМПУЛЬСИВНОСТЬ...САМОКОНТРОЛЬ.,
independ = Ш7...НЕЗАВИСИМОСТЬ...СОГЛАСИЕ.)
У нас много пропущенных значений, поэтому просто оставим полныеdata = data[complete.cases(data), ]Фишка мультиклассификации в пакете mlr в том, что целевая переменная принимается в one hot encoding и притом еще тип переменной - logical :)
data[, c('day', 'everyweek', 'never', 'moreweek')] = one_hot(as.data.table(data$feedback))
data[, c('day', 'everyweek', 'never', 'moreweek')] =lapply(data[, c('day', 'everyweek', 'never', 'moreweek')], as.logical)
Это выглядит так
head(data[, c('day', 'everyweek', 'never', 'moreweek')])
day everyweek never moreweek
1 FALSE FALSE TRUE FALSE
2 FALSE FALSE TRUE FALSE
15 FALSE TRUE FALSE FALSE
16 FALSE TRUE FALSE FALSE
24 FALSE FALSE TRUE FALSE
25 FALSE TRUE FALSE FALSE
Далее приводим данные в необходимый формат mlr
data = data[, -1]
labels = c('day', 'everyweek', 'never', 'moreweek')
data.task = makeMultilabelTask(id = "multi", data = data[, -1], target = labels)
Обратите внимание, я беру data[, -1], потому что первая переменная эта наша "сырая" целевая переменная,с ней наш прогноз будет очень точным:)
Хотя он и так будет точным, потому что я не буду разбивать на трейн / тест сплит, думаю, вы это сделаете сами.
Далее указываем учителя
mrf = makeLearner("multilabel.randomForestSRC", predict.type = "prob")
getParamSet("multilabel.randomForestSRC")
Обращаю внимание, ниже я даю ссылку на источник, там вы можете увидеть вариант, как можно решать задачи мультиклассификации в mlr с помощью учителя бинарного классификатора.
Задаем параметры (я даю мало, вы же сами выбираете то, что вам нра.
mrf_param = makeParamSet(
makeIntegerParam("ntree",lower = 50, upper = 500),
makeIntegerParam("mtry", lower = 1, upper = 6)
)
Далее обычные операции подготовки данных
rancontrol = makeTuneControlRandom(maxit = 7L) # мы ограничимся 5 итерациями, вы пробуйте больше от 50
set_cv = makeResampleDesc("CV", stratify = FALSE, iters = 3)
А здесь важно обратить вминание, что stratify = TRUE не работает, и методы кросс валидации мультиклассификации в mlr отличаются от просто классификации и регрессии (см. ниже по ссылке).Тренируем
rf_tune = tuneParams(learner = mrf, resampling = set_cv, task = data.task,
par.set = mrf_param, control = rancontrol)
mrf.tree = setHyperPars(mrf, par.vals = rf_tune$x)
mrforest = mlr::train(mrf.tree, data.task)
getLearnerModel(mrforest)
Далее я хочу получить метрики качества модели, но поскольку я не разбивал на трейн и тест сеты, то использую полный набор данных.
mrfmodel = predict(mrforest, data.task)
head(mrfmodel$data)
getMultilabelBinaryPerformances(mrfmodel, measures = list(acc, mmce, auc))
acc.test.mean mmce.test.mean auc.test.mean
day 0.9486964 0.05130362 0.9882470
everyweek 0.9276703 0.07232969 0.9821936
never 0.9520606 0.04793944 0.9902949
moreweek 0.9091674 0.09083263 0.9787145
Вам главный вопрос: почему мы получили такие хорошие метрики качества? Нарисую ROC AUC для вариант day - площадь под кривой того, что руководитель будет давать обратную связь каждый день
ROCRpred1 = prediction(mrfmodel$data$prob.day , mrfmodel$data$truth.day ) # as.integer(data$day) as.numeric(ROCR::performance(ROCRpred1, "auc")@y.values) ROCRperf1 = ROCR::performance(ROCRpred1, "tpr", "fpr") par(cex.axis=1, cex=1) plot(ROCRperf1 , colorize=TRUE, lwd = 6, print.cutoffs.at=seq(0,1,by=0.1), text.adj=c(-0.2,1.7), cex.axis=2) grid(nx=NULL, ny=NULL, lwd = 3) lines(x = c(0,1), y = c(0,1), lty=5)Картинка такая
__________________________________________________________
shapiro.test()
cor.test(x, y, method = 'pearson')R может нам выдать что-то типа
Pearson's product-moment correlation data: x and y t = 4.5451, df = 79, p-value = 1.95e-05 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.2630944 0.6127261 sample estimates: cor 0.455289p-value у нас 1.95e-05, это меньше 0, 05, что говорит о наличии связи, а сила связи 0, 455289, что говорит о средней силе связи (связь считается сильной, если показатель связи > 0, 5, слабой, если < 0, 1)
cor.test(x, y, method = ' spearman')
table()Или с помощью
prop.table()Если мы хотим посмотреть распределение в %. Отвечаем на вопрос о наличии связи с помощью критерия Хи-квадрат.
chisq.test()И на вопрос о силе связи мы отвечаем с помощью критерия Крамера
library(lsr) cramersV()
"devtools" "HRanalytics"Установить их можно так, как указано по ссылке Пакет R HR analytics для HR аналитиков. В пакете HRanalytics собраны часть результатов нашего исследования факторов эффективности и текучести персонала, и, если кто-то помнит, респонденты могли пройти тесты Лаборатории Гуманитарные Технологии:
0:28:38 0:27:40 0:28:51 0:21:30Т.е. у нас сразу встает задача перевести данный формат в удобоваримый. Давайте по порядку.
library(dplyr)э library(HRanalytics) data(survey) str(survey)В файле survey куча ненужных нам переменных, поэтому мы формируем отдельный файл в помощью команды select пакета dplyr
tests1 =select(survey, overall=Ш1...ОБЩИЙ.БАЛЛ., time=ПРОДОЛЖИТЕЛЬНОСТЬ ) tests1 = tests1[complete.cases(tests1$time),]Вторая команда нужна, чтобы убрать много пустых строк: у нас далеко не все респонденты проходили тесты. Итого у нас 1 208 строк, что достаточно для исследования. Теперь самое интересное: давайте переведем наш формат времени в минуты.
obj = strsplit(as.character(tests1$time), ":") # делим текст на объекты по запятой obj[1:5] # посмотрим, что представляет из себя объект tests1$часы = as.integer(lapply(obj, '[', 1)) # превращаем последовательно объект в часы, минуты и секунды tests1$минуты = as.integer(lapply(obj, '[', 2)) tests1$секунды = as.integer(lapply(obj, '[', 3)) tests1$времятеста = tests1$часы*60+tests1$минуты+tests1$секунды/60 # а теперь сведем в одну переменную tests1$времятеста = round(tests1$времятеста, 1) # ограничим одним знаком # ну и для марафета ограничимся одним знаком после запятойИ посмотрим, что из себя время прохождения теста представляет
summary(tests1$времятеста) Min. 1st Qu. Median Mean 3rd Qu. Max. 3.60 22.37 25.45 35.47 28.90 1092.00Диаграмма
hist(tests1$времятеста)
tests2 = subset(tests1, tests1$времятеста <= 60) hist(tests2$времятеста) summary(tests2$времятеста) Min. 1st Qu. Median Mean 3rd Qu. Max. 3.60 22.20 25.40 25.97 28.60 59.10Эти результаты кажутся более осмысленным, хотя стоило бы еще как-то обрезать слева: ну нереально пройти батарею тестов за 3.6 минуты, согласны?
shapiro.test(tests2$времятеста) Shapiro-Wilk normality test data: tests2$времятеста W = 0.92625, p-value < 2.2e-16 > shapiro.test(tests2$overall) Shapiro-Wilk normality test data: tests2$overall W = 0.99435, p-value = 0.0002018P-value менее 0, 05 говорит нам о том, что распределение обоих переменных распределено ненормально, значит для проверки гипотезы о наличии связи между временем прохождения теста и интеллектом мы должны применять не коэффициент корреляции Пирсона, а непараметрический - Спирмена.
hist(tests2$overall)А на диаграмме распределение ответов по шкале общих способностей теста КТО кажется вполне нормальным
cor.test(tests2$времятеста, tests2$overall, method = c("spearman") )
Spearman's rank correlation rho
data: tests2$времятеста and tests2$overall
S = 317090000, p-value = 8.552e-09
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.166812
p-value = 8.552e-09 и означает, что связь есть, -0.166812 связь не очень тесная, но она есть и она отрицательная. Последнее выглядит логичным: чем больше время прохождения теста, тем меньше интеллект.
> tests3 = subset(tests2, tests2$времятеста >= 15)
> cor.test(tests3$времятеста, tests3$overall, method = c("spearman") )
Spearman's rank correlation rho
data: tests3$времятеста and tests3$overall
S = 310710000, p-value = 2.575e-13
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.2130846
P-value и коэффициент корреляции показывают более высокую связь, чем в предыдущем случае. Таким образом, мы можем говорить о наличии связи между временем прохождения теста и интеллектом (шкалой общих способностей теста КТО).