# R/curveFunctions.R In bdots: Bootstrapped Differences of Time Series

#### Documented in logistic

```#' Logistic curve function for nlme
#'
#' Logistic function used in fitting nlme curve for observations
#'
#' @param dat subject data to be used
#' @param y outcome variable
#' @param time time variable
#' @param params \code{NULL} unless user wants to specify starting parameters for gnls
#' @param ... just in case
#'
#' @details \code{y ~ mini + (peak - mini) / (1 + exp(4 * slope * (cross - (time)) / (peak - mini)))}
#' @export
logistic <- function(dat, y, time, params = NULL, ...) {

logisticPars <- function(dat, y, time, ...) {
time <- dat[[time]]
y <- dat[[y]]

## Remove cases with zero variance
if (var(y) == 0) {
return(NULL)
}

mini <- min(y)
peak <- max(y)
r <- (peak - mini)
cross <- time[which.min(abs(0.5*r - y))]

# slope
q75 <- .75 * r + mini
q25 <- .25 * r + mini
time75 <- time[which.min(abs(q75 - y))]
time25 <- time[which.min(abs(q25 - y))]
slope <- (q75 - q25) / (time75 - time25)

return(c(mini = mini, peak = peak, slope = slope, cross = cross))
}

if (is.null(params)) {
params <- logisticPars(dat, y, time)
} else {
if (length(params) != 4) stop("logistic requires 4 parameters be specified for refitting")
if (!all(names(params) %in% c("mini", "peak", "slope", "cross"))) {
stop("logistic parameters for refitting must be correctly labeled")
}
}
## Return NA list if var(y) is 0
if (is.null(params)) {
return(NULL)
}
y <- str2lang(y)
time <- str2lang(time)
ff <- bquote(.(y) ~ mini + (peak - mini) / (1 + exp(4 * slope * (cross - (.(time))) / (peak - mini))))
attr(ff, "parnames") <- names(params)
return(list(formula = ff, params = params))
}

## -----------

#' Double Gauss curve function for nlme
#'
#' Double Gauss function used in fitting nlme curve for observations
#'
#' @param dat subject data to be used
#' @param y outcome variable, character vector
#' @param time time variable, character vector
#' @param concave Boolean
#' @param params \code{NULL} unless user wants to specify starting parameters for gnls
#' @param ... just in case
#'
#' @details User should only have to worry about setting concavity
#' of this function
#'
#' \code{y ~ (time < mu) * (exp(-1 * (time - mu) ^ 2
#' / (2 * sig1 ^ 2)) * (ht - base1) + base1)
#' + (mu <= time) * (exp(-1 * (time - mu) ^ 2
#'                          / (2 * sig2 ^ 2)) * (ht - base2) + base2)}
#' @export
doubleGauss <- function(dat, y, time, params = NULL, concave = TRUE, ...) {

if (is.null(params)) {
params <- dgaussPars(dat, y, time, concave)
} else {
if (length(params) != 6) stop("doubleGauss requires 6 parameters be specified for refitting")
if (!all(names(params) %in% c("mu", "ht", "sig1", "sig2", "base1", "base2"))) {
stop("doubleGauss parameters for refitting must be correctly labeled")
}
}
## Return NA list if var(y) is 0
if (is.null(params)) {
return(NULL)
}

y <- str2lang(y)
time <- str2lang(time)
ff <- bquote(.(y) ~ (.(time) < mu) * (exp(-1 * (.(time) - mu) ^ 2
/ (2 * sig1 ^ 2)) * (ht - base1) + base1)
+ (mu <= .(time)) * (exp(-1 * (.(time) - mu) ^ 2
/ (2 * sig2 ^ 2)) * (ht - base2) + base2))
attr(ff, "parnames") <- names(params)
return(list(formula = ff, params = params))
}

## Take this out of doubleGauss function proper
dgaussPars <- function(dat, y, time, conc) {
time <- dat[[time]]
y <- dat[[y]]

## Remove cases with zero variance
if (var(y) == 0) {
return(NULL)
}

mu <- ifelse(conc, time[which.max(y)], time[which.min(y)])
ht <- ifelse(conc, max(y), min(y))
base1 <- ifelse(conc, min(y[time < mu]), max(y[time < mu]))
base2 <- ifelse(conc, min(y[time > mu]), max(y[time > mu]))

## A little more involved
y1 <- y - base1
y1 <- rev(y1[time <= mu])
time1 <- rev(time[time <= mu])
totalY1 <- sum(y1)
sigma1 <- mu - time1[which.min(abs((pnorm(1) - pnorm(-1)) * totalY1 - cumsum(y1)))]

y2 <- y - base2
y2 <- y2[time >= mu]
time2 <- time[time >= mu]
totalY2 <- sum(y2)
sigma2 <- time2[which.min(abs((pnorm(1) - pnorm(-1)) * totalY2 - cumsum(y2)))] - mu

return(c(mu = mu, ht = ht, sig1 = sigma1, sig2 = sigma2,
base1 = base1, base2 = base2))
}

## -----------

#' Polynomial curve function for nlme
#'
#' Polynomial function used in fitting nlme curve for observations
#'
#' @param dat subject data to be used
#' @param y outcome variable
#' @param time time variable
#' @param degree degree of polynomial
#' @param raw Boolean, use raw polynomials?
#' @param params \code{NULL} unless user wants to specify starting parameters for gnls
#' @param ... just in case
#'
#' @details It's recommended that one uses raw polynomials for this function for
#' numerical stability. As inference is not performed on the parameters themselves,
#' this should have minimial consequences
#'
#' @details \code{y ~ mini + (peak - mini) / (1 + exp(4 * slope * (cross - (time)) / (peak - mini)))}
#' @export
polynomial <- function(dat, y, time, degree, raw = TRUE, params = NULL, ...) {

polyPars <- function(dat, y, time, degree, raw) {
## Remove cases with zero variance
if (var(dat[[y]]) == 0) {
return(NULL)
}
pp <- lm(dat[[y]] ~ poly(dat[[time]], degree = degree, raw = raw))
setNames(coef(pp), c(paste0("beta", seq(degree + 1L))))
}

if (is.null(params)) {
params <- polyPars(dat, y, time, degree, raw)
} else {
if (length(params) != degree + 1L) stop("poly requires that number of params matches degree + 1 (for intercept)")
if (!all(names(params) %in% paste0("beta", seq(degree + 1L)))) {
stop("polynomial parameters must be labeled beta1, ..., beta[degree + 1]")
}
}
## Return NA list if var(y) is 0
if (is.null(params)) {
return(NULL)
}

time_names <- paste0("I(Time^", seq(degree + 1L) - 1L , ")")

ff <- paste(names(params), time_names, sep = "*", collapse = "+")
ff <- str2lang(ff)
y <- str2lang(y)
ff <- bquote(.(y) ~ .(ff))
attr(ff, "parnames") <- names(params)
return(list(formula = ff, params = params))
}

#' Linear curve function
#'
#' Linear function used in fitting nlme curve for observations
#'
#' @param dat subject data to be used
#' @param y outcome variable
#' @param time time variable
#' @param params \code{NULL} unless user wants to specify starting parameters for gnls
#' @param ... just in case
#'
#' @details Don't use this function please
#'
#' @details \code{y ~ slope*time + intercept}
#' @export
linear <- function(dat, y, time, params = NULL, ...) {
linearPars <- function(dat, y, time, ...) {
time <- dat[[time]]
y <- dat[[y]]

## Remove cases with zero variance
if (var(y) == 0) {
return(NULL)
}
bb <- min(y)
mm <- (max(y) - bb) / max(time)

return(c(intercept = bb, slope = mm))
}

if (is.null(params)) {
params <- linearPars(dat, y, time)
} else {
if (length(params) != 2) stop("linear requires 2 parameters be specified for refitting")
if (!all(names(params) %in% c("intercept", "slope"))) {
stop("linear parameters for refitting must be correctly labeled")
}
}
## Return NA list if var(y) is 0
if (is.null(params)) {
return(NULL)
}
y <- str2lang(y)
time <- str2lang(time)
ff <- bquote(.(y) ~ slope * .(time) + intercept)
attr(ff, "parnames") <- names(params)
return(list(formula = ff, params = params))
}

#' Exponential curve function
#'
#' Exponential function used in fitting nlme curve for observations
#'
#' @param dat subject data to be used
#' @param y outcome variable
#' @param time time variable
#' @param params \code{NULL} unless user wants to specify starting parameters for gnls
#' @param ... just in case
#'
#' @details Remove any values of zero, or jitter, before using with bdotsFit
#'
#' @details \code{y ~ x_0 exp(k beta)}
#' @export
expCurve <- function(dat, y, time, params = NULL, ...) {
estExpPars <- function(dat, y, time) {
tt <- lm(log(dat[[y]]) ~ dat[[time]])
x0 <- exp(coef(tt)[1])
k <- coef(tt)[2]
names(x0) <- names(k) <- NULL
return(c(x0 = x0, k = k))
}

if (any(dat[[y]] <= 0)) {
message("Subjects with values <= 0 will not be fit")
return(NULL)
}

if (is.null(params)) {
params <- estExpPars(dat, y, time)
} else {
# put checks here
}
y <- str2lang(y)
time <- str2lang(time)
ff <- bquote(.(y) ~ x0 * exp(.(time) * k))
attr(ff, "parnames") <- names(params)
return(list(formula = ff, params = params))
}

#' ## -----------
#'
#' #' Logistic curve function for nlme
#' #'
#' #' DELETE THIS FUNCTION BEFORE GIT PUSH
#' #'
#' #' @param dat subject data to be used
#' #' @param y outcome variable
#' #' @param time time variable
#' #' @param params \code{NULL} unless user wants to specify starting parameters for gnls
#' #' @param ... just in case
#' #'
#' #' @details \code{y ~ mini + (peak - mini) / (1 + exp(4 * slope * (cross - (time)) / (peak - mini)))}
#' #' @export
#' logistic2 <- function(dat, y, time, params = NULL, ...) {
#'
#'   logisticPars <- function(dat, y, time, ...) {
#'     dat[, looks := mean(looks), by = starttime]
#'     dat <- unique(dat)
#'     time <- dat[[time]]
#'     y <- dat[[y]]
#'
#'     ## Remove cases with zero variance
#'     if (var(y) == 0) {
#'       return(NULL)
#'     }
#'
#'     mini <- min(y)
#'     peak <- max(y)
#'     r <- (peak - mini)
#'     cross <- time[which.min(abs(0.5*r - y))]
#'
#'     # slope
#'     q75 <- .75 * r + mini
#'     q25 <- .25 * r + mini
#'     time75 <- time[which.min(abs(q75 - y))]
#'     time25 <- time[which.min(abs(q25 - y))]
#'     slope <- (q75 - q25) / (time75 - time25)
#'
#'     return(c(mini = mini, peak = peak, slope = slope, cross = cross))
#'   }
#'
#'   if (is.null(params)) {
#'     params <- logisticPars(dat, y, time)
#'   } else {
#'     if (length(params) != 4) stop("logistic requires 4 parameters be specified for refitting")
#'     if (!all(names(params) %in% c("mini", "peak", "slope", "cross"))) {
#'       stop("logistic parameters for refitting must be correctly labeled")
#'     }
#'   }
#'   ## Return NA list if var(y) is 0
#'   if (is.null(params)) {
#'     return(NULL)
#'   }
#'   y <- str2lang(y)
#'   time <- str2lang(time)
#'   ff <- bquote(.(y) ~ mini + (peak - mini) / (1 + exp(4 * slope * (cross - (.(time))) / (peak - mini))))
#'   attr(ff, "parnames") <- names(params)
#'   return(list(formula = ff, params = params))
#' }
```

## Try the bdots package in your browser

Any scripts or data that you put into this service are public.

bdots documentation built on Jan. 7, 2023, 1:18 a.m.