Nothing
#' 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))
#' }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.