地震分析(2/4):分类变量探索性分析

5月1日,二千零一十九
By

(本文首次发表于 R编程-数据科学+,并对 188bet appR博主

    类别

    1. 基础统计学

    标签

    1. 数据可视化
    2. 探索性分析
    3. R程序设计

    这是第二部分我们的邮政系列金宝搏网址关于在30天的特定时间窗内报告地震和类似事件的公开数据集的探索性分析。在以下内容中,我们将分析数据集的分类变量。分类变量可以是有限的,通常固定一些可能的值。因子变量是分类变量,可以是数值变量,也可以是字符串变量。r将分类变量存储到一个因子中。他们的分析可能需要不同于定量变量的统计工具。

    Packages

    I am going to take advantage of the following packages.

    suppresspackagestartupmessages(library(ggplot2))suppresspackagestartupmessages(library(dplyr))suppresspackagestartupmessages(library(hmisc))suppresspackagestartupmessages(library(lubridate))suppresspackagestartupmessages(library(vcd))suppresspackagestartupmessages(library(vclutra))suppresspackagestartupmessages(lib

    Packages versions are herein listed.

    包<-c(“ggplot2”,“DPLYR”“HMISC”“卢布里特”“VCD”“VCDEXTRA”,“gmodels”)版本<-lapply(软件包,packageversion)版本c<-do.call(c,版本)data.frame(packages=packages,version=as.character(version_c))。
    ##包装版本1 ggplot2 3.1.0 2 dplyr 0.8.0.1 3 hmisc 4.2.0 4 lubridate 1.7.4 5 vcd 1.4.4 6 vcdextra 0.7.1 7 gmodels 2.18.1

    在Windows-10上运行以下R语言版本。

    R.版本
    ####平台x86_64-w64-mingw32 arch x86_64 os mingw32 system x86_64,MingW32状态uage r version.string r version 3.5.2(2018-12-20)昵称蛋壳冰屋

    获取数据

    As shown in the first post,我们从下载地震数据集开始分析地震.usgs.gov站点,特别是过去30天的数据集风格。请注意,这样的地震数据集是逐日更新的,覆盖了数据收集的最后30天。此外,它不是最新的可用数据集,几周前我收集的。如果这样的数据集还没有出现在我们的工作区中,我们下载并保存它以加载到地震局部变量

    if(“all_week.csv%”,in%dir(“.”))==false)url<-“https://severse.usgs.gov/severs/feed/v1.0/summary/all_month.csv”下载.file(url=url,destfile = "all_week.csv")}quakes <- read.csv("all_month.csv",header=TRUE,SEP=,,stringsasFactors=false)地震$time<-ymd_hms(地震$time)地震$updated<-ymd_hms(地震$updated)地震$magtype<-as.factor(地震$magtype)地震$net<-as.factor(地震$net)地震$type<-as.factor(地震$type)地震$status<-as.factor(地震$status)地震$locationsource<-as.factor(地震$locationsource)地震$magsource<-as.factor(地震$locationsource)S$源)

    探索性分析——分类变量

    如果分类变量的类是因素.

    (因素变量<-名称class)==“因子”))
    ##[1]“magtype”“net”“type”“status”“[5]“locationSource”“magsource”
    长度(因子变量)
    α〔1〕6

    hmisc包中的describe()函数也可以用于分类变量。

    describe(quakes[,factor_vars])
    地震,因子变量6变量8407观测n缺少明显8407 0 10 35;值mb lg md mh ml mw mwr mWW \35; \ \35\\\35\\\35;\\\\\\\\\35;\\\\\\\35\\\\\\\\35\\\\\\\\\\;AK(2469,0.294)CI(1344)0.160)高压(253)0.030)ismpkansas(8,0.001)LD(4)0)MB(157)0.019)NC(1435)0.171)nm(28)0.003),nn(604,0.072)PR(427)0.051)硒(15)0.002)美国(897)0.107)乌苏(588)0.070)UW(178)0.021)化学爆炸(2,0)地震(8232,0.979)爆炸(58,0.007)冰震(16)0.002)其他事件(3)0)quarry blast (95,## 0.011),岩爆(1)0.000)频率1691 6716比例0.201 0.799--------------位置源8407 0 15值AK CI高压ISMP LD MB NC NM NN OK频率2470 1344 253 8 4 157 1435 28 604 6比例0.294 0.160 0.030 0.0010.000 0.000 0.019 0.171 0.003 0.003 0.072 0.001 \35\\35\\\\\35;\\35;\.000 0.000 0 0.0 0 0.0 0.0 0 0.0 0 0.171 0.003 0.003 0.0.0 0 0.0 0.0.0.0 0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.8407 015值AK-CI HV ISMP LD MB NC NM NN OK频率2480 1344 253 8 4 157 1435 28 604 5比例0.295 0.160 0.030 0.001 0.000 0.019 0.171 0.003 0.072 0.001 Vvalue pr se us uu uw频率427 15 881 588 178比例0.051 0.002 0.105 0.070 0.021

    We notice from the magType description that two records have a null string magType.然后我们将它们替换为NA值。

    地震$magtype[地震$magtype=“”]<-na

    To understand relationship or dependencies among categorical variables,我们利用各种类型的表和图形方法。此外,如果两个主要变量之间的关系对于所考虑的分层变量的所有级别都相同或不同,则可以包含分层变量,以突出这一点。

    应急表被称为单向的仅涉及一个分类变量时的味道。他们说双向的当涉及两个分类变量时,诸如此类(N-way

    例如,这是magtype变量的单向应急表。

    (tbl<-表(地震$magtype))
    ####mb mb_lg md mh ml mun mw mwr mww__0 604 47 2423 14 5203 2 4 19 89

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,AES(X=MaGyType)fill=magtype))+geom_bar(stat='count')+theme(axis.text.x=element_text(angle=45,hjust = 1)) + guides(fill=FALSE)

    基于事件注册网络的单向事件应急表。

    表(地震净值)
    ####AK CI高压ismpkansas ld mb 2469 1344 253 8 4 157 NC NM nn pr se us 1435 28 604 427 15 897 uu uw 588 178

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,aes(x = net,fill=net))+geom_bar(stat='count')+辅助线(fill=false)

    基于事件类型的单向应急表。

    表(地震类型)
    ## ## chemical explosion         earthquake          explosion ##                  2               8232                 58 ##          ice quake        other event       quarry blast ##                 16                  3                 95 ##         rock burst ##                  1

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,AES(x=型)fill=type))+geom_bar(stat='count')+theme(axis.text.x=element_text(angle=45,hjust = 1)) + guides(fill = FALSE)

    基于事件状态的单向应急表。

    表(地震$状态)
    ####自动审查1691 6716

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,AES(x=状态,fill=status))+geom_bar(stat='count')+辅助线(fill=false)

    基于位置源的单向事件应急表。

    表(地震$locationSource)
    ####AK CI高压ISMP LD MB NC NM NN OK PR SE US UU UW 2470 1344 253 8 4 157 1435 28 604 6 427 15 890 588 178

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,aes(x=位置源,fill=locationSource))+geom_bar(stat='count')+辅助线(fill=false)

    基于最初编写此事件报告大小的网络的事件单向应急表。

    表(地震$magsource)
    ####AK CI高压ISMP LD MB NC NM NN OK PR SE US UU UW 2480 1344 253 8 4 157 1435 28 604 5 427 15 881 588 178

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,aes(x=磁源,fill = magSource)) + geom_bar(stat='count') + guides(fill=FALSE)

    基于网络作为数据贡献者和MAGTYPE的双向应急表,即用于计算事件首选大小的方法或算法。

    表(地震$net,地震震级
    ####mb-mb_-lg-md-mh-mw-mw-mw ak 0 1 0 0 0 2467 0 0 0 0 0 1 ci 0 0 0 0 0 3 1339 2 0 0 0 0 hv 0 0 0 0 145 0 108 0 0 0 0 0 0 0 0 0 ismpkansas 0 0 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0##0 0 0 0 0 4 0 0 0 0 0 mb 0 0 0 0 6 0 151 0 0 0 0 0 0 nc 0 0 0 1423 0 7 0 3 0 0 nm 0 0 0 28 0 0 0 0 0 0 0 0 0 0 nn 0 0 0 0 0 0 0 604 0 0 0 0 0 0 0 0 pr 00 0 425 0 2 0 0 0 0 0 SE 0 0 0 0 0 15 0 0 0 0 0 0 0 0 0 0 US 0 603 47 0 0 140 0 0 0 19 88 UU 0 0 0 365 0 222 0 1 0 0 UW 0 0 0 0 16 11 151 0 0 0 0 0 0 0 0 0

    图中显示了与条形图相同的图形表示。

    ggplot(数据=地震,AES(X=Net),fill=magtype))+geom_bar(stat='count')

    从这张双向应急表开始:

    (tbl<-表(地震$net,地震震级)
    ####自动审查AK 1105 1364 CI 29 1315 HV 90 163 ISMPKANSAS 0 8 LD 0 4 MB 0 157 NC 460 975 NM 0 28 NN 6598 PR 0 427 SE 0 15 US 0 897 UU 0 588 UW 1 177

    它对应的行比例表如图所示。每一行显示给定网络的自动/回顾地震的比例。每行总计一行。

    支柱表(TBL)1)
    ####自动审查35 \35; \35\35\35\\\35\35\\\35\\\35\\\ \\35; \\ \\35\\\\\35\\\\\\\\NC 0.320557491 0.679442509牛米0.0000000001.000000000 nn 0.009933775 0.990066225 pr 0.000000000 1.000000000 se 0.000000000 1.000000000 us 0.000000000 1.000000000 uU 0.000000000 1.000000000 uW 0.005617978 0.994382022

    列比例表。每一行显示了给定特定状态下地震网络的比例(自动/审查)。每列总计一个。

    支柱表(TBL)2)
    ##             ##                 automatic     reviewed##   ak         0.6534594914 0.2030970816##   ci         0.0171496156 0.1958010721##   hv         0.0532229450 0.0242703990##   ismpkansas 0.0000000000 0.0011911852##   ld         0.0000000000 0.0005955926##   mb         0.0000000000 0.0233770101##   nc         0.2720283856 0.1451756998##   nM 0.0000000000000000.000.0041691483 \35\35\35;\\35;\\\35\\\35; \\\35; \\\\\\\35\\\\\\\\\\\\\\\\\\\

    让我们将条形图作为表示基于网络的事件计数的图,同时,用不同的填充颜色突出显示其状态。

    ggplot(数据=地震,aes(x = net,fill=status))+geom_bar(stat='count')

    我们希望提供更好的图形表示,在不同的地位比例可以更好地感知。然后我们构建一个数据框架,以频率形式收集我们的网络+地震状态信息。从它开始,我们将显示得到的尖晶石图。

    tbl_df<-as.data.frame(tbl)colnames(tbl_df)<-c(“net”,“地位”“freq”)tbl_df$net<-系数(tbl_df$net)tbl_df
    ##网络状态频率\\35\353535353535\35\\35\\\35\\35\35\\35\\\35\35\35\35\35\\\35\35\\35 10 pr自动0 11自动0 35353535353535353535353535353535\\\35\\\\35\35\\35\35\\\\35\35\35\\35\\35\ 21 NC审查975 22纳米审查28 23 NN审查598 24 PR审查427 25 SE审查15 26 US审查897 27 UU审查588 28 UW审查177

    在范围[0,1]内,spineplot使网络状态的计数差异更加明显。

    (xtabs-res<-xtabs(freq~net+状态,data = tbl_df))
    ##状态网络自动审查AK 1105 1364 CI 29 1315 HV 90 163 ISMPKANSAS 0 8 LD 0 4 MB 0 157 NC 460 975 NM 0 28 NN6 598 PR 0 427 SE 0 15 US 0 897 UU 0 588 UW 1 177
    尖晶石图

    函数的作用是:使用公式样式的输入创建数据的交叉表。此外,applying the summary() to the xtabs() result,我们得到了所有因素独立性的卡方检验,同时指示箱数和尺寸。

    键入“净”计数<-地震%>%分组依据(状态,net)%>%dplyr::summarse(freq=n())xtabs-res<-xtabs(freq~net+状态,data=type_net_count)摘要(xtabs_res)
    ##调用:xtabs(公式=freq~net+状态,data = type_net_count)## Number of cases in table: 8407 ## Number of factors: 2 ## Test for independence of all factors:##  Chisq = 2082.2,df = 13,p值=0卡方近似值可能不正确

    用稍微不同的方法得到的相同的尖晶石图。

    尖晶石图

    页边距表是总结分类数据的另一种方法。

    (mt<-利润表(tbl,2∶1)
    ####AK CI HV ISMPKANSAS LD MB NC NM NN PR SE自动1105 29 90 0 0 0 0 460 0 6 0 0审查1364 1315 163 8 4 157 975 28 598 427 15 US UU UW自动0 1审查897 588 177

    可以通过addMargins()函数添加每行和每列的总和。

    addmargins(mt)
    ####AK CI高压Ismpkansas LD MB NC NM NN PR SE自动1105 29 90 0 0 0 460 0 6 0审查1364 1315 163 8 4 157 975 28 598 427 15和2469 1344 253 8 4 157 1435 28 604 427 15美国UU UW SUM自动0 0 1 1691审查897 588 177 6716 SUM 897 588 178 8407

    The Cross Tabulation result is shown.

    mt<-利润表(表(地震$status,地震类型:2:1)十字工作台(mt,prop.chisq = FALSE,支柱=真,PR= r=真,format=“spss”)格式
    ## ##    Cell Contents## |-------------------------|## |                   Count |## |             Row Percent |## |          Column Percent |## |           Total Percent |## |-------------------------|## ## Total Observations in Table:  8407 ## ##                    |  ##                    | automatic  |  reviewed  | Row Total | ## -------------————————————————————————————————————————\35; 0 2 \\\\.000%.000% \\\      \\\\  35;--|##地震6544 8232 \\\\124; \\\ \\ \\  \124; \124; 8232 \\\\\\\\\\\\\\\\\\\——爆炸3 55 58|##| 5.172%;94.828%;0.690%\35\\\\\\\\\\\\\\\\\\\\ \\\\\\\\\\\\\\\\\\\\\\\%0.000%100.000%|0.190%\35\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\.000%0.036%0.000%|0.045%\\\\\35; \\\\  0.000%1.415%||##|0.000% 0.000%0.012%----—————————————————————————————————————————————列合计1691 6716 8407 20.114%79.886%20.114%79.886%--|####

    借助VCD包中的结构,我们展示了一个三向应急表。

    可构造(类型~status+magtype,数据=地震
    ##                   type chemical explosion earthquake explosion ice quake other event quarry blast rock burst## status    magType                                                                                           ## automatic                               0          0         0         0           0            0          0##           mb                            0          1         0         0           0            0          0##           mb_lg                         0          0         0         0           0            0          0##           md                            0        515         0         0           0            0          0##           mh                            0          0         0         0           0            0          0##           ml                            0       1172         3         0           0            0          0##           mun                           0          0         0         0           0            0          0##           mw                            0          0         0         0           0            0          0##           mwr                           0          0         0         0           0            0          0##           mww                           0          0         0         0           0            0          0## reviewed                                0          0         0         0           0            0          0##           mb                            0        603         0         0           0            0          0##           mb_lg                         0         47         0         0           0            0          0##           md                            0       1893         2         0           0           13          0##           mh                            0         14         0         0           0            0          0##           ml                            0       3873        53        16           3           82          1##           mun                           2          0         0         0           0            0          0##           mw                            0          4         0         0           0            0          0##           mwr                           0         19         0         0           0            0          0##           mww                           0         89         0         0           0            0          0

    地理区域对比深度类别地震分析

    Based on ref.[5]定义我们创建新的因素变量。

    地震<-地震%>%滤波器(类型=“地震”)地震<-地震[完整案例(地震[,C(“纬度”)“经度”,“深度”)],]发射球<-as.factor(ifelse(地震$Latitude>=0,“北方”,“南”))电磁波系数(ifelse(地震$经度>=0,"east",“西部”))地震$地区<-ifelse(地震$纬度>0&地震$经度>0,"NE",ifelse(地震$纬度>0&地震$经度<0,“西北”,ifelse(earthquakes$latitude < 0 & earthquakes$longitude > 0,“SE”“sw”))地震$深度_型<-ifelse(地震$深度<=70,“浅”ifelse(地震$depth<=300,“中间”“深”))地震$region<-factor(地震$region)地震$depth_type<-factor(地震$depth_type)

    将显示结果表。

    (tbl<-表(地震$region,地震$深度_型)
    ####深中浅东北9 65 197西北0 526 7109东南10 40 109西南31 45 91
    支柱表(TBL)1)
    ####深中浅东北0.03321033 0.23985240 0.72693727西北0.00000000 0.06889325 0.93110675东南0.06289308 0.25157233 0.68553459西南0.18562874 0.26946108 0.54491018
    支柱表(TBL)2)
    ####深中浅东北0.18000000 0.09615385 0.02624567西北0.00000000 0.77810651 0.94710898东南0.20000000 0.05917160 0.01452172西南0.620000000 0.06656805 0.01212363

    We wonder if the region plays a role for the depth type classification.我们运行下面的卡方测试来解决这个问题。

    Chisq.试验(地震$region,地震$深度类型,Simulate.p.Value=真)
    ## ##  Pearson's Chi-squared test with simulated p-value (based on 2000##  replicates)## ## data:  earthquakes$region and earthquakes$depth_type## X-squared = 1322.4,DF=Na,p值=0.0004998

    报告的P值突出显示,具有统计学意义。

    交叉表如下。

    mt<-利润表(tbl,2:1)十字工作台(mt,prop.chisq = FALSE,支柱=真,PR= r=真,format=“spss”)格式
    ## ##    Cell Contents## |-------------------------|## |                   Count |## |             Row Percent |## |          Column Percent |## |           Total Percent |## |-------------------------|## ## Total Observations in Table:  8232 ## ##              |  ##              |       NE  |       NW  |       SE  |       SW  | Row Total | ## ---————————————————————————————————————————————————————————————————\35\35深9 0 10 31 \\\\\\\\\\————————————————————————————————————————————————————————————————\\\\\\\\\\\\%18.563%0.109%0.000%| 0.121% |    0.377% |           | ## -------------|-----------|-----------|-----------|-----------|-----------|## intermediate |       65  |      526  |       40  |       45  |      676  | ##              |    9.615% |   77.811% |    5.917% |    6.657% |    8.212% | ##              |   23.985% |    6.889% |   25.157% |   26.946% |           | ##              |    0.790% |    6.390% |    0.486% |    0.547% |           | ## -------------|-----------|-----------|-----------|-----------|-----------|##      shallow |      197  |     7109  |      109  |       91  |     7506  | ##              |    2.625% |   94.711% |    1.452% |    1.212% |   91.181% | ##              |   72.694% |   93.111% |   68.553% |   54.491% |           | ##              |    2.393% |   86.358% |    1.324% |    1.105% |           | ## -------------|-----------|-----------|-----------|-----------|-----------|## Column Total |      271  |     7635  |      159  |      167  |     8232  | ##              |    3.292% |   92.748% |    1.931% |    2.029% |           | ##————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————————--|####

    尖晶石图能够突出比例上比条形图更好的差异。下面我们展示了所得到的尖晶石图。

    类型_区域_计数<-地震%>%分组_按(深度_类型,region)%>%dplyr::summarse(freq=n())xtabs-res<-xtabs(freq~depth_type+region,data=type_region_count)尖晶石图(Xtabs_res)

    Under independence hypothesis,我们可以计算预期频率如下。

    (exp<-独立性表(mt))
    ####东北-西北-东南-西南-深1.646016 46.37391 0.9657434 1.014334中间22.254130 626.97522 13.0568513 13.713800-浅247.099854 6961.65087 144.9774052 152.271866

    我们现在以频率的形式提供应急表。

    tbl_df <- as.data.frame(tbl)colnames(tbl_df) <- c("region",“德斯梯型”“freq”)tbl_df$region<-factor(tbl_df$region)tbl_df$depth_type<-factor(tbl_df$depth_type)tbl_df
    ##区域深度\\\\\35\\35\\35\35\\\35\\\35\\35\\\\\\\35\\\\\\\\\\35; 10 NW浅7109_11 SE SHallow 109 12 SW浅91

    关联模式可以通过vcDextra包中的cmhtest()或筛网图来显示。

    cmhtest(频率~区域+深度_类型,数据=地震)
    ##Cochran Mantel Haenszel地区深度\35\35\35\35\35\35\35\35\35\\\\\\35\35\\35;\\\35\\\\\\\\\35\\\\\\\\6789E-282号

    筛网图突出显示的频率不同于由独立性决定的预期频率。每个矩形的面积总是与预期的频率成正比,但是,观察到的频率由每个矩形中的平方数表示。

    筛子(频率~区域+深度类型,数据=地震,阴影=真,传说=真,labeling=标记_值)

    此外,阴影的颜色和强度以及皮尔逊残差的传说突出了与预期的偏差。

    关联图在前景中偏离了独立性,在这个意义上,每个盒子的面积与(观察到的-预期的)频率成比例。观察到的>预期频率高于基线代表独立性的细胞,而含有低于预期频率的细胞则低于预期频率。

    assoc(~ region + depth_type,数据=地震,阴影=真,传说=真,labeling=标记_值)

    镶嵌图可以应用于N向列联表。对于如图所示的双向工作台,马赛克显示就像一个分组条形图,其中,钢筋的高度(或宽度)表示一个变量的相对频率,而每根钢筋的截面宽度(或高度)表示给定第一个变量的第二个变量的条件频率。

    马赛克(~区域+深度类型,数据=地震,阴影=真,传说=真,labeling=标记剩余量)

    Please note that in above mosaic plot,西北地区深部地震的频率等于零,没有明确的图形表示。There are several labeling option for the mosaic plot,参考表5.1。〔4〕。在我们的例子中,我们选择用Pearson残差值标记细胞。皮尔逊残差的计算如图所示。

    残差<(mt-exp)/sqrt(exp)圆(残差,1)
    ##东北-西北-东南-西南-深5.7-6.8 9.2 29.8-中间9.1-4.0 7.5 8.4-浅3.2 1.8-3.0-5.0

    在哪里?mt是利润表,EXP是基于独立假设的预期值,如前面计算的那样。

    如果您想更详细地了解分类数据探索性分析,阅读参考文献。[4]作为技术内容和示例的宝贵来源。

    如果你有任何问题,请在下面发表评论。

    工具书类

    1. 地震数据集
    2. EathQuake数据集术语
    3. 介绍统计学与R,第二版,P.Dalgaard斯普林格
    4. 离散数据分析与R,M.友好的DMeyerCRC出版社
    5. 确定地震深度

    相关帖子

    1. 地震分析(1/4):定量变量探索性分析
    2. How to calculate the correlation coefficients for more than two variables
    3. R中的六西格玛DMAIC系列——第5部分
    4. Dow Jones Stock Market Index (2/4): Trade Volume Exploratory Analysis
    5. 道琼斯股票市场指数(1/4):对数回报探索性分析

    对客人发帖感兴趣?我们很乐意与我们的社区分享您的代码和想法。

    留下评论作者,please follow the link and comment on their blog: R编程-数据科学+.

    188bet appR博客提供 每日电子邮件更新金宝搏网址 R新闻与 教程关于以下主题: 数据科学大数据, r作业,可视化(可视化) ggplot2箱形图地图动画)程序设计(程序) 演播室SweaveLaTeXSQL日食吉特哈多普刮网)统计 回归主成分分析时间序列交易还有更多…



    如果你走这么远,why not 订阅更新 从站点?选择您的口味: 电子邮件推特1188bet app,或 脸谱网

    Comments are closed.

    搜索R-Blo188bet appggers


    Sponsors

    千万不要错过更新!
    订阅R-Bloggers188bet app接收
    最新R帖子的电子邮件。
    (您将不再看到此消息。)

    单击此处关闭(此弹出窗口将不再出现)