1
- # # ----echo = FALSE-------------------------------------------------------------
1
+ # # ----echo = FALSE------------------------------------------------------
2
2
NOT_CRAN <- identical(tolower(Sys.getenv(" NOT_CRAN" )), " true" )
3
3
knitr :: opts_chunk $ set(
4
4
collapse = TRUE ,
@@ -7,10 +7,11 @@ knitr::opts_chunk$set(
7
7
eval = NOT_CRAN
8
8
)
9
9
10
- # # ----setup, include=FALSE-----------------------------------------------------
10
+
11
+ # # ----setup, include=FALSE----------------------------------------------
11
12
knitr :: opts_chunk $ set(
12
13
echo = TRUE ,
13
- dpi = 150 ,
14
+ dpi = 100 ,
14
15
fig.asp = 0.8 ,
15
16
fig.width = 6 ,
16
17
out.width = " 60%" ,
@@ -19,45 +20,54 @@ library(mvgam)
19
20
library(ggplot2 )
20
21
theme_set(theme_bw(base_size = 12 , base_family = ' serif' ))
21
22
22
- # # -----------------------------------------------------------------------------
23
+
24
+ # # ----------------------------------------------------------------------
23
25
simdat <- sim_mvgam(n_series = 4 , T = 24 , prop_missing = 0.2 )
24
26
head(simdat $ data_train , 16 )
25
27
26
- # # -----------------------------------------------------------------------------
28
+
29
+ # # ----------------------------------------------------------------------
27
30
class(simdat $ data_train $ series )
28
31
levels(simdat $ data_train $ series )
29
32
30
- # # -----------------------------------------------------------------------------
33
+
34
+ # # ----------------------------------------------------------------------
31
35
all(levels(simdat $ data_train $ series ) %in% unique(simdat $ data_train $ series ))
32
36
33
- # # -----------------------------------------------------------------------------
37
+
38
+ # # ----------------------------------------------------------------------
34
39
summary(glm(y ~ series + time ,
35
40
data = simdat $ data_train ,
36
41
family = poisson()))
37
42
38
- # # -----------------------------------------------------------------------------
43
+
44
+ # # ----------------------------------------------------------------------
39
45
summary(gam(y ~ series + s(time , by = series ),
40
46
data = simdat $ data_train ,
41
47
family = poisson()))
42
48
43
- # # -----------------------------------------------------------------------------
49
+
50
+ # # ----------------------------------------------------------------------
44
51
gauss_dat <- data.frame (outcome = rnorm(10 ),
45
52
series = factor (' series1' ,
46
53
levels = ' series1' ),
47
54
time = 1 : 10 )
48
55
gauss_dat
49
56
50
- # # -----------------------------------------------------------------------------
57
+
58
+ # # ----------------------------------------------------------------------
51
59
gam(outcome ~ time ,
52
60
family = betar(),
53
61
data = gauss_dat )
54
62
55
- # # ----error=TRUE---------------------------------------------------------------
63
+
64
+ # # ----error=TRUE--------------------------------------------------------
56
65
mvgam(outcome ~ time ,
57
66
family = betar(),
58
67
data = gauss_dat )
59
68
60
- # # -----------------------------------------------------------------------------
69
+
70
+ # # ----------------------------------------------------------------------
61
71
# A function to ensure all timepoints within a sequence are identical
62
72
all_times_avail = function (time , min_time , max_time ){
63
73
identical(as.numeric(sort(time )),
@@ -81,18 +91,21 @@ if(any(checked_times$all_there == FALSE)){
81
91
cat(' All series have observations at all timepoints :)' )
82
92
}
83
93
84
- # # -----------------------------------------------------------------------------
94
+
95
+ # # ----------------------------------------------------------------------
85
96
bad_times <- data.frame (time = seq(1 , 16 , by = 2 ),
86
97
series = factor (' series_1' ),
87
98
outcome = rnorm(8 ))
88
99
bad_times
89
100
90
- # # ----error = TRUE-------------------------------------------------------------
101
+
102
+ # # ----error = TRUE------------------------------------------------------
91
103
get_mvgam_priors(outcome ~ 1 ,
92
104
data = bad_times ,
93
105
family = gaussian())
94
106
95
- # # -----------------------------------------------------------------------------
107
+
108
+ # # ----------------------------------------------------------------------
96
109
bad_times %> %
97
110
dplyr :: right_join(expand.grid(time = seq(min(bad_times $ time ),
98
111
max(bad_times $ time )),
@@ -101,12 +114,14 @@ bad_times %>%
101
114
dplyr :: arrange(time ) - > good_times
102
115
good_times
103
116
104
- # # ----error = TRUE-------------------------------------------------------------
117
+
118
+ # # ----error = TRUE------------------------------------------------------
105
119
get_mvgam_priors(outcome ~ 1 ,
106
120
data = good_times ,
107
121
family = gaussian())
108
122
109
- # # -----------------------------------------------------------------------------
123
+
124
+ # # ----------------------------------------------------------------------
110
125
bad_levels <- data.frame (time = 1 : 8 ,
111
126
series = factor (' series_1' ,
112
127
levels = c(' series_1' ,
@@ -115,85 +130,100 @@ bad_levels <- data.frame(time = 1:8,
115
130
116
131
levels(bad_levels $ series )
117
132
118
- # # ----error = TRUE-------------------------------------------------------------
133
+
134
+ # # ----error = TRUE------------------------------------------------------
119
135
get_mvgam_priors(outcome ~ 1 ,
120
136
data = bad_levels ,
121
137
family = gaussian())
122
138
123
- # # -----------------------------------------------------------------------------
139
+
140
+ # # ----------------------------------------------------------------------
124
141
setdiff(levels(bad_levels $ series ), unique(bad_levels $ series ))
125
142
126
- # # -----------------------------------------------------------------------------
143
+
144
+ # # ----------------------------------------------------------------------
127
145
bad_levels %> %
128
146
dplyr :: mutate(series = droplevels(series )) - > good_levels
129
147
levels(good_levels $ series )
130
148
131
- # # ----error = TRUE-------------------------------------------------------------
149
+
150
+ # # ----error = TRUE------------------------------------------------------
132
151
get_mvgam_priors(outcome ~ 1 ,
133
152
data = good_levels ,
134
153
family = gaussian())
135
154
136
- # # -----------------------------------------------------------------------------
155
+
156
+ # # ----------------------------------------------------------------------
137
157
miss_dat <- data.frame (outcome = rnorm(10 ),
138
158
cov = c(NA , rnorm(9 )),
139
159
series = factor (' series1' ,
140
160
levels = ' series1' ),
141
161
time = 1 : 10 )
142
162
miss_dat
143
163
144
- # # ----error = TRUE-------------------------------------------------------------
164
+
165
+ # # ----error = TRUE------------------------------------------------------
145
166
get_mvgam_priors(outcome ~ cov ,
146
167
data = miss_dat ,
147
168
family = gaussian())
148
169
149
- # # -----------------------------------------------------------------------------
170
+
171
+ # # ----------------------------------------------------------------------
150
172
miss_dat <- list (outcome = rnorm(10 ),
151
173
series = factor (' series1' ,
152
174
levels = ' series1' ),
153
175
time = 1 : 10 )
154
176
miss_dat $ cov <- matrix (rnorm(50 ), ncol = 5 , nrow = 10 )
155
177
miss_dat $ cov [2 ,3 ] <- NA
156
178
157
- # # ----error=TRUE---------------------------------------------------------------
179
+
180
+ # # ----error=TRUE--------------------------------------------------------
158
181
get_mvgam_priors(outcome ~ cov ,
159
182
data = miss_dat ,
160
183
family = gaussian())
161
184
162
- # # ----fig.alt = "Plotting time series features for GAM models in mvgam"--------
185
+
186
+ # # ----fig.alt = "Plotting time series features for GAM models in mvgam"----
163
187
plot_mvgam_series(data = simdat $ data_train ,
164
188
y = ' y' ,
165
189
series = ' all' )
166
190
167
- # # ----fig.alt = "Plotting time series features for GAM models in mvgam"--------
191
+
192
+ # # ----fig.alt = "Plotting time series features for GAM models in mvgam"----
168
193
plot_mvgam_series(data = simdat $ data_train ,
169
194
y = ' y' ,
170
195
series = 1 )
171
196
172
- # # ----fig.alt = "Plotting time series features for GAM models in mvgam"--------
197
+
198
+ # # ----fig.alt = "Plotting time series features for GAM models in mvgam"----
173
199
plot_mvgam_series(data = simdat $ data_train ,
174
200
newdata = simdat $ data_test ,
175
201
y = ' y' ,
176
202
series = 1 )
177
203
178
- # # -----------------------------------------------------------------------------
204
+
205
+ # # ----------------------------------------------------------------------
179
206
data(" all_neon_tick_data" )
180
207
str(dplyr :: ungroup(all_neon_tick_data ))
181
208
182
- # # -----------------------------------------------------------------------------
209
+
210
+ # # ----------------------------------------------------------------------
183
211
plotIDs <- c(' SCBI_013' ,' SCBI_002' ,
184
212
' SERC_001' ,' SERC_005' ,
185
213
' SERC_006' ,' SERC_012' ,
186
214
' BLAN_012' ,' BLAN_005' )
187
215
188
- # # -----------------------------------------------------------------------------
216
+
217
+ # # ----------------------------------------------------------------------
189
218
model_dat <- all_neon_tick_data %> %
190
219
dplyr :: ungroup() %> %
191
220
dplyr :: mutate(target = ixodes_scapularis ) %> %
192
221
dplyr :: filter(plotID %in% plotIDs ) %> %
193
222
dplyr :: select(Year , epiWeek , plotID , target ) %> %
194
223
dplyr :: mutate(epiWeek = as.numeric(epiWeek ))
195
224
196
- # # -----------------------------------------------------------------------------
225
+
226
+ # # ----------------------------------------------------------------------
197
227
model_dat %> %
198
228
# Create all possible combos of plotID, Year and epiWeek;
199
229
# missing outcomes will be filled in as NA
@@ -207,7 +237,8 @@ model_dat %>%
207
237
dplyr :: select(siteID , plotID ) %> %
208
238
dplyr :: distinct()) - > model_dat
209
239
210
- # # -----------------------------------------------------------------------------
240
+
241
+ # # ----------------------------------------------------------------------
211
242
model_dat %> %
212
243
dplyr :: mutate(series = plotID ,
213
244
y = target ) %> %
@@ -216,33 +247,39 @@ model_dat %>%
216
247
dplyr :: select(- target , - plotID ) %> %
217
248
dplyr :: arrange(Year , epiWeek , series ) - > model_dat
218
249
219
- # # -----------------------------------------------------------------------------
250
+
251
+ # # ----------------------------------------------------------------------
220
252
model_dat %> %
221
253
dplyr :: ungroup() %> %
222
254
dplyr :: group_by(series ) %> %
223
255
dplyr :: arrange(Year , epiWeek ) %> %
224
256
dplyr :: mutate(time = seq(1 , dplyr :: n())) %> %
225
257
dplyr :: ungroup() - > model_dat
226
258
227
- # # -----------------------------------------------------------------------------
259
+
260
+ # # ----------------------------------------------------------------------
228
261
levels(model_dat $ series )
229
262
230
- # # ----error=TRUE---------------------------------------------------------------
263
+
264
+ # # ----error=TRUE--------------------------------------------------------
231
265
get_mvgam_priors(y ~ 1 ,
232
266
data = model_dat ,
233
267
family = poisson())
234
268
235
- # # -----------------------------------------------------------------------------
269
+
270
+ # # ----------------------------------------------------------------------
236
271
testmod <- mvgam(y ~ s(epiWeek , by = series , bs = ' cc' ) +
237
272
s(series , bs = ' re' ),
238
273
trend_model = ' AR1' ,
239
274
data = model_dat ,
240
275
backend = ' cmdstanr' ,
241
276
run_model = FALSE )
242
277
243
- # # -----------------------------------------------------------------------------
278
+
279
+ # # ----------------------------------------------------------------------
244
280
str(testmod $ model_data )
245
281
246
- # # -----------------------------------------------------------------------------
282
+
283
+ # # ----------------------------------------------------------------------
247
284
code(testmod )
248
285
0 commit comments