-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path6final_models.R
140 lines (107 loc) · 4.01 KB
/
6final_models.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
# Script to fit models on data
library(MASS)
library(stsm)
# Load libraries ----
library(rstanarm)
library(brms) # for models
library(bayesplot)
library(ggplot2)
library(dplyr)
library(tidybayes)
library(modelr)
#library(tidyverse)
library(caret)
library(readr)
library(purrr)
library(parallel)
##########################################################################
# FITTING STANDARD NEGATIVE BINOMIAL MODEL
##########################################################################
negbinner <- function(x, theta = 1.5, n = 60){
# create training and test sets
# set seed
set.seed(456)
trainIndex <- round(0.8*length(x$count))
# Create the data sets
fireTrain <- x[1:trainIndex,]
fireTest <- x[(trainIndex+1):n,]
# Fit model on training set
glmNB <- MASS::glm.nb(count ~ max_temp +
rainfall, data = fireTrain,
link = "log")
# Predict on training set
predictions_train <- predict(glmNB,
newdata = fireTrain, type = "response")
# Predict on testing set
predictions_test <- predict(glmNB,
newdata = fireTest, type = "response")
# get the rmse
train_rmse <- caret::RMSE(round(predictions_train),fireTrain$count)
test_rmse <- caret::RMSE(round(predictions_test),fireTest$count)
# get the MASE
test_mase <- Metrics::mase(actual = fireTest$count,
predicted = round(predictions_test))
# Get the bias
test_bias <- Metrics::bias(actual = fireTest$count,
predicted = round(predictions_test))
# Dispersion parameter
cbind(rmse_train = train_rmse, rmse_test = test_rmse, mase_test = test_mase,
bias_test = test_bias,
theta = theta, n = n)
}
#########################################################################
##########################################################################
# FITTING BAYESIAN NEGATIVE BINOMIAL MODEL
##########################################################################
stanbinner <- function(x, theta = 1.5, n = 60){
# create training and test sets
# set seed
set.seed(456)
# Add time by month index
x$time <- rep(1:12, length.out = n)
trainIndex <- round(0.8*length(x$count))
# Create the datasets
fireTrain <- x[1:trainIndex,]
fireTest <- x[(trainIndex+1):n,]
# Add prior means
get_prior_means <- function(x){
library(dplyr)
x %>% group_by(time) %>%
summarize(count_mean = mean(count)) %>%
mutate(costhet = cos((2*12*pi)/time),
sinthet = sin((2*12*pi)/time))%>%
data.frame()
}
p_means = get_prior_means(fireTrain)
# Join means to train data
fireTrain2 <- fireTrain %>%
inner_join(p_means, by = "time")
# Join means to test data
fireTest2 <- fireTest %>%
inner_join(p_means, by = "time")
# Fit model on training set
stanNB <- rstanarm::stan_glm.nb(count ~ max_temp+
rainfall + costhet + sinthet,
data = fireTrain2,
link = "log")
# Predict on training set
predictions_train <- predict(stanNB,
newdata = fireTrain2, type = "response")
# Predict on testing set
predictions_test <- predict(stanNB,
newdata = fireTest2, type = "response")
# get the rmse
train_rmse <- caret::RMSE(round(predictions_train),fireTrain2$count)
test_rmse <- caret::RMSE(round(predictions_test),fireTest2$count)
# get the MASE
test_mase <- Metrics::mase(actual = fireTest2$count,
predicted = round(predictions_test))
# Get the bias
test_bias <- Metrics::bias(actual = fireTest2$count,
predicted = round(predictions_test))
# Dispersion parameter
cbind(rmse_train = train_rmse, rmse_test = test_rmse, mase_test = test_mase,
bias_test = test_bias,
theta = theta, n = n)
}
##########################################################################