MSHA.Rmd

Analysis of MSHA Inspections, Violations, and Accident Data
========================================================

_An open and ongoing analysis of the data provided by the [US Mine Safety & Health Administration](http://www.msha.gov/OpenGovernmentData/OGIMSHA.asp). See [github.com/mwfrost/MSHA](https://github.com/mwfrost/MSHA) for details._

```{r}
require(reshape)
require(plyr)
require(ggplot2)
require(lubridate)
require(xtable)
require(scales)
require(RJSONIO)

setwd('~/Code/MSHA/')
source('./lib/calendarHeat.r')
Sys.setlocale(locale="C")
```
## Load Data

### Load Mine Data
```{r}
mdat <- read.table('./data/msha_source/Mines.TXT', header=T, sep="|", fill=T, as.is=c(1:59),quote="\"")
t(head(mdat,3))
mdat_wv <- subset(mdat, STATE == 'WV' & COAL_METAL_IND == 'C')

```

*rmongodb* approach

Write each mine to a document in a mongodb collection

```{r mongo.write.mines}
require(rmongodb)

data.frame.to.mongo <- function(df, db, col, fields=names(df)) {
  mongo <- mongo.create()
  if (mongo.is.connected(mongo)) {
    for(i in 1:nrow(df)) {
        b <- mongo.bson.from.list(df[i,fields])
        #print(b)
        #print('----------------')
        if (!is.null(b)){
          mongo.insert(mongo, paste(db, col, sep='.'), b)
          }
        else {
          print('NULL BSON object')
          print(t(df[i,]))
              }
      }  
  }
}



db <- 'msha'
col <- 'mines'
fields <- setdiff(names(mdat), 'DIRECTIONS_TO_MINE')
# data.frame.to.mongo(mdat, db, col, fields)

```

```{r rmongodb.2}
  mongo <- mongo.create()
buf <- mongo.bson.buffer.create()
        mongo.bson.buffer.append(buf, "MINE_ID", 1519137)
        query <- mongo.bson.from.buffer(buf)

mongo.find(mongo, 'msha', query)

```

PITA.
Try the RMongo library instead.

Some sample Mongo aggregations on the Mines collection
Based on these examples: 

http://docs.mongodb.org/manual/reference/sql-aggregation-comparison/

http://docs.mongodb.org/manual/reference/sql-comparison/

Example mongodb search: `db.mines.find({MINE_ID : 1519137})`

```{r rmongo}
require(RMongo)

mongo <- mongoDbConnect(dbName="msha", host="localhost",port='27017')

output <- dbAggregate(mongo, "mines",
                       c('{ 
                         $match: { NO_EMPLOYEES: {$gt: 0} } 
                         }'
                        ,'{ 
                         $group: { 
                         _id: "$PRIMARY_SIC", 
                         employees : {$sum:"$NO_EMPLOYEES"}
                         } 
                         }'
                      )
                      )

print(output)

output <- lapply(output, fromJSON)

output <- do.call(rbind, output)
output <- data.frame(output)

system.time(
  dbAggregate(mongo, "mines",
                       c('{ $match: { NO_EMPLOYEES: {$gt: 0} } }'
                        ,'{ $group: { _id: "$PRIMARY_SIC", employees : {$sum:"$NO_EMPLOYEES"}} }'
                      )
                      )
)

```


```{r ddply.employees}
 system.time(
  ddply(mdat[mdat$NO_EMPLOYEES > 0,], .(PRIMARY_SIC), summarize, employees=sum(NO_EMPLOYEES))
           )
```


```{r rmongo.sample.queries}
dbGetQuery(mongo, 'mines', '{STATE: "WV", "COAL_METAL_IND" :"C"}')

# Robinson Run
dbGetQuery(mongo, 'mines', '{MINE_ID: 4601318}', '{MINE_ID:1 }')

```

```{r mongo.functions}
function minebyid(id) {
    mine <- db.mines.findOne({_id: id})
   # mine.inspections <- db.comments.find({postId: id}).toArray()
    return mine
}


```


### Inspections
To periodically rebuild this WV set of records, run `./01-extract-wv.R`
```{r}
idat <- read.csv("./data/wv_idat.csv")
idat <- merge(idat, mdat_wv[,c('MINE_ID','CURRENT_MINE_NAME','CURRENT_CONTROLLER_NAME')])
t(head(idat,3))

```


```{r mongo.inspections}
# add all the Robinson Run inspections
rr.i <- idat[idat$MINE_ID == 4601318,]
#data.frame.to.mongo(rr.i, 'msha', 'inspections')
#  data.frame.to.mongo(idat, 'msha', 'inspections')
system.time(
  dbAggregate(mongo, "inspections",
                       c('{ $match: { CURRENT_CONTROLLER_NAME: "CONSOL Energy Inc" } }'
                        ,'{ $group: { _id: "$CURRENT_MINE_NAME", inspections : {$sum:1}} }'
                      )
                      )
)

system.time(
  ddply(
    subset(idat, CURRENT_CONTROLLER_NAME == 'CONSOL Energy Inc'), 
    .(CURRENT_MINE_NAME) , nrow
    )
)

```


### Accidents
```{r}
adat <- read.table('./data/msha_source/Accidents.TXT', header=T, sep="|", fill=T, quote="\"",comment.char = "")

data.frame.to.mongo(adat, 'msha', 'accidents')

# WV is FIPS 54
adat <- subset(adat, FIPS_STATE_CD==54)
adat <- merge(adat, mdat_wv[,c('MINE_ID','CURRENT_MINE_NAME','CURRENT_CONTROLLER_NAME')])
adat$accident.key <- seq(1:nrow(adat))
t(head(adat,3))
```

### Violations

```{r}
vdat <- read.csv("./data/wv_vdat.csv")
vdat <- merge(vdat, mdat_wv[,c('MINE_ID','CURRENT_MINE_NAME','CURRENT_CONTROLLER_NAME')])
t(head(vdat,3))
```


Reduce mdat_wv to only mines that are not abandoned
```{r}
mdat_wv <- mdat_wv[grep('Abandoned.*',mdat_wv$CURRENT_MINE_STATUS, invert=TRUE),]
nrow(mdat_wv)
```
Reduce inspection and violation tables accordingly
```{r}
idat <- subset(idat, MINE_ID %in% mdat_wv$MINE_ID)
vdat <- subset(vdat, MINE_ID %in% mdat_wv$MINE_ID)
nrow(idat)
nrow(vdat)
```
## Create a time series by mine and month

Add real dates to the various data.frames

```{r}

mdat_wv$current.status.dt <- mdy(mdat_wv$CURRENT_STATUS_DT)
idat$inspection.begin.dt <- mdy(idat$INSPECTION_BEGIN_DT)
vdat$violation.occur.dt <- mdy(vdat$VIOLATION_OCCUR_DT)
adat$accident.dt <- mdy(adat$ACCIDENT_DT)

names(adat)
names(idat)
names(vdat)
names(mdat)
```

Count violations by day.

```{r}
vxday <- ddply(vdat, .(violation.occur.dt), nrow)
names(vxday) <- c('date','vcount')
```

Make a calendar heat map

```{r}
calendarHeat(subset(vxday, year(date)==2011)$date , subset(vxday, year(date)==2011)$vcount)

```


Normalize inspections, violations, and accidents into "event" objects

```{r}
events <- 
  rbind(
  data.frame(event='inspection', MINE_ID = idat$MINE_ID, event.day=idat$inspection.begin.dt, event.key=idat$EVENT_NO)
  ,
  data.frame(event='violation', MINE_ID = vdat$MINE_ID, event.day=vdat$violation.occur.dt,event.key=vdat$EVENT_NO) 
  ,
  data.frame(event='accident', MINE_ID = adat$MINE_ID, event.day=adat$accident.dt, event.key=adat$accident.key) 
  )

```
edays, short for event days, is a temporary hack


```{r}
edays<-cast(events, event.day+MINE_ID~event, value='event.key', fun.aggregate='min')

eventsRR <- subset(edays, MINE_ID == 4601318)

eventsRR$day.decimal <- decimal_date(eventsRR$event.day)

eventsRR$lastinsp <- ifelse(eventsRR$inspection > 0 & eventsRR$inspection < Inf, eventsRR$day.decimal , cummax(ifelse(eventsRR$inspection > 0 & eventsRR$inspection < Inf, eventsRR$day.decimal, 0))) 

eventsRR$lastviol <- ifelse(eventsRR$violation > 0 & eventsRR$violation < Inf, eventsRR$day.decimal , cummax(ifelse(eventsRR$violation > 0 & eventsRR$violation < Inf, eventsRR$day.decimal, 0))) 


eventsRR$lastacc <- ifelse(eventsRR$accident > 0 & eventsRR$accident < Inf, eventsRR$day.decimal , cummax(ifelse(eventsRR$accident > 0 & eventsRR$accident < Inf, eventsRR$day.decimal, 0))) 


eventsRR$sinceinsp <- eventsRR$day.decimal - eventsRR$lastinsp 
eventsRR$sinceviol <- eventsRR$day.decimal - eventsRR$lastviol
eventsRR$sinceacc <- eventsRR$day.decimal - eventsRR$lastacc

eventsRR$year <- floor(eventsRR$day.decimal)

ggplot(subset(eventsRR, year> 2000)) + 
  geom_line(aes(x=day.decimal - year, y=sinceinsp),color='blue') +
#  geom_line(aes(x=day.decimal - year, y=sinceviol),color='black') +
  geom_line(aes(x=day.decimal - year, y=sinceacc),color='red') +
 facet_wrap('year')

```