探索性数据分析week4

最新推荐文章于 2024-04-06 21:43:56 发布

林思

最新推荐文章于 2024-04-06 21:43:56 发布

阅读量845

点赞数

分类专栏： datascience 探索性数据分析

本文链接：https://blog.csdn.net/u014596936/article/details/38293197

版权

datascience 同时被 2 个专栏收录

25 篇文章 0 订阅

订阅专栏

探索性数据分析

4 篇文章 0 订阅

订阅专栏

本week主要为两case

其一为聚类

其二为空气污染案例

数据：2012和1999年pm2.5的值，以观察两者之间有否不同，从而论证近年来空气洁净计划的成果

step 1：先大概瞄一眼数据

因此决定读取数据方法为

pm0<- read.table("RD_501_88101_1999-0.txt", comment.char ="#", header =FALSE, sep ="|", na.strings ="")

head(pm0)
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14
1 RD I 1 27 1 88101 1 7 105 120 19990103 00:00 NA AS
2 RD I 1 27 1 88101 1 7 105 120 19990106 00:00 NA AS
3 RD I 1 27 1 88101 1 7 105 120 19990109 00:00 NA AS
4 RD I 1 27 1 88101 1 7 105 120 19990112 00:00 8.841
5 RD I 1 27 1 88101 1 7 105 120 19990115 00:00 14.920
6 RD I 1 27 1 88101 1 7 105 120 19990118 00:00 3.878
V15 V16 V17 V18 V19 V20 V21 V22 V23 V24 V25 V26 V27 V28
1 3 NA NA NA NA NA NA NA NA NA NA NA NA
2 3 NA NA NA NA NA NA NA NA NA NA NA NA
3 3 NA NA NA NA NA NA NA NA NA NA NA NA
4 3 NA NA NA NA NA NA NA NA NA NA NA NA
5 3 NA NA NA NA NA NA NA NA NA NA NA NA
6 3 NA NA NA NA NA NA NA NA NA NA NA NA

cnames<-strsplit(cnames,"|", fixed =TRUE) 此处fixed=T表示准确匹配

cnames<-strsplit(cnames,"|",fixed=T)
> cnames
[[1]]
[1] "# RD"
[2] "Action Code"
[3] "State Code"
[4] "County Code"
[5] "Site ID"
[6] "Parameter"
[7] "POC"
[8] "Sample Duration"
[9] "Unit"
[10] "Method"
[11] "Date"
[12] "Start Time"
[13] "Sample Value"
[14] "Null Data Code"
[15] "Sampling Frequency"
[16] "Monitor Protocol (MP) ID"
[17] "Qualifier - 1"
[18] "Qualifier - 2"
[19] "Qualifier - 3"
[20] "Qualifier - 4"
[21] "Qualifier - 5"
[22] "Qualifier - 6"
[23] "Qualifier - 7"
[24] "Qualifier - 8"
[25] "Qualifier - 9"
[26] "Qualifier - 10"
[27] "Alternate Method Detectable Limit"
[28] "Uncertainty"

names(pm0)<-make.names(cnames[[1]]) 将名字中间的空格去掉，换之以.以符合名字的命名规则

head(pm0)
X..RD Action.Code State.Code County.Code Site.ID Parameter
1 RD I 1 27 1 88101
2 RD I 1 27 1 88101
3 RD I 1 27 1 88101
4 RD I 1 27 1 88101
5 RD I 1 27 1 88101
6 RD I 1 27 1 88101
POC Sample.Duration Unit Method Date Start.Time
1 1 7 105 120 19990103 00:00
2 1 7 105 120 19990106 00:00
3 1 7 105 120 19990109 00:00
4 1 7 105 120 19990112 00:00
5 1 7 105 120 19990115 00:00
6 1 7 105 120 19990118 00:00
Sample.Value Null.Data.Code Sampling.Frequency
1 NA AS 3
2 NA AS 3
3 NA AS 3
4 8.841 3
5 14.920 3
6 3.878 3
Monitor.Protocol..MP..ID Qualifier...1 Qualifier...2
1 NA NA
2 NA NA
3 NA NA
4 NA NA
5 NA NA
6 NA NA
Qualifier...3 Qualifier...4 Qualifier...5 Qualifier...6
1 NA NA NA NA
2 NA NA NA NA
3 NA NA NA NA
4 NA NA NA NA
5 NA NA NA NA
6 NA NA NA NA
Qualifier...7 Qualifier...8 Qualifier...9 Qualifier...10
1 NA NA NA NA
2 NA NA NA NA
3 NA NA NA NA
4 NA NA NA NA
5 NA NA NA NA
6 NA NA NA NA
Alternate.Method.Detectable.Limit Uncertainty
1 NA NA
2 NA NA
3 NA NA
4 NA NA
5 NA NA
6 NA NA

x0 <- pm0$Sample.Value
> class(x0)
[1] "numeric"
> str(x0)
num [1:117421] NA NA NA 8.84 14.92 ...
> summary(x0)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 7.20 11.50 13.74 17.90 157.10 13217
> sum(is.na(x0))
[1] 13217
> mean(is.na(x0))
[1] 0.1125608

可见出现了miss data

 
 
  
  pm1 
  
  <- read.table
  
  (
  
  "RD_501_88101_2012-0.txt"
  
  , comment.char 
  
  = 
  
  "#"
  
  , header 
  
  = 
  
  FALSE
  
  , sep 
  
  = 
  
  "|"
  
  , na.strings 
  
  = 
  
  ""
  
  , nrow 
  
  = 
  
  1304290
  
  )
 
 
 
 
  
  names
  
  (pm1
  
  ) 
  
  <- 
  
  make.names
  
  (cnames
  
  [[
  
  1
  
  ]])
 
 
 
 
  
  head
  
  (pm1
  
  )
 
 
 
 
  
  dim
  
  (pm1
  
  )
 
 
 
 
  
  x1 
  
  <- pm1
  
  $Sample.Value
 
 
 
 
  
  class
  
  (x1
  
  )

summary(x1)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-10.00 4.10 8.00 9.79 13.00 985.00 78313
> summary(x0)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 7.20 11.50 13.74 17.90 157.10 13217

mean(is.na(x1))## Are missing values important here?

[1] 0.06004263

boxplot(x1,x2)

boxplot(log10(x0), log10(x1))

Warning messages:
1: In boxplot.default(log10(x0), log10(x1)) : NaNs produced
2: In bplt(at[i], wid = width[i], stats = z$stats[, i], out = z$out[z$group == :
Outlier (-Inf) in boxplot 1 is not drawn
3: In bplt(at[i], wid = width[i], stats = z$stats[, i], out = z$out[z$group == :
Outlier (-Inf) in boxplot 2 is not drawn

产生了一堆warning，原因在于原始数据中存在负值

pm2.5为负值，想想有些小奇怪

 
 
  
  negative 
  
  <- x1 
  
  < 
  
  0
 
 
 
 
  
  sum
  
  (negative
  
  , na.rm 
  
  = 
  
  T
  
  )
 
 
 
 
  
  mean
  
  (negative
  
  , na.rm 
  
  = 
  
  T
  
  )

共有24463个负值，平均为0.01995388

来瞄一下负值出现的原因

 
 
  
  dates 
  
  <- pm1
  
  $Date
 
 
 
 
  
  str
  
  (dates
  
  )
 
 
 
 
  
  dates 
  
  <- 
  
  as.Date
  
  (
  
  as.character
  
  (dates
  
  ), 
  
  "%Y%m%d"
  
  )
 
 
 
 
  
  str
  
  (dates
  
  )
 
 
 
 
  
  hist
  
  (dates
  
  , 
  
  "month"
  
  )  
  
  ## Check what's going on in months 1--6

从图中可得出一些结论

好似冬天的miss要多于夏天的

好吧，老师说与其瞄全国的中位数，何不就锁定某个位置的某一个监控器呢，或者锁定某个州呢

 
 
  
  ## Find a monitor for New York State that exists in both datasets
 
 
 
 
  
  site0 
  
  <- 
  
  unique
  
  (
  
  subset
  
  (pm0
  
  , State.Code 
  
  == 
  
  36
  
  , 
  
  c
  
  (County.Code
  
  , Site.ID
  
  )))
 
 
 
 
  
  site1 
  
  <- 
  
  unique
  
  (
  
  subset
  
  (pm1
  
  , State.Code 
  
  == 
  
  36
  
  , 
  
  c
  
  (County.Code
  
  , Site.ID
  
  )))
 
 
 
 
  
  site0 
  
  <- 
  
  paste
  
  (site0
  
  [,
  
  1
  
  ], site0
  
  [,
  
  2
  
  ], sep 
  
  = 
  
  "."
  
  )
 
 
 
 
  
  site1 
  
  <- 
  
  paste
  
  (site1
  
  [,
  
  1
  
  ], site1
  
  [,
  
  2
  
  ], sep 
  
  = 
  
  "."
  
  )
 
 
 
 
  
  str
  
  (site0
  
  )
 
 
 
 
  
  str
  
  (site1
  
  )
 
 
 
 
  
  both 
  
  <- 
  
  intersect
  
  (site0
  
  , site1
  
  )
 
 
 
 
  
  print
  
  (both
  
  )

 
 
  
  ## Find how many observations available at each monitor
 
 
 
 
  
  pm0
  
  $county.site 
  
  <- 
  
  with
  
  (pm0
  
  , 
  
  paste
  
  (County.Code
  
  , Site.ID
  
  , sep 
  
  = 
  
  "."
  
  ))
 
 
 
 
  
  pm1
  
  $county.site 
  
  <- 
  
  with
  
  (pm1
  
  , 
  
  paste
  
  (County.Code
  
  , Site.ID
  
  , sep 
  
  = 
  
  "."
  
  ))
 
 
 
 
  
  cnt0 
  
  <- 
  
  subset
  
  (pm0
  
  , State.Code 
  
  == 
  
  36 
  
  & county.site 
  
  %in% both
  
  )
 
 
 
 
  
  cnt1 
  
  <- 
  
  subset
  
  (pm1
  
  , State.Code 
  
  == 
  
  36 
  
  & county.site 
  
  %in% both
  
  )
 
 
 
 
  
  sapply
  
  (
  
  split
  
  (cnt0
  
  , cnt0
  
  $county.site
  
  ), 
  
  nrow
  
  )
 
 
 
 
  
  sapply
  
  (
  
  split
  
  (cnt1
  
  , cnt1
  
  $county.site
  
  ), 
  
  nrow
  
  )
 
 
 
 
  
  

 
 
 
 
  
  ## Choose county 63 and side ID 2008
 
 
 
 
  
  pm1sub 
  
  <- 
  
  subset
  
  (pm1
  
  , State.Code 
  
  == 
  
  36 
  
  & County.Code 
  
  == 
  
  63 
  
  & Site.ID 
  
  == 
  
  2008
  
  )
 
 
 
 
  
  pm0sub 
  
  <- 
  
  subset
  
  (pm0
  
  , State.Code 
  
  == 
  
  36 
  
  & County.Code 
  
  == 
  
  63 
  
  & Site.ID 
  
  == 
  
  2008
  
  )
 
 
 
 
  
  dim
  
  (pm1sub
  
  )
 
 
 
 
  
  dim
  
  (pm0sub
  
  )
 
 
 
 
  
  

 
 
 
 
  
  ## Plot data for 2012
 
 
 
 
  
  dates1 
  
  <- pm1sub
  
  $Date
 
 
 
 
  
  x1sub 
  
  <- pm1sub
  
  $Sample.Value
 
 
 
 
  
  plot
  
  (dates1
  
  , x1sub
  
  )
 
 
 
 
  
  dates1 
  
  <- 
  
  as.Date
  
  (
  
  as.character
  
  (dates1
  
  ), 
  
  "%Y%m%d"
  
  )
 
 
 
 
  
  str
  
  (dates1
  
  )
 
 
 
 
  
  plot
  
  (dates1
  
  , x1sub
  
  )
 
 
 
 
  
  

 
 
 
 
  
  ## Plot data for 1999
 
 
 
 
  
  dates0 
  
  <- pm0sub
  
  $Date
 
 
 
 
  
  dates0 
  
  <- 
  
  as.Date
  
  (
  
  as.character
  
  (dates0
  
  ), 
  
  "%Y%m%d"
  
  )
 
 
 
 
  
  x0sub 
  
  <- pm0sub
  
  $Sample.Value
 
 
 
 
  
  plot
  
  (dates0
  
  , x0sub
  
  )
 
 
 
 
  
  

 
 
 
 
  
  ## Plot data for both years in same panel
 
 
 
 
  
  par
  
  (mfrow 
  
  = 
  
  c
  
  (
  
  1
  
  , 
  
  2
  
  ), mar 
  
  = 
  
  c
  
  (
  
  4
  
  , 
  
  4
  
  , 
  
  2
  
  , 
  
  1
  
  ))
 
 
 
 
  
  plot
  
  (dates0
  
  , x0sub
  
  , pch 
  
  = 
  
  20
  
  )
 
 
 
 
  
  abline
  
  (h 
  
  = median
  
  (x0sub
  
  , na.rm 
  
  = 
  
  T
  
  ))
 
 
 
 
  
  plot
  
  (dates1
  
  , x1sub
  
  , pch 
  
  = 
  
  20
  
  )  
  
  ## Whoa! Different ranges
 
 
 
 
  
  abline
  
  (h 
  
  = median
  
  (x1sub
  
  , na.rm 
  
  = 
  
  T
  
  ))
 
 
 
 
  
  

 
 
 
 
  
  ## Find global range
 
 
 
 
  
  rng 
  
  <- 
  
  range
  
  (x0sub
  
  , x1sub
  
  , na.rm 
  
  = 
  
  T
  
  )
 
 
 
 
  
  rng
 
 
 
 
  
  par
  
  (mfrow 
  
  = 
  
  c
  
  (
  
  1
  
  , 
  
  2
  
  ), mar 
  
  = 
  
  c
  
  (
  
  4
  
  , 
  
  4
  
  , 
  
  2
  
  , 
  
  1
  
  ))
 
 
 
 
  
  plot
  
  (dates0
  
  , x0sub
  
  , pch 
  
  = 
  
  20
  
  , ylim 
  
  = rng
  
  )
 
 
 
 
  
  abline
  
  (h 
  
  = median
  
  (x0sub
  
  , na.rm 
  
  = 
  
  T
  
  ))
 
 
 
 
  
  plot
  
  (dates1
  
  , x1sub
  
  , pch 
  
  = 
  
  20
  
  , ylim 
  
  = rng
  
  )
 
 
 
 
  
  abline
  
  (h 
  
  = median
  
  (x1sub
  
  , na.rm 
  
  = 
  
  T
  
  ))
 
 
 
 
  
  

 
 
 
 
  
  ## Show state-wide means and make a plot showing trend
 
 
 
 
  
  head
  
  (pm0
  
  )
 
 
 
 
  
  mn0 
  
  <- 
  
  with
  
  (pm0
  
  , 
  
  tapply
  
  (Sample.Value
  
  , State.Code
  
  , 
  
  mean
  
  , na.rm 
  
  = 
  
  T
  
  ))
 
 
 
 
  
  str
  
  (mn0
  
  )
 
 
 
 
  
  summary
  
  (mn0
  
  )
 
 
 
 
  
  mn1 
  
  <- 
  
  with
  
  (pm1
  
  , 
  
  tapply
  
  (Sample.Value
  
  , State.Code
  
  , 
  
  mean
  
  , na.rm 
  
  = 
  
  T
  
  ))
 
 
 
 
  
  str
  
  (mn1
  
  )
 
 
 
 
  
  

 
 
 
 
  
  ## Make separate data frames for states / years
 
 
 
 
  
  d0 
  
  <- 
  
  data.frame
  
  (state 
  
  = 
  
  names
  
  (mn0
  
  ), mean 
  
  = mn0
  
  )
 
 
 
 
  
  d1 
  
  <- 
  
  data.frame
  
  (state 
  
  = 
  
  names
  
  (mn1
  
  ), mean 
  
  = mn1
  
  )
 
 
 
 
  
  mrg 
  
  <- 
  
  merge
  
  (d0
  
  , d1
  
  , by 
  
  = 
  
  "state"
  
  )
 
 
 
 
  
  dim
  
  (mrg
  
  )
 
 
 
 
  
  head
  
  (mrg
  
  )
 
 
 
 
  
  

 
 
 
 
  
  ## Connect lines
 
 
 
 
  
  par
  
  (mfrow 
  
  = 
  
  c
  
  (
  
  1
  
  , 
  
  1
  
  ))
 
 
 
 
  
  with
  
  (mrg
  
  , plot
  
  (
  
  rep
  
  (
  
  1
  
  , 
  
  52
  
  ), mrg
  
  [, 
  
  2
  
  ], xlim 
  
  = 
  
  c
  
  (
  
  .5
  
  , 
  
  2.5
  
  )))
 
 
 
 
  
  with
  
  (mrg
  
  , points
  
  (
  
  rep
  
  (
  
  2
  
  , 
  
  52
  
  ), mrg
  
  [, 
  
  3
  
  ]))
 
 
 
 
  
  segments
  
  (
  
  rep
  
  (
  
  1
  
  , 
  
  52
  
  ), mrg
  
  [, 
  
  2
  
  ], 
  
  rep
  
  (
  
  2
  
  , 
  
  52
  
  ), mrg
  
  [, 
  
  3
  
  ])