R中Web抓取和解析数据 | 研究H-1b数据（2）

本文概述

• 探索性数据分析：基本工资
• 分类EDA：案例状态
• EDA时间：归档时间
• 奖励：使用多个功能
• 总结

srcmini的博客”数据可以帮助你的H-1B签证申请”向你介绍了多年来对H-1B数据的分析结果。现在, 是时候开始自己了！ Ted Kwartler将通过一系列R教程指导你完成此过程。

在我们的上一篇文章中摘录, 我想向你展示如何在探索性数据分析(EDA)工作中添加更多视觉效果, 并介绍一些新功能。上一次你学习了如何通过网络抓取180万H1B签证申请, 更正变量类型, 检查一些行并创建基本可视化效果, 例如barplot和scatterplot。

在这篇文章中, 我假设你能够使用fread加载数据, 将变量更正为日期, 将薪水更改为数字, 并创建新的数据功能(如年和月)。

首先, 添加一些可视化库。

library(ggplot2)
library(ggthemes)

探索性数据分析：基本工资

在执行EDA工作时, 摘要功能非常有用。在这种情况下, 你将summary()函数应用于数字矢量$ base_salary。除了数值之外, 你还可以将summary()函数应用于字符串和因子以获得适当的描述性统计信息。你甚至可以将其应用于整个数据帧以获取按列统计信息。这里summary()将打印H1B工资的最小值, 第一和第三四分位数, 中位数, 均值和最大值。

summary(h1b.data$base_salary)

在摘要打印输出中, 应注意最大值为7, 278, 872, 788。这绝对是一个错误！查找具有最大值的行的一种方法是使用which.max()。 which.max()函数确定给定向量中具有第一个最大值的行$ base_salary。

在下面的代码中, 使用此行号为数据框建立索引, 以便打印整个记录。

h1b.data[which.max(h1b.data$base_salary)]

那是一位高薪的IBM软件工程师。毫不奇怪, H1B案件状态为”撤回”。删除这样的异常值的一种快速方法是简单地将一个负号与which.max语句一起传递。仅出于教育目的, 你可以在除去最高薪水之前和之后使用nrow验证数据框中的行数。

nrow(h1b.data)
h1b.data<-h1b.data[-which.max(h1b.data$base_salary)]
nrow(h1b.data)

删除单行很繁琐, 而且通常不是删除异常值的最佳方法。相反, 你可以绘制薪水以直观地了解数据分布。在此代码中, 最里面的功能是排序。 Sort()接受薪水向量, 默认情况下, sort()将按升序排列值。下一个最里面的函数是tail()函数。在此示例中, tail()将采用排序后的向量并将其子集化为最后50个值。由于顺序是递增的, tail()将对50个最大的薪水进行排序。最后, 基于R的绘图功能将自动创建50个薪水的点状图。

plot(tail(sort(h1b.data$base_salary), 50))

即使取消了最高薪水, 最近五位的薪水也从1000万美元到7.5亿美元不等。

在下一部分中, 你将学习如何使用ifelse语句创建任意截止值。这是一次删除离群值的有效方法, 而不是一次。

ifelse语句的工作方式如下：如果某些逻辑声明为TRUE, 则执行一个动作, 否则(else)执行另一个动作。

在此代码中, 你将使用ifelse重写$ basesalary向量。然后, 在向量上调用summary()以显示ifesle语句如何更改输出。

具体来说, 代码遵循以下逻辑：如果$ base_salary大于$ 250, 000, 则将其更改为NA, 否则使用原始的$ base_salary。

h1b.data$base_salary<-ifelse(h1b.data$base_salary>250000, NA, h1b.data$base_salary)
summary(h1b.data$base_salary)

请注意, 摘要输出现在将7, 137 NA添加到其他描述性统计信息中。此外, 最大值现在为250000。

对于连续数据, 你可以使用核密度图探索分布。核密度图可估计样本的总体分布。在这里, 你首先使用ggplot创建基本的可视层。

首先传入h1b.data, 然后使用列名base_salary设置X美感。接下来, 添加geom_density图层以创建内核密度图。在参数内, 以下代码指定填充颜色, “深色”和” NA”作为轮廓色。在最后一层中, 你使用ggtheme的gdocs主题并删除Y轴文本。

ggplot(h1b.data, aes(base_salary)) +
  geom_density(fill='darkred', color='NA') + 
  theme_gdocs() + 
  theme(axis.text.y=element_blank())

你会看到H1B薪水大部分在$ 50K以上, 但在$ 100K以下。

分类EDA：案例状态

现在, 你已经对探索数字值有了基本的了解, 可以过渡到分类值。

分类值表示唯一的类属性。例如, 从每天的7个”级别”开始, 一周中的某天可以是星期一, 星期二, 星期三等。在R中, 分类值可能会错误地变成字符串, 这意味着它们即使在谈论”星期二”的多天时也都是唯一的。相反, 两个类别的星期四共享相同的属性, 因为它们都共享”星期四”级别。

在本节中, 你将查看$ case_status向量。 H1B签证数据在此列中重复, 使人们认为它们是申请状态的因素, 而不是唯一的字符串。要检查这些信息, 请使用as.factor()确保向量是一个因子。这会将类数据更改为分类值。接下来, 使用table()来汇总信息：

h1b.data$case_status<-as.factor(h1b.data$case_status)
table(h1b.data$case_status)

再一次, 值不是平均分配的。 180万行中的1.7m已”认证”。另外95k是” WITHDRAWN”或” DENIED”, 而其余因子水平总计小于20。

下面的代码演示了如何删除异常值。首先, 你创建一个小的向量滴, 其中包含要删除的因子水平。每个因子水平由竖线” |”分隔操作员。管道是Enter键上方的一条直线。在逻辑表达式中, 符号表示”或”。

因此, 放置对象说的是”质量无效或待定…或未分配或拒绝”。

drops<-c('INVALIDATED|PENDING QUALITY AND COMPLIANCE REVIEW - UNASSIGNED|REJECTED')

接下来, 如果行之一包含drop的级别, 则使用grepl()声明TRUE。 grepl()命令是一个古老的UNIX命令, 用于搜索模式。如果找到该模式, 则输出为TRUE, 否则为FALSE。

在这种情况下, 你将从drop中传递文本作为搜索模式。在下一个参数中, 将grepl()应用于$ case_status向量。

请注意, 以grepl表达式开头的感叹号：此符号表示否定。

放在一起, 当grepl返回TRUE时, 求反。

h1b.data<-h1b.data[!grepl(drops, h1b.data$case_status), ]

即使删除R, R也会保留未使用的电平。在这种情况下, h1b.data $ case_status上的调用表仍将显示为” REJECTED”(带零)。删除未使用级别的一种方法是在向量​​上调用factor。

你可以再次使用table()重新检查提示：

table(h1b.data$case_status)
h1b.data$case_status<-factor(h1b.data$case_status)
table(h1b.data$case_status)

可视化新因子分布的一种方法是使用数据可视化库ggplot。

在第一层中, 你传入数据框并使用aes设置美观。这将创建视觉的基本status.plot。下一层使用geom_bar添加几何条形, 然后使用coord_flip翻转x, y轴, 最后使用theme_gdocs()添加主题信息。 GDocs主题是预定义调色板的快捷方式。在最后一层, 你可以通过删除Y轴标题和图例来调整GDocs主题。

status.plot <- ggplot(h1b.data, aes(factor(case_status), fill=factor(case_status)))
status.plot + 
   geom_bar()+ 
   coord_flip()+ 
   theme_gdocs()+
   theme(axis.title.y=element_blank(), legend.position="none")

该图比较了H1B病例状态因素。

薪资和现况：按要素EDA的数字

探索数据的另一种方法是查看变量之间的相互作用。如果你正在进行预测建模, 并且可以通过因变量检查输入, 则这特别有用。

在本节中, 你将按H1B案例状态浏览工资分配。

在此代码中, 你正在构建箱线图。箱形图是查看人口分布的紧凑方法。

下图在左侧垂直放置一个正态分布, 在右侧放置一个箱形图。在方框图中, 人口中位数表示为方框内的黑线。该线将人口分为第二和第三四分位数。框外延伸的”胡须”代表第一和第四四分位数。最后, 圆点是为了表示分布的真实离群值。

这是与箱形图比较的正态分布, 用于比较。

你可以使用ggplot构造箱形图并传递H1B数据。 geom_boxplot()函数是创建箱形图分布的图层。在此函数中, 你可以使用诸如” CERTIFIED”之类的值将X轴的美感设置为case_status。在es中, Y轴是base_salary, 最后一个参数只是对不同因素进行颜色编码。

下一部分代码说明了第三种方法, 用于消除异常值。以前, 你使用which.max()和ifelse来删除异常值。

在ggplot可视图中, 你还可以通过限制轴来省略记录。

在这里, 你正在使用ylim并传入0, 100000。这意味着Y轴只会扩展到100, 000美元, 而ggplot会忽略其他任何轴。最后, 和以前一样, 你添加theme_gdocs色彩并通过删除X轴文本来调整主题。

ggplot(h1b.data) +
  geom_boxplot(aes(factor(case_status), base_salary, fill=as.factor(case_status))) +
  ylim(0, 100000) + 
  theme_gdocs() + 
  scale_fill_gdocs() + 
  theme(axis.text.x=element_blank())

方框图显示拒绝H1B申请的薪水较低。看来, 年薪在75, 000美元以上的人更有可能获得认证。

EDA时间：归档时间

探索这些数据的另一种方法是按时间。假设你想知道何时申请H1B签证, 你将需要专注于CERTIFIED申请。使用子集, 传递整个数据帧, 然后传递逻辑条件。

case_status等于” CERTIFIED”的任何行都将保留在case.month中。

case.month<-subset(h1b.data, h1b.data$case_status=='CERTIFIED')

现在, 你将需要沿着两个维度(月份和年份)描述认证的应用程序。

你将使用tapply()传递case.month $ case $ status, 然后将其他两个维放入列表中。最后一个参数是要应用于数组的函数。在这种情况下, length将逐年返回经过认证的H1B申请的数量。在第二行中, 调用case.month以显示值。

case.month<-tapply(case.month$case_status, list(case.month$submit_month, case.month$submit_yr), length)
case.month

如果使用data.table的melt()函数对数据进行整形, 则绘制该数据将更加容易。此代码会将月份的年份和值分为三列。你可以使用head()检查新形状。

case.month<-melt(case.month)
head(case.month)

(可选)你可以使用match()中的常数month.abb将月份更改为数值。在这里, 你将Var1从月份的缩写重写为数字。如果执行此操作, 则会丢失绘图中的月份标签, 但这些行按时间顺序排列, 而不是按月份名称的字母顺序排列：

case.month$Var1<-match(case.month$Var1, month.abb)

接下来, 为时间探索创建时间线。

首先, 传递case.month对象。接下来, 添加代表月份的X轴, 然后按年份分组和颜色。然后像以前一样添加ggtheme美学。调整主题以删除x标题, 最后最后一层指定X轴应具有表示月份的1到12的整数。

ggplot(case.month) + 
  geom_line(aes(Var1, value, group=as.factor(Var2), color=as.factor(Var2))) + 
  theme_gdocs() + 
  scale_color_gdocs() + 
  theme(axis.title.x=element_blank())+
  scale_x_continuous(breaks=seq(1, 12))

3月(第3个月)是获得H1B应用认证的高峰。

既然你知道自己的申请很有可能在三月份获得认证, 那么你需要知道处理过程需要多长时间。

这意味着你需要计算Submit_date和start_date之间的差异。你可以使用difftime()计算日期向量之间的差异。首先, 你传递的是commit_date, 然后是start_date。最后一个参数声明时间单位。

submit.start<- difftime(h1b.data$submit_date , h1b.data$start_date , units = c("days"))

掌握了每个H1B应用程序的时差后, 就可以打印均值和中值。根据这些数据, 你应该在三月前75-90天左右提交申请, 以提高机会！

mean(submit.start)
median(submit.start)

奖励：使用多个功能

我将展示使用此数据的最后一种方法是检查顶级雇主, 职位和工作地点, 以便你可以集中精力进行工作。

要快速收集多个功能的频率值, 可以编写一个函数而不是冗余代码。

对于此分析, 你将创建列名称的向量。

h1b.features <- c('employer', 'job_title', 'location')

接下来, 编写一个名为freq.get的函数。将向量名称作为x传递时, freq.get使用table()为该列构造一个提示。接下来, 对向量进行排序, 应用名称并构建数据帧。最后, 数据帧从上到下重新排序。

freq.get<-function(x){
  f.table<-table(h1b.data[[x]])
  f.table<-sort(f.table, decreasing=TRUE)
  f.names<-names(f.table)
  f.df<-data.frame(names=f.names, freq=f.table)
  f.df<-f.df[order(f.df$freq, decreasing = T), ]
  return(f.df)
}

而不是将此功能应用于每列, 而是使用pbapply包中的pblapply()。传入h1b.feature, 然后传入freq.get。输出是3个频率数据帧的列表。你在第二行中使用名称添加列表名称。

all.freq <- pblapply(h1b.features, freq.get)
names(all.freq)<-h1b.features

现在, 你可以轻松查看每个列中最频繁的值。请记住, 你可以在每一列上单独完成此操作, 但是在广泛的数据集中编写函数可以节省时间, 并且不易出错。

head(all.freq$employer, 12)
head(all.freq$location, 12)
head(all.freq$job_title, 12)

总结

希望你发现这些EDA帖子有趣。在此过程中, 你学习了如何通过Web抓取数据, 更正变量类以及最后获得基本功能以获得描述性统计信息。作为视觉学习者, 我在进行EDA时总是寻找有趣而有见地的可视化效果, 因此希望这些帖子能添加到你的可视化技能中。

如果你愿意的话, 可以使用此数据创建一个很酷的H1B可视化并将其发送给我@tkwartler！

请继续关注”使用R探索H-1B数据”教程系列的另一部分！同时, 请考虑查看我们的data.table备忘单, R数据框架教程, R中的数据分析, data.table Way课程或R的dplyr课程中的Data Manipulation。