Nothing
#' @importFrom rlang .data
fit_lin <- function(time, y) {
if (missing(y)) {
stop(gettextf("Missing 'y' vector"))
}
{
if (missing(time)) {
stop(gettextf("Missing 'time' vector"))
}
}
for (i in 1:length(y)) {
if (y[i] > 1) {
stop(gettextf("values must between 0 and 1)"))
}
if (y[i] >= 1) {
y[i] <- 0.999
gettextf("values = '1' converted to '0.999'")
}
if (y[i] == 0) {
y[i] <- 0.0001
gettextf("values = '0' converted to '0.0001'")
}
}
epi <- data.frame(time, y)
model =
value =
v0 =
v0_se=
v0_ci_lwr =
v0_ci_upr =
CCC =
best_model =
linear =
linearized =
name =
r =
r_se =
r_ci_lwr =
y0 =
y0_ci_lwr =
y0_ci_upr =
y0 =
df =
y0_se =
r =
r_se =
CCC =
r_ci_lwr =
y0_ci_lwr =
y0_ci_upr =
model =
time =
y = NULL
epi$exponit <- log(epi$y)
epi$monit <- log(1 / (1 - epi$y))
epi$logit <- log(epi$y / (1 - epi$y))
epi$gompit <- -log(-log(epi$y))
# colnames(epi) <- c("time", "y", "exponit", "monit", "logit", "gompit")
# Suppress summarise info
options(dplyr.summarise.inform = FALSE)
z <- epi %>%
tidyr::gather(3:6, key = "model", value = "value") %>%
dplyr::group_by(model) %>%
dplyr::summarise(
r = summary(stats::lm(value ~ time))$coefficients[2, 1],
r_se = summary(stats::lm(value ~ time))$coefficients[2, 2],
r_ci_lwr = as.numeric(as.matrix(stats::confint(stats::lm(value ~ time)))[2, 1]),
r_ci_upr = as.numeric(as.matrix(stats::confint(stats::lm(value ~ time)))[2, 2]),
v0 = as.numeric(as.matrix((stats::lm(value ~ time))$coefficients)[1, 1]),
v0_se = summary(stats::lm(value ~ time))$coefficients[1, 2],
v0_ci_lwr = as.numeric(as.matrix(stats::confint(stats::lm(value ~ time)))[1, 1]),
v0_ci_upr = as.numeric(as.matrix(stats::confint(stats::lm(value ~ time)))[1, 2]),
r_squared = summary(stats::lm(value ~ time))$r.squared,
RSE = summary(stats::lm(value ~ time))$sigma,
CCC = DescTools::CCC(stats::lm(value ~ time)$fitted.values, value)$rho.c$est
) %>%
dplyr::mutate(
y0 = dplyr::case_when(
model == "exponit" ~ exp(v0),
model == "gompit" ~ exp(-exp(-v0)),
model == "logit" ~ (1 / (1 + exp(-v0))),
model == "monit" ~ 1 - exp(-v0)
),
y0_ci_lwr = dplyr::case_when(
model == "exponit" ~ exp(v0_ci_lwr),
model == "gompit" ~ exp(-exp(-v0_ci_lwr)),
model == "logit" ~ (1 / (1 + exp(-v0_ci_lwr))),
model == "monit" ~ 1 - exp(-v0_ci_lwr)
),
y0_ci_upr = dplyr::case_when(
model == "exponit" ~ exp(v0_ci_upr),
model == "gompit" ~ exp(-exp(-v0_ci_upr)),
model == "logit" ~ (1 / (1 + exp(-v0_ci_upr))),
model == "monit" ~ 1 - exp(-v0_ci_upr)
),
model = dplyr::case_when(
model == "exponit" ~ "Exponential",
model == "gompit" ~ "Gompertz",
model == "logit" ~ "Logistic",
model == "monit" ~ "Monomolecular"
)
) %>%
dplyr::select(-v0_ci_lwr, -v0_ci_upr) %>%
dplyr::arrange(-CCC) %>%
dplyr::mutate(best_model = 1:4) %>%
dplyr::select(best_model, 1:13)
# z
epi2 = epi %>%
dplyr::group_by(time) %>%
dplyr::mutate(rep = 1:dplyr::n())
predicted <- epi2 %>%
dplyr::mutate(Exponential = dplyr::filter(z, model == "Exponential")$y0 * exp(dplyr::filter(z, model == "Exponential")$r * time)) %>%
dplyr::mutate(Monomolecular = 1 - (1 - dplyr::filter(z, model == "Monomolecular")$y0) * exp(-dplyr::filter(z, model == "Monomolecular")$r * time)) %>%
dplyr::mutate(Logistic = 1 / (1 + ((1 - dplyr::filter(z, model == "Logistic")$y0) / dplyr::filter(z, model == "Logistic")$y0) * exp(-dplyr::filter(z, model == "Logistic")$r * time))) %>%
dplyr::mutate(Gompertz = exp(log(dplyr::filter(z, model == "Gompertz")$y0) * exp(-dplyr::filter(z, model == "Gompertz")$r * time))) %>%
tidyr::pivot_longer(c(3:6,8:11)) %>%
dplyr::mutate(type = dplyr::case_when(name %in% c("exponit","monit","logit","gompit")~"linearized",
name %in% c("Exponential","Monomolecular","Logistic","Gompertz")~"predicted")) %>%
dplyr::mutate(model = dplyr::case_when(name == "exponit"~ "Exponential",
name == "monit" ~ "Monomolecular",
name == "logit" ~ "Logistic",
name == "gompit" ~ "Gompertz",
name %in% c("Exponential","Monomolecular","Logistic","Gompertz")~name )) %>%
dplyr::select(-name) %>%
tidyr::pivot_wider(values_from = "value",
names_from = c("type"),
id_cols = c("time","y","model","rep")) %>%
dplyr::select(-rep) %>%
dplyr::mutate(residual = y - predicted)
za <- z %>%
dplyr::mutate(Estimate = r, Std.error = r_se, Lower = r_ci_lwr, Upper = r_ci_lwr)
zb <- z %>%
dplyr::mutate(Estimate = y0, Linearized = v0, lin.SE = v0_se, Lower = y0_ci_lwr, Upper = y0_ci_upr)
z1 <- as.matrix(z[, c("CCC", "r_squared", "RSE")])
rownames(z1) <- as.matrix(za[, "model"])
z2 <- as.matrix(za[, c("Estimate", "Std.error", "Lower", "Upper")])
rownames(z2) <- as.matrix(za[, "model"])
z3 <- as.matrix(zb[, c("Estimate","Linearized","lin.SE", "Lower", "Upper")])
rownames(z3) <- as.matrix(zb[, "model"])
a <- list(
header = "Results",
Stats = round(z1, 4),
`Infection rate` = z2,
`Initial inoculum` = z3,
data = as.data.frame(predicted),
stats_all = z
)
class(a) <- "fit_lin"
return(a)
}
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.