annotate dennis_gam_initial_functions.R @ 2:0e7066603eab draft default tip

Uploaded
author mnhn65mo
date Mon, 06 Aug 2018 04:45:49 -0400
parents 5b126f770671
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
1 ### R-Script Adapted from script provided by the CEH, UK BY: Reto Schmucki [ reto.schmucki@mail.mcgill.ca]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
2 ### DATE: 14 July 2014 function to run two stage model in DENNIS et al. 2013
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
3
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
4
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
5 .onAttach <- function(libname, pkgname) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
6 packageStartupMessage(" The regionalGAM package that is no longer maintained, \n use the new rbms package instead. \n
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
7 devtools::install_github(\"RetoSchmucki/rbms\", force=TRUE)")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
8 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
9
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
10
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
11 #' year_day_func Function
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
12 #' This function generate the full sequence of days, months and include the observation to that file.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
13 #' @param sp_data A data.frame with your observation.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
14 #' @keywords year days
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
15 #' @export
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
16 #' @author Reto Schmucki
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
17 #' @examples
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
18 #' year_day_func()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
19
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
20
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
21 # FUNCTIONS
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
22
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
23 year_day_func = function(sp_data) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
24
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
25 year <- unique(sp_data$YEAR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
26
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
27 origin.d <- paste(year, "01-01", sep = "-")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
28 if ((year%%4 == 0) & ((year%%100 != 0) | (year%%400 == 0))) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
29 nday <- 366
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
30 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
31 nday <- 365
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
32 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
33
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
34 date.serie <- as.POSIXlt(seq(as.Date(origin.d), length = nday, by = "day"), format = "%Y-%m-%d")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
35
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
36 dayno <- as.numeric(julian(date.serie, origin = as.Date(origin.d)) + 1)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
37 month <- as.numeric(strftime(date.serie, format = "%m"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
38 week <- as.numeric(strftime(date.serie, format = "%W"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
39 week_day <- as.numeric(strftime(date.serie, format = "%u"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
40 day <- as.numeric(strftime(date.serie, format = "%d"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
41
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
42 site_list <- sp_data[!duplicated(sp_data$SITE), c("SITE")]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
43
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
44 all_day_site <- data.frame(SPECIES = sp_data$SPECIES[1], SITE = rep(site_list, rep(nday, length(site_list))),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
45 YEAR = sp_data$YEAR[1], MONTH = month, WEEK = week, DAY = day, DAY_WEEK = week_day, DAYNO = dayno,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
46 COUNT = NA)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
47
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
48 count_index <- match(paste(sp_data$SITE, sp_data$DAYNO, sep = "_"), paste(all_day_site$SITE, all_day_site$DAYNO,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
49 sep = "_"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
50 all_day_site$COUNT[count_index] <- sp_data$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
51 site_count_length <- aggregate(sp_data$COUNT, by = list(sp_data$SITE), function(x) list(1:length(x)))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
52 names(site_count_length$x) <- as.character(site_count_length$Group.1)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
53 site_countno <- utils::stack(site_count_length$x)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
54 all_day_site$COUNTNO <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
55 all_day_site$COUNTNO[count_index] <- site_countno$values # add count number to ease extraction of single count
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
56
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
57 # Add zero to close observation season two weeks before and after the first and last
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
58 first_obs <- min(all_day_site$DAYNO[!is.na(all_day_site$COUNT)])
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
59 last_obs <- max(all_day_site$DAYNO[!is.na(all_day_site$COUNT)])
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
60
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
61 closing_season <- c((first_obs - 11):(first_obs - 7), (last_obs + 7):(last_obs + 11))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
62
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
63 # If closing season is before day 1 or day 365, simply set the first and last 5 days to 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
64 if (min(closing_season) < 1)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
65 closing_season[1:5] <- c(1:5)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
66 if (max(closing_season) > nday)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
67 closing_season[6:10] <- c((nday - 4):nday)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
68
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
69 all_day_site$COUNT[all_day_site$DAYNO %in% closing_season] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
70 all_day_site$ANCHOR <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
71 all_day_site$ANCHOR[all_day_site$DAYNO %in% closing_season] <- 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
72
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
73 all_day_site <- all_day_site[order(all_day_site$SITE, all_day_site$DAYNO), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
74
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
75 return(all_day_site)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
76 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
77
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
78
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
79 #' trap_area Function
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
80 #'
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
81 #' This function compute the area under the curve using the trapezoid method.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
82 #' @param x A vector or a two columns matrix.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
83 #' @param y A vector, Default is NULL
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
84 #' @keywords trapezoid
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
85 #' @export
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
86 #' @examples
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
87 #' trap_area()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
88
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
89
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
90 trap_area = function(x, y = NULL) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
91 # If y is null and x has multiple columns then set y to x[,2] and x to x[,1]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
92 if (is.null(y)) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
93 if (length(dim(x)) == 2) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
94 y = x[, 2]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
95 x = x[, 1]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
96 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
97 stop("ERROR: need to either specifiy both x and y or supply a two column data.frame/matrix to x")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
98 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
99 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
100
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
101 # Check x and y are same length
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
102 if (length(x) != length(y)) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
103 stop("ERROR: x and y need to be the same length")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
104 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
105
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
106 # Need to exclude any pairs that are NA for either x or y
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
107 rm_inds = which(is.na(x) | is.na(y))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
108 if (length(rm_inds) > 0) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
109 x = x[-rm_inds]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
110 y = y[-rm_inds]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
111 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
112
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
113 # Determine values of trapezoids under curve Get inds
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
114 inds = 1:(length(x) - 1)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
115 # Determine area using trapezoidal rule Area = ( (b1 + b2)/2 ) * h where b1 and b2 are lengths of bases
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
116 # (the parallel sides) and h is the height (the perpendicular distance between two bases)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
117 areas = ((y[inds] + y[inds + 1])/2) * diff(x)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
118
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
119 # total area is sum of all trapezoid areas
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
120 tot_area = sum(areas)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
121
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
122 # Return total area
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
123 return(tot_area)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
124 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
125
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
126
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
127 #' trap_index Function
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
128 #'
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
129 #' This function compute the area under the curve (Abundance Index) across species, sites and years
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
130 #' @param sp_data A data.frame containing species count data generated from the year_day_func()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
131 #' @param y A vector, Default is NULL
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
132 #' @keywords Abundance index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
133 #' @export
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
134 #' @examples
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
135 #' trap_index()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
136
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
137
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
138
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
139 trap_index = function(sp_data, data_col = "IMP", time_col = "DAYNO", by_col = c("SPECIES", "SITE", "YEAR")) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
140
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
141 # Build output data.frame
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
142 out_obj = unique(sp_data[, by_col])
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
143 # Set row.names to be equal to collapsing of output rows (will be unique, you need them to make uploading
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
144 # values back to data.frame will be easier)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
145 row.names(out_obj) = apply(out_obj, 1, paste, collapse = "_")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
146
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
147 # Using this row.names from out_obj above as index in by function to loop through values all unique combs
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
148 # of by_cols and fit trap_area to data
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
149 ind_dat = by(sp_data[, c(time_col, data_col)], apply(sp_data[, by_col], 1, paste, collapse = "_"), trap_area)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
150
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
151 # Add this data to output object
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
152 out_obj[names(ind_dat), "SINDEX"] = round(ind_dat/7, 1)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
153
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
154 # Set row.names to defaults
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
155 row.names(out_obj) = NULL
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
156
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
157 # Return output object
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
158 return(out_obj)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
159 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
160
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
161
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
162 #' flight_curve Function
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
163 #' This function compute the flight curve across and years
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
164 #' @param your_dataset A data.frame containing original species count data. The data format is a csv (comma "," separated) with 6 columns with headers, namely "species","transect_id","visit_year","visit_month","visit_day","sp_count"
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
165 #' @keywords standardize flight curve
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
166 #' @export
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
167 #' @examples
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
168 #' flight_curve()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
169
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
170
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
171 flight_curve <- function(your_dataset) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
172
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
173 if("mgcv" %in% installed.packages() == "FALSE") {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
174 print("mgcv package is not installed.")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
175 x <- readline("Do you want to install it? Y/N")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
176 if (x == 'Y') {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
177 install.packages("mgcv")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
178 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
179 if (x == 'N') {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
180 stop("flight curve can not be computed without the mgcv package, sorry")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
181 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
182 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
183 your_dataset$DAYNO <- strptime(paste(your_dataset$DAY, your_dataset$MONTH,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
184 your_dataset$YEAR, sep = "/"), "%d/%m/%Y")$yday + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
185 dataset <- your_dataset[, c("SPECIES", "SITE", "YEAR", "MONTH",
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
186 "DAY", "DAYNO", "COUNT")]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
187 sample_year <- unique(dataset$YEAR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
188 sample_year <- sample_year[order(sample_year)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
189 if (length(sample_year) >1 ) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
190 for (y in sample_year) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
191 dataset_y <- dataset[dataset$YEAR == y, ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
192 nsite <- length(unique(dataset_y$SITE))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
193 # Determine missing days and add to dataset
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
194 sp_data_all <- year_day_func(dataset_y)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
195 if (nsite > 200) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
196 sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
197 200, replace = F)]), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
198 sp_data_all <- sp_data_all
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
199 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
200 sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
201 print(paste("Fitting the GAM for",as.character(sp_data_all$SPECIES[1]),"and year",y,"with",length(unique(sp_data_all$SITE)),"sites :",Sys.time()))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
202 if(length(unique(sp_data_all$SITE))>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
203 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
204 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
205 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
206 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
207 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
208 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
209 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
210 # Give a second try if the GAM does not converge the first time
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
211 if (class(gam_obj_site)[1] == "try-error") {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
212 # Determine missing days and add to dataset
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
213 sp_data_all <- year_day_func(dataset_y)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
214 if (nsite > 200) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
215 sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
216 200, replace = F)]), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
217 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
218 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
219 sp_data_all <- sp_data_all
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
220 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
221 sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
222 print(paste("Fitting the GAM for",sp_data_all$SPECIES[1],"at year", y,"with",length(unique(sp_data_all$SITE)),"sites :",Sys.time(),"second try"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
223 if(length(unique(sp_data_all$SITE))>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
224 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
225 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
226 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
227 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
228 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
229 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
230 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
231 if (class(gam_obj_site)[1] == "try-error") {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
232 print(paste("Error in fitting the flight period for",sp_data_all$SPECIES[1],"at year", y,"no convergence after two trial"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
233 sp_data_all[, "FITTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
234 sp_data_all[, "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
235 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
236 sp_data_all[, "NM"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
237 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
238 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
239 # Generate a list of values for all days from the additive model and use
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
240 # these value to fill the missing observations
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
241 sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
242 c("trimDAYNO", "SITE")], type = "response")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
243 # force zeros at the beginning end end of the year
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
244 sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
245 sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
246 # if infinite number in predict replace with NA.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
247 if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
248 print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
249 sp_data_all[, "FITTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
250 sp_data_all[, "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
251 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
252 sp_data_all[, "NM"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
253 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
254 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
255 sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
256 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
257 # Define the flight curve from the fitted values and append them over
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
258 # years (this is one flight curve per year for all site)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
259 site_sums <- aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
260 FUN = sum)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
261 # Rename sum column
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
262 names(site_sums)[names(site_sums) == "x"] <- "SITE_YR_FSUM"
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
263 # Add data to sp_data data.frame (ensure merge does not sort the data!)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
264 sp_data_all = merge(sp_data_all, site_sums, by <- c("SITE"),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
265 all = TRUE, sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
266 # Calculate normalized values
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
267 sp_data_all[, "NM"] <- sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
268 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
269 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
270 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
271 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
272 # Generate a list of values for all days from the additive model and use
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
273 # these value to fill the missing observations
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
274 sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
275 c("trimDAYNO", "SITE")], type = "response")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
276 # force zeros at the beginning end end of the year
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
277 sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
278 sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
279 # if infinite number in predict replace with NA.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
280 if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
281 print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
282 sp_data_all[, "FITTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
283 sp_data_all[, "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
284 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
285 sp_data_all[, "NM"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
286 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
287 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
288 sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
289 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
290 # Define the flight curve from the fitted values and append them over
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
291 # years (this is one flight curve per year for all site)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
292 site_sums = aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
293 FUN = sum)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
294 # Rename sum column
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
295 names(site_sums)[names(site_sums) == "x"] = "SITE_YR_FSUM"
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
296 # Add data to sp_data data.frame (ensure merge does not sort the data!)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
297 sp_data_all = merge(sp_data_all, site_sums, by = c("SITE"), all = TRUE,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
298 sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
299 # Calculate normalized values
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
300 sp_data_all[, "NM"] = sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
301 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
302 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
303 sp_data_filled <- sp_data_all
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
304 flight_curve <- data.frame(species = sp_data_filled$SPECIES, year = sp_data_filled$YEAR,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
305 week = sp_data_filled$WEEK, DAYNO = sp_data_filled$DAYNO, DAYNO_adj = sp_data_filled$trimDAYNO,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
306 nm = sp_data_filled$NM)[!duplicated(paste(sp_data_filled$YEAR,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
307 sp_data_filled$DAYNO, sep = "_")), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
308 flight_curve <- flight_curve[order(flight_curve$DAYNO), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
309 # bind if exist else create
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
310 if (is.na(flight_curve$nm[1])) next()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
311 if ("flight_pheno" %in% ls()) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
312 flight_pheno <- rbind(flight_pheno, flight_curve)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
313 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
314 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
315 flight_pheno <- flight_curve
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
316 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
317 } # end of year loop
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
318 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
319 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
320 y <- unique(dataset$YEAR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
321 dataset_y <- dataset[dataset$YEAR == y, ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
322 nsite <- length(unique(dataset_y$SITE))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
323 # Determine missing days and add to dataset
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
324 sp_data_all <- year_day_func(dataset_y)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
325 if (nsite > 200) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
326 sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
327 200, replace = F)]), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
328 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
329 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
330 sp_data_all <- sp_data_all
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
331 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
332 sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
333 print(paste("Fitting the GAM for",sp_data_all$SPECIES[1],"at year", y,":",Sys.time()))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
334 if(length(unique(sp_data_all$SITE))>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
335 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
336 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
337 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
338 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
339 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
340 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
341 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
342 # Give a second try if the GAM does not converge the first time
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
343 if (class(gam_obj_site)[1] == "try-error") {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
344 # Determine missing days and add to dataset
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
345 sp_data_all <- year_day_func(dataset_y)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
346 if (nsite > 200) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
347 sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
348 200, replace = F)]), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
349 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
350 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
351 sp_data_all <- sp_data_all
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
352 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
353 sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
354 print(paste("Fitting the GAM for",sp_data_all$SPECIES[1],"at year", y,"with",length(unique(sp_data_all$SITE)),"sites :",Sys.time(),"second try"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
355 if(length(unique(sp_data_all$SITE))>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
356 gam_obj_site <- try(mgcv::bam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) - 1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
357 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
358 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
359 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
360 gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") -1,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
361 data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
362 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
363 if (class(gam_obj_site)[1] == "try-error") {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
364 print(paste("Error in fitting the flight period for",sp_data_all$SPECIES[1],"at year", y,"no convergence after two trial"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
365 sp_data_all[, "FITTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
366 sp_data_all[, "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
367 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
368 sp_data_all[, "NM"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
369 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
370 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
371 # Generate a list of values for all days from the additive model and use
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
372 # these value to fill the missing observations
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
373 sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
374 c("trimDAYNO", "SITE")], type = "response")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
375 # force zeros at the beginning end end of the year
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
376 sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
377 sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
378 # if infinite number in predict replace with NA.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
379 if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
380 print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
381 sp_data_all[, "FITTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
382 sp_data_all[, "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
383 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
384 sp_data_all[, "NM"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
385 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
386 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
387 sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
388 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
389 # Define the flight curve from the fitted values and append them over
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
390 # years (this is one flight curve per year for all site)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
391 site_sums <- aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
392 FUN = sum)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
393 # Rename sum column
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
394 names(site_sums)[names(site_sums) == "x"] <- "SITE_YR_FSUM"
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
395 # Add data to sp_data data.frame (ensure merge does not sort the data!)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
396 sp_data_all = merge(sp_data_all, site_sums, by <- c("SITE"),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
397 all = TRUE, sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
398 # Calculate normalized values
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
399 sp_data_all[, "NM"] <- sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
400 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
401 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
402 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
403 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
404 # Generate a list of values for all days from the additive model and use
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
405 # these value to fill the missing observations
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
406 sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
407 c("trimDAYNO", "SITE")], type = "response")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
408 # force zeros at the beginning end end of the year
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
409 sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
410 sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
411 # if infinite number in predict replace with NA.
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
412 if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
413 print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
414 sp_data_all[, "FITTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
415 sp_data_all[, "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
416 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
417 sp_data_all[, "NM"] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
418 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
419 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
420 sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
421 sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
422 # Define the flight curve from the fitted values and append them over
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
423 # years (this is one flight curve per year for all site)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
424 site_sums = aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
425 FUN = sum)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
426 # Rename sum column
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
427 names(site_sums)[names(site_sums) == "x"] = "SITE_YR_FSUM"
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
428 # Add data to sp_data data.frame (ensure merge does not sort the data!)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
429 sp_data_all = merge(sp_data_all, site_sums, by = c("SITE"), all = TRUE,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
430 sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
431 # Calculate normalized values
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
432 sp_data_all[, "NM"] = sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
433 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
434 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
435 sp_data_filled <- sp_data_all
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
436 flight_curve <- data.frame(species = sp_data_filled$SPECIES, year = sp_data_filled$YEAR,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
437 week = sp_data_filled$WEEK, DAYNO = sp_data_filled$DAYNO, DAYNO_adj = sp_data_filled$trimDAYNO,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
438 nm = sp_data_filled$NM)[!duplicated(paste(sp_data_filled$YEAR,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
439 sp_data_filled$DAYNO, sep = "_")), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
440 flight_curve <- flight_curve[order(flight_curve$DAYNO), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
441 if (is.na(flight_curve$nm[1])){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
442 flight_pheno <- data.frame()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
443 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
444 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
445 # bind if exist else create
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
446 if ("flight_pheno" %in% ls()) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
447 flight_pheno <- rbind(flight_pheno, flight_curve)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
448 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
449 else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
450 flight_pheno <- flight_curve
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
451 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
452 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
453 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
454 return(flight_pheno)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
455 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
456
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
457
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
458 #' abundance_index Function
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
459 #'
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
460 #' This function compute the Abundance Index across sites and years from your dataset and the regional flight curve
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
461 #' @param your_dataset A data.frame containing original species count data. The data format is a csv (comma "," separated) with 6 columns with headers, namely "species","transect_id","visit_year","visit_month","visit_day","sp_count"
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
462 #' @param flight_pheno A data.frame for the regional flight curve computed with the function flight_curve
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
463 #' @keywords standardize flight curve
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
464 #' @export
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
465 #' @examples
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
466 #' abundance_index()
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
467
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
468 abundance_index <- function(your_dataset,flight_pheno) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
469
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
470 your_dataset$DAYNO <- strptime(paste(your_dataset$DAY, your_dataset$MONTH,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
471 your_dataset$YEAR, sep = "/"), "%d/%m/%Y")$yday + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
472
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
473 dataset <- your_dataset[, c("SPECIES", "SITE", "YEAR", "MONTH",
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
474 "DAY", "DAYNO", "COUNT")]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
475
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
476 sample_year <- unique(dataset$YEAR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
477 sample_year <- sample_year[order(sample_year)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
478
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
479
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
480 if (length(sample_year)>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
481
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
482 for (y in sample_year) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
483
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
484 year_pheno <- flight_pheno[flight_pheno$year == y, ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
485
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
486 dataset_y <- dataset[dataset$YEAR == y, ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
487
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
488 sp_data_site <- year_day_func(dataset_y)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
489 sp_data_site$trimDAYNO <- sp_data_site$DAYNO - min(sp_data_site$DAYNO) + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
490
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
491 sp_data_site <- merge(sp_data_site, year_pheno[, c("DAYNO", "nm")],
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
492 by = c("DAYNO"), all.x = TRUE, sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
493
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
494 # compute proportion of the flight curve sampled due to missing visits
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
495 pro_missing_count <- data.frame(SITE = sp_data_site$SITE, WEEK = sp_data_site$WEEK,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
496 NM = sp_data_site$nm, COUNT = sp_data_site$COUNT, ANCHOR = sp_data_site$ANCHOR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
497 pro_missing_count$site_week <- paste(pro_missing_count$SITE, pro_missing_count$WEEK,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
498 sep = "_")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
499 siteweeknocount <- aggregate(pro_missing_count$COUNT, by = list(pro_missing_count$site_week),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
500 function(x) sum(!is.na(x)) == 0)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
501 pro_missing_count <- pro_missing_count[pro_missing_count$site_week %in%
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
502 siteweeknocount$Group.1[siteweeknocount$x == TRUE], ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
503 pro_count_agg <- aggregate(pro_missing_count$NM, by = list(pro_missing_count$SITE),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
504 function(x) 1 - sum(x, na.rm = T))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
505 names(pro_count_agg) <- c("SITE", "PROP_PHENO_SAMPLED")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
506
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
507 # remove samples outside the monitoring window
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
508 sp_data_site$COUNT[sp_data_site$nm==0] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
509
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
510 # Compute the regional GAM index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
511
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
512 if(length(unique(sp_data_site$SITE))>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
513 glm_obj_site <- glm(COUNT ~ factor(SITE) + offset(log(nm)) - 1, data = sp_data_site,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
514 family = quasipoisson(link = "log"), control = list(maxit = 100))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
515 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
516 glm_obj_site <- glm(COUNT ~ offset(log(nm)) - 1, data = sp_data_site,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
517 family = quasipoisson(link = "log"), control = list(maxit = 100))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
518 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
519
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
520 sp_data_site[, "FITTED"] <- predict.glm(glm_obj_site, newdata = sp_data_site,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
521 type = "response")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
522 sp_data_site[, "COUNT_IMPUTED"] <- sp_data_site$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
523 sp_data_site[is.na(sp_data_site$COUNT), "COUNT_IMPUTED"] <- sp_data_site$FITTED[is.na(sp_data_site$COUNT)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
524
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
525 ## add fitted value for missing mid-week data
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
526 sp_data_site <- sp_data_site[!paste(sp_data_site$DAY_WEEK, sp_data_site$COUNT) %in%
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
527 c("1 NA", "2 NA", "3 NA", "5 NA", "6 NA", "7 NA"), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
528
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
529 ## remove all added mid-week values for weeks with real count
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
530 ## (observation)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
531 sp_data_site$site_week <- paste(sp_data_site$SITE, sp_data_site$WEEK,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
532 sep = "_")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
533 siteweekcount <- aggregate(sp_data_site$COUNT, by = list(sp_data_site$site_week),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
534 function(x) sum(!is.na(x)) > 0)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
535 sp_data_site <- sp_data_site[!(is.na(sp_data_site$COUNT) & (sp_data_site$site_week %in%
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
536 siteweekcount$Group.1[siteweekcount$x == TRUE])), names(sp_data_site) !=
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
537 "site_week"]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
538
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
539 ## Compute the regional GAM index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
540 print(paste("Compute index for",sp_data_site$SPECIES[1],"at year", y,"for",length(unique(sp_data_site$SITE)),"sites:",Sys.time()))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
541 regional_gam_index <- trap_index(sp_data_site, data_col = "COUNT_IMPUTED",
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
542 time_col = "DAYNO", by_col = c("SPECIES", "SITE", "YEAR"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
543
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
544 cumu_index <- merge(regional_gam_index, pro_count_agg, by = c("SITE"),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
545 all.x = TRUE, sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
546 names(cumu_index) <- c("SITE", "SPECIES", "YEAR", "regional_gam", "prop_pheno_sampled")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
547
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
548 cumu_index <- cumu_index[order(cumu_index$SITE), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
549
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
550 # bind if exist else create
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
551 if ("cumullated_indices" %in% ls()) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
552 cumullated_indices <- rbind(cumullated_indices, cumu_index)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
553 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
554 cumullated_indices <- cumu_index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
555 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
556
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
557 } # end of year loop
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
558
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
559 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
560
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
561 y <- unique(dataset$YEAR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
562 year_pheno <- flight_pheno[flight_pheno$year == y, ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
563
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
564 dataset_y <- dataset[dataset$YEAR == y, ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
565
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
566 sp_data_site <- year_day_func(dataset_y)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
567 sp_data_site$trimDAYNO <- sp_data_site$DAYNO - min(sp_data_site$DAYNO) + 1
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
568
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
569 sp_data_site <- merge(sp_data_site, year_pheno[, c("DAYNO", "nm")],
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
570 by = c("DAYNO"), all.x = TRUE, sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
571
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
572 # compute proportion of the flight curve sampled due to missing visits
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
573 pro_missing_count <- data.frame(SITE = sp_data_site$SITE, WEEK = sp_data_site$WEEK,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
574 NM = sp_data_site$nm, COUNT = sp_data_site$COUNT, ANCHOR = sp_data_site$ANCHOR)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
575 pro_missing_count$site_week <- paste(pro_missing_count$SITE, pro_missing_count$WEEK,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
576 sep = "_")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
577 siteweeknocount <- aggregate(pro_missing_count$COUNT, by = list(pro_missing_count$site_week),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
578 function(x) sum(!is.na(x)) == 0)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
579 pro_missing_count <- pro_missing_count[pro_missing_count$site_week %in%
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
580 siteweeknocount$Group.1[siteweeknocount$x == TRUE], ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
581 pro_count_agg <- aggregate(pro_missing_count$NM, by = list(pro_missing_count$SITE),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
582 function(x) 1 - sum(x, na.rm = T))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
583 names(pro_count_agg) <- c("SITE", "PROP_PHENO_SAMPLED")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
584
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
585 # remove samples outside the monitoring window
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
586 sp_data_site$COUNT[sp_data_site$nm==0] <- NA
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
587
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
588 # Compute the regional GAM index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
589 if(length(unique(sp_data_site$SITE))>1){
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
590 glm_obj_site <- glm(COUNT ~ factor(SITE) + offset(log(nm)) - 1, data = sp_data_site,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
591 family = quasipoisson(link = "log"), control = list(maxit = 100))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
592 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
593 glm_obj_site <- glm(COUNT ~ offset(log(nm)) - 1, data = sp_data_site,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
594 family = quasipoisson(link = "log"), control = list(maxit = 100))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
595 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
596
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
597 sp_data_site[, "FITTED"] <- predict.glm(glm_obj_site, newdata = sp_data_site,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
598 type = "response")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
599 sp_data_site[, "COUNT_IMPUTED"] <- sp_data_site$COUNT
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
600 sp_data_site[is.na(sp_data_site$COUNT), "COUNT_IMPUTED"] <- sp_data_site$FITTED[is.na(sp_data_site$COUNT)]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
601
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
602 # add fitted value for missing mid-week data
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
603 sp_data_site <- sp_data_site[!paste(sp_data_site$DAY_WEEK, sp_data_site$COUNT) %in%
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
604 c("1 NA", "2 NA", "3 NA", "5 NA", "6 NA", "7 NA"), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
605
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
606 # remove all added mid-week values for weeks with real count
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
607 # (observation)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
608 sp_data_site$site_week <- paste(sp_data_site$SITE, sp_data_site$WEEK,
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
609 sep = "_")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
610 siteweekcount <- aggregate(sp_data_site$COUNT, by = list(sp_data_site$site_week),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
611 function(x) sum(!is.na(x)) > 0)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
612 sp_data_site <- sp_data_site[!(is.na(sp_data_site$COUNT) & (sp_data_site$site_week %in%
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
613 siteweekcount$Group.1[siteweekcount$x == TRUE])), names(sp_data_site) !=
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
614 "site_week"]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
615
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
616 # Compute the regional gam index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
617 print(paste("Compute index for",sp_data_site$SPECIES[1],"at year", y,"for",length(unique(sp_data_site$SITE)),"sites:",Sys.time()))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
618 regional_gam_index <- trap_index(sp_data_site, data_col = "COUNT_IMPUTED",
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
619 time_col = "DAYNO", by_col = c("SPECIES", "SITE", "YEAR"))
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
620
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
621 cumu_index <- merge(regional_gam_index, pro_count_agg, by = c("SITE"),
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
622 all.x = TRUE, sort = FALSE)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
623 names(cumu_index) <- c("SITE", "SPECIES", "YEAR", "regional_gam", "prop_pheno_sampled")
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
624
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
625 cumu_index <- cumu_index[order(cumu_index$SITE), ]
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
626
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
627 # bind if exist else create
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
628 if ("cumullated_indices" %in% ls()) {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
629 cumullated_indices <- rbind(cumullated_indices, cumu_index)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
630 } else {
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
631 cumullated_indices <- cumu_index
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
632 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
633
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
634 }
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
635
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
636 return(cumullated_indices)
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
637
5b126f770671 Uploaded
mnhn65mo
parents:
diff changeset
638 }