Nothing
lee <- function (conc, time, points = 3, method = c("ols", "lad", "hub", "npr"), lt = TRUE){
lad <- function(y, x) {
resid <- Inf
for (i in 1:length(y)) {
for (j in 1:length(y)) {
if ((i != j) & (x[j] != x[i])) {
slope <- (y[j] - y[i])/(x[j] - x[i])
intct <- y[j] - slope * x[j]
absresid <- abs(y - (intct + slope * x))
if (sum(absresid) < resid) {
mad <- median(absresid)
resid <- sum(absresid)
k <- slope
d <- intct
}
}
}
}
return(list(k = as.double(k), d = as.double(d), resid = as.double(resid),
mad = as.double(mad)))
}
hub <- function(y, x, mad = lad(y = y, x = x)$mad) {
hubloss <- function(kd) {
absresid <- abs(y - kd[[1]] * x - kd[[2]])
khuber <- 2.2245 * mad
sum(ifelse(absresid < khuber, absresid * absresid,
khuber * (2 * absresid - khuber)))
}
start <- as.vector(c(lm(y ~ x)$coef[2], lm(y ~ x)$coef[1]))
res <- optim(start, hubloss, method = "Nelder-Mead",
control = c(reltol = 1e-09))
return(list(k = as.double(res$par[1]), d = as.double(res$par[2]),
resid = as.double(res$value)))
}
npr <- function(y, x) {
weighted.median <- function(w, x) {
data <- data.frame(x, w, stringsAsFactors = TRUE)
data <- data[order(data$x), ]
i <- 1
while (sum(data$w[1:i]) <= 0.5) {
i <- i + 1
}
ifelse(sum(data$w[1:i - 1]) == 0.5, return((data$x[i -
1] + data$x[i])/2), return(data$x[i]))
}
total <- 0
for (i in 1:(length(y) - 1)) {
for (j in (i + 1):length(y)) {
total <- total + abs(x[i] - x[j])
}
}
l <- 1
b <- array(1:(length(y) * (length(y) - 1)/2))
w <- array(1:(length(y) * (length(y) - 1)/2))
for (i in 1:(length(y) - 1)) {
for (j in (i + 1):length(y)) {
b[l] <- ((y[i] - y[j])/(x[i] - x[j]))
w[l] <- abs(x[i] - x[j])/total
l <- l + 1
}
}
data <- subset(data.frame(w = as.vector(w), b = as.vector(b)),
b != Inf & b != -Inf)
k <- weighted.median(w = data$w, x = data$b)
d <- median(y - k * x)
e <- y - k * x - d
resid <- sum((rank(e) - 1/2 * (length(y) + 1)) * e)
return(list(k = as.double(k), d = as.double(d), resid = as.double(resid)))
}
ols <- function(y, x) {
res <- lm(y ~ x)
return(list(k = as.double(res$coef[2]), d = as.double(res$coef[1]),
resid = as.double(sum(res$resid * res$resid))))
}
method = match.arg(method)
if (!is.vector(time) || !is.vector(conc)) {
stop("argument time and/or conc invalid")
}
if (length(time) != length(conc)) {
stop("time and conc differ in length")
}
if (!is.logical(lt)) {
stop("argument lt invalid")
}
if (points < 2) {
stop("not enough points in terminal phase")
}
data <- na.omit(data.frame(conc, time, stringsAsFactors = TRUE))
if (any(data$conc <= 0)) {
data$conc[data$conc <= 0] <- NA
warning("concentration below or equal to zero were omitted")
data <- na.omit(data)
}
if (nrow(data) < 4) {
stop("a minimum of 4 observations are required")
}
data <- data[order(data$time), ]
n <- nrow(data)
conc <- log10(data$conc)
time <- data$time
switch(method, lad = {
model <- lad(y = conc, x = time)
}, ols = {
model <- ols(y = conc, x = time)
}, hub = {
model <- hub(y = conc, x = time, mad = lad(y = conc,
x = time)$mad)
}, npr = {
model <- npr(y = conc, x = time)
}, )
final.term.model <- model
final.init.model <- model
resid <- model$resid
final.chgpt <- NA
if (model$k >= 0) {
resid <- Inf
final.term.model$k <- NA
final.init.model$k <- NA
}
if (points > n - 2) {
stop("not enough points for inital phase")
}
for (i in 2:(n - points)) {
init.conc <- conc[1:i]
init.time <- time[1:i]
term.conc <- conc[(i + 1):n]
term.time <- time[(i + 1):n]
switch(method, lad = {
init.model <- lad(y = init.conc, x = init.time)
term.model <- lad(y = term.conc, x = term.time)
}, ols = {
init.model <- ols(y = init.conc, x = init.time)
term.model <- ols(y = term.conc, x = term.time)
}, hub = {
init.model <- hub(y = init.conc, x = init.time, mad = lad(y = init.conc,
x = init.time)$mad)
term.model <- hub(y = term.conc, x = term.time, mad = lad(y = term.conc,
x = term.time)$mad)
}, npr = {
init.model <- npr(y = init.conc, x = init.time)
term.model <- npr(y = term.conc, x = term.time)
}, )
lower <- data$time[i]
upper <- data$time[i + 1]
chgpt <- (init.model$d - term.model$d)/(term.model$k -
init.model$k)
if (chgpt > lower & chgpt < upper) {
if (!lt) {
if (sum(term.model$resid, init.model$resid) <
resid) {
final.init.model <- init.model
final.term.model <- term.model
final.chgpt <- as.double(chgpt)
resid <- sum(term.model$resid, init.model$resid)
}
}
if (lt) {
if (init.model$k <= term.model$k & term.model$k <
0 & sum(term.model$resid, init.model$resid) <
resid) {
final.init.model <- init.model
final.term.model <- term.model
final.chgpt <- as.double(chgpt)
resid <- sum(term.model$resid, init.model$resid)
}
}
}
}
init.hl <- -log10(2)/final.init.model$k
term.hl <- -log10(2)/final.term.model$k
if (is.na(init.hl) | is.na(term.hl)) {
warning("No model evaluated")
}
parms <- data.frame(initial = as.double(c(init.hl, final.init.model$k,
final.init.model$d)), terminal = as.double(c(term.hl, final.term.model$k,
final.term.model$d)), stringsAsFactors = TRUE)
rownames(parms) <- c("halflife", "slope", "intercept")
res <- list(parms = parms, chgpt = as.double(final.chgpt),
conc = 10^conc, time = time, method = "lee")
class(res) <- "halflife"
return(res)
}
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.