-
Notifications
You must be signed in to change notification settings - Fork 60
/
g.analyse.perday.R
617 lines (607 loc) · 35.1 KB
/
g.analyse.perday.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
g.analyse.perday = function(ndays, firstmidnighti, time, nfeatures,
midnightsi, metashort, averageday,
doiglevels, nfulldays,lastmidnight, ws3, ws2, qcheck,
fname, idloc, sensor.location, wdayname, tooshort, includedaycrit,
doquan, quantiletype, doilevels, domvpa,
mvpanames, wdaycode, ID,
deviceSerialNumber, ExtFunColsi, myfun, desiredtz = "",
params_247 = c(), params_phyact = c(),
...) {
#get input variables
input = list(...)
if (length(input) > 0 ||
length(params_247) == 0 || length(params_phyact) == 0) {
# Extract and check parameters if user provides more arguments than just the parameter arguments,
# or if params_[...] aren't specified (so need to be filled with defaults).
# So, inside GGIR this will not be used, but it is used when g.analyse is used on its own
# as if it was still the old g.analyse function
params = extract_params(params_247 = params_247,
params_phyact = params_phyact,
input = input,
params2check = c("247", "phyact")) # load default parameters
params_247 = params$params_247
params_phyact = params$params_phyact
rm(params)
}
startatmidnight = endatmidnight = 0
if (nfulldays >= 1) {
if (firstmidnighti == 1) { #if measurement starts at midnight
ndays = ndays - 1
startatmidnight = 1
}
if (lastmidnight == time[length(time)] & nrow(metashort) < ((60/ws3) * 1440)) { #if measurement ends at midnight
ndays = ndays - 1
endatmidnight = 1
}
}
daysummary = matrix("",ceiling(ndays),nfeatures)
ds_names = rep("",nfeatures)
windowsummary = ws_names = c()
#===============================================
# Features per day (based on on single variables)
qwindow_names = qwindowbackup = params_247[["qwindow"]]
qwindow_actlog = FALSE
if (is.data.frame(params_247[["qwindow"]]) == TRUE) {
qwindow_actlog = TRUE
}
unique_dates_recording = unique(as.Date(iso8601chartime2POSIX(time[c(seq(1, length(time),
by = (3600/ws2) * 12),
length(time))], tz = desiredtz)))
for (di in 1:ndays) { #run through days
params_247[["qwindow"]] = qwindowbackup
if (is.data.frame(params_247[["qwindow"]]) == TRUE) {
currentdate = which(params_247[["qwindow"]]$date == unique_dates_recording[di])
if (length(currentdate) > 0) {
qwindow_times = as.character(unlist(params_247[["qwindow"]]$qwindow_times[currentdate]))
qwindow_names = as.character(unlist(params_247[["qwindow"]]$qwindow_names[currentdate]))
params_247[["qwindow"]] = qwindow_values_backup = unlist(params_247[["qwindow"]]$qwindow_values[currentdate])
qwindow_names = qwindow_names[!is.na(params_247[["qwindow"]])]
qwindow_times = qwindow_times[!is.na(params_247[["qwindow"]])]
}
params_247[["qwindow"]] = params_247[["qwindow"]][!is.na(params_247[["qwindow"]])]
if (length(params_247[["qwindow"]]) == 1) params_247[["qwindow"]] = c()
if (length(params_247[["qwindow"]]) == 0 | length(currentdate) == 0) {
qwindow_times = c("00:00","24:00")
qwindow_names = c("daystart","dayend")
params_247[["qwindow"]] = qwindow_values_backup = c(0, 24)
}
} else {
if (length(params_247[["qwindow"]]) == 0) {
params_247[["qwindow"]] = c(0, 24)
}
qwindow_values_backup = params_247[["qwindow"]]
qwindow_times = as.character(params_247[["qwindow"]])
qwindow_names = as.character(params_247[["qwindow"]])
}
# extract day from matrix D and qcheck
if (startatmidnight == 1 & endatmidnight == 1) {
qqq1 = midnightsi[di] * (ws2/ws3) #a day starts at 00:00
qqq2 = (midnightsi[(di + 1)]*(ws2/ws3)) - 1
} else if (startatmidnight == 1 & endatmidnight == 0) {
if (di < floor(ndays)) { #applying floor because when day there is day saving time then ndays is not an integer
qqq1 = (midnightsi[di] - 1) * (ws2 / ws3) + 1
qqq2 = ((midnightsi[(di + 1)] - 1) * (ws2 / ws3)) #remove +1 because I got some extra seconds per day
} else if (di == floor(ndays)) {
qqq1 = (midnightsi[di] - 1) * (ws2 / ws3) + 1
qqq2 = nrow(metashort)
}
} else if (startatmidnight == 0 & endatmidnight == 0) {
if (di == 1) {
qqq1 = 1
qqq2 = ((midnightsi[di] - 1) * (ws2 / ws3))
} else if (di > 1 & di < floor(ndays)) {
qqq1 = (midnightsi[(di - 1)] - 1) * (ws2 / ws3) + 1 # a day starts at 00:00
qqq2 = ((midnightsi[di] - 1) * (ws2 / ws3))
} else if (di == floor(ndays)) {
qqq1 = (midnightsi[(di - 1)] - 1)*(ws2 / ws3) + 1 # a day starts at 00:00
qqq2 = nrow(metashort)
}
} else if (startatmidnight == 0 & endatmidnight == 1) {
if (di == 1) {
qqq1 = 1
qqq2 = (midnightsi[di] * (ws2 / ws3)) - 1
} else if (di > 1 & di <= floor(ndays)) {
qqq1 = midnightsi[(di - 1)] * (ws2 / ws3) # a day starts at 00:00
qqq2 = (midnightsi[di] * (ws2 / ws3)) - 1
}
}
if (qqq2 > nrow(metashort)) qqq2 = nrow(metashort)
vari = as.matrix(metashort[qqq1:qqq2, ])
val = qcheck[qqq1:qqq2]
nvalidhours_qwindow = rep(0, length(params_247[["qwindow"]]) - 1)
nhours_qwindow = rep(0, length(params_247[["qwindow"]]) - 1)
# Ignore qwindow values that are not possible for this day
LENVAL_hours = length(val) / (60 * (60 / ws3)) #11.2
if (length(which(round(LENVAL_hours) %in% 23:25 == TRUE)) == 0) {
if (di == 1) {
# Following 8 lines were turned off on June 5 2018 because first and last day are imputed,
# but turned on again on 14 March 2019 because when the first day is half the relative start
# of the windows must be different. Yes, the half days are imputed, but we are still only
# interested in the non-imputed part. This is probably were the confusion came from.
hours2delta = 24 - LENVAL_hours
qw_select = which(params_247[["qwindow"]] > hours2delta)
# if (qw_select[1] > 1) qw_select = c(qw_select[1] - 1, qw_select)
# params_247[["qwindow"]] = params_247[["qwindow"]][qw_select]
# # qwindow_names = qwindow_names[qw_select]
qwindowindices = params_247[["qwindow"]] - hours2delta # - LENVAL_hours # because 1 is now different
if (length(which(qwindowindices < 0)) > 0) qwindowindices[which(qwindowindices < 0)] = 0
} else if (di == ndays) {
qwindowindices = params_247[["qwindow"]]
}
} else {
hours2delta = 0
qwindowindices = params_247[["qwindow"]]
}
deltaLengthQwindow = 0
if (length(params_247[["qwindow"]]) < 2) params_247[["qwindow"]] = c()
if (length(params_247[["qwindow"]]) > 0) {
if (length(qwindow_names) == 1) {
warning("Argument to qwindow is invalid, requires a vector of at least length 2")
}
if (length(qwindow_names) == 2) {
if (params_247[["qwindow"]][1] != 0 | params_247[["qwindow"]][2] != 24) {
if((qwindowindices[2] * 60 * (60 / ws3)) <= length(val)) {
valq = val[((qwindowindices[1] * 60 * (60 / ws3)) + 1):(qwindowindices[2] * 60 * (60 / ws3))]
} else {
valq = val[((qwindowindices[1] * 60 * (60 / ws3)) + 1):length(val)]
}
nvalidhours_qwindow = length(which(valq == 0)) / (3600 / ws3)
nhours_qwindow = length(valq) / (3600 / ws3) #valid hours per day (or half a day)
}
} else if (length(qwindow_names) > 2) {
deltaLengthQwindow = length(qwindow_names) - length(qwindowindices)
for (qwi in 1:(length(qwindowindices) - 1)) {
startindex = qwindowindices[qwi] * 60 * (60/ws3)
endindex = qwindowindices[qwi + 1] * 60 * (60/ws3)
if (startindex == endindex) {
# inexistent qwindow (may occure in the first day if recording started later than first qwindow)
valq = c()
} else if (startindex <= length(val) & endindex <= length(val)) {
valq = val[(startindex + 1):endindex]
} else if (startindex <= length(val) & endindex >= length(val)) {
valq = val[(startindex + 1):length(val)]
} else if (startindex >= length(val) & endindex >= length(val)) {
valq = c()
}
if (length(valq) > 0) {
nvalidhours_qwindow[qwi + deltaLengthQwindow] = length(which(valq == 0)) / (3600 / ws3)
nhours_qwindow[qwi + deltaLengthQwindow] = length(valq) / (3600 / ws3) #valid hours per day (or half a day)
} else {
nvalidhours_qwindow[qwi + deltaLengthQwindow] = ""
nhours_qwindow[qwi + deltaLengthQwindow] = "" #valid hours per day (or half a day)
}
}
}
}
val = as.numeric(val)
nvalidhours = length(which(val == 0)) / (3600 / ws3) #valid hours per day (or half a day)
nhours = length(val) / (3600 / ws3) #valid hours per day (or half a day)
#start collecting information about the day
fi = 1
check_daysummary_size = function(daysummary, ds_names, fi) {
if (fi > ncol(daysummary) - 1000) {
expand = fi - (ncol(daysummary) - 1000)
daysummary = cbind(daysummary, matrix("", nrow(daysummary), expand + 1000))
ds_names = c(ds_names, rep("", expand + 1000))
}
invisible(list(daysummary = daysummary, ds_names = ds_names))
}
daysummary[di,fi] = ID
idremember = daysummary[di,fi]
ds_names[fi] = "ID"; fi = fi + 1
daysummary[di,fi] = fname
ds_names[fi] = "filename"; fi = fi + 1
calendardate = unlist(strsplit(as.character(vari[1,1])," "))[1]
daysummary[di,fi] = calendardate
daysummary[di,(fi + 1)] = sensor.location
daysummary[di,(fi + 2)] = nvalidhours
daysummary[di,(fi + 3)] = nhours
ds_names[fi:(fi + 3)] = c("calendar_date","bodylocation","N valid hours","N hours")
fi = fi + 4
if (length(params_247[["qwindow"]]) > 0) {
if (length(params_247[["qwindow"]]) > 2 | params_247[["qwindow"]][1] != 0 | params_247[["qwindow"]][2] != 24) {
for (qwi in 1:(length(qwindow_names) - 1)) {
tmp_name = c(paste0("N_valid_hours_", qwindow_names[qwi], "-", qwindow_names[qwi + 1], "hr"),
paste0("N_hours_", qwindow_names[qwi], "-", qwindow_names[qwi + 1], "hr"))
matchi1 = which(ds_names %in% tmp_name[1] == TRUE)
matchi2 = which(ds_names %in% tmp_name[2] == TRUE)
if (length(matchi1) == 1 & length(matchi2) == 1) {
ds_names[matchi1] == tmp_name[1]
ds_names[matchi2] == tmp_name[2]
if (length(which(qwindow_values_backup[qwi] %in% params_247[["qwindow"]] == TRUE)) > 0) {
daysummary[di,matchi1] = nvalidhours_qwindow[qwi]
daysummary[di,matchi2] = nhours_qwindow[qwi]
} else {
daysummary[di,matchi1] = "" # note that we only consider qwindows when there is data for the entire window, otherwise duration is 0
daysummary[di,matchi2] = ""
}
fi = fi + 2
} else { # variable does not exists, so add it,
ds_names = c(ds_names,tmp_name)
daysummary = cbind(daysummary, matrix("", nrow(daysummary), 2))
newncol = ncol(daysummary)
if (length(which(qwindow_values_backup[qwi] %in% params_247[["qwindow"]] == TRUE)) > 0) {
daysummary[di, (newncol - 1):newncol] = c(nvalidhours_qwindow[qwi], nhours_qwindow[qwi])
} else {
daysummary[di, (newncol - 1):newncol] = c("", "")
}
}
}
}
}
if (di > 1 & length(which(ds_names == "weekday")) > 0) { #qwindow_actlog == TRUE &
fi = which(ds_names == "weekday")
}
#--------------------------------------
weekdays = c("Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday")
weekdays = rep(weekdays, 104) # hardcoded maximum number of weeks of 104, not ideal
if (di == 1) {
daysummary[di,fi] = wdayname
} else {
daysummary[di,fi] = weekdays[wdaycode + (di - 1)]
}
daysummary[di, (fi + 1)] = di #day number relative to start of measurement
ds_names[fi:(fi + 1)] = c("weekday","measurementday")
fi = fi + 2
if (qwindow_actlog == TRUE | length(params_247[["qwindow"]]) > 2) { # add column with all timestamps in single character
ds_names[fi] = "qwindow_timestamps"
daysummary[di, fi] = paste0(qwindow_times, collapse = "_")
fi = fi + 1
ds_names[fi] = "qwindow_names"
daysummary[di,fi] = paste0(qwindow_names, collapse = "_")
fi = fi + 1
}
correct_fi = function(di, ds_names, fi, varname) {
# function to reset variable index to allign with previous day
if (di > 1) {
varexists = which(ds_names == varname)
if (length(varexists) > 0) {
fi = varexists
} else {
fi = which(ds_names == "")[1]
}
}
return(fi)
}
if (tooshort == 0) {
if (nvalidhours >= includedaycrit[1]) {
#--------------------------------------
# Features per day and per segment of the day
# guided by qwindow, which is a vector of indices to slice up the day
keepindex_46 = keepindex_48 = matrix(NA, length(2:ncol(metashort)), 2)
# generate objects to help with selecting the slices and to give names to the respective variables
anwi_t0 = 1 # analysis window time 0
anwi_t1 = nrow(as.matrix(vari)) # analysis window time 1
anwi_nameindices = "_0-24hr"
anwi_index = 1
if (length(params_247[["qwindow"]]) > 0) {
if (length(params_247[["qwindow"]]) == 2 & params_247[["qwindow"]][1] == 0 & params_247[["qwindow"]][2] == 24) {
} else {
oddqwindow = 1:(length(qwindowindices) - 1)
evenqwindow = 2:length(qwindowindices)
anwi_t0 = c(anwi_t0, ((qwindowindices[oddqwindow] * 60 * (60 / ws3)) + 1))
anwi_t1 = c(anwi_t1, (qwindowindices[evenqwindow] * 60 * (60 / ws3)))
for (iin in 1:length(oddqwindow)) {
anwi_nameindices = c(anwi_nameindices,paste0("_", qwindow_names[oddqwindow[iin]], "-",
qwindow_names[evenqwindow[iin]], "hr"))
}
}
} else {
anwi_t0 = c(anwi_t0, ((0 * 60 * (60 / ws3)) + 1))
anwi_t1 = c(anwi_t1,(24 * 60 * (60 / ws3)))
anwi_nameindices = c(anwi_nameindices, "")
}
fi_remember = fi
for (anwi_index in 1:length(anwi_t0)) {
if (anwi_index != 1 & di == 1) {
# increase value of fi to leave enough space for the variables to be calculated in second day of measurement
shift = (deltaLengthQwindow * (fi - fi_remember))
fi = fi + shift
fi_remember = fi
}
new = check_daysummary_size(daysummary, ds_names, fi)
daysummary = new$daysummary
ds_names = new$ds_names
L5M5window_name = anwi_nameindices[anwi_index]
anwindices = anwi_t0[anwi_index]:anwi_t1[anwi_index] # indices of varnum corresponding to a segment
if (length(anwindices) > 0 & all(diff(anwindices) > 0)) { # negative diff(anwindices) may occur in the first day if a qwindow is not within the recorded hours
minames = colnames(metashort)
for (mi in 2:ncol(metashort)) { #run through metrics (for features based on single metrics)
#=======================================
# Motivation on the code below:
# Standardise the number of hours in a day as an attempt to create a
# fair comparison between days in terms of day length. For example, to
# compare time spent in intensity levels or MVPA if days are not of
# equal length, such as when a recording starts in
# the middle of the day.
# To do this we impute the missing part based on the average day
# (literally called averageday in the code).
# Note: We only do this here in part 2 and not in part 5, and it has
# been this way since the early days of GGIR.
# In part 2, GGIR aims to use as much data as possible to provide
# estimates of behaviour on each recording day, even for the
# incomplete first and last recording day. As a result, it is
# important to account for imbalance in day length, which we do below.
# In part 5, however, GGIR forces the user to only work with complete
# days and by that the day length is less of a problem and not accounted for.
NRV = length(which(is.na(as.numeric(as.matrix(vari[,mi]))) == FALSE))
# Note: vari equals the imputed time series (metahsort) data from one day
varnum = as.numeric(as.matrix(vari[,mi])) # Note: varnum is one column of vari
deltaLength = NRV - length(averageday[, (mi - 1)])
if (deltaLength < 0) {
# Less than 24 hours: Append data from averageday
if (di == 1) {
# On first day of recording append the average day to the start
varnum = c(averageday[1:abs(deltaLength), (mi - 1)], varnum)
# readjust anwi indices in case that varnum has been imputed
if (max(anwi_t1) < length(varnum)) { # since GGIR always calculates full window, max(anwi_t1) should always equals length(varnum)
anwi_t0 = anwi_t0 + abs(deltaLength)
anwi_t1 = anwi_t1 + abs(deltaLength)
# then, we reset the minimum anwi_t0 to 1 to consider the imputed varnum
# anwi_t0[which(anwi_t0 == min(anwi_t0))] = 1
anwi_t0[1] = 1
}
} else {
# When it is not the first day of recording
if (NRV == 23) { # day has 23 hours (assuming DST)
# Append data after 2nd hour
startMissingHour = 2 * 60 * (60/ws3) + 1
enMissingHour = 3 * 60 * (60/ws3)
varnum = c(varnum[1:(startMissingHour - 1)], averageday[startMissingHour:enMissingHour, (mi - 1)],
varnum[startMissingHour, length(varnum)])
} else { # day has less than 24 hours for another reason
# Append the average day to the end
a56 = length(averageday[,(mi - 1)]) - abs(deltaLength) + 1
a57 = length(averageday[, (mi - 1)])
varnum = c(varnum,averageday[a56:a57, (mi - 1)])
}
}
} else if (deltaLength > 0) { # 25 hour days, assuming DST
# If there is more than 24 hours in a day then clock must
# have moved backward, remove the double hour.
startDoubleHour = 2 * 60 * (60/ws3) + 1
endDoubleHour = 3 * 60 * (60/ws3)
if (length(varnum) > endDoubleHour) {
varnum = varnum[-c(startDoubleHour:endDoubleHour)]
}
}
if (anwi_index != 1) {
if (length(anwindices) > 0) {
if (max(anwindices) > length(varnum)) {
anwindices = anwindices[which(anwindices <= length(varnum))]
}
varnum = as.numeric(varnum[anwindices]) #cut short varnum to match day segment of interest
} else {
varnum = c()
}
}
gUnitMetric = length(grep(x = colnames(metashort)[mi], pattern = "BrondCount|ZCX|ZCY|ZCZ|NeishabouriCount|ExtAct", invert = TRUE)) > 0
UnitReScale = ifelse(test = gUnitMetric, yes = 1000, no = 1)
# Starting filling output matrix daysummary with variables per day segment and full day.
if (minames[mi] %in% c("ENMO","LFENMO", "BFEN", "EN", "HFEN", "HFENplus", "MAD", "ENMOa",
"ZCX", "ZCY", "ZCZ", "BrondCount_x", "BrondCount_y",
"BrondCount_z", "NeishabouriCount_x", "NeishabouriCount_y",
"NeishabouriCount_z", "NeishabouriCount_vm", "ExtAct")) {
collectfi = c()
for (winhr_value in params_247[["winhr"]]) { # Variable (column) names
# We are first defining location of variable names, before calculating
# variables
ML5colna = c(paste0("L",winhr_value,"hr"), paste0("L",winhr_value),
paste0("M",winhr_value,"hr"), paste0("M",winhr_value))
if (length(params_247[["iglevels"]]) > 0 & length(params_247[["MX.ig.min.dur"]]) == 1) { # intensity gradient (as described by Alex Rowlands 2018)
if (winhr_value >= params_247[["MX.ig.min.dur"]]) {
for (li in 1:2) { # do twice, once for LX and once for MX
varnameig = paste0(paste0(ifelse(li == 1, yes = "L", no = "M"),winhr_value,"_"), c("ig_gradient","ig_intercept","ig_rsquared"))
ML5colna = c(ML5colna, varnameig)
}
}
}
if (length(params_247[["qM5L5"]]) > 0) {
ML5colna = c(ML5colna,
paste0("L", winhr_value, "_q", round(params_247[["qM5L5"]] * 100)),
paste0("M", winhr_value, "_q", round(params_247[["qM5L5"]] * 100)))
}
# add metric name and timewindow name
fi = correct_fi(di, ds_names, fi, varname = paste0(ML5colna[1],"_", colnames(metashort)[mi], "_mg",
L5M5window_name))
ds_names[fi:(fi - 1 + length(ML5colna))] = paste0(ML5colna, "_", colnames(metashort)[mi], "_mg",
L5M5window_name)
collectfi = c(collectfi, fi)
fi = fi + length(ML5colna)
}
if (length(varnum) > ((60 / ws3) * 60 * min(params_247[["winhr"]]) * 1.2)) { # Calculate values
for (winhr_value in params_247[["winhr"]]) {
exfi = collectfi[which(params_247[["winhr"]] == winhr_value)]
if (length(varnum) > (60 / ws3) * 60 * winhr_value * 1.2) { # Calculate values
# Time window for L5 & M5 analysis
t0_LFMF = 0 #start
# L5M5window[2] #end in 24 hour clock hours (if a value higher than 24 is chosen, it will take early hours of previous day to complete the 5 hour window
t1_LFMF = length(varnum) / (60 * (60 / ws3)) + (winhr_value - (params_247[["M5L5res"]] / 60))
ML5 = g.getM5L5(varnum, ws3, t0_LFMF, t1_LFMF, params_247[["M5L5res"]], winhr_value, qM5L5 = params_247[["qM5L5"]],
iglevels = params_247[["iglevels"]], MX.ig.min.dur = params_247[["MX.ig.min.dur"]])
ML5colna = colnames(ML5)
ML5 = as.numeric(ML5)
if (anwi_index > 1) {
L5M5shift = qwindow_values_backup[anwi_index - 1] #qwindow
} else {
L5M5shift = 0
}
if (length(ML5) > 3) {
ML5 = as.numeric(ML5)
ML5[c(1, 3)] = ML5[c(1, 3)] + L5M5shift
daysummary[di, (exfi):(exfi + length(ML5) - 1)] = ML5
} else {
daysummary[di, (exfi):(exfi + length(ML5) - 1)] = ""
}
} else {
daysummary[di,(exfi):(exfi + (4 * length(params_247[["winhr"]])) - 1 + (length(params_247[["qM5L5"]]) * 2))] = ""
}
# (Below) 4 is the length of ML5 and then 2 extra variables for every qM5L5 value
# winhr is not considered because we are in the winhr loop:
exfi = exfi + 4 + (length(params_247[["qM5L5"]]) * 2)
}
} else {
daysummary[di,collectfi] = ""
}
if (anwindices[1] == 1 & length(anwindices) > 6*60*(60/ws3)) { # only derive if 1-6am falls within window
fi = correct_fi(di, ds_names, fi, varname = paste0("mean_", colnames(metashort)[mi], "_mg_1-6am"))
daysummary[di,fi] = mean(varnum[((1 * 60 * (60 / ws3)) + 1):(6 * 60 * (60 / ws3))]) * UnitReScale #from 1am to 6am
ds_names[fi] = paste0("mean_",colnames(metashort)[mi],"_mg_1-6am"); fi = fi + 1
}
if (anwi_nameindices[anwi_index] == "") { # for consistency with previous GGIR version
anwi_nameindices[anwi_index] = "_24hr"
}
cn_metashort = colnames(metashort)
fi = correct_fi(di, ds_names, fi, varname = paste0("mean_", cn_metashort[mi], "_mg",
anwi_nameindices[anwi_index]))
if (length(varnum) > 0) {
daysummary[di,fi] = mean(varnum) * UnitReScale
} else {
daysummary[di,fi] = ""
}
ds_names[fi] = paste0("mean_", cn_metashort[mi], "_mg", anwi_nameindices[anwi_index]); fi = fi + 1
if (anwi_nameindices[anwi_index] == "_24hr") {
anwi_nameindices[anwi_index] = ""
}
if (doquan == TRUE) {
q46 = c()
q46 = quantile(varnum, probs = params_247[["qlevels"]], na.rm = T, type = quantiletype) * UnitReScale
keepindex_46[mi - 1,] = c(fi,(fi + (length(params_247[["qlevels"]]) - 1)))
namesq46 = rep(0,length(rownames(as.matrix(q46))))
for (rq46i in 1:length(rownames(as.matrix(q46)))) {
tmp1 = rownames(as.matrix(q46))[rq46i]
tmp2 = as.character(unlist(strsplit(tmp1, "%")))
namesq46[rq46i] = paste0("p", tmp2, "_", colnames(metashort)[mi], "_mg",
anwi_nameindices[anwi_index])
}
fi = correct_fi(di, ds_names, fi, varname = namesq46[1])
if (length(varnum) > 0) {
daysummary[di, fi:(fi + (length(params_247[["qlevels"]]) - 1))] = q46
} else {
daysummary[di, fi:(fi + (length(params_247[["qlevels"]]) - 1))] = ""
}
ds_names[fi:(fi + (length(params_247[["qlevels"]]) - 1))] = namesq46
fi = fi + length(params_247[["qlevels"]])
}
if (doilevels == TRUE) {
q48 = c()
#times 1000 to convert to mg only if it g-unit metric
q47 = cut((varnum * UnitReScale), breaks = params_247[["ilevels"]], right = FALSE)
q47 = table(q47)
q48 = (as.numeric(q47) * ws3) / 60 #converting to minutes
keepindex_48[mi - 1,] = c(fi, (fi + (length(q48) - 1)))
namesq47 = rep(0, length(rownames(q47)))
for (rq47i in 1:length(rownames(q47))) {
namesq47[rq47i] = paste0(rownames(q47)[rq47i], "_", colnames(metashort)[mi], "_mg",
anwi_nameindices[anwi_index])
}
fi = correct_fi(di, ds_names, fi, varname = namesq47[1])
if (length(varnum) > 0) {
daysummary[di,fi:(fi + (length(q48) - 1))] = q48
} else {
daysummary[di, fi:(fi + (length(q48) - 1))] = ""
}
ds_names[fi:(fi + (length(q48) - 1))] = namesq47
fi = fi + length(q48)
}
if (doiglevels == TRUE) { # intensity gradient (as described by Alex Rowlands 2018)
q49 = c()
q50 = cut((varnum * UnitReScale), breaks = params_247[["iglevels"]], right = FALSE)
q50 = table(q50)
q49 = (as.numeric(q50) * ws3) / 60 #converting to minutes
x_ig = zoo::rollmean(params_247[["iglevels"]], k = 2)
y_ig = q49
igout = g.intensitygradient(x_ig, y_ig)
varnameig = paste0(c("ig_gradient", "ig_intercept", "ig_rsquared"),
paste0("_", colnames(metashort)[mi], anwi_nameindices[anwi_index]))
fi = correct_fi(di, ds_names, fi, varname = varnameig[1])
if (length(varnum) > 0) {
daysummary[di, fi:(fi + 2)] = as.vector(unlist(igout))
} else {
daysummary[di, fi:(fi + 2)] = rep("", 3)
}
ds_names[fi:(fi + 2)] = varnameig
fi = fi + 3
}
#=========================================
if (domvpa == TRUE) {
for (mvpai in 1:length(params_phyact[["mvpathreshold"]])) {
mvpa = rep(0,6)
if (length(varnum) < 100) {
mvpa[1:6] = 0
} else {
# METHOD 1: time spent above threhold based on 5 sec epoch
mvpa[1] = length(which(varnum * UnitReScale >= params_phyact[["mvpathreshold"]][mvpai])) / (60/ws3) #time spent MVPA in minutes
# METHOD 2: time spent above threshold based on 1minute epoch
varnum2 = cumsum(c(0, varnum))
select = seq(1, length(varnum2), by = 60/ws3)
varnum3 = diff(varnum2[round(select)]) / abs(diff(round(select)))
mvpa[2] = length(which(varnum3 * UnitReScale >= params_phyact[["mvpathreshold"]][mvpai])) #time spent MVPA in minutes
# METHOD 3: time spent above threshold based on 5minute epoch
select = seq(1, length(varnum2), by = 300/ws3)
varnum3 = diff(varnum2[round(select)]) / abs(diff(round(select)))
mvpa[3] = length(which(varnum3 * UnitReScale >= params_phyact[["mvpathreshold"]][mvpai])) * 5 #time spent MVPA in minutes
# METHOD 4: time spent above threshold
boutduration = params_phyact[["mvpadur"]][1] * (60/ws3) # per minute
rr1 = matrix(0, length(varnum), 1)
p = which(varnum * UnitReScale >= params_phyact[["mvpathreshold"]][mvpai]); rr1[p] = 1
getboutout = g.getbout(x = rr1, boutduration = boutduration,
boutcriter = params_phyact[["boutcriter"]], ws3 = ws3)
mvpa[4] = length(which(getboutout == 1)) / (60/ws3) #time spent MVPA in minutes
# METHOD 5: time spent above threshold 5 minutes
boutduration = params_phyact[["mvpadur"]][2] * (60/ws3) #per five minutes
rr1 = matrix(0, length(varnum), 1)
p = which(varnum * UnitReScale >= params_phyact[["mvpathreshold"]][mvpai]); rr1[p] = 1
getboutout = g.getbout(x = rr1, boutduration = boutduration,
boutcriter = params_phyact[["boutcriter"]], ws3 = ws3)
mvpa[5] = length(which(getboutout == 1)) / (60/ws3) #time spent MVPA in minutes
# METHOD 6: time spent above threshold 10 minutes
boutduration = params_phyact[["mvpadur"]][3] * (60/ws3) # per ten minutes
rr1 = matrix(0,length(varnum),1)
p = which(varnum * UnitReScale >= params_phyact[["mvpathreshold"]][mvpai]); rr1[p] = 1
getboutout = g.getbout(x = rr1, boutduration = boutduration,
boutcriter = params_phyact[["boutcriter"]], ws3 = ws3)
mvpa[6] = length(which(getboutout == 1)) / (60/ws3) #time spent MVPA in minutes
}
if (length(which(is.nan(mvpa) == TRUE)) > 0) mvpa[which(is.nan(mvpa) == TRUE)] = 0
mvpanames[,mvpai] = c( paste0("MVPA_E", ws3, "S_T", params_phyact[["mvpathreshold"]][mvpai]),
paste0("MVPA_E1M_T", params_phyact[["mvpathreshold"]][mvpai]),
paste0("MVPA_E5M_T", params_phyact[["mvpathreshold"]][mvpai]),
paste0("MVPA_E", ws3, "S_B", params_phyact[["mvpadur"]][1], "M", (params_phyact[["boutcriter"]] * 100), "%_T", params_phyact[["mvpathreshold"]][mvpai]),
paste0("MVPA_E", ws3, "S_B", params_phyact[["mvpadur"]][2], "M", (params_phyact[["boutcriter"]] * 100), "%_T", params_phyact[["mvpathreshold"]][mvpai]),
paste0("MVPA_E", ws3, "S_B", params_phyact[["mvpadur"]][3], "M", (params_phyact[["boutcriter"]] * 100), "%_T", params_phyact[["mvpathreshold"]][mvpai]))
for (fillds in 1:6) {
mvpavarname = paste0(mvpanames[fillds,mvpai], "_", colnames(metashort)[mi], anwi_nameindices[anwi_index])
fi = correct_fi(di, ds_names, fi, varname = mvpavarname)
daysummary[di,fi] = mvpa[fillds]
ds_names[fi] = mvpavarname; fi = fi + 1
}
}
}
}
if (mi %in% ExtFunColsi == TRUE) { # INSERT HERE VARIABLES DERIVED WITH EXTERNAL FUNCTION
rti = which(ExtFunColsi == mi)
if (myfun$reporttype[rti] == "event") { # For the event report type we take the sum
varnameevent = paste0(colnames(metashort)[mi], "_sum", anwi_nameindices[anwi_index])
fi = correct_fi(di, ds_names, fi, varname = varnameevent)
daysummary[di,fi] = sum(varnum)
ds_names[fi] = varnameevent; fi = fi + 1
} else if (myfun$reporttype[rti] == "scalar") { # For the scalar report type we take the mean
varnamescalar = paste0(colnames(metashort)[mi], "_mean", anwi_nameindices[anwi_index])
fi = correct_fi(di, ds_names, fi, varname = varnamescalar)
daysummary[di,fi] = mean(varnum)
ds_names[fi] = varnamescalar; fi = fi + 1
} else if (myfun$reporttype[rti] == "type") { # For type we calculate time spent in each class
# Not implemented yet
}
}
}
}
}
}
rm(val); rm(vari)
}
}
invisible(list(daysummary = daysummary, ds_names = ds_names,
windowsummary = windowsummary, ws_names = ws_names))
}