
########################## Kõik aastad NOAA yritused

rm(list = ls())

library(plyr)
library(ggmap)
library(ggplot2)
library(weathermetrics)
library(corrplot)
library(ROCR)
library(Amelia)
library(caret)
library(e1071)
library(gains)
library(mboost)
library(ROSE)


## Home
setwd("C:/Users/Vladimir/Dropbox/master rescue/thesis dev/sheets")
## Work
setwd("C:/Users/vovvie/Dropbox/master rescue/thesis dev/sheets")








###################------------------------ 2010
###----
## Andmete prep
d2010 = read.csv("test2010.csv", header=T, sep=";", fill=T) # 37883
data10 <- d2010[!is.na(d2010$Väljakutse.nr),]
data10 = subset(data10, !duplicated(Väljakutse.nr))

data10 = subset(data10, select=c("Sündmuse.liik..SOS.", "VALLA_NIMI", "Registreerimise.aeg"))
names(data10) <- c("Sündmus", "Vald", "Aeg")


###################------------------------ 2011
###----
## Andmete prep
d2011 = read.csv("test2011.csv", header=T, sep=";", fill=T) # 38445
data11 <- d2011[!is.na(d2011$Väljakutse.nr),]
data11 = subset(data11, !duplicated(Väljakutse.nr))

data11 = subset(data11, select=c("Sündmuse.liik..SOS.", "VALLA_NIMI", "Registreerimise.aeg"))
names(data11) <- c("Sündmus", "Vald", "Aeg")


###################------------------------ 2012
###----
## Andmete prep
d2012 = read.csv("test2012.csv", header=T, sep=";", fill=T) # 34581
data12 <- d2012[!is.na(d2012$Väljakutse.nr),]
data12 = subset(data12, !duplicated(Väljakutse.nr))

data12 = subset(data12, select=c("Sündmuse.liik..SOS.", "VALLA_NIMI", "Registreerimise.aeg"))
names(data12) <- c("Sündmus", "Vald", "Aeg")


###################------------------------ 2013
###----
## Andmete prep
d2013 = read.csv("test2013.csv", header=T, sep=";", fill=T) # 42612
data13 <- d2013[!is.na(d2013$VÄLJAKUTSE_NR),]
data13 = subset(data13, !duplicated(VÄLJAKUTSE_NR))

data13 = subset(data13, select=c("SÜNDMUS", "Vald", "REGISTREERIMISE_AEG"))
names(data13) <- c("Sündmus", "Vald", "Aeg")


###################------------------------ Combining all years into one data frame

allyears <- rbind(data11, data12, data13) # 67904

##### yritused

yritused = read.csv("koik_yritused.csv", header=T, sep=";", fill=T)
yri = subset(yritused, Aasta == '2011' | Aasta == '2012' | Aasta == '2013')
yri = subset(yri, Linn == "Tartu")
ddkoikyri = yri[c("Kuupäev")]

ddkoikyri$Kuupäev <- format(as.Date(ddkoikyri$Kuupäev, format="%d.%m.%Y"), "%Y.%m.%d")

ddkoikyri <- unique(count(ddkoikyri, "Kuupäev"))
names(ddkoikyri) <- c("Day", "pidu")

##### noaa

ilmastik = read.csv("noaa_ilmastik.csv", header=T, sep=",", fill=T) # 37883
log_ilmastik = ilmastik
log_ilmastik$PRCP[log_ilmastik$PRCP < -9000] <- NA
log_ilmastik$SNWD[log_ilmastik$SNWD < -9000] <- NA
log_ilmastik$TAVG[log_ilmastik$TAVG < -9000] <- NA
log_ilmastik$TMAX[log_ilmastik$TMAX < -9000] <- NA
log_ilmastik$TMIN[log_ilmastik$TMIN < -9000] <- NA

# next stuff for ilmastik
log_ilmastik = subset(log_ilmastik, STATION_NAME == 'TARTU EN')
log_ilmastik$Y.M = format(as.Date(as.character(log_ilmastik$DATE), format="%Y%m%d"), "%Y-%m")
log_ilmastik$YEAR = format(as.Date(as.character(log_ilmastik$DATE), format="%Y%m%d"), "%Y")
log_ilmastik$Day = format(as.Date(as.character(log_ilmastik$DATE), format="%Y%m%d"), "%Y.%m.%d")
log_ilmastik$PRCP <- as.numeric(log_ilmastik$PRCP)
log_ilmastik = subset(log_ilmastik, YEAR == '2011' | YEAR == '2012' | YEAR == '2013')
summary(log_ilmastik)

findMeanPrecipitation <- function(arg1){
  ddd = subset(log_ilmastik, Y.M == arg1 & !is.na(PRCP))
  return(mean(ddd$PRCP))
}

findMeanTemp <- function(arg1){
  ddd = subset(log_ilmastik, Y.M == arg1 & !is.na(TAVG))
  return(mean(ddd$TAVG))
}

findMeanMax <- function(arg1){
  ddd = subset(log_ilmastik, Y.M == arg1 & !is.na(TMAX))
  return(mean(ddd$TMAX))
}

log_ilmastik$PRCP[is.na(log_ilmastik$PRCP)] <- lapply(log_ilmastik$Y.M[is.na(log_ilmastik$PRCP)], findMeanPrecipitation)
log_ilmastik$TAVG[is.na(log_ilmastik$TAVG)] <- lapply(log_ilmastik$Y.M[is.na(log_ilmastik$TAVG)], findMeanTemp)
log_ilmastik$TMAX[is.na(log_ilmastik$TMAX)] <- lapply(log_ilmastik$Y.M[is.na(log_ilmastik$TMAX)], findMeanMax)
log_ilmastik$PRCP <- as.numeric(log_ilmastik$PRCP)
log_ilmastik$TAVG <- as.numeric(log_ilmastik$TAVG)
log_ilmastik$TMAX <- as.numeric(log_ilmastik$TMAX)

ilmastik <- subset(log_ilmastik, select=c("Day", "PRCP", "TAVG", "TMAX"))


##### Sündmuste prep
syndmused <- subset(allyears, grepl("liiklus", tolower(Sündmus)))
#syndmused <- subset(allyears, grepl("loodus", tolower(Sündmus)))
"
syndmused <- subset(allyears, grepl('insult', tolower(Sündmus)) | grepl('pea', tolower(Sündmus)) 
| grepl('inimene', tolower(Sündmus)) | grepl('nõrkus', tolower(Sündmus))
| grepl('valu', tolower(Sündmus)) | grepl('kõht', tolower(Sündmus))
| grepl('hingam', tolower(Sündmus)) | grepl('kael', tolower(Sündmus))
| grepl('teadvus', tolower(Sündmus)) | grepl('vere', tolower(Sündmus))
| grepl('jäse', tolower(Sündmus)) | grepl('kõrgus', tolower(Sündmus))
| grepl('iiveldus', tolower(Sündmus)) | grepl('lamab', tolower(Sündmus)) 
| grepl('torke', tolower(Sündmus)) | grepl('rütm', tolower(Sündmus)))
"
#syndmused <- subset(allyears, grepl("tulekahju", tolower(Sündmus)))
#syndmused_d11 <- subset(data11, grepl('inimene', tolower(Sündmus)))
#syndmused <- allyears

## Dates
dd <- subset(syndmused, select=c("Vald", "Aeg"))

dd$Vald <- tolower(dd$Vald)
dd <- subset(dd, Vald == 'tartu')

dd$Paev <- format(as.Date(dd$Aeg, format="%d.%m.%Y"), "%Y.%m.%d")

dd <- unique(count(dd, "Paev")) # 145
names(dd) <- c("Day", "paa")

##### Merge

andmed = merge(dd, ddkoikyri, by = c("Day"), all=T)
andmed = merge(andmed, ilmastik, by = c("Day"), all=T)

andmed$TAVG <- fahrenheit.to.celsius(andmed$TAVG)
andmed$TMAX <- fahrenheit.to.celsius(andmed$TMAX)
andmed$PRCP <- as.numeric(andmed$PRCP)

andmed$pidu[is.na(andmed$pidu)] <- 0
andmed$rescue.event <- ifelse(is.na(andmed$paa), 0, 1)
summary(andmed)

## 


newdata = subset(andmed, Day >= '2011.01.01' & Day <= '2013.12.31')
newdata = newdata[,c(3:7)]
#newdata$rescue.event <- as.factor(newdata$rescue.event)
names(newdata) <- c("Public.events", "Precipitation", "Avg.temperature", "Max.temperature", "Rescue.event.occurs")

summary(newdata)

## splitting the data

train <- newdata[1:776,]
test <- newdata[777:1096,]
table(newdata$Rescue.event.occurs)

## ovun.sample

data_balanced = ovun.sample(Rescue.event.occurs ~ ., data=train, method='under', N=500)$data
table(data_balanced$Rescue.event.occurs)

## logistic regression

logr <- glm(Rescue.event.occurs ~ Public.events, data=train, family=binomial(link='logit'))
summary(logr)

## variable selection

null = glm(Rescue.event.occurs ~ 1, data=train, family=binomial(link='logit'))
full = glm(Rescue.event.occurs ~ .:., data=train, family=binomial(link='logit'))
step(null, scope=list(lower=null, upper=full, direction='forward'))

## ANOVA

anova(logr, test="Chisq")

## boosting

lm.boosted = glmboost(Rescue.event.occurs ~ ., data=train)
summary(lm.boosted)

## ROC

predicted <- predict(logr, newdata=test, type="response")
prediction <- prediction(predicted, test$Rescue.event.occurs)
prediction.performance <- performance(prediction, measure="tpr", x.measure="fpr")
plot(prediction.performance, main='ROC curve, occurrence of nature caused rescue events in Tallinn')

auc <- performance(prediction, measure="auc")
auc <- auc@y.values[[1]]
auc

## confusion matrix

pred = predict(logr, test, type='response')
blah = data.frame(pred)
blah$test <- ifelse(blah$pred >= 0.5, 1, 0)

confusionMatrix(blah$test, test$Rescue.event.occurs)


## Gain

nrow(subset(test, Rescue.event.occurs == 1))

predicted <- predict(logr, newdata=test, type="response")
prediction <- prediction(predicted, test$Rescue.event.occurs)
gain <- performance(prediction, "tpr", "rpp")
plot(gain, main="Gain chart, fires in Tartu")
lines(c(0,1), c(0,1), col='red')
legend('bottomright', legend=c('Ideal model', 'Our model', 'No model'), col=c('green', 'black', 'red'), lty=1, cex=1)

# nature Tallinn (52)
#lines(seq(0, 1, 0.05), c(0, 0.31, 0.615, 0.923, rep(1, each=17)), col='green')

# road Tallinn (113)
#lines(seq(0, 1, 0.05), c(0, 0.14, 0.28, 0.42, 0.566, 0.71, 0.85, 0.99, rep(1, each=13)), col='green')

# nature Tartu (18)
#lines(seq(0, 1, 0.05), c(0, 0.89, rep(1, each=19)), col='green')

# road Tartu (58)
#lines(seq(0, 1, 0.05), c(0, 0.276, 0.552, 0.828, rep(1, each=17)), col='green')

# fire Tartu (220)
#lines(seq(0, 1, 0.05), c(0, 0.07, 0.145, 0.218, 0.29, 0.36, 0.436, 0.51, 0.581, 0.654, 0.73,  0.8, 0.87, 0.94, rep(1, each=7)), col='green')

# fire Kuressaare (31)
#lines(seq(0, 1, 0.05), c(0, 0.516, rep(1, each=19)), col='green')


