From bdb684cace04c12788149a044676c440804e4455 Mon Sep 17 00:00:00 2001 From: JamesStagge Date: Tue, 12 Sep 2017 11:40:35 -0600 Subject: [PATCH] Fixed fitting and objects for MF models --- NAMESPACE | 1 + R/model_fit.R | 131 +++++++++++++++++++++++++++++++++++++++++----- R/paleo_objects.R | 9 ++-- man/mf_fit.Rd | 21 ++++++++ man/paleo.fit.Rd | 2 +- 5 files changed, 147 insertions(+), 17 deletions(-) create mode 100644 man/mf_fit.Rd diff --git a/NAMESPACE b/NAMESPACE index 2e0fbaa..4fbdd77 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -11,6 +11,7 @@ export(fit_model) export(fit_norm_dist) export(flow_reconstr) export(flow_to_norm) +export(mf_fit) export(norm_to_flow) export(paleo.fit) export(paleo.norm) diff --git a/R/model_fit.R b/R/model_fit.R index 9a7fec8..8b0a68c 100755 --- a/R/model_fit.R +++ b/R/model_fit.R @@ -13,21 +13,13 @@ ### Generic wrapper for fitting model fit_model <- function(method, regmethod="lm", reconst_data, annual_norm=NULL, monthly_norm=NULL, pred_ts=NULL, reg_eq=NULL, monthly_obs=NULL,...){ - if(method == "mf"){ - - } - else if (method == "ap" | method == "apr") { - ### Check that annual reconstructed timeseries and normalization are from same data - expect_that(annual_norm$prefix, is_identical_to(paste0(reconst_data$site_prefix, "_annual"))) - - ### Process object + ################################################ + ### Process Reconstructed time series + ################################################# flow_ts <- reconst_data$ts ### Test if the reconstruction uses calendar years are_cal_years <- is.null(reconst_data$wy_first_month) - ### Calculate normalized annual flow - flow_ts$annual_norm <- flow_to_norm(flow_series=reconst_data, dist_object=annual_norm) - ### Create a monthly time series of all combinations of months and year ### For calendar years, need the range of years and 12 months if (are_cal_years){ @@ -40,9 +32,47 @@ fit_model <- function(method, regmethod="lm", reconst_data, annual_norm=NULL, mo monthly_ts <- expand.grid(month = seq(1,12), year = seq(year_range[1], year_range[2])) monthly_ts$water_year <- usgs_wateryear(monthly_ts$year, monthly_ts$month, reconst_data$wy_first_month) } - + ### Set an order column to allow easier resorting at the end - monthly_ts$t <- seq(1,dim(monthly_ts)[1]) + monthly_ts$t <- seq(1,dim(monthly_ts)[1]) + + ################################################ + ### Begin to fit model + ################################################# + if(method == "mf"){ + + ### Merge annual flows and month + monthly_ts <- merge(monthly_ts, flow_ts, all.x=TRUE) + ### Rename flow column to annual recon + names(monthly_ts)[which(names(monthly_ts) == "flow")] <- "annual_recon" + + ### Re-sort to obtain time series + monthly_ts <- monthly_ts[with(monthly_ts, order(t)), ] + + ### Process observed flow + obs_ts <- monthly_obs$ts + + ### Test if the reconstruction uses calendar years + are_cal_years <- is.null(reconst_data$wy_first_month) + + ### Run the MF model + if(are_cal_years == TRUE) { + mf_prop <- mf_fit(obs_ts=obs_ts, wy_first_month=1) + } else { + mf_prop <- mf_fit(obs_ts=obs_ts, wy_first_month=reconst_data$wy_first_month) + } + + ### Prepare to return fit object + x <- paleo.fit(method=method, reconst_data=reconst_data, mf_prop=mf_prop) + x$reconst_data$ts <- monthly_ts + x$reconst_data$time_scale <- "monthly" + } + else if (method == "ap" | method == "apr") { + ### Check that annual reconstructed timeseries and normalization are from same data + expect_that(annual_norm$prefix, is_identical_to(paste0(reconst_data$site_prefix, "_annual"))) + + ### Calculate normalized annual flow + flow_ts$annual_norm <- flow_to_norm(flow_series=reconst_data, dist_object=annual_norm) ### Merge annual flows and month monthly_ts <- merge(monthly_ts, flow_ts, all.x=TRUE) @@ -121,8 +151,20 @@ flow_reconstr <- function(recon_model, post_proc=FALSE, interpolate=FALSE){ site_prefix <- recon_model$reconst_data$site_prefix if(method == "mf"){ + recon_ts <- recon_model$reconst_data + + ### Merge with MF proportion + mf_ts <- merge(recon_ts$ts, recon_model$mf_prop, by="month", all.x=TRUE) + ### Multiply monthly fraction by annual reconstruction * 12 + mf_ts$flow_est <- mf_ts$prop_mean * mf_ts$annual_recon*12 + ### Re-sort to obtain time series + mf_ts <- mf_ts[with(mf_ts, order(t)), ] + recon_ts$ts <- mf_ts + ### Create object to return + monthly_rec_ts <- paleo.ts(ts=recon_ts, time_scale="monthly", site_prefix=site_prefix, wy_first_month=recon_model$first_month_wy, method=method) + } else if (method == "ap" | method == "apr") { if (method == "ap") { @@ -173,6 +215,69 @@ return(monthly_rec_ts) + +#' Generic No documentation +#' +#' This function needs documentation. +#' +#' @param data A dataframe +#' @param write_folder Folder in which to save results +#' @param write_file Name of file to be saved +#' +#' @return p Plot +#' +#' +#' @export +mf_fit <- function(obs_ts, wy_first_month=1) { +require(data.table) + +### Add a water year column +obs_ts$water_year <- usgs_wateryear(obs_ts$year, obs_ts$month, wy_first_month) + +################################################ +### Calculate Annual Flow +################################################# +obs_ts <- data.table(obs_ts) + +### Calculate flow by water year +flow_obs_annual <- obs_ts[,list(annual_sum=sum(flow), annual_mean=mean(flow)),by="water_year"] +### Merge back with observed flow +obs_ts <- merge(obs_ts,flow_obs_annual, by=c("water_year"), all.x = T) + + +################################################ +### Calculate Monthly Fraction +################################################# +mf_ts <- obs_ts + +### Only consider those water years with 12 months of data, so first sum months with flows greater than 0. Negative flows will confuse the proportion method. +months_by_wy <- mf_ts[,list(count_months=sum(flow>0, na.rm=TRUE)), by=water_year] + +### Determine full water years and subset null_dt to only full years +full_wys <- months_by_wy$water_year[months_by_wy$count_months==12] +mf_ts <- mf_ts[mf_ts$water_year %in% full_wys,] + +### Calculate monthly proportion +mf_ts$prop_month <- mf_ts$flow/mf_ts$annual_sum + +### Summarize null_dt by calculating the mean and median proportion across all years +null_prop <- mf_ts[,list(prop_mean=mean(prop_month, na.rm=TRUE)), by=month] + +### Rescale so that months always add up to 1 +null_prop$prop_mean <- null_prop$prop_mean/sum(null_prop$prop_mean) +null_prop <- data.frame(null_prop) + +### Re-sort by months +null_prop <- null_prop[order(null_prop$month),] + +### Return monthly proportion +return(null_prop) + +} + + + + #' Generic No documentation diff --git a/R/paleo_objects.R b/R/paleo_objects.R index bb1ef38..3960b38 100755 --- a/R/paleo_objects.R +++ b/R/paleo_objects.R @@ -127,11 +127,11 @@ plot.paleo.norm <- function(x, ...) { #' #' #' @export -paleo.fit <- function(method, reconst_data, annual_norm=NULL, monthly_norm=NULL, first_month_wy=10, reg_eq=NULL, pred_ts=NULL){ +paleo.fit <- function(method, reconst_data, annual_norm=NULL, monthly_norm=NULL, first_month_wy=10, reg_eq=NULL, pred_ts=NULL, mf_prop=NULL){ x <- list(method=method, reconst_data=reconst_data, first_month_wy=first_month_wy) if(method == "mf"){ - model_list <- list() + model_list <- list(mf_prop = mf_prop) } else if (method == "ap") { model_list <- list(annual_norm=annual_norm, monthly_norm=monthly_norm) } else if (method == "apr") { @@ -158,7 +158,10 @@ paleo.fit <- function(method, reconst_data, annual_norm=NULL, monthly_norm=NULL, #' #' @export coef.paleo.fit <- function(fit, ...) { - if (fit$method == "ap") { + if (fit$method == "mf") { + coef_mat <- matrix(fit$mf_prop$prop_mean,1,12) + rownames(coef_mat) <- "mf_fraction" + } else if (fit$method == "ap") { coef_mat <- matrix(rep(1,12),1, 12) rownames(coef_mat) <- "annual_norm" } else if (fit$method == "apr") { diff --git a/man/mf_fit.Rd b/man/mf_fit.Rd new file mode 100644 index 0000000..d52611d --- /dev/null +++ b/man/mf_fit.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/model_fit.R +\name{mf_fit} +\alias{mf_fit} +\title{Generic No documentation} +\usage{ +mf_fit(obs_ts, wy_first_month = 1) +} +\arguments{ +\item{data}{A dataframe} + +\item{write_folder}{Folder in which to save results} + +\item{write_file}{Name of file to be saved} +} +\value{ +p Plot +} +\description{ +This function needs documentation. +} diff --git a/man/paleo.fit.Rd b/man/paleo.fit.Rd index c152442..fd8888b 100644 --- a/man/paleo.fit.Rd +++ b/man/paleo.fit.Rd @@ -5,7 +5,7 @@ \title{Generic No documentation} \usage{ paleo.fit(method, reconst_data, annual_norm = NULL, monthly_norm = NULL, - first_month_wy = 10, reg_eq = NULL, pred_ts = NULL) + first_month_wy = 10, reg_eq = NULL, pred_ts = NULL, mf_prop = NULL) } \arguments{ \item{data}{A dataframe}