R/growth_model.r

In ShellChron: Builds Chronologies from Oxygen Isotope Profiles in Shells

#' Function that models a d18O curve through SST and GR sinusoids
#' 
#' The core function of the ShellChron growth model. Uses growth
#' rate and SST (Sea Surface Temperature) sinusoids to model d18O
#' data to be matched with the input. In the ShellChron modeling 
#' routine, this function is optimized using the SCEUA algorithm
#' and applied on sliding windows through the dataset to estimate
#' the age of each datapoint
#' @param pars List of parameters for temperature and growth rate sinusoids
#' \code{pars <- c(T_amp, T_pha, T_av, G_amp, G_pha, G_av, G_skw)}
#' @param T_per Period of SST sinusoid (in days; default = 365)
#' @param G_per Period of growth rate sinusoid (in days; default = 365)
#' @param years Number of years to be modeled (default = 1)
#' @param t_int Time interval (in days; default = 1)
#' @param transfer_function Transfer function used to convert d18Oc to temperature
#' data.
#' @param d18Ow Either a single value (constant d18Ow) or a vector of length
#' equal to the period in SST data (365 days by default) containing information
#' about seasonality in d18Ow. Defaults to constant d18Ow of 0 permille VSMOW
#' (the modern mean ocean value)
#' @param Dsam Vector of \code{D} values serving as input (keep unit consistent
#' throughout model)
#' @param Osam Vector of \code{d18Oc} values serving as input (in permille VPDB)
#' @param t_maxtemp Timing of the warmest day of the year (in julian day; 
#' default = 182.5, or May 26th halfway through the year)
#' @param plot Should results of modeling be plotted? \code{TRUE/FALSE}
#' @param MC Number of Monte Carlo simulations to apply for error propagation
#' Default = 1000
#' @param D_err OPTIONAL: Vector containing errors on \code{Dsam}
#' @param O_err OPTIONAL: Vector containing errors on \code{Osam}
#' @param return String indicating whether to return just the Sum of Squared
#' Residuals ("SSR") or a matrix containing the results of the model and the
#' propagated uncertainties (if applicable)
#' @return Depending on the value of the "return" parameter either a single
#' value representing the Sum of Squared Residuals ("SSR") as a measure for
#' the closeness of the match between modeled d18O and input values, or a
#' matrix containing the full result of the modeling including propagated
#' uncertainties if applicable.
#' @references package dependencies: ggplot2 3.2.1
#' function dependencies: temperature_curve, d18O_model, growth_rate_curve,
#' mc_err_orth
#' 
#'   \doi{10.1016/j.palaeo.2017.09.034}
#' @examples
#' # Set parameters
#' G_amp <- 20
#' G_per <- 365
#' G_pha <- 100
#' G_av <- 15
#' G_skw <- 70
#' T_amp <- 20
#' T_per <- 365
#' T_pha <- 150
#' T_av <- 15
#' pars <- c(T_amp, T_pha, T_av, G_amp, G_pha, G_av, G_skw)
#' d18Ow <- 0
#' # Create dummy data
#' Dsam <- seq(1, 40, 1)
#' Osam <- sin((2 * pi * (seq(1, 40, 1) - 8 + 30 / 4)) / 30)
#' # Test returning residual sum of squares for optimization
#' SSR <- growth_model(pars, T_per, G_per, Dsam = Dsam, Osam = Osam,
#'     return = "SSR")
#' # Test returning full model result
#' resmat <- growth_model(pars, T_per, G_per, Dsam = Dsam, Osam = Osam,
#'     return = "result")
#' @export
growth_model <- function(pars, # Growth model function to optimize using sceua
    T_per = 365, # Period of SST sinusoid (default = 365 days)
    G_per = 365, # Period of growth sinusoid (default = 365 days)
    years = 1, # Default year = 1
    t_int = 1, # Time interval, default = 1 day
    transfer_function = "KimONeil97", # Transfer function used to convert d18Oc values of the material into temperature (default = "KimONeil97"; Kim and O'Neil, 1997)
    d18Ow = "default", # d18Ow vector, default = constant at 0 permille VSMOW
    Dsam, # Depth data
    Osam, # d18O data
    t_maxtemp = 182.5, # Time (day) at which maximum temperature is reached. Default is 182.5 (exactly halfway through the year, or 1st of June)
    plot = FALSE, # Plot results?
    MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (error on D and d18O measurements not considered)
    D_err = NULL, # Optional: include standard deviation of uncertainty on depth measurements, default = empty (no error incorporated in the calculation, same effect as MC = 0)
    O_err = NULL, # Optional: include standard deviation of uncertainty on d18O measurements, default = empty (no error incorporated in the calculation, same effect as MC = 0)
    return = "SSR" # What to return? Default = the total sum of squares to be minimized by the SCEUA algorithm
    ){

    # Extract variables
    T_amp <- pars[1]
    T_pha = pars[2]
    T_av = pars[3]

    G_amp = pars[4]
    G_pha = pars[5]
    G_av = pars[6]
    G_skw = pars[7]

    T_par <- c(T_amp, T_per, T_pha, T_av)
    G_par <- c(G_amp, G_per, G_pha, G_av, G_skw)

    # Run functions for SST and GR
    SST <- temperature_curve(T_par, years, t_int)
    d18Oc <- d18O_model(SST, d18Ow, transfer_function)
    if(length(d18Oc) == 1){ # Catch errors in d18O model
        return("ERROR: Supplied transfer function is not recognized")
    }
    GR <- growth_rate_curve(G_par, years, t_int)

    # Create distance (D) vector from GR
    D <- cumsum(GR[,2]) + Dsam[1]

    # Compare model result with data
    Omod <- approx( # Interpolate modeled d18Oc values to positions along the record.
        x = D,
        y = d18Oc[,2],
        xout = Dsam,
        method = "linear",
        rule = 2
    )

    residuals <- Osam - Omod$y # Calculate residuals

    if(return == "SSR"){
        SSR <- sum(residuals^2) # Calculate sum of squared residuals
        return(SSR)
    }else{
        # Compare model result with data
        t <- approx( # Interpolate modeled time values to positions along the record.
            x = D,
            y = d18Oc[,1],
            xout = Dsam,
            method = "linear",
            rule = 2
        )

        gr <- approx( # Interpolate modeled growth rate values to positions along the record.
            x = D,
            y = GR[,2],
            xout = Dsam,
            method = "linear",
            rule = 2
        )

        Tmod <- approx( # Interpolate modeled temperature values to positions along the record.
            x = D,
            y = SST[,2],
            xout = Dsam,
            method = "linear",
            rule = 2
        )

        TY <- (t$y - T_pha + t_maxtemp) %% 365 # Calculate time of year relative to the first T maximum in the record (T_pha) and the assumed day of the year when it occurs (default = 182.5)

        # Prepare export matrix
        resmat <- cbind(Dsam, Osam, Omod$y, residuals, TY, gr$y, Tmod$y) # Calculate matrix of end results
        colnames(resmat) <- c("Dsam", "Osam", "Omod", "Residuals", "Time_of_year", "Instantaneous_growth_rate", "Modeled_temperature")
        
        # Optional: Incorporate measurement errors on D and d18O
        if(MC > 0){
            print("Propagating measurement uncertainties")
            if(is.null(D_err)){ # If no error on D is given, but error analysis is asked (MC > 0), then D_err is zero for all samples
                D_err <- rep(0, length(Dsam))
            }
            if(is.null(O_err)){ # If no error on d18Oc is given, but error analysis is asked (MC > 0), then O_err is zero for all samples
                O_err = rep(0, length(Osam))
            }
            # Propagate the combined effect of error on D and d18Oc on the model fit
            Proj_err <- mc_err_orth(Dsam, D_err, Osam, O_err, D, d18Oc, MC) # Combine error on D and d18Oc through orthogonal projection on the modeled D-d18Oc curve
            D_err_comb <- Proj_err$X_err # Isolate uncertainty in X-domain for further processing

            Drange <- cbind((Dsam - D_err_comb) %% D[length(D)], (Dsam + D_err_comb) %% D[length(D)]) # Find upper and lower boundaries of D error (1 SD)
            Drange_pos <- apply(abs(outer(Drange, D, FUN = "-")), c(1, 2), which.min) # Find positions of D ranges

            # Approximate SDs of d18Oc, t, gr and Tmod from range of D values within +/- 1 SD
            Omod_SD <- vector(length = length(Drange_pos[,1]))
            t_SD <- vector(length = length(Drange_pos[,1]))
            gr_SD <- vector(length = length(Drange_pos[,1]))
            Tmod_SD <- vector(length = length(Drange_pos[,1]))
            for(j in 1:length(Drange_pos[,1])){ # Loop through all positions in the window and propagate the D error onto the results
                if(Drange_pos[j, 2] < Drange_pos[j, 1]){ # If the D range contains the year edge, incorporate the outer edges of the year.
                # Find the range of modeled values contained in the 1 SD range of D error and approximate SD of modeled values from in and max values
                    Omod_SD[j] <- (max(d18Oc[c(Drange_pos[j, 1]:length(d18Oc[, 1]), 1:Drange_pos[j, 2]), 2]) - min(d18Oc[c(Drange_pos[j, 1]:length(d18Oc[, 1]), 1:Drange_pos[j, 2]), 2])) / 2
                    t_SD[j] <- length(c(Drange_pos[j, 1]:length(d18Oc[, 1]), 1:Drange_pos[j, 2])) / 2
                    gr_SD[j] <- (max(GR[c(Drange_pos[j, 1]:length(GR[, 1]), 1:Drange_pos[j, 2]), 2]) - min(GR[c(Drange_pos[j, 1]:length(GR[, 1]), 1:Drange_pos[j, 2]), 2])) / 2
                    Tmod_SD[j] <- (max(SST[c(Drange_pos[j, 1]:length(SST[, 1]), 1:Drange_pos[j, 2]), 2]) - min(SST[c(Drange_pos[j, 1]:length(SST[, 1]), 1:Drange_pos[j, 2]), 2])) / 2
                }else{ # If the D range is contained within one year, range of d18Oc values is located between the lower and upper boundaries of D (+/- 1 SD)
                    Omod_SD[j] <- (max(d18Oc[Drange_pos[j, 1]:Drange_pos[j, 2], 2]) - min(d18Oc[Drange_pos[j, 1]:Drange_pos[j, 2], 2])) / 2
                    t_SD[j] <- length(Drange_pos[j, 1]:Drange_pos[j, 2]) / 2
                    gr_SD[j] <- (max(GR[Drange_pos[j, 1]:Drange_pos[j, 2], 2]) - min(GR[Drange_pos[j, 1]:Drange_pos[j, 2], 2])) / 2
                    Tmod_SD[j] <- (max(SST[Drange_pos[j, 1]:Drange_pos[j, 2], 2]) - min(SST[Drange_pos[j, 1]:Drange_pos[j, 2], 2])) / 2
                }
            }
            resmat <- cbind(resmat, Omod_SD, t_SD, gr_SD, Tmod_SD)
        }else{
            resmat <- cbind(resmat, matrix(rep(NA, 4 * length(resmat[,1])), ncol = 4)) # Add NA-filled columns if no MC error propagation is done
        }
        colnames(resmat) <- c("Dsam",
            "Osam",
            "Omod",
            "Residuals",
            "Time_of_year",
            "Instantaneous_growth_rate",
            "Modeled_temperature",
            "Omod_SD",
            "Time_of_Year_SD",
            "Instantaneous_growth_rate_SD",
            "Modeled_temperature_SD")

        if(plot == TRUE){
            plotdf <- as.data.frame(cbind(resmat, D_err, O_err))
            dev.new()
            ggplot2::ggplot(plotdf, ggplot2::aes(x = Dsam, y = Osam)) +
                ggplot2::geom_point() + # Plot original data
                ggplot2::geom_errorbar(ggplot2::aes(ymin = Osam - O_err, # Add error bars on measurement (1 SD)
                    ymax = Osam + O_err),
                    width = 100,
                    col = "black") +
                ggplot2::geom_errorbarh(ggplot2::aes(xmin = Dsam - D_err,
                    xmax = Dsam + D_err),
                    height = 0.05,
                    col = "black") +
                ggplot2::geom_point(ggplot2::aes(Dsam, Omod), col="red") + # Plot modeled d18Oc on top of data
                ggplot2::geom_errorbar(ggplot2::aes(ymin = Omod - Omod_SD, # Add error bars on model result (1 SD)
                    ymax = Omod + Omod_SD),
                    width = 100,
                    col = "red") +
                ggplot2::geom_errorbarh(ggplot2::aes(Dsam,
                    Omod,
                    xmin = Dsam - D_err,
                    xmax = Dsam + D_err),
                    height = 0.05,
                    col = "red")
        }
        return(resmat)
    }
}