R/rMR.R

#'@export
#'@import graphics
plotRaw <-
    function(data, DO.var.name, time.var.name = "std.time",
             start.time = data$x[1],
             end.time = data$x[length(data$x)], ...){

        data$x <- eval(parse(text = paste("data$", time.var.name, sep = "")))
        data <- data[data$x >= start.time & data$x <= end.time,]
        data$y <- eval(parse(text = paste("data$", DO.var.name, sep = "")))

        plot(x=data$x, y=data$y, ...)
    }



#'@export
#'@import stats
sumsq <-
    function(x){
        mx<-mean(x)
        sq.dev<-((x-mx)^2)
        sumsquares<-sum(sq.dev)
        return(sumsquares)
    }


#'@export
#'@import stats
tot.rss <-
    function(data, break.pt, xvar, yvar){
        data$x <- eval(parse(text = paste("data$", xvar, sep="")))
        data$y <- eval(parse(text = paste("data$", yvar, sep="")))
        d1 <- data[data$x >= break.pt,]
        m1 <- lm(d1$y~d1$x)
        
        d2 <- data[data$x < break.pt,]
        m2 <- lm(d2$y~d2$x)
        
        trss <- (sumsq(m1$residuals)+sumsq(m2$residuals))
        return(trss)
    }

#'@export

#'@import stats
#'@import biglm
background.resp <-
    function(data, DO.var.name, time.var.name = "std.time",
             start.time, end.time, col.vec = c("black","red"), ...){
       
        orig <- "1970-01-01 00:00:00 UTC"
        
        data$y <- eval(parse(text = paste("data$", DO.var.name, sep = "")))
        data$x <- eval(parse(text = paste("data$", time.var.name, sep = "")))

        if(class(start.time) != class(end.time)) {
            stop ("start time and end time must be of same atomic class")
        }
        if(is.character(start.time)==T){
            data <- data[data$x >= as.POSIXct(start.time,
                                              origin=orig) &
                             data$x<= as.POSIXct(end.time,
                                                 origin=orig),]  
        }else if (class(start.time)[1]==("POSIXct") &
                      class(start.time)[2]==("POSIXt")){
        
            data <- data[data$x >= start.time &
                             data$x <= end.time] 
        }
        
        #calculate lm for background rate
        m1 <- biglm(y ~ x, data)
        
        MR1 <- coef(m1)[2]
        
        if (MR1>=0){warning("slope control 1 negative")}
        x<-data$x
        y<-data$y
        
        # Too keep point edges from getting obscured and hazy
        dots<-list(...)
        if(exists("dots$cex")==TRUE){
            cex.val<-0.7
        }else{cex.val<-dots$cex}
        
        if(exists("dots$pch")==TRUE){
            pch.val<-1
        }else{pch.val<-dots$pch}
        
        plot(x, y, type="n", ...)
        points(x=x, y=y, cex=cex.val, pch=pch.val, col=col.vec[1])
        abline(coefficients(m1), col = col.vec[2],...)
        
        return(summary(m1))
    }

#'@export
Barom.Press <-
    function(elevation.m, units = "atm"){
        if(units == "atm") {fact <-1}
        else if(units == "kpa") {
            fact <- 101.325
        }else if(units == "mmHg") {
            fact <- 760
        }else{
            stop("invalid pressure units, must be
             'atm', 'kpa', or 'mmHg'")
        }
        P<-exp(-9.80665 * 0.0289644 * elevation.m / 
                   (8.31447 * 288.15))*fact
        return(P)
    }


#'@export
#'
DO.saturation <-
    function(DO.mgl, temp.C, elevation.m = NULL,
             bar.press = NULL, bar.units = NULL,
             salinity = 0 , salinity.units = "pp.thou"){
        if(is.null(bar.press) == TRUE){
            DO.sat<- DO.mgl / Eq.Ox.conc(temp.C, elevation.m = elevation.m,
                                         bar.press = NULL, bar.units=NULL,
                                         salinity = salinity,
                                         salinity.units = salinity.units)   
        }else{
            DO.sat<- DO.mgl / Eq.Ox.conc(temp.C, bar.press = bar.press,
                                         bar.units = bar.units,
                                         salinity = salinity,
                                         salinity.units = salinity.units) 
        }
        
        return(DO.sat)
    }

#'@export
DO.unit.convert <-
    function(x, DO.units.in, DO.units.out, 
             bar.units.in, bar.press, temp.C,
             bar.units.out = "mmHg",
             salinity = 0 , salinity.units = "pp.thou"){
        if(bar.units.in == bar.units.out){
            bar.press <- bar.press
        }else if(bar.units.in == "atm"){
            if(bar.units.out == "mmHg"){
                bar.press <- bar.press * 760
            }else if(bar.units.out == "kpa"){
                bar.press <- bar.press * 101.32501
            }else{
                stop("invalid 'bar.units.out' -- must be 'atm', 'mmHg', 'kpa'")
            }
        }else if(bar.units.in == "mmHg"){
            if(bar.units.out == "atm"){
                bar.press <- bar.press / 760
            }else if(bar.units.out == "kpa"){
                bar.press <- bar.press * 101.32501 / 760
            }else{
                stop("invalid 'bar.units.out' -- must be 'atm', 'mmHg', 'kpa'")
            }
        }else if(bar.units.in == "kpa"){
            if(bar.units.out == "atm"){
                bar.press <- bar.press / 101.32501
            }
            if(bar.units.out == "mmHg"){
                bar.press <- bar.press * 760 / 101.32501
            }
        }
        
        
        if (DO.units.in == "pct"){
            DO.pct <- x / 100
        }else{
            eq.o2.in <- Eq.Ox.conc(temp.C, bar.units = bar.units.out, 
                                   bar.press =bar.press,
                                   out.DO.meas = DO.units.in,
                                   salinity = salinity,
                                   salinity.units = salinity.units)
            DO.pct <- x / eq.o2.in
        }
        
        if (DO.units.out == "pct"){
            DO.conc <- DO.pct * 100
        }else{
            eq.o2.out <- Eq.Ox.conc(temp.C, bar.units = bar.units.out,
                                    bar.press = bar.press,
                                    out.DO.meas = DO.units.out,
                                    salinity = salinity,
                                    salinity.units = salinity.units)
            DO.conc <- DO.pct * eq.o2.out
        }
        
        return(DO.conc)
    }

#'@export
Eq.Ox.conc <-
    function(temp.C, elevation.m = NULL,
             bar.press = NULL, bar.units = NULL,
             out.DO.meas = "mg/L", salinity = 0,
             salinity.units = "pp.thou"){
        tk <- 273.15 + temp.C
        if(is.null(elevation.m) == FALSE &&
           is.null(bar.units) == FALSE){
            stop("'bar.units' must be NULL if 'elevation.m' is assigned a value. ")
        }
        
        if( out.DO.meas == "PP"){
            
            if(is.null(bar.press) == FALSE && is.null(elevation.m) == TRUE){
                bar.press <- bar.press
            } else if(is.null(bar.press) == TRUE &&
                      is.null(elevation.m) == FALSE){
                bar.press <- Barom.Press (elevation.m,
                                          units = bar.units)
            }else{
                stop("EITHER 'elevation.m' or 'bar.press' must be assigned
                     a value. The other argument must be NULL.")
            }

                DO <- bar.press*0.20946
            
            
            }else if (out.DO.meas == "mg/L"){
                
                if(is.null(bar.press) == FALSE &&
                   is.null(elevation.m) == TRUE){
                    if(bar.units == "atm"){
                        bar.press <- bar.press
                    }else if(bar.units == "kpa"){
                        bar.press <- bar.press / 101.325
                    }else if(bar.units == "mmHg"){
                        bar.press <- bar.press / 760
                    }else{
                        stop("invalid pressure units, must be
                             'atm', 'kpa', or 'mmHg'")
                    }
                    } else if(is.null(bar.press) == TRUE &&
                              is.null(elevation.m) == FALSE){
                        bar.press <- Barom.Press (elevation.m, units = "atm")
                    }else{
                        stop("EITHER 'elevation.m' or 'barom.press' must be assigned
                             a value. The other argument must be NULL.")
                    }
                ## Benson and Krause eqns, USGS 2011 ##
                
                A1 <- -139.34411
                A2 <- 1.575701e5
                A3 <- 6.642308e7
                A4 <- 1.243800e10
                A5 <- 8.621949e11

                DO <- exp(A1 + (A2/tk) -
                                  (A3/(tk^2)) +
                                  (A4/(tk^3)) -
                                  (A5/(tk^4)))
                
            }else{
                stop("must specify 'out.DO.meas' as 'mg/L' or 'PP'")
            }
                # salinity factor #
                if(salinity.units == "uS"){
                    sal <- salinity*5.572e-4 + (salinity^2)*2.02e-9 
                }else if(salinity.units == "pp.thou"){
                    sal <- salinity
                }else{
                    stop("salinity.units must be 'uS' or 'pp.thou'")
                }
                
                Fs <- exp(-sal*(0.017674 - (10.754/tk) + (2140.7/(tk^2))))
                
                ## Pressure factor determination ##
                theta <- 0.000975 - 
                    temp.C*1.426e-5 +
                    (temp.C^2)*6.436e-8
                
                u <- exp(11.8571 - (3840.70/tk) - (216961/(tk^2)))
                
                # pressure factor #
                Fp <- ((bar.press - u)*(1-(theta*bar.press))) /
                      ((1-u)*(1-theta))
                              
                Cp <- DO * Fs * Fp

        return(Cp)
            }

#'@export

#'@import stats
#'@import graphics
#'@import biglm
get.pcrit <-
    function(data, DO.var.name, MR.var.name = NULL, Pcrit.below,
             time.interval, time.var = NULL,
             start.time, stop.time, time.units = "sec",
             Pcrit.type = "both",
             col.vec = c("black", "gray60", "red", "blue"),...){
        
        data$DO<- eval(parse(text = paste("data$", DO.var.name,
                                              sep = "")))
            
        
        ## set time denominator based on specified time.units ##
        
        if(is.null(time.var)==FALSE){
            if(time.units == "sec"){
                t.denom <- 1
            }else if(time.units == "min"){
                t.denom <- 60
            }else if(time.units == "hr"){
                t.denom <- 3600
            }
            data$time <- eval(parse(text =
                                   paste("data$", time.var, sep = "")))
            
            data <- data[data$time >= start.time
                         & data$time <= stop.time,]
            
            data$time <- as.numeric(data$time) - min(as.numeric(data$time)) 
        }

        
        
        if(any(is.na(data$DO)==TRUE)){
            warning("DO variable contains missing values")
        }
         
        
        if(is.null(MR.var.name) == TRUE){
            
            if(is.null(time.interval) == TRUE){
                stop("if using DO to calculate rates & Pcrit,
                     'time.interval' must be specified")
            }
            
            calc.MRs <- data.frame()
            
            i<-0
            
            while(i + time.interval < length(data$time)){
                the.interval <- data[data$time >= i &
                                         data$time < i + time.interval,]
                m <- lm(DO ~ time, data = the.interval)
                MR <- m$coefficients[2]*(-1)
                DO <- mean(the.interval$DO, na.rm = TRUE)
                time <- round(the.interval$time[1])
                df.row <- t(c(DO, MR))
                calc.MRs <- rbind(calc.MRs, df.row)
                i <- i + time.interval
                
            }
            
            
            data <- calc.MRs
            names(data) <- c("DO", "MR")
            data$MR <- data$MR * t.denom
            
        }else if(is.null(MR.var.name) == FALSE){
            data$MR <- eval(parse(text = paste("data$", MR.var.name,
                                               sep = "")))
        }
        
        
        
        
        minimum.DO <- 
            as.numeric(min(data$DO <= Pcrit.below, na.rm=TRUE))
        
        zone.of.interest <- data[data$DO <= Pcrit.below,]
        zone.of.interest <- zone.of.interest[order(zone.of.interest$DO, 
                                                   decreasing=TRUE),]
        
        RSS.table<-data.frame()
        for(w in 2:length(zone.of.interest[,1])-1){
            if(zone.of.interest$DO[w] != minimum.DO){
                break.point <- zone.of.interest$DO[w]
                RSS <- tot.rss(data = data,
                               break.pt = break.point,
                               xvar = "DO", yvar = "MR")
                RSS.row <- cbind(RSS, zone.of.interest$DO[w])
                RSS.table <- rbind(RSS.table,RSS.row)    
            }
        }
        
        DO.idx <- RSS.table[RSS.table$RSS == min(RSS.table$RSS), 2]
        idx.rname <- row.names(RSS.table[RSS.table[,2] == DO.idx,])
        DO.idx.low <- RSS.table[row.names(RSS.table) == 
                                    (as.numeric(idx.rname)-1), 2]

        midpoint.approx <- mean(c(DO.idx, DO.idx.low))
        
        
        ##creating linear models##
        dat1 <- data[data$DO > DO.idx,]
        mod.1 <- lm(MR~DO, data=dat1)
        
        dat2 <- data[data$DO <= DO.idx,]
        mod.2 <- lm(MR~DO, data=dat2)
        
        sm1 <- summary(mod.1)
        sm2 <- summary(mod.2)
        
        adjr2.pre <- sm1$adj.r.squared
        adjr2.post <- sm2$adj.r.squared
        
        ## Plot: line intersect##
        
        plot(MR~DO, data, type= "n", ...)
      
        intersect<-(mod.2$coefficients[1] - mod.1$coefficients[1]) /
            (mod.1$coefficients[2] - mod.2$coefficients[2])
        names(intersect)<-NULL
        if(is.na(mod.2$coefficients[2])==T){
            abline(mod.1$coefficients[1], 0)
        }else{abline(coef = mod.1$coefficients, col = col.vec[2], ...)}
        abline(coef = mod.2$coefficients, col = col.vec[2], ...)

        
        if (Pcrit.type == "lm" | Pcrit.type == "both"){
            abline(v = intersect, col = col.vec[3], lty=2, ...)            
        }
        if (Pcrit.type == "midpoint" | Pcrit.type == "both"){
            abline(v = midpoint.approx, col = col.vec[4], lty=3, ...)            
        }
       # Too keep point edges from getting obscured and hazy
        dots<-list(...)
        if(exists("dots$cex")==TRUE){
            cex.val<-0.7
        }else{cex.val<-dots$cex}
        
        if(exists("dots$pch")==TRUE){
            pch.val<-1
        }else{pch.val<-dots$pch}
            
        points(x = c(dat1$DO, dat2$DO), y = c(dat1$MR, dat2$MR), 
                   cex = cex.val, pch = pch.val, col = col.vec[1])

        
        dat.pre<-data[data$DO>=(2*intersect),]
        
        P.crit<-as.data.frame(cbind(intersect, midpoint.approx,
                                    adjr2.pre, adjr2.post))
        names(P.crit)<-c("Pcrit.lm", "Pcrit.midpoint",
                         "Adj.r2.above", "Adj.r2.below")

        above.Pc <- list(mod.1)
        names(above.Pc) <- "above"
        below.Pc <- list(mod.2)
        names(below.Pc) <- "below"
        
        Pc<-c(P.crit, above.Pc, below.Pc)
        return(Pc)
        
    }

#'@export

#'@import stats
#'@import graphics
#'@import biglm
MR.loops <-
    function(data, DO.var.name, time.var.name = "std.time",
             in.DO.meas = "mg/L", out.DO.meas = "mg/L",
             start.idx, stop.idx, system.vol = 1,
             background.consumption = 0,
             background.indices = NULL,
             temp.C, elevation.m = NULL,
             bar.press = NULL, bar.units = "atm",
             PP.units, time.units = "sec",
             col.vec = c("black", "red"),...){
        ## format the time vectors into POSIX ##
        orig = "1970-01-01 00:00:00 UTC"
        start.idx <- as.POSIXct((start.idx), origin = orig)
        stop.idx <- as.POSIXct((stop.idx), origin = orig)
        
        ## make sure num of start and stop indices agree ##
        if (length(start.idx) != length(stop.idx)){
            stop ("number of start times not equal
                  to number of stop times")
        }
        
        
        ## set time denominator based on specified time.units ##
        if(time.units == "sec"){
            t.denom <- 1
        }else if(time.units == "min"){
            t.denom <- 60
        }else if(time.units == "hr"){
            t.denom <- 3600
        }
        
        
        ## set response variable ##
        data$y <- eval(parse(text = paste("data$", DO.var.name, sep = "")))
        ## set time variable ##
        data$x <- eval(parse(text = paste("data$", time.var.name, sep = "")))
        
        ## set background DO consumption rate ##
        if(is.null(background.indices) == TRUE){
            bgd.slope.int <- background.consumption
            bgd.slope.slope <- 0
            
        }else if(is.null(background.indices) == FALSE) {
            if(length(background.indices) < 2){
                stop("background.indices must be NULL or have a length >= 2")
            }else{
                ## if there are multiple background calibrations ##
                bgd.mod <- lm(background.consumption ~ 
                                  as.POSIXct(background.indices))
                bgd.slope.int <- bgd.mod$coefficients[1]
                bgd.slope.slope <- bgd.mod$coefficients[2]
            }
            
        }
        
        if(is.null(bar.press) == FALSE &&
           is.null(elevation.m) == FALSE){
            stop("Either 'bar.press' or 'elevation.m' should be NULL")
        }
        
        ## barometric pressure ##
        if(is.null(bar.press) == FALSE){
            if(is.character(bar.press) == TRUE){
                bar.press <- eval(parse(
                    text = paste("data$",
                                 bar.press, sep = "")))
            }else if(is.numeric(bar.press) == TRUE){
                bar.press <- bar.press
            }else{
                stop("'bar.press' must be 'NULL', numeric, or
                     the col.name for barometric pressure")
            }
            }
        
        
        ## Temperature ##
        if (is.character(temp.C) == TRUE){
            temp.C <- eval(parse(
                text = paste("data$", temp.C, sep = "")))
        }else if(is.numeric(temp.C) == TRUE){
            temp.C <- temp.C
        }else{
            stop("invalid temp.C argument")
        }
        
        # DO sat conversions #
        if (in.DO.meas == "pct"){
            data$y <- (data$y /100) * 
                Eq.Ox.conc(temp.C = temp.C,
                           elevation.m = elevation.m,
                           bar.press = bar.press,
                           bar.units = bar.units,
                           out.DO.meas = out.DO.meas)
        }else if (in.DO.meas == "PP"){
            fraction <- data$y / 
                Eq.Ox.conc(temp.C = temp.C,
                           elevation.m = elevation.m,
                           bar.press = bar.press,
                           bar.units = bar.units,
                           out.DO.meas = "PP")
            data$y <- fraction *
                Eq.Ox.conc(temp.C = temp.C,
                           elevation.m = elevation.m,
                           bar.press = bar.press,
                           bar.units = bar.units,
                           out.DO.meas = out.DO.meas)
        }else if(in.DO.meas == "mg/L"){
            data$y <- data$y
        }else{
            stop("invalid 'in.DO.meas' argument:
                 must be 'pct', 'PP', or 'mg/L' ")
        }
        
        if(out.DO.meas == "mg/L"){
            data$y <- system.vol * data$y
            # Now data$y in units of mg, not mg/L #
        }else if(out.DO.meas == "PP"){
            # converting MR to PP units #
            data$y <- DO.unit.convert(x = data$y,
                                      DO.units.in = "mg/L",
                                      DO.units.out = "PP",
                                      bar.units.in = "atm",
                                      bar.units.out = PP.units)
        }else if(out.DO.meas == "pct"){
            data$y <- fraction * 100
        }
        
        ## the following loop adds adj.y as a variable.     ##
        ## it adjusts DO level by iincorporating background ##
        ## respiration rate.                                ##
        
        data$adj.y <- rep(0, length(data[,1]))
        tab.wid <- length(data[1,])
        data.new <- NULL
        
        for(i in 1:length(start.idx)){
            
            dsn <- data[data$x >= as.POSIXct(start.idx[i]) &
                            data$x <= as.POSIXct(stop.idx[i]), ]
            dsn$adj.y <- (dsn$y- (as.numeric(dsn$x) -
                                      as.numeric(start.idx[i]))*
                                     ((as.numeric(dsn$x) * bgd.slope.slope) +
                                          bgd.slope.int))
            
            data.new <- rbind(data.new, dsn)
        }
        
        data <- data.new[(data.new$x >= start.idx[1] - 600) &
                             data.new$x <= (tail(stop.idx,1) + 600),]
        x<-data$x
        y<-data$y
        
        plot(x, y, type="n",...)
        
        name.num<-as.character(c(1:length(start.idx)))
        ms<-list()
        MR.summary<-data.frame()
        
        dots<-list(...)
        if(exists("dots$cex")==TRUE){
            cex.val<-0.7
        }else{cex.val<-dots$cex}
        
        if(exists("dots$pch")==TRUE){
            pch.val<-1
        }else{pch.val<-dots$pch}
        
        for(i in 1:length(start.idx)){
            dat <- data.new[data.new$x >= start.idx[i]
                            & data.new$x <= stop.idx[i],]
            
            
            mk <- biglm(adj.y ~ x, data=dat)
            ms[[i]] <- mk
            
            points(dat$x, dat$adj.y, col = col.vec[1],
                   pch = pch.val, cex = cex.val)
            names(ms[[i]])<-paste(names(ms[[i]]), name.num[i], sep=".")
            abline(coef(ms[[i]]),
                   col = col.vec[2], ...)
            
            MR <- coef(mk)[2]*-1
            sds <- summary(mk)$mat[2,4]*sqrt(length(dat[,1]))
            rsquare <- summary(mk)$rsq
            mrrow <- t(c(MR, sds, rsquare))
            MR.summary <- rbind(MR.summary,mrrow)
        }
        
        names(MR.summary) <- c("MR", "sd.slope", "r.square")
        MR.summary[,c(1,2)] <- MR.summary[,c(1,2)] * t.denom
        
        
        ofthejedi <- list(MR.summary, ms)
        names(ofthejedi) <- c("MR.summary", "lm.details")
        return(ofthejedi)
        }

#'@export
#'@import utils
get.witrox.data <-
    function(data.name, lines.skip, delimit="tab", choose.names = F,
             chosen.names = NULL,
             format){
        if(delimit == "tab"){
            separate = "\t"
        }else if(delimit == "space"){
            separate = ""
        }else if(delimit == "comma"){
            separate = ","
        }
        
        if(choose.names==F){
            d<-read.table(data.name, sep=separate, skip=lines.skip,
                          header =T,check.names=F)
            invalid.names<-colnames(d)
            valid.names<-make.names(colnames(d))
            var.names<-NULL
            for(i in 1:length(invalid.names)){
                if(invalid.names[i] == valid.names[i]){
                    var.names[i] <- valid.names[i]
                }else{
                    split.name.period <- as.vector(
                        strsplit(valid.names[i], fixed = T,
                                 split = ".")[[1]])
                    split.name.period <- paste(split.name.period, sep="",
                                               collapse="")

                }
            }
        }else if(choose.names==T){
            d<-read.table(data.name, sep = separate, skip = lines.skip,
                          header = F)
            var.names<-chosen.names
        }
        
        colnames(d) <- var.names
        d[,1]<-as.character(d[,1]) 
        d$std.time <- as.POSIXct(strptime(d[,1], format = format),
                                 origin = "1970-01-01 00:00:00 UTC")
        return(d)
        
    }
tyler-l-moulton/rMR documentation built on May 3, 2019, 1:53 p.m.