Commit ⋅ weecology/portalcasting

#' @title Cast Portal Rodents Models
#' @title Forecast Portal Rodents Models
#'
#' @description Cast the Portal rodent population data using the (updated if needed) data and models in a portalcasting directory. \cr \cr
#' @description Forecast the Portal rodent population data using the (updated if needed) data and models in a portalcasting directory. \cr \cr
#'  \code{portalcast} wraps around \code{cast} to allow multiple runs of multiple models, with data preparation as needed between runs occurring via \code{prepare_data}. \cr \cr
#'  \code{cast} runs a single cast of multiple models across data sets.
#'

@@ -161,7 +161,7 @@

  # patch

  cast_tab$obs   <- NA
  cast_dataset   <- gsub("_interp", "", cast_tab$dataset)
  cast_dataset   <- gsub("dm_", "", gsub("_interp", "", cast_tab$dataset))
  ucast_dataset  <- unique(cast_dataset)
  ncast_datasets <- length(ucast_dataset)


@@ -0,0 +1,175 @@

1	+	# these functions are to be moved to the ewRL package
2	+	# the functionality needs to be split out better again
3	+	# see also the code in portalr
4	+
5	+	#' @title Download the Portal Predictions Repository Archive
6	+	#'
7	+	#' @description Downloads a specific \code{version} of the Portal Predictions repository from either GitHub or Zenodo (based on \code{source}) into the \code{<main>/raw} sub.
8	+	#'
9	+	#' @param main \code{character} value defining the main component of the portalcasting directory tree.
10	+	#'
11	+	#' @param resources_sub \code{character} value defining the resources subdirectory of the portalcasting directory tree.
12	+	#'
13	+	#' @param version \code{character} version of the data to download. Default \code{"latest"} downloads the most recent (by date published). \code{NULL} means no download.
14	+	#'
15	+	#' @param source \code{character} indicator of the source for the download. Either \code{"github"} (default) or \code{"github"}.
16	+	#'
17	+	#' @param pause Positive \code{integer} or integer \code{numeric} seconds for pausing during steps around unzipping that require time delay.
18	+	#'
19	+	#' @param timeout Positive \code{integer} or integer \code{numeric} seconds for timeout on downloads. Temporarily overrides the \code{"timeout"} option in \code{\link[base]{options}}.
20	+	#'
21	+	#' @param quiet \code{logical} indicator if progress messages should be quieted.
22	+	#'
23	+	#' @param verbose \code{logical} indicator if detailed messages should be printed.
24	+	#'
25	+	#' @note There are two calls to \code{link[base]{Sys.sleep}} for \code{pause} seconds each to allow for the file unzipping, copying, and such to catch up.
26	+	#'
27	+	#' @return \code{NULL}, \code{\link[base]{invisible}}-ly.
28	+	#'
29	+	#' @examples
30	+	#' \donttest{
31	+	#'
32	+	#' create_dir(main = "./portalcasting")
33	+	#' download_archive(main = "./portalcasting")
34	+	#' }
35	+	#'
36	+	#' @name download archive
37	+	#'
38	+	#' @export
39	+	#'
40	+	download_archive <- function(main = ".",
41	+	resources_sub = "resources",
42	+	version = "latest",
43	+	source = "github",
44	+	quiet = FALSE,
45	+	verbose = FALSE,
46	+	pause = 30,
47	+	timeout = getOption("timeout")) {
48	+
49	+
50	+	return_if_null(version)
51	+
52	+	timeout_backup <- getOption("timeout")
53	+	on.exit(options(timeout = timeout_backup))
54	+	options(timeout = timeout)
55	+
56	+	version <- tolower(version)
57	+	source <- tolower(source)
58	+
59	+	if (source == "zenodo") {
60	+
61	+	base_url <-
62	+
63	+	got <- GET(url = "https://zenodo.org/api/records/",
64	+	query = list(q = "conceptrecid:833438",
65	+	size = 9999,
66	+	all_versions = "true"))
67	+
68	+	stop_for_status(x = got,
69	+	task = paste0("locate Zenodo concept record"))
70	+
71	+	contents <- content(x = got)
72	+
73	+	metadata <- lapply(FUN = getElement,
74	+	X = contents,
75	+	name = "metadata")
76	+	versions <- sapply(FUN = getElement,
77	+	X = metadata,
78	+	name = "version")
79	+	pub_date <- sapply(FUN = getElement,
80	+	X = metadata,
81	+	name = "publication_date")
82	+
83	+	selected <- ifelse(test = version == "latest",
84	+	yes = which.max(as.Date(pub_date)),
85	+	no = which(versions == version))
86	+
87	+	if (length(selected) == 0){
88	+
89	+	stop("Failed to locate version `", version, "` on Zenodo")
90	+
91	+	}
92	+
93	+	zipball_url <- contents[[selected]]$files[[1]]$links$download
94	+	version <- ifelse(test = version == "latest",
95	+	yes = metadata[[selected]]$version,
96	+	no = version)
97	+
98	+	} else if (source == "github") {
99	+
100	+	url <- ifelse(version == "latest",
101	+	"https://api.github.com/repos/weecology/portalPredictions/releases/latest",
102	+	paste0("https://api.github.com/repos/weecology/portalPredictions/releases/tags/", version))
103	+
104	+	got <- GET(url = url)
105	+
106	+	stop_for_status(x = got,
107	+	task = paste0("locate version `", version, "` on GitHub"))
108	+
109	+	zipball_url <- content(got)$zipball_url
110	+
111	+	version <- ifelse(version == "latest", content(got)$tag_name, version)
112	+
113	+	} else {
114	+
115	+	stop("`source` must be either 'zenodo' or 'github'")
116	+
117	+	}
118	+
119	+	messageq("Downloading archive version `", version, "` ...", quiet = quiet)
120	+
121	+	temp <- file.path(tempdir(), "portalPredictions.zip")
122	+	final <- file.path(main, resources_sub, "portalPredictions")
123	+
124	+	result <- tryCatch(
125	+	expr = download.file(url = zipball_url,
126	+	destfile = temp,
127	+	quiet = !verbose,
128	+	mode = "wb"),
129	+	error = function(x){NA})
130	+
131	+	if (is.na(result)) {
132	+
133	+	warning("Archive version `", version, "` could not be downloaded")
134	+	return(invisible())
135	+
136	+	}
137	+
138	+
139	+	if (file.exists(final)) {
140	+
141	+	old_files <- list.files(path = final,
142	+	full.names = TRUE,
143	+	all.files = TRUE,
144	+	recursive = TRUE,
145	+	include.dirs = FALSE)
146	+
147	+	file.remove(old_files)
148	+
149	+	unlink(x = final,
150	+	recursive = TRUE)
151	+
152	+	}
153	+
154	+	folder_name <- unzip(temp, list = TRUE)$Name[1]
155	+
156	+	temp_unzip <- file.path(main, resources_sub, folder_name)
157	+
158	+	unzip(temp, exdir = file.path(main, resources_sub))
159	+
160	+	Sys.sleep(pause)
161	+
162	+	dir.create(final)
163	+
164	+	file.copy(list.files(temp_unzip, full.names = TRUE),
165	+	final,
166	+	recursive = TRUE)
167	+
168	+	Sys.sleep(pause)
169	+
170	+	unlink(temp_unzip, recursive = TRUE)
171	+	file.remove(temp)
172	+
173	+	invisible()
174	+
175	+	}

@@ -47,7 +47,6 @@

         recursive    = TRUE,
         showWarnings = FALSE)


  write_directory_configuration(main     = main, 
                                settings = settings, 
                                quiet    = quiet)

@@ -1,174 +1,3 @@

#' @title Download the Portal Predictions Repository Archive
#'
#' @description Downloads a specific \code{version} of the Portal Predictions repository from either GitHub or Zenodo (based on \code{source}) into the \code{<main>/raw} sub.
#'
#' @param main \code{character} value defining the main component of the portalcasting directory tree. 
#'
#' @param resources_sub \code{character} value defining the resources subdirectory of the portalcasting directory tree. 
#'
#' @param version \code{character} version of the data to download. Default \code{"latest"} downloads the most recent (by date published). \code{NULL} means no download. 
#'
#' @param source \code{character} indicator of the source for the download. Either \code{"github"} (default) or \code{"github"}.
#'
#' @param pause Positive \code{integer} or integer \code{numeric} seconds for pausing during steps around unzipping that require time delay. 
#'
#' @param timeout Positive \code{integer} or integer \code{numeric} seconds for timeout on downloads. Temporarily overrides the \code{"timeout"} option in \code{\link[base]{options}}.
#'
#' @param quiet \code{logical} indicator if progress messages should be quieted.
#'
#' @param verbose \code{logical} indicator if detailed messages should be printed.
#'
#' @note There are two calls to \code{link[base]{Sys.sleep}} for \code{pause} seconds each to allow for the file unzipping, copying, and such to catch up.
#'
#' @return \code{NULL}, \code{\link[base]{invisible}}-ly.
#'
#' @examples
#'  \donttest{
#'
#'   create_dir(main = "./portalcasting")
#'   download_archive(main = "./portalcasting")
#'  } 
#'
#' @name download archive
#'
#' @export
#'
download_archive <- function(main          = ".",
                             resources_sub = "resources",
                             version       = "latest", 
                             source        = "github",
                             quiet         = FALSE,
                             verbose       = FALSE,
                             pause         = 30,
                             timeout       = getOption("timeout")) {


  return_if_null(version)

  timeout_backup <- getOption("timeout")
  on.exit(options(timeout = timeout_backup))
  options(timeout = timeout) 

  version <- tolower(version)
  source  <- tolower(source)

  if (source == "zenodo") {

    base_url <-  

    got <- GET(url   = "https://zenodo.org/api/records/", 
               query = list(q            = "conceptrecid:833438",
                            size         = 9999, 
                            all_versions = "true"))

    stop_for_status(x    = got,
                    task = paste0("locate Zenodo concept record"))

    contents <- content(x = got)    

    metadata <- lapply(FUN  = getElement, 
                       X    = contents, 
                       name = "metadata")
    versions <- sapply(FUN  = getElement, 
                       X    = metadata, 
                       name = "version")
    pub_date <- sapply(FUN  = getElement, 
                       X    = metadata, 
                       name = "publication_date")

    selected <- ifelse(test = version == "latest",
                       yes  = which.max(as.Date(pub_date)),
                       no   = which(versions == version))

    if (length(selected) == 0){

      stop("Failed to locate version `", version, "` on Zenodo")
   
    }
    
    zipball_url <- contents[[selected]]$files[[1]]$links$download     
    version <- ifelse(test = version == "latest", 
                      yes  = metadata[[selected]]$version,
                      no   = version)

  } else if (source == "github") {

    url <- ifelse(version == "latest", 
                  "https://api.github.com/repos/weecology/portalPredictions/releases/latest",
                  paste0("https://api.github.com/repos/weecology/portalPredictions/releases/tags/", version))

    got <- GET(url = url)

    stop_for_status(x    = got, 
                    task = paste0("locate version `", version, "` on GitHub"))

    zipball_url <- content(got)$zipball_url      
 
    version <- ifelse(version == "latest", content(got)$tag_name, version)

  } else {

    stop("`source` must be either 'zenodo' or 'github'")

  }
  
  messageq("Downloading archive version `", version, "` ...", quiet = quiet)

  temp  <- file.path(tempdir(), "portalPredictions.zip")
  final <- file.path(main, resources_sub, "portalPredictions")

  result <- tryCatch(
              expr  = download.file(url      = zipball_url, 
                                    destfile = temp, 
                                    quiet    = !verbose, 
                                    mode     = "wb"),
              error = function(x){NA})

  if (is.na(result)) {

    warning("Archive version `", version, "` could not be downloaded")
    return(invisible())

  }


  if (file.exists(final)) {

    old_files <- list.files(path         = final,
                            full.names   = TRUE,
                            all.files    = TRUE,
                            recursive    = TRUE,
                            include.dirs = FALSE)

    file.remove(old_files)

    unlink(x         = final, 
           recursive = TRUE)

  }

  folder_name <- unzip(temp, list = TRUE)$Name[1]

  temp_unzip <- file.path(main, resources_sub, folder_name)

  unzip(temp, exdir = file.path(main, resources_sub))

  Sys.sleep(pause)

  dir.create(final)

  file.copy(list.files(temp_unzip, full.names = TRUE), 
            final, 
            recursive = TRUE)

  Sys.sleep(pause)

  unlink(temp_unzip, recursive = TRUE)
  file.remove(temp)

  invisible()

}

#' @title Download Climate Forecasts
#'

@@ -2,30 +2,45 @@

#'
#' @description Most users will not want or need to change the directory folders and file names, but it is helpful to have them be flexible for certain circumstances, and this function gathers them into a list for pipeline functionality.
#'
#' @param directory_configuration_file \code{character} value of the path to the directory config YAML.
#' @param files \code{list} of \code{character} names of standard files, see \code{\link{directory_files}}.
#'
#' @param moons_file \code{character} name of the file for saving the moons data.
#' @param subdirectories \code{list} of \code{character} names of standard subdirectories, see \code{\link{directory_subdirectories}}.
#'
#' @param covariates_file,historical_covariates_file,forecast_covariates_file \code{character} name of the files for saving the covariate data (combined historical and forecast) and each component.
#' @param resources \code{list} of \code{list}s for standard resources, see \code{\link{directory_resources}}.
#'
#' @param metadata_file \code{character} name of the file for saving the forecast metadata.
#' @param directory_configuration \code{character} name for the directory configuration YAML.
#'
#' @param dataset_controls_file \code{character} name of the file for saving the rodent dataset controls.
#' @param moons \code{character} name for the lunar data csv.
#'
#' @param model_controls_file \code{character} name of the file for saving the model controls.
#' @param covariates \code{character} name for the combined historical and forecast covariates csv.
#'
#' @param cast_evaluations_file \code{character} name of the file for saving the cast evaluations.
#' @param historical_covariates \code{character} name for the historical covariates csv.
#'
#' @param subdirectories \code{character} vector of the subdirectory names. Default includes \code{raw}, \code{data}, \code{models}, \code{fits}, and \code{casts}. 
#' @param forecast_covariates \code{character} name for the forecast covariates csv.
#'
#' @param PortalData \code{list} with \code{source} and \code{version} elements that are \code{character} values for the source and version of the Portal Data to download. Default values retrieve the latest data from github. \cr \cr
#'                   See \code{\link[portalr]{download_observations}}.
#' @param metadata \code{character} name for the Forecast metadata YAML.
#'
#' @param portalPredictions \code{list} with \code{source} and \code{version} elements that are \code{character} values for the source and version of the archive to download. Default values point to github, but \code{version = NULL} indicates no download. \cr \cr 
#'                          See \code{\link{download_archive}}.
#' @param dataset_controls \code{character} name for the YAML of dataset control list(s).
#'
#' @param climate_forecast \code{list} with \code{source} and \code{version} elements that are \code{character} values for the source and version of the climate forecasts to download. Default values retrieve the current day's forecast from the Northwest Knowledge Network's North American Multi-Model Ensemble (NMME) climate forecasts. \cr \cr 
#'                         See \code{\link{download_climate_forecasts}}.
#' @param model_controls \code{character} name for the YAML of model controls list(s).
#'
#' @param cast_evaluations \code{character} name for the forecast evaluations csv.
#'
#' @param resources \code{character} name for the resources subdirectory.
#'
#' @param data \code{character} name for the data subdirectory.
#'
#' @param models \code{character} name for the models subdirectory.
#'
#' @param fits \code{character} name for the fits subdirectory.
#'
#' @param forecasts \code{character} name for the forecasts subdirectory.
#'
#' @param PortalData \code{list} of \code{source} and \code{version} elements of \code{character} values for the Portal Data download. Default values retrieve the latest data from github
#'
#' @param portalPredictions \code{list} of \code{source} and \code{version} elements of \code{character} values for the archive download. Default values point to github, but \code{version = NULL} indicates no download.
#'
#' @param climate_forecast \code{list} of \code{source}, \code{version}, and \code{data} elements of \code{character} values for the climate forecasts download. Default values retrieve the current day's forecast of min, mean, and max temperature and precipitation from the Northwest Knowledge Network's North American Multi-Model Ensemble (NMME) climate forecasts.
#'
#' @param save \code{logical} indicator controlling if the output should be saved out.
#'

@@ -35,7 +50,7 @@

#'
#' @param download_timeout Positive \code{integer} or integer \code{numeric} seconds for timeout on downloads. Temporarily overrides the \code{"timeout"} option in \code{\link[base]{options}}.
#'
#' @return Named \code{list} of settings for the directory.
#' @return Named \code{list} of settings for the directory (for \code{directory_settings}) or \code{list} of settings components (for \code{directory_files}, \code{directory_subdirectories}, and \code{directory_resources}).
#'
#' @name directory settings
#'

@@ -45,37 +60,17 @@

#'
#' @export
#'
directory_settings <- function (directory_configuration_file = "directory_configuration.yaml",
                                moons_file                   = "moon_dates.csv",
                                covariates_file              = "covariates.csv",
                                historical_covariates_file   = "historical_covariates.csv",
                                forecast_covariates_file     = "forecast_covariates.csv",
                                dataset_controls_file        = "dataset_controls.yaml", 
                                model_controls_file          = "model_controls.yaml",
                                cast_evaluations_file        = "cast_evaluations.csv",
                                metadata_file                = "metadata.yaml",
                                subdirectories               = list("forecasts" = "forecasts", "fits" = "fits", "models" = "models", "resources" = "resources", "data" = "data"),
                                PortalData                   = list(source = "github", version = "latest"),
                                portalPredictions            = list(source = "github", version = NULL),
                                climate_forecast             = list(source = "NMME", version = as.character(Sys.Date()), data = c("tasmin", "tasmean", "tasmax", "pr")),
                                save                         = TRUE,
                                overwrite                    = TRUE, 
                                unzip_pause                  = 30,
                                download_timeout             = getOption("timeout")) {

  list(files            = list(directory_configuration = directory_configuration_file,
                               moons                   = moons_file,
                               covariates              = covariates_file,
                               historical_covariates   = historical_covariates_file,
                               forecast_covariates     = forecast_covariates_file,
                               dataset_controls        = dataset_controls_file,
                               model_controls          = model_controls_file,
                               cast_evaluations        = cast_evaluations_file,
                               metadata                = metadata_file),
directory_settings <- function (files             = directory_files( ),
                                subdirectories    = directory_subdirectories( ),
                                resources         = directory_resources( ),
                                save              = TRUE,
                                overwrite         = TRUE, 
                                unzip_pause       = 30,
                                download_timeout  = getOption("timeout")) {

  list(files            = files,
       subdirectories   = subdirectories,
       resources        = list(PortalData              = PortalData,
                               portalPredictions       = portalPredictions,
                               climate_forecast        = climate_forecast),
       resources        = resources,
       repository       = "portalPredictions",
       save             = save, 
       overwrite        = overwrite, 

@@ -84,16 +79,81 @@

84	79
85	80		}
86	81
	82	+	#' @rdname directory-settings
	83	+	#'
	84	+	#' @export
	85	+	#'
	86	+	directory_files <- function (directory_configuration = "directory_configuration.yaml",
	87	+	moons = "moon_dates.csv",
	88	+	covariates = "covariates.csv",
	89	+	historical_covariates = "historical_covariates.csv",
	90	+	forecast_covariates = "forecast_covariates.csv",
	91	+	dataset_controls = "dataset_controls.yaml",
	92	+	model_controls = "model_controls.yaml",
	93	+	cast_evaluations = "cast_evaluations.csv",
	94	+	metadata = "metadata.yaml") {
	95	+
	96	+	list(directory_configuration = directory_configuration,
	97	+	moons = moons,
	98	+	covariates = covariates,
	99	+	historical_covariates = historical_covariates,
	100	+	forecast_covariates = forecast_covariates,
	101	+	dataset_controls = dataset_controls,
	102	+	model_controls = model_controls,
	103	+	cast_evaluations = cast_evaluations,
	104	+	metadata = metadata)
	105	+
	106	+	}
	107	+
87	108
88	109		#' @rdname directory-settings
89	110		#'
90	111		#' @export
91	112		#'
92		-	production_settings <- function ( ) {
	113	+	directory_subdirectories <- function (forecasts = "forecasts",
	114	+	fits = "fits",
	115	+	models = "models",
	116	+	resources = "resources",
	117	+	data = "data") {
	118	+
	119	+	list(forecasts = forecasts,
	120	+	fits = fits,
	121	+	models = models,
	122	+	resources = resources,
	123	+	data = data)
	124	+
	125	+	}
	126	+
	127	+	#' @rdname directory-settings
	128	+	#'
	129	+	#' @export
	130	+	#'
	131	+	directory_resources <- function (PortalData = list(source = "github",
	132	+	version = "latest"),
	133	+	portalPredictions = list(source = "github",
	134	+	version = NULL),
	135	+	climate_forecast = list(source = "NMME",
	136	+	version = as.character(Sys.Date()),
	137	+	data = c("tasmin", "tasmean", "tasmax", "pr"))) {
	138	+
	139	+	list(PortalData = PortalData,
	140	+	portalPredictions = portalPredictions,
	141	+	climate_forecast = climate_forecast)
	142	+
	143	+	}
	144	+
	145	+
	146	+	#' @rdname directory-settings
	147	+	#'
	148	+	#' @export
	149	+	#'
	150	+	production_settings <- function (download_timeout = max(getOption("timeout"), 600)) {
	151	+
	152	+	resources <- directory_resources(portalPredictions = list(source = "github",
	153	+	version = "latest"))
93	154
94		-	directory_settings(portalPredictions = list(source = "github",
95		-	version = "latest"),
96		-	download_timeout = 600)
	155	+	directory_settings(resources = resources,
	156	+	download_timeout = download_timeout)
97	157
98	158		}
99	159

@@ -74,15 +74,27 @@

  weather_data <- weather(level = "newmoon", fill = TRUE, horizon = 28, path = file.path(main, settings$subdirectories$resources))
  ndvi_data    <- ndvi(level = "newmoon", fill = TRUE, path = file.path(main, settings$subdirectories$resources))

  out                     <- weather_data
  out$source              <- "historic"
  out$ndvi                <- NA
  moon_match              <- match(ndvi_data$newmoonnumber, out$newmoonnumber)
  out$ndvi[moon_match]    <- ndvi_data$ndvi
  out$ndvi_12_month       <- as.numeric(filter(out$ndvi, rep(1, 12), sides = 1))
  out$warm_precip_3_month <- as.numeric(filter(out$warm_precip, rep(1, 3), sides = 1))
  rodents_table <- read_rodents_table(main     = main,
                                      dataset = "controls", 
                                      settings = settings) 

  out                              <- weather_data
  out$source                       <- "historic"
  out$ndvi                         <- NA
  moon_match                       <- match(ndvi_data$newmoonnumber, out$newmoonnumber)
  out$ndvi[moon_match]             <- ndvi_data$ndvi
  out$ndvi_12_month                <- as.numeric(filter(out$ndvi, rep(1, 12), sides = 1))
  out$warm_precip_3_month          <- as.numeric(filter(out$warm_precip, rep(1, 3), sides = 1))
  out$spectab_controls             <- NA
  for (i in 1:nrow(out)) {
    spot <- which(rodents_table$newmoonnumber == out$newmoonnumber[i])
    if (length(spot) == 1) {
      out$spectab_controls[i] <- rodents_table$DS[spot] 
    }
  }
  out$spectab_controls <- as.numeric(round(pmax(na.interp(out$spectab_controls), 0) ))

  cols_to_keep <- c("newmoonnumber", "date", "mintemp", "maxtemp", "meantemp", "precipitation", "warm_precip", "ndvi", "warm_precip_3_month", "ndvi_12_month", "source")
  cols_to_keep <- c("newmoonnumber", "date", "mintemp", "maxtemp", "meantemp", "precipitation", "warm_precip", "ndvi", "warm_precip_3_month", "ndvi_12_month", "spectab_controls", "source")

  write_data(x         = out[ , cols_to_keep], 
             main      = main, 

@@ -242,7 +254,10 @@

  ndvi_12_months        <- filter(c(historical_covariates$ndvi, ndvis), rep(1, 12), sides = 1)
  warm_precip_3_months  <- filter(c(historical_covariates$warm_precip, warm_precips), rep(1, 3), sides = 1)

  past <- list(past_obs = 1, past_mean = 12)
  spectab_controls_mod  <- tsglm(historical_covariates$spectab_controls, model = past, distr = "poisson", link = "log")

  spectab_controls_cast <- predict(spectab_controls_mod, nhist_time)$pred

  hist_climate_forecasts <- data.frame(newmoonnumber       = hist_time, 
                                       date                = moons$newmoondate[match(hist_time, moons$newmoonnumber)],

@@ -254,6 +269,7 @@

                                       ndvi                = ndvis,
                                       warm_precip_3_month = warm_precip_3_months[(length(warm_precip_3_months) - nhist_time + 1):length(warm_precip_3_months)], 
                                       ndvi_12_month       = ndvi_12_months[(length(ndvi_12_months) - nhist_time + 1):length(ndvi_12_months)], 
                                       spectab_controls    = spectab_controls_cast,
                                       source              = "forecast")



@@ -1,9 +1,582 @@

	1	+	#' @rdname prefab_model_functions
	2	+	#'
	3	+	#' @export
	4	+	#'
	5	+	jags_logistic_competition <- function (main = ".",
	6	+	dataset = "controls",
	7	+	settings = directory_settings(),
	8	+	control_runjags = runjags_control(),
	9	+	quiet = FALSE,
	10	+	verbose = FALSE) {
	11	+
	12	+
	13	+	dataset <- tolower(dataset)
	14	+	messageq(" -jags_logistic_competition for ", dataset, quiet = quiet)
	15	+
	16	+	monitor <- c("mu", "sigma", "r_int", "log_K_int", "log_K_slope")
	17	+
	18	+
	19	+	inits <- function (data = NULL) {
	20	+
	21	+	rngs <- c("base::Wichmann-Hill", "base::Marsaglia-Multicarry", "base::Super-Duper", "base::Mersenne-Twister")
	22	+
	23	+	log_mean_past_count <- log(mean(data$past_count, na.rm = TRUE))
	24	+	log_max_past_count <- log(max(data$past_count, na.rm = TRUE))
	25	+
	26	+	sd_log_past_count <- sd(log(data$past_count), na.rm = TRUE)
	27	+
	28	+	function (chain = chain) {
	29	+
	30	+	mu <- rnorm(1, log_mean_past_count, 0.1)
	31	+	log_K_int <- rnorm(1, log_max_past_count, 0.01)
	32	+
	33	+	log_K_int <- max(c(mu, log_K_int)) + 0.1
	34	+
	35	+	list(.RNG.name = sample(rngs, 1),
	36	+	.RNG.seed = sample(1:1e+06, 1),
	37	+	mu = mu,
	38	+	sigma = runif(1, 0, 0.0005),
	39	+	r_int = rnorm(1, 0, 0.01),
	40	+	log_K_int = log_K_int,
	41	+	log_K_slope = rnorm(1, 0, 0.01))
	42	+
	43	+	}
	44	+
	45	+	}
	46	+
	47	+	jags_model <- "model {
	48	+
	49	+	# priors
	50	+
	51	+	mu ~ dnorm(log_mean_past_count, 5)
	52	+	sigma ~ dunif(0, 0.001)
	53	+	tau <- pow(sigma, -2)
	54	+	r_int ~ dnorm(0, 5)
	55	+	log_K_int ~ dnorm(log_max_past_count, 5)
	56	+	log_K_slope ~ dnorm(0, 1)
	57	+
	58	+	# initial state
	59	+
	60	+	log_X[1] <- mu
	61	+	X[1] <- exp(log_X[1])
	62	+	count[1] ~ dpois(X[1])
	63	+
	64	+	# exampand parameters
	65	+
	66	+	for (i in 1:N) {
	67	+
	68	+	r[i] <- r_int
	69	+	K[i] <- exp(log_K_int + log_K_slope * spectab[i])
	70	+
	71	+	}
	72	+
	73	+
	74	+	# through time
	75	+
	76	+	for(i in 2:N) {
	77	+
	78	+	# Process model
	79	+
	80	+	pred_X[i] <- X[i-1] * exp(r[i] * (1 - (X[i - 1] / K[i])))
	81	+	pred_log_X[i] <- log(pred_X[i])
	82	+	log_X[i] ~ dnorm(pred_log_X[i], tau)
	83	+	X[i] <- exp(log_X[i])
	84	+
	85	+	# observation model
	86	+
	87	+	count[i] ~ dpois(X[i])
	88	+
	89	+	}
	90	+
	91	+	}"
	92	+
	93	+
	94	+
	95	+
	96	+	model_name <- "jags_logistic_covariates"
	97	+
	98	+	runjags.options(silent.jags = control_runjags$silent_jags,
	99	+	silent.runjags = control_runjags$silent_jags)
	100	+
	101	+	rodents_table <- read_rodents_table(main = main,
	102	+	dataset = dataset,
	103	+	settings = settings)
	104	+	covariates <- read_covariates(main = main,
	105	+	settings = settings)
	106	+	metadata <- read_metadata(main = main,
	107	+	settings = settings)
	108	+
	109	+	dataset_controls <- metadata$controls_r[[dataset]]
	110	+	start_moon <- metadata$time$start_moon
	111	+	end_moon <- metadata$time$end_moon
	112	+	true_count_lead <- length(metadata$time$rodent_cast_moons)
	113	+	confidence_level <- metadata$confidence_level
	114	+	dataset_controls <- metadata$dataset_controls[[dataset]]
	115	+
	116	+	species <- species_from_table(rodents_tab = rodents_table,
	117	+	total = TRUE,
	118	+	nadot = TRUE)
	119	+	temp_species <- read_model_controls(main = main, settings = settings)$jags_logistic_competition$species
	120	+	if (temp_species == "all") {
	121	+	species <- species
	122	+	} else {
	123	+	species <- species[species %in% temp_species]
	124	+	}
	125	+	nspecies <- length(species)
	126	+	mods <- named_null_list(species)
	127	+	casts <- named_null_list(species)
	128	+	cast_tab <- data.frame()
	129	+
	130	+	for (i in 1:nspecies) {
	131	+
	132	+	s <- species[i]
	133	+	ss <- gsub("NA.", "NA", s)
	134	+	messageq(" -", ss, quiet = !verbose)
	135	+
	136	+	moon_in <- which(rodents_table$newmoonnumber >= start_moon & rodents_table$newmoonnumber <= end_moon)
	137	+	past_moon_in <- which(rodents_table$newmoonnumber < start_moon)
	138	+	moon <- rodents_table[moon_in, "newmoonnumber"]
	139	+	moon <- c(moon, metadata$time$rodent_cast_moons)
	140	+	past_moon <- rodents_table[past_moon_in, "newmoonnumber"]
	141	+
	142	+	ntraps <- rodents_table[moon_in, "ntraps"]
	143	+	na_traps <- which(is.na(ntraps) == TRUE)
	144	+	ntraps[na_traps] <- 0
	145	+	cast_ntraps <- rep(max(ntraps), true_count_lead)
	146	+	ntraps <- c(ntraps, cast_ntraps)
	147	+	past_ntraps <- rodents_table[past_moon_in, "ntraps"]
	148	+
	149	+	species_in <- which(colnames(rodents_table) == s)
	150	+	count <- rodents_table[moon_in, species_in]
	151	+
	152	+	if (sum(count, na.rm = TRUE) == 0) {
	153	+
	154	+	next()
	155	+
	156	+	}
	157	+
	158	+	cast_count <- rep(NA, true_count_lead)
	159	+	count <- c(count, cast_count)
	160	+	past_count <- rodents_table[past_moon_in, species_in]
	161	+	no_count <- which(is.na(past_count) == TRUE)
	162	+
	163	+	if (length(no_count) > 0) {
	164	+
	165	+	past_moon <- past_moon[-no_count]
	166	+	past_count <- past_count[-no_count]
	167	+	past_ntraps <- past_ntraps[-no_count]
	168	+
	169	+	}
	170	+
	171	+	log_mean_past_count <- log(mean(past_count, na.rm = TRUE))
	172	+	log_max_past_count <- log(max(past_count, na.rm = TRUE))
	173	+
	174	+	warm_rain_three_months <- scale(covariates$warm_precip_3_month[covariates$newmoonnumber >= start_moon & covariates$newmoonnumber <= max(metadata$time$covariate_cast_moons)])
	175	+	ndvi_twelve_months <- scale(covariates$ndvi_12_month[covariates$newmoonnumber >= start_moon & covariates$newmoonnumber <= max(metadata$time$covariate_cast_moons)])
	176	+	spectab_controls <- scale(covariates$spectab_controls[covariates$newmoonnumber >= start_moon & covariates$newmoonnumber <= max(metadata$time$covariate_cast_moons)])
	177	+
	178	+
	179	+	data <- list(count = count,
	180	+	ntraps = ntraps,
	181	+	N = length(count),
	182	+	moon = moon,
	183	+	log_mean_past_count = log_mean_past_count,
	184	+	log_max_past_count = log_max_past_count,
	185	+	past_count = past_count,
	186	+	spectab = spectab_controls[ , 1])
	187	+
	188	+
	189	+	obs_pred_times <- metadata$time$rodent_cast_moons
	190	+	obs_pred_times_spot <- obs_pred_times - metadata$time$start_moon
	191	+
	192	+	train_text <- paste0("train: ", metadata$start_moon, " to ", metadata$end_moon)
	193	+	test_text <- paste0("test: ", metadata$end_moon + 1, " to ", metadata$end_moon + metadata$lead)
	194	+	model_text <- paste0("model: ", model_name, "; ", train_text, "; ", test_text)
	195	+
	196	+	if (control_runjags$cast_obs) {
	197	+
	198	+	obs_pred <- paste0("X[", obs_pred_times_spot, "]")
	199	+	monitor <- c(monitor, obs_pred)
	200	+
	201	+	}
	202	+
	203	+	mods[[i]] <- run.jags(model = jags_model,
	204	+	monitor = monitor,
	205	+	inits = inits(data),
	206	+	data = data,
	207	+	n.chains = control_runjags$nchains,
	208	+	adapt = control_runjags$adapt,
	209	+	burnin = control_runjags$burnin,
	210	+	sample = control_runjags$sample,
	211	+	thin = control_runjags$thin,
	212	+	modules = control_runjags$modules,
	213	+	method = control_runjags$method,
	214	+	factories = control_runjags$factories,
	215	+	mutate = control_runjags$mutate,
	216	+	summarise = FALSE,
	217	+	plots = FALSE)
	218	+
	219	+	if (control_runjags$cast_obs) {
	220	+
	221	+	nchains <- control_runjags$nchains
	222	+	vals <- mods[[i]]$mcmc[[1]]
	223	+
	224	+	if (nchains > 1) {
	225	+
	226	+	for (j in 2:nchains) {
	227	+
	228	+	vals <- rbind(vals, mods[[i]]$mcmc[[j]])
	229	+
	230	+	}
	231	+
	232	+	}
	233	+
	234	+	pred_cols <- grep("X", colnames(vals))
	235	+	vals <- vals[ , pred_cols]
	236	+	point_forecast <- round(apply(vals, 2, mean), 3)
	237	+	HPD <- HPDinterval(as.mcmc(vals))
	238	+	lower_cl <- round(HPD[ , "lower"], 3)
	239	+	upper_cl <- round(HPD[ , "upper"], 3)
	240	+	casts_i <- data.frame(Point.Forecast = point_forecast,
	241	+	lower_cl = lower_cl,
	242	+	upper_cl = upper_cl,
	243	+	moon = obs_pred_times)
	244	+
	245	+	colnames(casts_i)[2:3] <- paste0(c("Lo.", "Hi."), confidence_level * 100)
	246	+	rownames(casts_i) <- NULL
	247	+	casts[[i]] <- casts_i
	248	+
	249	+	cast_tab_i <- data.frame(cast_date = metadata$time$cast_date,
	250	+	cast_month = metadata$time$rodent_cast_months,
	251	+	cast_year = metadata$time$rodent_cast_years,
	252	+	moon = metadata$time$rodent_cast_moons,
	253	+	currency = dataset_controls$args$output,
	254	+	model = model_name,
	255	+	dataset = dataset,
	256	+	species = ss,
	257	+	estimate = point_forecast,
	258	+	lower_pi = lower_cl,
	259	+	upper_pi = upper_cl,
	260	+	cast_group = metadata$cast_group,
	261	+	confidence_level = metadata$confidence_level,
	262	+	start_moon = metadata$time$start_moon,
	263	+	end_moon = metadata$time$end_moon)
	264	+
	265	+	cast_tab <- rbind(cast_tab, cast_tab_i)
	266	+
	267	+	}
	268	+
	269	+	}
	270	+
	271	+	metadata <- update_list(metadata,
	272	+	models = "jags_logistic_covariates",
	273	+	datasets = dataset,
	274	+	dataset_controls = dataset_controls)
	275	+
	276	+	list(metadata = metadata,
	277	+	cast_tab = cast_tab,
	278	+	model_fits = mods,
	279	+	model_casts = casts)
	280	+
	281	+
	282	+	}
	283	+
	284	+
	285	+	#' @rdname prefab_model_functions
	286	+	#'
	287	+	#' @export
	288	+	#'
	289	+	jags_logistic_competition_covariates <- function (main = ".",
	290	+	dataset = "controls",
	291	+	settings = directory_settings(),
	292	+	control_runjags = runjags_control(),
	293	+	quiet = FALSE,
	294	+	verbose = FALSE) {
	295	+
	296	+
	297	+	dataset <- tolower(dataset)
	298	+	messageq(" -jags_logistic_competition_covariates for ", dataset, quiet = quiet)
	299	+
	300	+	monitor <- c("mu", "sigma", "r_int", "r_slope", "log_K_int", "log_K_slope_NDVI", "log_K_slope_DS")
	301	+
	302	+
	303	+	inits <- function (data = NULL) {
	304	+
	305	+	rngs <- c("base::Wichmann-Hill", "base::Marsaglia-Multicarry", "base::Super-Duper", "base::Mersenne-Twister")
	306	+
	307	+	log_mean_past_count <- log(mean(data$past_count, na.rm = TRUE))
	308	+	log_max_past_count <- log(max(data$past_count, na.rm = TRUE))
	309	+
	310	+	sd_log_past_count <- sd(log(data$past_count), na.rm = TRUE)
	311	+
	312	+	function (chain = chain) {
	313	+
	314	+	mu <- rnorm(1, log_mean_past_count, 0.1)
	315	+	log_K_int <- rnorm(1, log_max_past_count, 0.01)
	316	+
	317	+	log_K_int <- max(c(mu, log_K_int)) + 0.1
	318	+
	319	+	list(.RNG.name = sample(rngs, 1),
	320	+	.RNG.seed = sample(1:1e+06, 1),
	321	+	mu = mu,
	322	+	sigma = runif(1, 0, 0.0005),
	323	+	r_int = rnorm(1, 0, 0.01),
	324	+	r_slope = rnorm(1, 0, 0.1),
	325	+	log_K_int = log_K_int,
	326	+	log_K_slope_NDVI = rnorm(1, 0, 0.01),
	327	+	log_K_slope_DS = rnorm(1, 0, 0.01))
	328	+
	329	+	}
	330	+
	331	+	}
	332	+
	333	+	jags_model <- "model {
	334	+
	335	+	# priors
	336	+
	337	+	mu ~ dnorm(log_mean_past_count, 5)
	338	+	sigma ~ dunif(0, 0.001)
	339	+	tau <- pow(sigma, -1/2)
	340	+	r_int ~ dnorm(0, 5)
	341	+	r_slope ~ dnorm(0, 1)
	342	+	log_K_int ~ dnorm(log_max_past_count, 5)
	343	+	log_K_slope_NDVI ~ dnorm(0, 4)
	344	+	log_K_slope_DS ~ dnorm(0, 4)
	345	+
	346	+	# initial state
	347	+
	348	+	log_X[1] <- mu
	349	+	X[1] <- exp(log_X[1])
	350	+	count[1] ~ dpois(X[1])
	351	+
	352	+	# exampand parameters
	353	+
	354	+	for (i in 1:N) {
	355	+
	356	+	r[i] <- r_int + r_slope * warm_rain_three_months[i]
	357	+	K[i] <- exp(log_K_int + log_K_slope_DS * spectab[i] + log_K_slope_NDVI * ndvi_twelve_months[i])
	358	+
	359	+	}
	360	+
	361	+
	362	+	# through time
	363	+
	364	+	for(i in 2:N) {
	365	+
	366	+	# Process model
	367	+
	368	+	pred_X[i] <- X[i-1] * exp(r[i] * (1 - (X[i - 1] / K[i])))
	369	+	pred_log_X[i] <- log(pred_X[i])
	370	+	log_X[i] ~ dnorm(pred_log_X[i], tau)
	371	+	X[i] <- exp(log_X[i])
	372	+
	373	+	# observation model
	374	+
	375	+	count[i] ~ dpois(X[i])
	376	+
	377	+	}
	378	+
	379	+	}"
	380	+
	381	+
	382	+
	383	+
	384	+	model_name <- "jags_logistic_competition_covariates"
	385	+
	386	+	runjags.options(silent.jags = control_runjags$silent_jags,
	387	+	silent.runjags = control_runjags$silent_jags)
	388	+
	389	+	rodents_table <- read_rodents_table(main = main,
	390	+	dataset = dataset,
	391	+	settings = settings)
	392	+	covariates <- read_covariates(main = main,
	393	+	settings = settings)
	394	+	metadata <- read_metadata(main = main,
	395	+	settings = settings)
	396	+
	397	+	dataset_controls <- metadata$controls_r[[dataset]]
	398	+	start_moon <- metadata$time$start_moon
	399	+	end_moon <- metadata$time$end_moon
	400	+	true_count_lead <- length(metadata$time$rodent_cast_moons)
	401	+	confidence_level <- metadata$confidence_level
	402	+	dataset_controls <- metadata$dataset_controls[[dataset]]
	403	+
	404	+	species <- species_from_table(rodents_tab = rodents_table,
	405	+	total = TRUE,
	406	+	nadot = TRUE)
	407	+	temp_species <- read_model_controls(main = main, settings = settings)$jags_logistic_competition_covariates$species
	408	+	if (temp_species == "all") {
	409	+	species <- species
	410	+	} else {
	411	+	species <- species[species %in% temp_species]
	412	+	}
	413	+	nspecies <- length(species)
	414	+	mods <- named_null_list(species)
	415	+	casts <- named_null_list(species)
	416	+	cast_tab <- data.frame()
	417	+
	418	+	for (i in 1:nspecies) {
	419	+
	420	+	s <- species[i]
	421	+	ss <- gsub("NA.", "NA", s)
	422	+	messageq(" -", ss, quiet = !verbose)
	423	+
	424	+	moon_in <- which(rodents_table$newmoonnumber >= start_moon & rodents_table$newmoonnumber <= end_moon)
	425	+	past_moon_in <- which(rodents_table$newmoonnumber < start_moon)
	426	+	moon <- rodents_table[moon_in, "newmoonnumber"]
	427	+	moon <- c(moon, metadata$time$rodent_cast_moons)
	428	+	past_moon <- rodents_table[past_moon_in, "newmoonnumber"]
	429	+
	430	+	ntraps <- rodents_table[moon_in, "ntraps"]
	431	+	na_traps <- which(is.na(ntraps) == TRUE)
	432	+	ntraps[na_traps] <- 0
	433	+	cast_ntraps <- rep(max(ntraps), true_count_lead)
	434	+	ntraps <- c(ntraps, cast_ntraps)
	435	+	past_ntraps <- rodents_table[past_moon_in, "ntraps"]
	436	+
	437	+	species_in <- which(colnames(rodents_table) == s)
	438	+	count <- rodents_table[moon_in, species_in]
	439	+
	440	+	if (sum(count, na.rm = TRUE) == 0) {
	441	+
	442	+	next()
	443	+
	444	+	}
	445	+
	446	+	cast_count <- rep(NA, true_count_lead)
	447	+	count <- c(count, cast_count)
	448	+	past_count <- rodents_table[past_moon_in, species_in]
	449	+	no_count <- which(is.na(past_count) == TRUE)
	450	+
	451	+	if (length(no_count) > 0) {
	452	+
	453	+	past_moon <- past_moon[-no_count]
	454	+	past_count <- past_count[-no_count]
	455	+	past_ntraps <- past_ntraps[-no_count]
	456	+
	457	+	}
	458	+
	459	+	log_mean_past_count <- log(mean(past_count, na.rm = TRUE))
	460	+	log_max_past_count <- log(max(past_count, na.rm = TRUE))
	461	+
	462	+	warm_rain_three_months <- scale(covariates$warm_precip_3_month[covariates$newmoonnumber >= start_moon & covariates$newmoonnumber <= max(metadata$time$covariate_cast_moons)])
	463	+	ndvi_twelve_months <- scale(covariates$ndvi_12_month[covariates$newmoonnumber >= start_moon & covariates$newmoonnumber <= max(metadata$time$covariate_cast_moons)])
	464	+	spectab_controls <- scale(covariates$spectab_controls[covariates$newmoonnumber >= start_moon & covariates$newmoonnumber <= max(metadata$time$covariate_cast_moons)])
	465	+
	466	+
	467	+	data <- list(count = count,
	468	+	ntraps = ntraps,
	469	+	N = length(count),
	470	+	moon = moon,
	471	+	log_mean_past_count = log_mean_past_count,
	472	+	log_max_past_count = log_max_past_count,
	473	+	past_count = past_count,
	474	+	warm_rain_three_months = warm_rain_three_months[ , 1],
	475	+	ndvi_twelve_months = ndvi_twelve_months[ , 1],
	476	+	spectab = spectab_controls[ , 1])
	477	+
	478	+
	479	+	obs_pred_times <- metadata$time$rodent_cast_moons
	480	+	obs_pred_times_spot <- obs_pred_times - metadata$time$start_moon
	481	+
	482	+	train_text <- paste0("train: ", metadata$start_moon, " to ", metadata$end_moon)
	483	+	test_text <- paste0("test: ", metadata$end_moon + 1, " to ", metadata$end_moon + metadata$lead)
	484	+	model_text <- paste0("model: ", model_name, "; ", train_text, "; ", test_text)
	485	+
	486	+	if (control_runjags$cast_obs) {
	487	+
	488	+	obs_pred <- paste0("X[", obs_pred_times_spot, "]")
	489	+	monitor <- c(monitor, obs_pred)
	490	+
	491	+	}
	492	+
	493	+	mods[[i]] <- run.jags(model = jags_model,
	494	+	monitor = monitor,
	495	+	inits = inits(data),
	496	+	data = data,
	497	+	n.chains = control_runjags$nchains,
	498	+	adapt = control_runjags$adapt,
	499	+	burnin = control_runjags$burnin,
	500	+	sample = control_runjags$sample,
	501	+	thin = control_runjags$thin,
	502	+	modules = control_runjags$modules,
	503	+	method = control_runjags$method,
	504	+	factories = control_runjags$factories,
	505	+	mutate = control_runjags$mutate,
	506	+	summarise = FALSE,
	507	+	plots = FALSE)
	508	+
	509	+	if (control_runjags$cast_obs) {
	510	+
	511	+	nchains <- control_runjags$nchains
	512	+	vals <- mods[[i]]$mcmc[[1]]
	513	+
	514	+	if (nchains > 1) {
	515	+
	516	+	for (j in 2:nchains) {
	517	+
	518	+	vals <- rbind(vals, mods[[i]]$mcmc[[j]])
	519	+
	520	+	}
	521	+
	522	+	}
	523	+
	524	+	pred_cols <- grep("X", colnames(vals))
	525	+	vals <- vals[ , pred_cols]
	526	+	point_forecast <- round(apply(vals, 2, mean), 3)
	527	+	HPD <- HPDinterval(as.mcmc(vals))
	528	+	lower_cl <- round(HPD[ , "lower"], 3)
	529	+	upper_cl <- round(HPD[ , "upper"], 3)
	530	+	casts_i <- data.frame(Point.Forecast = point_forecast,
	531	+	lower_cl = lower_cl,
	532	+	upper_cl = upper_cl,
	533	+	moon = obs_pred_times)
	534	+
	535	+	colnames(casts_i)[2:3] <- paste0(c("Lo.", "Hi."), confidence_level * 100)
	536	+	rownames(casts_i) <- NULL
	537	+	casts[[i]] <- casts_i
	538	+
	539	+	cast_tab_i <- data.frame(cast_date = metadata$time$cast_date,
	540	+	cast_month = metadata$time$rodent_cast_months,
	541	+	cast_year = metadata$time$rodent_cast_years,
	542	+	moon = metadata$time$rodent_cast_moons,
	543	+	currency = dataset_controls$args$output,
	544	+	model = model_name,
	545	+	dataset = dataset,
	546	+	species = ss,
	547	+	estimate = point_forecast,
	548	+	lower_pi = lower_cl,
	549	+	upper_pi = upper_cl,
	550	+	cast_group = metadata$cast_group,
	551	+	confidence_level = metadata$confidence_level,
	552	+	start_moon = metadata$time$start_moon,
	553	+	end_moon = metadata$time$end_moon)
	554	+
	555	+	cast_tab <- rbind(cast_tab, cast_tab_i)
	556	+
	557	+	}
	558	+
	559	+	}
	560	+
	561	+	metadata <- update_list(metadata,
	562	+	models = "jags_logistic_covariates",
	563	+	datasets = dataset,
	564	+	dataset_controls = dataset_controls)
	565	+
	566	+	list(metadata = metadata,
	567	+	cast_tab = cast_tab,
	568	+	model_fits = mods,
	569	+	model_casts = casts)
	570	+
	571	+
	572	+	}
	573	+
1	574		#' @rdname prefab_model_functions
2	575		#'
3	576		#' @export
4	577		#'
5	578		jags_logistic_covariates <- function (main = ".",
6		-	dataset = "dm_controls",
	579	+	dataset = "controls",
7	580		settings = directory_settings(),
8	581		control_runjags = runjags_control(),
9	582		quiet = FALSE,

@@ -51,7 +624,7 @@


    mu          ~  dnorm(log_mean_past_count, 5)
    sigma       ~  dunif(0, 1) 
    tau         <- pow(sigma, -1/2)
    tau         <- pow(sigma, -2)
    r_int       ~  dnorm(0, 5)
    r_slope     ~  dnorm(0, 1)
    log_K_int   ~  dnorm(log_max_past_count, 5)

@@ -118,6 +691,13 @@

  species <- species_from_table(rodents_tab = rodents_table, 
                                total       = TRUE, 
                                nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$jags_logistic_covariates$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }

  nspecies <- length(species)
  mods     <- named_null_list(species)
  casts    <- named_null_list(species)

@@ -282,7 +862,7 @@

#' @export
#'
jags_logistic <- function (main            = ".", 
                           dataset         = "dm_controls",  
                           dataset         = "controls",  
                           settings        = directory_settings(), 
                           control_runjags = runjags_control(), 
                           quiet           = FALSE, 

@@ -328,7 +908,7 @@


    mu    ~  dnorm(log_mean_past_count, 5)
    sigma ~  dunif(0, 1) 
    tau   <- pow(sigma, -1/2)
    tau   <- pow(sigma, -2)
    r     ~  dnorm(0, 5)
    log_K ~  dnorm(log_max_past_count, 5)
    K     <- exp(log_K) 

@@ -383,6 +963,13 @@

  species <- species_from_table(rodents_tab = rodents_table, 
                                total       = TRUE, 
                                nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$jags_logistic$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }

  nspecies <- length(species)
  mods     <- named_null_list(species)
  casts    <- named_null_list(species)

@@ -543,7 +1130,7 @@

#' @export
#'
jags_RW <- function (main            = ".", 
                     dataset         = "dm_controls",  
                     dataset         = "controls",  
                     settings        = directory_settings(),
                     control_runjags = runjags_control(), 
                     quiet           = FALSE, 

@@ -577,7 +1164,7 @@


    mu    ~  dnorm(log_mean_past_count, 5)
    sigma ~  dunif(0, 1) 
    tau   <- pow(sigma, -1/2)
    tau   <- pow(sigma, -2)

    # initial state


@@ -627,6 +1214,12 @@

  species <- species_from_table(rodents_tab = rodents_table, 
                                total       = TRUE, 
                                nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$jags_RW$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }
  nspecies <- length(species)
  mods     <- named_null_list(species)
  casts    <- named_null_list(species)

@@ -823,7 +1416,7 @@

                             burnin      = 1e4, 
                             sample      = 1e4, 
                             thin        = 1, 
                             modules     = "", 
                             modules     = "glm", 
                             method      = "interruptible", 
                             factories   = "", 
                             mutate      = NA, 

@@ -52,6 +52,8 @@

#'  \code{jags_RW} fits a log-scale density random walk with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference. \cr \cr
#'  \code{jags_logistic} fits a log-scale density logistic-growth (r-K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference. 
#'  \code{jags_logistic_covariates} fits a log-scale density logistic-growth (r-K) based on covariates (warm rain influencing r, NDVI influencing K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference. 
#'  \code{jags_logistic_competition} fits a log-scale density logistic-growth (r-K) based on covariates (competitor density influencing K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference. 
#'  \code{jags_logistic_competition_covariates} fits a log-scale density logistic-growth (r-K) based on covariates (warm rain influencing r, NDVI and competitor density influencing K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference. 
#'
#' @details 
#'  \code{AutoArima} \cr

@@ -100,6 +102,8 @@

#'   jags_RW()
#'   jags_logistic()
#'   jags_logistic_covariates()
#'   jags_logistic_competition()
#'   jags_logistic_competition_covariates()
#'  }
#'
#' @name prefab_model_functions

@@ -126,6 +130,14 @@

  species       <- species_from_table(rodents_tab = rodents_table, 
                                      total       = TRUE, 
                                      nadot       = TRUE)

  temp_species  <- read_model_controls(main = main, settings = settings)$AutoArima$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }

  nspecies      <- length(species)

  metadata <- read_metadata(main     = main,

@@ -219,6 +231,12 @@

  species       <- species_from_table(rodents_tab = rodents_table, 
                                      total       = TRUE, 
                                      nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$NaiveArima$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }
  nspecies      <- length(species)

  metadata <- read_metadata(main     = main,

@@ -311,6 +329,12 @@

  species       <- species_from_table(rodents_tab = rodents_table, 
                                      total       = TRUE, 
                                      nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$ESSS$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }
  nspecies      <- length(species)

  metadata <- read_metadata(main     = main,

@@ -402,6 +426,12 @@

  species       <- species_from_table(rodents_tab = rodents_table, 
                                      total       = TRUE, 
                                      nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$nbGARCH$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }
  nspecies      <- length(species)

  metadata <- read_metadata(main     = main,

@@ -506,6 +536,12 @@

  species       <- species_from_table(rodents_tab = rodents_table, 
                                      total       = TRUE, 
                                      nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$nbsGARCH$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }
  nspecies      <- length(species)

  metadata <- read_metadata(main     = main,

@@ -633,6 +669,12 @@

  species       <- species_from_table(rodents_tab = rodents_table, 
                                      total       = TRUE, 
                                      nadot       = TRUE)
  temp_species  <- read_model_controls(main = main, settings = settings)$pevGARCH$species
  if (temp_species == "all") {
    species <- species
  } else {
    species <- species[species %in% temp_species]
  }
  nspecies      <- length(species)

  metadata <- read_metadata(main     = main,

@@ -91,7 +91,7 @@

  # patch

  # patch
  cast_tab$dataset <- gsub("_interp", "", cast_tab$dataset)
  cast_tab$dataset <- gsub("dm_", "", gsub("_interp", "", cast_tab$dataset))
  # patch

  cast_ids                <- ifnull(cast_ids, unique(cast_tab$cast_id))

@@ -761,7 +761,7 @@


    obs           <- read_rodents_table(main           = main, 
                                        settings       = settings,
                                        dataset = gsub("_interp", "", casts_meta$dataset))
                                        dataset = gsub("dm_", "", gsub("_interp", "", casts_meta$dataset)))
    colnames(obs) <- gsub("\\.", "", colnames(obs))
    sp_col        <- is_sp_col(obs, nadot = TRUE, total = TRUE)
    species       <- ifnull(species, colnames(obs)[sp_col])

@@ -949,20 +949,11 @@

                          start_moon  = 217, 
                          quiet       = FALSE) {

  model  <- ifnull(model, "Ensemble")
  model2 <- model

  if (!is.null(model) && tolower(model) == "ensemble") {

    model2 <- NULL

  }

  casts_meta <- select_casts(main      = main, 
                             settings  = settings, 
                             cast_ids  = cast_id,
                             end_moons = end_moon, 
                             models    = model2, 
                             models    = model, 
                             datasets  = dataset, 
                             quiet     = quiet)


@@ -986,10 +977,18 @@


  sp_col  <- is_sp_col(obs, nadot = TRUE, total = TRUE)
  species <- ifnull(species, colnames(obs)[sp_col])
  if (length(species) > 1) {
    if ("total" %in% species) {
      species <- "total"
    } else {
      species <- species[1]
    }

  }

  if (!all(species %in% colnames(obs))) {

    stop("cast not available for requested species")
    stop("observations not available for requested species")

  }


@@ -1012,7 +1011,12 @@


  }

  species   <- ifnull(species, unique(preds$species))
  if (!(species %in% preds$species)) {

    stop("cast not available for requested species")

  }

  match_sp  <- (preds$species %in% species)
  colnames  <- c("moon", "estimate", "lower_pi", "upper_pi")
  match_col <- (colnames(preds) %in% colnames)

@@ -1135,7 +1139,4 @@

1135	1139
1136	1140		mtext(text = title, side = 3, cex = 1.25, line = 0.5, at = 217, adj = 0)
1137	1141
1138		-	}
1139		-
1140		-
1141		-
	1142	+	}

Files		•	•	•	Coverage
R	3,428	3,194	0	234	93.17%
Project Totals (22 files)	3,428	3,194	0	234	93.17%

View CI Build Download

coverage:
  precision: 2
  round: down
  status:
    patch:
      default:
        target: 50%
    project:
      default:
        target: 50%

comment: false

Read the docs

Interact now

The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.

SVG

The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.

SVG

Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.

SVG

Learn more about the coverage sunburst graph

627	1214		species <- species_from_table(rodents_tab = rodents_table,
628	1215		total = TRUE,
629	1216		nadot = TRUE)
	1217	+	temp_species <- read_model_controls(main = main, settings = settings)$jags_RW$species
	1218	+	if (temp_species == "all") {
	1219	+	species <- species
	1220	+	} else {
	1221	+	species <- species[species %in% temp_species]
	1222	+	}
630	1223		nspecies <- length(species)
631	1224		mods <- named_null_list(species)
632	1225		casts <- named_null_list(species)

823	1416		burnin = 1e4,
824	1417		sample = 1e4,
825	1418		thin = 1,
826		-	modules = "",
	1419	+	modules = "glm",
827	1420		method = "interruptible",
828	1421		factories = "",
829	1422		mutate = NA,

633	669		species <- species_from_table(rodents_tab = rodents_table,
634	670		total = TRUE,
635	671		nadot = TRUE)
	672	+	temp_species <- read_model_controls(main = main, settings = settings)$pevGARCH$species
	673	+	if (temp_species == "all") {
	674	+	species <- species
	675	+	} else {
	676	+	species <- species[species %in% temp_species]
	677	+	}
636	678		nspecies <- length(species)
637	679
638	680		metadata <- read_metadata(main = main,

761	761
762	762		obs <- read_rodents_table(main = main,
763	763		settings = settings,
764		-	dataset = gsub("_interp", "", casts_meta$dataset))
	764	+	dataset = gsub("dm_", "", gsub("_interp", "", casts_meta$dataset)))
765	765		colnames(obs) <- gsub("\\.", "", colnames(obs))
766	766		sp_col <- is_sp_col(obs, nadot = TRUE, total = TRUE)
767	767		species <- ifnull(species, colnames(obs)[sp_col])

949	949		start_moon = 217,
950	950		quiet = FALSE) {
951	951
952		-	model <- ifnull(model, "Ensemble")
953		-	model2 <- model
954		-
955		-	if (!is.null(model) && tolower(model) == "ensemble") {
956		-
957		-	model2 <- NULL
958		-
959		-	}
960		-
961	952		casts_meta <- select_casts(main = main,
962	953		settings = settings,
963	954		cast_ids = cast_id,
964	955		end_moons = end_moon,
965		-	models = model2,
	956	+	models = model,
966	957		datasets = dataset,
967	958		quiet = quiet)
968	959

1		-	#' @title Cast Portal Rodents Models
	1	+	#' @title Forecast Portal Rodents Models
2	2		#'
3		-	#' @description Cast the Portal rodent population data using the (updated if needed) data and models in a portalcasting directory. \cr \cr
	3	+	#' @description Forecast the Portal rodent population data using the (updated if needed) data and models in a portalcasting directory. \cr \cr
4	4		#' \code{portalcast} wraps around \code{cast} to allow multiple runs of multiple models, with data preparation as needed between runs occurring via \code{prepare_data}. \cr \cr
5	5		#' \code{cast} runs a single cast of multiple models across data sets.
6	6		#'

161	161		# patch
162	162
163	163		cast_tab$obs <- NA
164		-	cast_dataset <- gsub("_interp", "", cast_tab$dataset)
	164	+	cast_dataset <- gsub("dm_", "", gsub("_interp", "", cast_tab$dataset))
165	165		ucast_dataset <- unique(cast_dataset)
166	166		ncast_datasets <- length(ucast_dataset)
167	167

47	47		recursive = TRUE,
48	48		showWarnings = FALSE)
49	49
50		-
51	50		write_directory_configuration(main = main,
52	51		settings = settings,
53	52		quiet = quiet)

2	2		#'
3	3		#' @description Most users will not want or need to change the directory folders and file names, but it is helpful to have them be flexible for certain circumstances, and this function gathers them into a list for pipeline functionality.
4	4		#'
5		-	#' @param directory_configuration_file \code{character} value of the path to the directory config YAML.
	5	+	#' @param files \code{list} of \code{character} names of standard files, see \code{\link{directory_files}}.
6	6		#'
7		-	#' @param moons_file \code{character} name of the file for saving the moons data.
	7	+	#' @param subdirectories \code{list} of \code{character} names of standard subdirectories, see \code{\link{directory_subdirectories}}.
8	8		#'
9		-	#' @param covariates_file,historical_covariates_file,forecast_covariates_file \code{character} name of the files for saving the covariate data (combined historical and forecast) and each component.
	9	+	#' @param resources \code{list} of \code{list}s for standard resources, see \code{\link{directory_resources}}.
10	10		#'
11		-	#' @param metadata_file \code{character} name of the file for saving the forecast metadata.
	11	+	#' @param directory_configuration \code{character} name for the directory configuration YAML.
12	12		#'
13		-	#' @param dataset_controls_file \code{character} name of the file for saving the rodent dataset controls.
	13	+	#' @param moons \code{character} name for the lunar data csv.
14	14		#'
15		-	#' @param model_controls_file \code{character} name of the file for saving the model controls.
	15	+	#' @param covariates \code{character} name for the combined historical and forecast covariates csv.
16	16		#'
17		-	#' @param cast_evaluations_file \code{character} name of the file for saving the cast evaluations.
	17	+	#' @param historical_covariates \code{character} name for the historical covariates csv.
18	18		#'
19		-	#' @param subdirectories \code{character} vector of the subdirectory names. Default includes \code{raw}, \code{data}, \code{models}, \code{fits}, and \code{casts}.
	19	+	#' @param forecast_covariates \code{character} name for the forecast covariates csv.
20	20		#'
21		-	#' @param PortalData \code{list} with \code{source} and \code{version} elements that are \code{character} values for the source and version of the Portal Data to download. Default values retrieve the latest data from github. \cr \cr
22		-	#' See \code{\link[portalr]{download_observations}}.
	21	+	#' @param metadata \code{character} name for the Forecast metadata YAML.
23	22		#'
24		-	#' @param portalPredictions \code{list} with \code{source} and \code{version} elements that are \code{character} values for the source and version of the archive to download. Default values point to github, but \code{version = NULL} indicates no download. \cr \cr
25		-	#' See \code{\link{download_archive}}.
	23	+	#' @param dataset_controls \code{character} name for the YAML of dataset control list(s).
26	24		#'
27		-	#' @param climate_forecast \code{list} with \code{source} and \code{version} elements that are \code{character} values for the source and version of the climate forecasts to download. Default values retrieve the current day's forecast from the Northwest Knowledge Network's North American Multi-Model Ensemble (NMME) climate forecasts. \cr \cr
28		-	#' See \code{\link{download_climate_forecasts}}.
	25	+	#' @param model_controls \code{character} name for the YAML of model controls list(s).
	26	+	#'
	27	+	#' @param cast_evaluations \code{character} name for the forecast evaluations csv.
	28	+	#'
	29	+	#' @param resources \code{character} name for the resources subdirectory.
	30	+	#'
	31	+	#' @param data \code{character} name for the data subdirectory.
	32	+	#'
	33	+	#' @param models \code{character} name for the models subdirectory.
	34	+	#'
	35	+	#' @param fits \code{character} name for the fits subdirectory.
	36	+	#'
	37	+	#' @param forecasts \code{character} name for the forecasts subdirectory.
	38	+	#'
	39	+	#' @param PortalData \code{list} of \code{source} and \code{version} elements of \code{character} values for the Portal Data download. Default values retrieve the latest data from github
	40	+	#'
	41	+	#' @param portalPredictions \code{list} of \code{source} and \code{version} elements of \code{character} values for the archive download. Default values point to github, but \code{version = NULL} indicates no download.
	42	+	#'
	43	+	#' @param climate_forecast \code{list} of \code{source}, \code{version}, and \code{data} elements of \code{character} values for the climate forecasts download. Default values retrieve the current day's forecast of min, mean, and max temperature and precipitation from the Northwest Knowledge Network's North American Multi-Model Ensemble (NMME) climate forecasts.
29	44		#'
30	45		#' @param save \code{logical} indicator controlling if the output should be saved out.
31	46		#'

35	50		#'
36	51		#' @param download_timeout Positive \code{integer} or integer \code{numeric} seconds for timeout on downloads. Temporarily overrides the \code{"timeout"} option in \code{\link[base]{options}}.
37	52		#'
38		-	#' @return Named \code{list} of settings for the directory.
	53	+	#' @return Named \code{list} of settings for the directory (for \code{directory_settings}) or \code{list} of settings components (for \code{directory_files}, \code{directory_subdirectories}, and \code{directory_resources}).
39	54		#'
40	55		#' @name directory settings
41	56		#'

45	60		#'
46	61		#' @export
47	62		#'
48		-	directory_settings <- function (directory_configuration_file = "directory_configuration.yaml",
49		-	moons_file = "moon_dates.csv",
50		-	covariates_file = "covariates.csv",
51		-	historical_covariates_file = "historical_covariates.csv",
52		-	forecast_covariates_file = "forecast_covariates.csv",
53		-	dataset_controls_file = "dataset_controls.yaml",
54		-	model_controls_file = "model_controls.yaml",
55		-	cast_evaluations_file = "cast_evaluations.csv",
56		-	metadata_file = "metadata.yaml",
57		-	subdirectories = list("forecasts" = "forecasts", "fits" = "fits", "models" = "models", "resources" = "resources", "data" = "data"),
58		-	PortalData = list(source = "github", version = "latest"),
59		-	portalPredictions = list(source = "github", version = NULL),
60		-	climate_forecast = list(source = "NMME", version = as.character(Sys.Date()), data = c("tasmin", "tasmean", "tasmax", "pr")),
61		-	save = TRUE,
62		-	overwrite = TRUE,
63		-	unzip_pause = 30,
64		-	download_timeout = getOption("timeout")) {
65		-
66		-	list(files = list(directory_configuration = directory_configuration_file,
67		-	moons = moons_file,
68		-	covariates = covariates_file,
69		-	historical_covariates = historical_covariates_file,
70		-	forecast_covariates = forecast_covariates_file,
71		-	dataset_controls = dataset_controls_file,
72		-	model_controls = model_controls_file,
73		-	cast_evaluations = cast_evaluations_file,
74		-	metadata = metadata_file),
	63	+	directory_settings <- function (files = directory_files( ),
	64	+	subdirectories = directory_subdirectories( ),
	65	+	resources = directory_resources( ),
	66	+	save = TRUE,
	67	+	overwrite = TRUE,
	68	+	unzip_pause = 30,
	69	+	download_timeout = getOption("timeout")) {
	70	+
	71	+	list(files = files,
75	72		subdirectories = subdirectories,
76		-	resources = list(PortalData = PortalData,
77		-	portalPredictions = portalPredictions,
78		-	climate_forecast = climate_forecast),
	73	+	resources = resources,
79	74		repository = "portalPredictions",
80	75		save = save,
81	76		overwrite = overwrite,

74	74		weather_data <- weather(level = "newmoon", fill = TRUE, horizon = 28, path = file.path(main, settings$subdirectories$resources))
75	75		ndvi_data <- ndvi(level = "newmoon", fill = TRUE, path = file.path(main, settings$subdirectories$resources))
76	76
77		-	out <- weather_data
78		-	out$source <- "historic"
79		-	out$ndvi <- NA
80		-	moon_match <- match(ndvi_data$newmoonnumber, out$newmoonnumber)
81		-	out$ndvi[moon_match] <- ndvi_data$ndvi
82		-	out$ndvi_12_month <- as.numeric(filter(out$ndvi, rep(1, 12), sides = 1))
83		-	out$warm_precip_3_month <- as.numeric(filter(out$warm_precip, rep(1, 3), sides = 1))
	77	+	rodents_table <- read_rodents_table(main = main,
	78	+	dataset = "controls",
	79	+	settings = settings)
	80	+
	81	+	out <- weather_data
	82	+	out$source <- "historic"
	83	+	out$ndvi <- NA
	84	+	moon_match <- match(ndvi_data$newmoonnumber, out$newmoonnumber)
	85	+	out$ndvi[moon_match] <- ndvi_data$ndvi
	86	+	out$ndvi_12_month <- as.numeric(filter(out$ndvi, rep(1, 12), sides = 1))
	87	+	out$warm_precip_3_month <- as.numeric(filter(out$warm_precip, rep(1, 3), sides = 1))
	88	+	out$spectab_controls <- NA
	89	+	for (i in 1:nrow(out)) {
	90	+	spot <- which(rodents_table$newmoonnumber == out$newmoonnumber[i])
	91	+	if (length(spot) == 1) {
	92	+	out$spectab_controls[i] <- rodents_table$DS[spot]
	93	+	}
	94	+	}
	95	+	out$spectab_controls <- as.numeric(round(pmax(na.interp(out$spectab_controls), 0) ))
84	96
85		-	cols_to_keep <- c("newmoonnumber", "date", "mintemp", "maxtemp", "meantemp", "precipitation", "warm_precip", "ndvi", "warm_precip_3_month", "ndvi_12_month", "source")
	97	+	cols_to_keep <- c("newmoonnumber", "date", "mintemp", "maxtemp", "meantemp", "precipitation", "warm_precip", "ndvi", "warm_precip_3_month", "ndvi_12_month", "spectab_controls", "source")
86	98
87	99		write_data(x = out[ , cols_to_keep],
88	100		main = main,

242	254		ndvi_12_months <- filter(c(historical_covariates$ndvi, ndvis), rep(1, 12), sides = 1)
243	255		warm_precip_3_months <- filter(c(historical_covariates$warm_precip, warm_precips), rep(1, 3), sides = 1)
244	256
	257	+	past <- list(past_obs = 1, past_mean = 12)
	258	+	spectab_controls_mod <- tsglm(historical_covariates$spectab_controls, model = past, distr = "poisson", link = "log")
245	259
	260	+	spectab_controls_cast <- predict(spectab_controls_mod, nhist_time)$pred
246	261
247	262		hist_climate_forecasts <- data.frame(newmoonnumber = hist_time,
248	263		date = moons$newmoondate[match(hist_time, moons$newmoonnumber)],

254	269		ndvi = ndvis,
255	270		warm_precip_3_month = warm_precip_3_months[(length(warm_precip_3_months) - nhist_time + 1):length(warm_precip_3_months)],
256	271		ndvi_12_month = ndvi_12_months[(length(ndvi_12_months) - nhist_time + 1):length(ndvi_12_months)],
	272	+	spectab_controls = spectab_controls_cast,
257	273		source = "forecast")
258	274
259	275

51	624
52	625		mu ~ dnorm(log_mean_past_count, 5)
53	626		sigma ~ dunif(0, 1)
54		-	tau <- pow(sigma, -1/2)
	627	+	tau <- pow(sigma, -2)
55	628		r_int ~ dnorm(0, 5)
56	629		r_slope ~ dnorm(0, 1)
57	630		log_K_int ~ dnorm(log_max_past_count, 5)

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

282	862		#' @export
283	863		#'
284	864		jags_logistic <- function (main = ".",
285		-	dataset = "dm_controls",
	865	+	dataset = "controls",
286	866		settings = directory_settings(),
287	867		control_runjags = runjags_control(),
288	868		quiet = FALSE,

543	1130		#' @export
544	1131		#'
545	1132		jags_RW <- function (main = ".",
546		-	dataset = "dm_controls",
	1133	+	dataset = "controls",
547	1134		settings = directory_settings(),
548	1135		control_runjags = runjags_control(),
549	1136		quiet = FALSE,

577	1164
578	1165		mu ~ dnorm(log_mean_past_count, 5)
579	1166		sigma ~ dunif(0, 1)
580		-	tau <- pow(sigma, -1/2)
	1167	+	tau <- pow(sigma, -2)
581	1168
582	1169		# initial state
583	1170

52	52		#' \code{jags_RW} fits a log-scale density random walk with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference. \cr \cr
53	53		#' \code{jags_logistic} fits a log-scale density logistic-growth (r-K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference.
54	54		#' \code{jags_logistic_covariates} fits a log-scale density logistic-growth (r-K) based on covariates (warm rain influencing r, NDVI influencing K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference.
	55	+	#' \code{jags_logistic_competition} fits a log-scale density logistic-growth (r-K) based on covariates (competitor density influencing K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference.
	56	+	#' \code{jags_logistic_competition_covariates} fits a log-scale density logistic-growth (r-K) based on covariates (warm rain influencing r, NDVI and competitor density influencing K) model with a Poisson observation process using JAGS (Just Another Gibbs Sampler; Plummer 2003) hierarchical Bayesian inference.
55	57		#'
56	58		#' @details
57	59		#' \code{AutoArima} \cr

100	102		#' jags_RW()
101	103		#' jags_logistic()
102	104		#' jags_logistic_covariates()
	105	+	#' jags_logistic_competition()
	106	+	#' jags_logistic_competition_covariates()
103	107		#' }
104	108		#'
105	109		#' @name prefab_model_functions

91	91		# patch
92	92
93	93		# patch
94		-	cast_tab$dataset <- gsub("_interp", "", cast_tab$dataset)
	94	+	cast_tab$dataset <- gsub("dm_", "", gsub("_interp", "", cast_tab$dataset))
95	95		# patch
96	96
97	97		cast_ids <- ifnull(cast_ids, unique(cast_tab$cast_id))

986	977
987	978		sp_col <- is_sp_col(obs, nadot = TRUE, total = TRUE)
988	979		species <- ifnull(species, colnames(obs)[sp_col])
	980	+	if (length(species) > 1) {
	981	+	if ("total" %in% species) {
	982	+	species <- "total"
	983	+	} else {
	984	+	species <- species[1]
	985	+	}
	986	+
	987	+	}
989	988
990	989		if (!all(species %in% colnames(obs))) {
991	990
992		-	stop("cast not available for requested species")
	991	+	stop("observations not available for requested species")
993	992
994	993		}
995	994

1012	1011
1013	1012		}
1014	1013
1015		-	species <- ifnull(species, unique(preds$species))
	1014	+	if (!(species %in% preds$species)) {
	1015	+
	1016	+	stop("cast not available for requested species")
	1017	+
	1018	+	}
	1019	+
1016	1020		match_sp <- (preds$species %in% species)
1017	1021		colnames <- c("moon", "estimate", "lower_pi", "upper_pi")
1018	1022		match_col <- (colnames(preds) %in% colnames)

1	coverage:
2	precision: 2
3	round: down
4	status:
5	patch:
6	default:
7	target: 50%
8	project:
9	default:
10	target: 50%
11
12	comment: false

Updating our web app