Executive Summary
How 'root length cm' responds as 'ferulic acid mM' increases, and where the midpoint sits.
The short answer
Ryegrass root length is cut in half at a ferulic acid concentration of 2.992 mM (95% CI 2.479–3.612 mM). This midpoint falls squarely within the tested dose range and is therefore interpolated from real observations rather than extrapolated.
The detail
The dose-response curve spans from 7.976 cm at low dose to 0.456 cm at high dose across 18 measurements at 6 distinct dose levels. The IC50 of 2.9923 mM has a confidence interval of 2.4791 to 3.6118 mM. The Hill slope is -2.809, indicating a steep transition zone. The curve explains 97.09% of the variation in root length. The replicate-based lack-of-fit test does not reject (p = 0.801), confirming the logistic shape is consistent with the observed scatter.
What this can't tell you
The tested doses ran from 0.9400 to 30.00 mM, so both the midpoint and both plateaus are measured rather than projected. However, the analysis assumes the relationship is causal only if ferulic acid was assigned experimentally; if doses were observed, the curve describes association, not causation.
Analysis Overview
Four-parameter logistic fit of 'root length cm' against 'ferulic acid mM' across 18 rows.
The short answer
Ryegrass root growth falls from about 8 cm at low ferulic acid exposure to about 0.5 cm at high exposure. The dose that cuts growth to half that span—the IC50—is 2.992 mM, with a 95% confidence range of 2.479 to 3.612 mM.
The detail
The fitted four-parameter logistic describes 'root length cm' as a function of log₁₀ 'ferulic acid mM' across 18 observations at 6 distinct dose levels. The low-dose plateau is 7.976 cm and the high-dose plateau is 0.4560 cm, for a total span of 7.520 cm. The IC50 of 2.992 mM falls squarely within the tested dose range (0.9400 to 30.00 mM), so it is interpolated from observed data rather than extrapolated. The Hill slope is −2.809, indicating a steep turnover: the useful dose window between little effect and most of the effect is narrow. The fitted model accounts for 97.09% of the variance in root length and leaves a residual standard error of 0.5657 cm.
What this can't tell you
The confidence interval (2.479 to 3.612 mM) reflects asymptotic Wald intervals, which run slightly optimistic for a nonlinear model; profile-likelihood intervals would be slightly wider. The model comparison shows that a simpler baseline-constrained three-parameter logistic has a lower AIC (34.34 vs. 36.05), and the nested test does not reject that constraint (F(1, 14) = 0.225, p = 0.643), so the fourth fitted parameter is not supported by the data and the four-parameter curve's precision may be overstated.
Data Quality
Row accounting, dose validity, and control handling.
The short answer
All 24 rows loaded without loss: 18 rows at positive ferulic acid doses were used to fit the curve, and 6 zero-dose control rows (mean root length 7.749 cm) were held out as the untreated baseline because zero has no logarithm.
The detail
No rows were removed for missing or invalid values. The 18 positive-dose rows span 6 distinct 'ferulic acid mM' levels. The 6 zero-dose control rows cannot be placed on a log₁₀ dose axis, so they are excluded from the curve fit and used instead to anchor the low-dose plateau (7.976 cm) against the untreated mean (7.749 cm). Both the dose and response columns were coerced to numbers under the rule that at least 95% of non-blank values must parse successfully. The response was checked for variation before curve fitting was attempted.
What this can't tell you
The preprocessing step does not reveal the replicate structure—how many observations sit at each dose level—though the residual diagnostics card shows the scatter. A finer-grained export listing dose level, replicate number, and individual root length would clarify whether replicates are balanced across doses and whether any single dose level drives the fit.
Fitted Dose-Response Curve
The 4PL curve drawn over the observed 'root length cm' values on a log10 'ferulic acid mM' axis.
The short answer
The fitted curve falls smoothly from 7.976 cm at low dose to 0.4560 cm at high dose, with observed points scattered evenly around it and no systematic drift at either end. The midpoint (IC50 = 2.992 mM) is marked by the vertical reference line.
The detail
The horizontal axis is log₁₀ 'ferulic acid mM', the scale on which the logistic curve is symmetric. The fitted curve runs from low-dose plateau 7.976 cm to high-dose plateau 0.4560 cm, turning over at log₁₀(2.992) = 0.476. The 6 zero-dose control rows (mean 7.749 cm) are drawn as a separate series one decade below the lowest tested dose—that position is a drawing convention for display, not a measured dose. Points scatter evenly around the fitted curve with no systematic arc or drift, supporting the logistic shape. Residuals at the low end and high end do not cluster systematically above or below the curve, which would signal a wrong functional form.
What this can't tell you
A visual scatter plot cannot quantify how much deviation from the curve is tolerable; that is what the lack-of-fit test answers (p = 0.801). The plot alone cannot reveal whether variance grows with dose (a funnel pattern), which would argue for a weighted fit; that requires the residual-diagnostics card.
Curve Parameters
The four fitted parameters with 95% confidence intervals.
| Parameter | Estimate | CI Low | CI High | Interpretation |
|---|---|---|---|---|
| Plateau at low 'ferulic acid mM' | 7.976 | 6.907 | 9.044 | The 'root length cm' the curve settles to as 'ferulic acid mM' approaches zero. The lowest tested doses reach this plateau, so it is observed rather than projected. |
| Plateau at high 'ferulic acid mM' | 0.456 | -0.0596 | 0.9716 | The 'root length cm' the curve saturates at as 'ferulic acid mM' grows large. The highest tested doses reach this plateau, so it is observed rather than projected. |
| IC50 (midpoint dose) | 2.992 | 2.479 | 3.612 | The 'ferulic acid mM' at which 'root length cm' is halfway between the two plateaus. It falls inside the tested range of 0.9400 to 30.00, so it is interpolated from observed doses. |
| Hill slope (steepness) | -2.809 | -4.09 | -1.528 | How sharply 'root length cm' turns over near the midpoint. A steeper slope means a narrower 'ferulic acid mM' window between little effect and most of the effect; the sign is negative because 'root length cm' falls with 'ferulic acid mM'. |
| Span (high plateau minus low plateau) | -7.519 | — | — | The total achievable change in 'root length cm' across the fitted curve, from 7.976 to 0.4560. |
The short answer
The IC50 is 2.9923 mM with a 95% confidence interval from 2.4791 to 3.6118 mM, a 1.457-fold range. Both plateaus (7.976 and 0.456 cm) were reached by the tested doses, so they are observed rather than projected.
The detail
The low-dose plateau estimate is 7.9755 cm (95% CI 6.907–9.044). The high-dose plateau is 0.456 cm (95% CI −0.0596 to 0.9716). The IC50 of 2.9923 mM is bracketed by tested doses from 0.9400 to 30.00 mM, so it is interpolated. The Hill slope is −2.8093 (95% CI −4.0902 to −1.5284), indicating a steep turn. These are Wald intervals from the nonlinear fit and are asymptotic, running slightly narrow compared with profile-likelihood intervals.
What this can't tell you
The IC50 interval is asymmetric around the estimate because it is estimated on the log₁₀ scale and back-transformed—this is correct for a dose but means the lower and upper bounds do not deviate equally in arithmetic space.
Model Comparison
4PL against a baseline-constrained 3PL and a plain log-linear model.
| Model | Parameters | Residual SE | Aic | Status | Aic Vs Best |
|---|---|---|---|---|---|
| 4PL (four-parameter logistic) | 4 | 0.5657 | 36.05 | fitted | 1.71 |
| 3PL (baseline held at the untreated control) | 3 | 0.5509 | 34.34 | fitted with the low-dose plateau held at the zero-dose control mean (7.749) | 0 |
| Log-linear (response on log10 dose) | 2 | 1.096 | 58.27 | fitted | 23.93 |
The short answer
The baseline-constrained three-parameter logistic has the lowest AIC (34.34) and is preferred over the four-parameter logistic (AIC 36.05). The nested test does not reject the constraint that the low-dose plateau equals the zero-dose control mean (F(1, 14) = 0.225, p = 0.643), so the simpler model is more defensible.
The detail
Three models were fitted to the same 18 rows: the 4PL (4 parameters, AIC 36.05, residual SE 0.5657), the baseline-constrained 3PL (3 parameters, AIC 34.34, residual SE 0.5509), and a log-linear model (2 parameters, AIC 58.27, residual SE 1.0964). The 3PL holds the low-dose plateau at the zero-dose control mean (7.749 cm) and estimates only the ceiling, midpoint, and slope. The F-test of the 4PL against the 3PL yields F(1, 14) = 0.225, p = 0.643: the fourth parameter (the free low-dose plateau) is not statistically earned. The log-linear model, included as a floor, has no floor, ceiling, or midpoint and cannot produce an IC50; it serves only to show that the curve structure is necessary.
What this can't tell you
On a sparse design with few dose levels or replicates, a four-parameter curve can fit noise rather than shape, and the AIC comparison is the main diagnostic. Here, the 3PL's lower AIC signals that the extra parameter in the 4PL is not supported, so the 4PL's reported IC50 interval is likely too narrow. Consider the 3PL estimates as more credible for this dataset.
Group Comparison
Per-group curves and the test of whether the midpoints differ.
The short answer
No grouping column was mapped, so one curve was fitted across all 18 observations and there are no groups to compare.
The detail
The analysis requires a grouping column (such as compound, channel, treatment arm, or cultivar) to fit separate curves per group and test whether their midpoints differ. With a single pooled curve across all rows, the comparison test cannot run.
What this can't tell you
If ryegrass cultivars, soil types, or experimental batches are recorded in the source data, mapping one as a grouping column would reveal whether the ferulic acid response differs among them. That would require a re-run with the grouping column specified in the analysis parameters.
Residual Diagnostics
Residuals against fitted values — the check the fit statistics cannot do.
The short answer
Residuals scatter evenly around zero with no systematic arc or funnel, and the lack-of-fit test does not reject the logistic shape (p = 0.801). Residuals are consistent with normality (Shapiro-Wilk p = 0.651).
The detail
The 18 residuals (observed minus predicted root length) have a standard error of 0.5657 cm. They cluster in a formless band around zero with no systematic drift—no arc that would signal the wrong functional form, no funnel that would indicate heteroscedasticity. The pure-error lack-of-fit test, which compares each dose level's mean deviation from the curve against its own replicate scatter, yields F(2, 12) = 0.226, p = 0.801: the curve shape is consistent with the observed replicate variation. The Shapiro-Wilk test for normality yields p = 0.651, so residuals do not deviate significantly from the normal distribution that the confidence intervals assume.
What this can't tell you
Residual diagnostics confirm that the logistic shape fits the observed scatter, but they do not resolve whether the four-parameter or three-parameter model is more credible—that is answered by the model-comparison card. The lack-of-fit test is robust only if replicates are present at each dose level; the residual plot shows this is the case here (three observations per dose level), so the test is valid.
Methods & Disclosure
Exactly how the curve was fitted, and what it can and cannot decide.
| Item | Detail |
|---|---|
| Model | Four-parameter logistic on log10 'ferulic acid mM': response = low plateau + (high plateau - low plateau) / (1 + 10^((log10(IC50) - log10(dose)) x hill)). |
| Estimation | Ordinary nonlinear least squares (base R nls) over 18 rows at 6 distinct positive 'ferulic acid mM' levels; residual degrees of freedom 14, residual standard error 0.5657. |
| Starting values | Derived from the data, not hard-coded: the plateaus start at the mean 'root length cm' at the lowest and highest tested dose, and normalising by them linearises the curve on the logit scale, so an ordinary least-squares line supplies the starting hill slope (2.125) and midpoint (3.660). |
| Confidence intervals | Wald intervals, estimate plus or minus 2.14 standard errors on 14 degrees of freedom. These are asymptotic: for a nonlinear model they are slightly optimistic compared with profile-likelihood intervals. |
| Midpoint interval | The midpoint is estimated on the log10 scale and back-transformed, so its interval (2.479 to 3.612) is asymmetric around 2.992 — which is the correct shape for a dose. |
| Zero-dose controls | 6 zero-dose control row(s) were found, mean 'root length cm' 7.749. They cannot sit on a log-dose axis, so they are excluded from the curve fit and drawn as their own series one decade below the lowest tested dose — that position is a drawing convention, not a measured dose. |
| Model comparison | The 4PL is compared against a three-parameter logistic (fitted with the low-dose plateau held at the zero-dose control mean (7.749)) and a log-linear model on AIC; the baseline-constrained three-parameter logistic has the lowest AIC here. |
| Lack of fit | Lack-of-fit F(2, 12) = 0.226, p = 0.801 |
| Convergence policy | If nls does not converge from the self-start, the bounded port algorithm and a grid of perturbed starts are tried. If all fail, the analysis stops and reports the failure with a diagnosis. A linear or log-linear fit is never substituted for the curve. |
| Causal standing | This is a fitted description of how 'root length cm' varies with observed 'ferulic acid mM'. Unless the doses were assigned experimentally, the curve is associated with dose and does not by itself establish that changing 'ferulic acid mM' causes the change in 'root length cm'. |
The short answer
The curve is a four-parameter logistic fitted by ordinary nonlinear least squares on log₁₀ ferulic acid mM, with starting values derived from the data. The confidence intervals are Wald intervals and therefore asymptotic—slightly optimistic for a nonlinear model.
The detail
The model is response = low plateau + (high plateau − low plateau) / (1 + 10^((log₁₀(IC50) − log₁₀(dose)) × hill)). Estimation used base R nls over 18 rows at 6 positive dose levels, with residual degrees of freedom 14 and residual standard error 0.5657. Starting values were derived from the data: plateaus from the extreme dose means, and hill slope (2.125) and midpoint (3.660) from a logit-linearised least-squares fit. Wald intervals are estimate ± 2.14 standard errors. Six zero-dose controls (mean 7.749 cm) were excluded from the curve fit and drawn separately. The 4PL was compared against a three-parameter logistic (with low plateau fixed at 7.749) and a log-linear model by AIC; the 3PL had the lowest AIC. Convergence policy: if nls fails, the bounded port algorithm and perturbed starts are tried; if all fail, the analysis stops and reports why rather than substituting a linear fit.
What this can't tell you
Unless ferulic acid was assigned experimentally, the curve describes association, not causation. Wald intervals are the optimistic end of the plausible range for a nonlinear model; profile-likelihood intervals would be slightly wider.
Dose-Response Curve Fitting — EC50 / IC50
How does response change as dose, spend, or exposure increases? The analysis fits a four-parameter logistic (4PL) curve by nonlinear least squares on log10 dose:
response = low plateau + (high plateau - low plateau) / (1 + 10 ^ ((log10(EC50) - log10(dose)) * hill))
and reports the four parameters with confidence intervals, the EC50 (or IC50 when the response falls) as the headline with its own interval, the fitted curve drawn over the observed points, residual diagnostics, and a comparison against a baseline-constrained 3PL and a plain log-linear model so a sparse dataset is not locked into an over-parameterised curve.
Why This Method?
A saturating response is not a straight line and not a log line: it has a floor, a ceiling, and a midpoint. The 4PL is the standard description of that shape, and its midpoint parameter — the EC50/IC50 — is the single number practitioners compare across compounds, channels, or campaigns.
What This Analysis Covers
- 4PL fit by base
nlswith data-derived self-starting values - Every parameter with a 95% interval; EC50/IC50 back-transformed from the
log-dose scale so its interval is properly asymmetric
- The fitted curve over the observed points on a log-dose axis, with
zero-dose controls shown separately and never fed to the log axis
- Residual diagnostics including a replicate-based lack-of-fit test
- 4PL vs 3PL vs log-linear model comparison
- Parallel curves per group with a formal test of whether the EC50s differ
Honesty
If nls cannot converge the module REFUSES and says why. It never falls back to a straight line and calls it a dose-response curve. If the tested doses do not bracket the midpoint, the EC50 is named as an extrapolation and the observed dose range is printed next to it.
Standard Library
Platform standard-library module (LAT-1441): runs on ANY dataset via the semantic mapping {dose, response, group}. All narrative is derived from the user's own column names and computed values.
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))Core Analysis Pipeline
Step 1: Resolve the mapped columns, humanized for every sentence
initial_rows <- nrow(df)
dose_h <- humanize_semantic("dose", col_map)
resp_h <- humanize_semantic("response", col_map)
has_group_col <- "group" %in% names(df)
group_h <- if (has_group_col) humanize_semantic("group", col_map) else "group"
if (!("dose" %in% names(df)) || !("response" %in% names(df))) {
stop(sprintf("Dose-response fitting needs both a dose column('%s') and a response column ('%s') mapped.",
dose_h, resp_h))
}Step 2: Coerce both to numeric under the 95% rule
coerce_num <- function(v, label_h, role) {
if (is.numeric(v)) return(as.numeric(v))
ch <- as.character(v)
non_blank <- !is.na(ch) & trimws(ch) != ""
conv <- suppressWarnings(as.numeric(ch))
if (sum(non_blank) == 0 ||
sum(!is.na(conv[non_blank])) < 0.95 * sum(non_blank)) {
stop(sprintf("The column '%s' was mapped as the %s but does not look numeric — fewer than 95%% of its values parse as numbers. Map a numeric column.",
label_h, role))
}
conv
}
dose_v <- coerce_num(df$dose, dose_h, "dose")
resp_v <- coerce_num(df$response, resp_h, "response")
grp_v <- if (has_group_col) {
g <- trimws(as.character(df$group))
g[is.na(g) | g == ""] <- "Missing"
g
} else rep("All observations", initial_rows)Step 3: Drop unusable rows and separate the zero-dose controls
ok <- !is.na(dose_v) & !is.na(resp_v)
n_drop_na <- sum(!ok)
dose_v <- dose_v[ok]; resp_v <- resp_v[ok]; grp_v <- grp_v[ok]
neg <- dose_v < 0
n_drop_nonpos <- sum(neg)
dose_v <- dose_v[!neg]; resp_v <- resp_v[!neg]; grp_v <- grp_v[!neg]
is_zero <- dose_v == 0
n_zero <- sum(is_zero)
zero_mean <- if (n_zero > 0) mean(resp_v[is_zero]) else NA_real_
zero_resp <- resp_v[is_zero]
zero_grp <- grp_v[is_zero]
d <- dose_v[!is_zero]; y <- resp_v[!is_zero]; g <- grp_v[!is_zero]
n <- length(d)
final_rows <- n + n_zero
rows_removed <- initial_rows - final_rowsStep 4: Hard guards, each naming the user's own columns
if (n < MIN_ROWS) {
stop(sprintf("Only %d rows have a positive '%s' and a usable '%s' value — a four-parameter dose-response curve needs at least %d. Zero-dose control rows (%d here) anchor the baseline but cannot be placed on a log-dose axis, so they do not count toward this minimum.",
n, dose_h, resp_h, MIN_ROWS, n_zero))
}
n_levels <- length(unique(d))
if (n_levels < MIN_LEVELS) {
stop(sprintf("'%s' has only %d distinct positive value(s) (%s). A four-parameter logistic curve has four unknowns and cannot be determined from fewer than %d distinct dose levels — the fit would be arbitrary rather than estimated.",
dose_h, n_levels,
paste(fmt_num(sort(unique(d))[1]), "to",
fmt_num(sort(unique(d))[n_levels])),
MIN_LEVELS))
}
if (!isTRUE(stats::var(y) > 0)) {
stop(sprintf("'%s' is constant — every value is identical — so there is no response to model against '%s'.",
resp_h, dose_h))
}
u <- log10(d)
dose_min <- min(d); dose_max <- max(d)Step 5: Refuse before fitting when the shape is not a dose-response
A 4PL describes a monotonic saturating curve. When the response does not move consistently with dose there is nothing for the curve to estimate, and reporting some fitted line as a dose-response would be a fabrication.
sp <- suppressWarnings(tryCatch(
stats::cor.test(u, y, method = "spearman"),
error = function(e) NULL))
sp_rho <- if (!is.null(sp)) unname(sp$estimate) else NA_real_
sp_p <- if (!is.null(sp)) sp$p.value else NA_real_
if (!is.finite(sp_rho) || !is.finite(sp_p) || abs(sp_rho) < 0.25 || sp_p >= 0.05) {
stop(sprintf("'%s' shows no consistent monotonic change across '%s' (Spearman rho %s, %s over %d rows and %d dose levels). A dose-response curve cannot be fitted to a response that does not rise or fall with dose — the shape may be flat, or it may peak in the middle of the range, which this method cannot represent. No curve is reported rather than a straight line dressed up as one.",
resp_h, dose_h, r3(sp_rho), fmt_pp(sp_p), n, n_levels))
}Step 6: Self-start, then fit the 4PL by nonlinear least squares
st <- dr_selfstart(u, y)
fit <- dr_fit4pl(u, y, st)
if (is.null(fit)) {
diag_bits <- character()
if (n_levels < 6) {
diag_bits <- c(diag_bits, sprintf("only %d distinct '%s' levels were tested, which is thin for a four-parameter curve",
n_levels, dose_h))
}
lev <- sort(unique(u))
mu_lev <- sapply(lev, function(v) mean(y[u == v]))
k <- length(mu_lev)
edge <- max(1L, floor(k / 4))
lowflat <- stats::sd(mu_lev[seq_len(edge + 1)])
highflat <- stats::sd(mu_lev[(k - edge):k])
spanobs <- abs(mu_lev[k] - mu_lev[1])
if (is.finite(lowflat) && is.finite(spanobs) && spanobs > 0 &&
lowflat > 0.25 * spanobs) {
diag_bits <- c(diag_bits, sprintf("the response is still moving at the lowest doses, so the lower plateau was never observed"))
}
if (is.finite(highflat) && is.finite(spanobs) && spanobs > 0 &&
highflat > 0.25 * spanobs) {
diag_bits <- c(diag_bits, sprintf("the response is still moving at the highest doses, so the upper plateau was never observed"))
}
if (abs(sp_rho) < 0.6) {
diag_bits <- c(diag_bits, sprintf("the dose-response trend is weak and noisy(Spearman rho %s)", r3(sp_rho)))
}
if (length(diag_bits) == 0) {
diag_bits <- "the residual surface has no stable minimum from any of the starting values tried"
}
stop(sprintf("The four-parameter dose-response curve for '%s' against '%s' did not converge. Likely cause: %s. Tested doses ran from %s to %s across %d levels and %d rows. No curve, EC50, or fitted parameter is reported — a straight-line fit is NOT substituted for the curve, because it would answer a different question.",
resp_h, dose_h,
paste(diag_bits, collapse = "; "),
fmt_num(dose_min), fmt_num(dose_max), n_levels, n))
}Step 7: Extract parameters. h is steepness; the DIRECTION of the
curve is the sign of (high plateau - low plateau), so a negative fitted h is the same curve with the plateaus swapped. Normalise to h > 0 so "low-dose plateau" always means what it says.
cf <- summary(fit)$coefficients
est <- cf[, 1]; se <- cf[, 2]
df_res <- stats::df.residual(fit)
tq <- stats::qt(0.975, max(1, df_res))
lo <- unname(est["lo"]); hi <- unname(est["hi"])
le50 <- unname(est["le50"]); h <- unname(est["h"])
lo_se <- unname(se["lo"]); hi_se <- unname(se["hi"])
le50_se <- unname(se["le50"]); h_se <- unname(se["h"])
if (is.finite(h) && h < 0) {
tmp <- lo; lo <- hi; hi <- tmp
tmp <- lo_se; lo_se <- hi_se; hi_se <- tmp
h <- -h
}
lo_ci <- c(lo - tq * lo_se, lo + tq * lo_se)
hi_ci <- c(hi - tq * hi_se, hi + tq * hi_se)
h_ci <- c(h - tq * h_se, h + tq * h_se)
ec50 <- 10^le50
ec50_lo <- 10^(le50 - tq * le50_se)
ec50_hi <- 10^(le50 + tq * le50_se)
increasing <- hi > lo
ec_label <- if (increasing) "EC50" else "IC50"
dir_word <- if (increasing) "rises" else "falls"
span_obs <- hi - loHill slope by the usual sign convention: negative for an inhibition curve.
hill_signed <- if (increasing) h else -h
hill_ci_signed <- if (increasing) h_ci else rev(-h_ci)
extrapolated <- ec50 < dose_min || ec50 > dose_maxStep 8: Were the plateaus actually observed, or are they projections?
lev <- sort(unique(u))
mu_lev <- sapply(lev, function(v) mean(y[u == v]))
k <- length(mu_lev)
tol_plateau <- 0.15 * abs(span_obs)
plateau_low_seen <- is.finite(tol_plateau) && tol_plateau > 0 &&
abs(mu_lev[1] - lo) <= tol_plateau
plateau_high_seen <- is.finite(tol_plateau) && tol_plateau > 0 &&
abs(mu_lev[k] - hi) <= tol_plateauStep 9: Residual diagnostics
fitted_v <- as.numeric(stats::fitted(fit))
resid_v <- as.numeric(stats::residuals(fit))
rss <- sum(resid_v^2)
tss <- sum((y - mean(y))^2)
rse <- sqrt(rss / max(1, df_res))
pseudo_r2 <- if (tss > 0) 1 - rss / tss else NA_real_Replicate-based lack-of-fit: pure error from repeated doses versus the remaining residual. This is the only honest test of curve shape when the design has replicates, and it is skipped (not faked) when it does not.
lof_p <- NA_real_; lof_note <- ""
df_pe <- n - n_levels
df_lof <- n_levels - 4
if (df_pe >= 1 && df_lof >= 1) {
ss_pe <- sum(sapply(unique(u), function(v) {
yy <- y[u == v]; sum((yy - mean(yy))^2)
}))
ss_lof <- rss - ss_pe
if (is.finite(ss_lof) && ss_lof > 0 && ss_pe > 0) {
f_lof <- (ss_lof / df_lof) / (ss_pe / df_pe)
lof_p <- stats::pf(f_lof, df_lof, df_pe, lower.tail = FALSE)
lof_note <- sprintf("Lack-of-fit F(%d, %d) = %s, %s", df_lof, df_pe,
r3(f_lof), fmt_pp(lof_p))
} else {
lof_note <- "The lack-of-fit test could not be computed from these replicates."
}
} else {
lof_note <- sprintf("No lack-of-fit test: it needs replicate measurements at repeated '%s' values and more than four distinct levels.", dose_h)
}
shapiro_p <- NA_real_
if (n >= 3 && n <= 5000) {
shapiro_p <- tryCatch(stats::shapiro.test(resid_v)$p.value,
error = function(e) NA_real_)
}Step 10: Competing models — a baseline-constrained 3PL and log-linear
The 3PL holds the zero-dose plateau at the UNTREATED CONTROL mean, so it is only defined when the data actually contain zero-dose rows. When they do not, it is reported as not applicable rather than anchored on a number invented for the occasion.
aic4 <- stats::AIC(fit)
fit3 <- NULL; aic3 <- NA_real_; rss3 <- NA_real_; df3 <- NA_real_
status3 <- ""
if (n_zero > 0 && is.finite(zero_mean)) {
b_fixed <- zero_mean
d3 <- data.frame(u = u, y = y)
fit3 <- tryCatch(
stats::nls(y ~ f4pl(u, b_fixed, hi, le50, h), data = d3,
start = list(hi = hi, le50 = le50, h = h),
algorithm = "port",
lower = c(hi = -Inf, le50 = min(u) - 3, h = 0.05),
upper = c(hi = Inf, le50 = max(u) + 3, h = 25),
control = stats::nls.control(maxiter = 200, warnOnly = FALSE)),
error = function(e) NULL, warning = function(w) NULL)
if (!is.null(fit3)) {
aic3 <- stats::AIC(fit3)
rss3 <- sum(stats::residuals(fit3)^2)
df3 <- stats::df.residual(fit3)
status3 <- sprintf("fitted with the low-dose plateau held at the zero-dose control mean(%s)",
fmt_num(b_fixed))
} else {
status3 <- "did not converge with the low-dose plateau held at the zero-dose control mean"
}
} else {
status3 <- sprintf("not applicable — the data contain no zero-dose '%s' rows to anchor an untreated baseline on",
dose_h)
}
fit_ll <- tryCatch(stats::lm(y ~ u), error = function(e) NULL)
aic_ll <- if (!is.null(fit_ll)) stats::AIC(fit_ll) else NA_real_
rss_ll <- if (!is.null(fit_ll)) sum(stats::residuals(fit_ll)^2) else NA_real_
r2_ll <- if (!is.null(fit_ll)) summary(fit_ll)$r.squared else NA_real_Nested F test, 4PL versus 3PL (the 3PL is the 4PL with one plateau fixed).
f_3pl_p <- NA_real_; f_3pl_note <- ""
if (!is.null(fit3) && is.finite(rss3) && is.finite(df3) && df3 > df_res) {
num <- (rss3 - rss) / (df3 - df_res)
den <- rss / df_res
if (is.finite(num) && is.finite(den) && den > 0 && num >= 0) {
f_stat <- num / den
f_3pl_p <- stats::pf(f_stat, df3 - df_res, df_res, lower.tail = FALSE)
f_3pl_note <- sprintf("F(%d, %d) = %s, %s", df3 - df_res, df_res,
r3(f_stat), fmt_pp(f_3pl_p))
}
}
aics <- c(fourpl = aic4, threepl = aic3, loglin = aic_ll)
aic_ok <- aics[is.finite(aics)]
best_model <- if (length(aic_ok) > 0) names(aic_ok)[which.min(aic_ok)] else "fourpl"
best_label <- switch(best_model,
fourpl = "the four-parameter logistic",
threepl = "the baseline-constrained three-parameter logistic",
loglin = "the plain log-linear model")
models_df <- data.frame(
model = c("4PL (four-parameter logistic)",
"3PL (baseline held at the untreated control)",
"Log-linear(response on log10 dose)"),
parameters = c(4L, 3L, 2L),
residual_se = round(c(rse,
if (!is.null(fit3)) sqrt(rss3 / max(1, df3)) else NA_real_,
if (!is.null(fit_ll)) summary(fit_ll)$sigma else NA_real_), 4),
aic = round(c(aic4, aic3, aic_ll), 2),
status = c("fitted", status3,
if (!is.null(fit_ll)) "fitted" else "did not fit"),
stringsAsFactors = FALSE
)
models_df$aic_vs_best <- round(models_df$aic - min(models_df$aic, na.rm = TRUE), 2)Step 11: Chart data — fitted curve drawn over the observed points.
The x axis is log10 dose because that is the scale the curve is symmetric on. Zero-dose controls have no log10 and are NOT silently placed at zero: they are drawn as their own series one decade below the lowest tested dose, and the prose says that position is a drawing convention.
set.seed(42)
obs_idx <- if (n > 900) sample(n, 900) else seq_len(n)
ctrl_x <- log10(dose_min) - 1
grid_u <- seq(min(u), max(u), length.out = 100)
if (has_group_col) {
keep_groups <- names(sort(table(g), decreasing = TRUE))
} else {
keep_groups <- "All observations"
}
chart_parts <- list()
chart_parts[[1]] <- data.frame(
plot_dose = round(u[obs_idx], 5),
response = round(y[obs_idx], 5),
series = if (has_group_col) paste0(g[obs_idx], " (observed)") else "Observed",
stringsAsFactors = FALSE
)
chart_parts[[2]] <- data.frame(
plot_dose = round(grid_u, 5),
response = round(f4pl(grid_u, lo, hi, le50, h), 5),
series = "Fitted 4PL curve",
stringsAsFactors = FALSE
)
if (n_zero > 0) {
zi <- if (n_zero > 150) sample(n_zero, 150) else seq_len(n_zero)
chart_parts[[3]] <- data.frame(
plot_dose = round(rep(ctrl_x, length(zi)), 5),
response = round(zero_resp[zi], 5),
series = "Zero-dose control(drawn one decade below the lowest dose)",
stringsAsFactors = FALSE
)
}Step 12: Per-group curves and a formal test of whether EC50s differ
n_groups <- 0L; group_test_p <- NA_real_; ec_ratio <- NA_real_
ec_ratio_label <- ""; group_note <- ""; group_excluded <- character(0)
group_df <- data.frame(
group = character(0), n = integer(0), dose_levels = integer(0),
ec50 = numeric(0), ec50_low = numeric(0), ec50_high = numeric(0),
hill = numeric(0), plateau_low = numeric(0), plateau_high = numeric(0),
stringsAsFactors = FALSE)
if (has_group_col) {
tab <- table(g)
cand <- names(sort(tab, decreasing = TRUE))
if (length(cand) > 6) {
group_excluded <- c(group_excluded, cand[7:length(cand)])
cand <- cand[1:6]
}
rows <- list(); rss_parts <- c(); n_parts <- c(); fit_groups <- character(0)
for (gg in cand) {
sel <- g == gg
ug <- u[sel]; yg <- y[sel]
if (length(ug) < MIN_ROWS || length(unique(ug)) < MIN_LEVELS ||
!isTRUE(stats::var(yg) > 0)) {
group_excluded <- c(group_excluded, gg); next
}
stg <- dr_selfstart(ug, yg)
fg <- dr_fit4pl(ug, yg, stg)
if (is.null(fg)) { group_excluded <- c(group_excluded, gg); next }
cfg <- summary(fg)$coefficients
eg <- cfg[, 1]; sg <- cfg[, 2]
dfg <- stats::df.residual(fg); tqg <- stats::qt(0.975, max(1, dfg))
glo <- unname(eg["lo"]); ghi <- unname(eg["hi"])
gle <- unname(eg["le50"]); gh <- unname(eg["h"])
gle_se <- unname(sg["le50"])
if (is.finite(gh) && gh < 0) { tmp <- glo; glo <- ghi; ghi <- tmp; gh <- -gh }
rows[[length(rows) + 1]] <- data.frame(
group = gg, n = length(ug), dose_levels = length(unique(ug)),
ec50 = round(10^gle, 5),
ec50_low = round(10^(gle - tqg * gle_se), 5),
ec50_high = round(10^(gle + tqg * gle_se), 5),
hill = round(if (ghi > glo) gh else -gh, 4),
plateau_low = round(glo, 4), plateau_high = round(ghi, 4),
stringsAsFactors = FALSE)
rss_parts <- c(rss_parts, sum(stats::residuals(fg)^2))
n_parts <- c(n_parts, length(ug))
fit_groups <- c(fit_groups, gg)Each group's own fitted curve joins the chart.
gug <- seq(min(ug), max(ug), length.out = 100)
chart_parts[[length(chart_parts) + 1]] <- data.frame(
plot_dose = round(gug, 5),
response = round(f4pl(gug, glo, ghi, gle, gh), 5),
series = paste0(gg, " (fitted)"),
stringsAsFactors = FALSE)
}
if (length(rows) > 0) {
group_df <- do.call(rbind, rows)
group_df <- group_df[order(group_df$ec50), , drop = FALSE]
rownames(group_df) <- NULL
n_groups <- nrow(group_df)
}
if (n_groups >= 2) {Full model = one 4PL per group (fitting them separately is exactly the pooled model with every parameter group-indexed, so the residual sums add). Reduced model = one shared midpoint, everything else free.
sel_all <- g %in% fit_groups
uu <- u[sel_all]; yy <- y[sel_all]
gi <- as.integer(factor(g[sel_all], levels = fit_groups))
K <- n_groups
n_all <- length(uu)
rss_full <- sum(rss_parts); df_full <- n_all - 4 * K
dd <- data.frame(u = uu, y = yy, gi = gi)
start_red <- list(
lo = group_df$plateau_low[match(fit_groups, group_df$group)],
hi = group_df$plateau_high[match(fit_groups, group_df$group)],
le50 = mean(log10(group_df$ec50)),
h = abs(group_df$hill[match(fit_groups, group_df$group)]))
fit_red <- tryCatch(
stats::nls(y ~ f4pl(u, lo[gi], hi[gi], le50, h[gi]), data = dd,
start = start_red,
control = stats::nls.control(maxiter = 300, warnOnly = FALSE)),
error = function(e) NULL, warning = function(w) NULL)
if (!is.null(fit_red) && df_full >= 1) {
rss_red <- sum(stats::residuals(fit_red)^2)
df_red <- n_all - (3 * K + 1)
num <- (rss_red - rss_full) / (df_red - df_full)
den <- rss_full / df_full
if (is.finite(num) && is.finite(den) && den > 0 && num >= 0) {
f_g <- num / den
group_test_p <- stats::pf(f_g, df_red - df_full, df_full,
lower.tail = FALSE)
group_note <- sprintf("Extra-sum-of-squares F(%d, %d) = %s, %s",
df_red - df_full, df_full, r3(f_g),
fmt_pp(group_test_p))
}
}
if (!is.finite(group_test_p)) {
group_note <- sprintf("The shared-midpoint comparison model did not converge, so no formal test of whether the %s midpoints differ is reported; compare the per-group intervals in the table instead.",
group_h)
}
ec_ratio <- max(group_df$ec50) / min(group_df$ec50)
ec_ratio_label <- sprintf("%s versus %s",
group_df$group[which.max(group_df$ec50)],
group_df$group[which.min(group_df$ec50)])
} else {
group_note <- sprintf("Fewer than two '%s' levels had enough data to carry their own curve, so no across-group comparison is reported.",
group_h)
}
} else {
group_note <- sprintf("No grouping column was mapped, so one curve was fitted across all %s rows.",
format(n, big.mark = ","))
}
curve_df <- do.call(rbind, chart_parts)
curve_df <- curve_df[order(curve_df$series, curve_df$plot_dose), , drop = FALSE]
rownames(curve_df) <- NULLStep 13: Tables
params_df <- data.frame(
parameter = c(sprintf("Plateau at low '%s'", dose_h),
sprintf("Plateau at high '%s'", dose_h),
sprintf("%s(midpoint dose)", ec_label),
"Hill slope(steepness)",
"Span(high plateau minus low plateau)"),
estimate = round(c(lo, hi, ec50, hill_signed, span_obs), 4),
ci_low = round(c(lo_ci[1], hi_ci[1], ec50_lo, hill_ci_signed[1], NA_real_), 4),
ci_high = round(c(lo_ci[2], hi_ci[2], ec50_hi, hill_ci_signed[2], NA_real_), 4),
interpretation = c(
sprintf("The '%s' the curve settles to as '%s' approaches zero. %s",
resp_h, dose_h,
if (plateau_low_seen)
"The lowest tested doses reach this plateau, so it is observed rather than projected."
else
sprintf("The lowest tested dose(%s) has not reached this plateau, so this value is a projection beyond the data.",
fmt_num(dose_min))),
sprintf("The '%s' the curve saturates at as '%s' grows large. %s",
resp_h, dose_h,
if (plateau_high_seen)
"The highest tested doses reach this plateau, so it is observed rather than projected."
else
sprintf("The highest tested dose(%s) has not reached this plateau, so this value is a projection beyond the data.",
fmt_num(dose_max))),
sprintf("The '%s' at which '%s' is halfway between the two plateaus. %s",
dose_h, resp_h,
if (extrapolated)
sprintf("It falls OUTSIDE the tested range of %s to %s, so it is an extrapolation, not a measured midpoint.",
fmt_num(dose_min), fmt_num(dose_max))
else
sprintf("It falls inside the tested range of %s to %s, so it is interpolated from observed doses.",
fmt_num(dose_min), fmt_num(dose_max))),
sprintf("How sharply '%s' turns over near the midpoint. A steeper slope means a narrower '%s' window between little effect and most of the effect; the sign is %s because '%s' %s with '%s'.",
resp_h, dose_h,
if (increasing) "positive" else "negative", resp_h, dir_word, dose_h),
sprintf("The total achievable change in '%s' across the fitted curve, from %s to %s.",
resp_h, fmt_num(lo), fmt_num(hi))),
stringsAsFactors = FALSE
)
resid_idx <- if (n > 900) obs_idx else seq_len(n)
residual_df <- data.frame(
fitted_value = round(fitted_v[resid_idx], 5),
residual = round(resid_v[resid_idx], 5),
stringsAsFactors = FALSE
)
residual_df <- residual_df[order(residual_df$fitted_value), , drop = FALSE]
rownames(residual_df) <- NULL
methods_df <- data.frame(
item = c("Model", "Estimation", "Starting values", "Confidence intervals",
"Midpoint interval", "Zero-dose controls", "Model comparison",
"Lack of fit", "Convergence policy", "Causal standing"),
detail = c(
sprintf("Four-parameter logistic on log10 '%s': response = low plateau + (high plateau - low plateau) / (1 + 10^((log10(%s) - log10(dose)) x hill)).",
dose_h, ec_label),
sprintf("Ordinary nonlinear least squares(base R nls) over %s rows at %d distinct positive '%s' levels; residual degrees of freedom %d, residual standard error %s.",
format(n, big.mark = ","), n_levels, dose_h, df_res, fmt_num(rse)),
sprintf("Derived from the data, not hard-coded: the plateaus start at the mean '%s' at the lowest and highest tested dose, and normalising by them linearises the curve on the logit scale, so an ordinary least-squares line supplies the starting hill slope (%s) and midpoint (%s).",
resp_h, r3(st$h), fmt_num(10^st$le50)),
sprintf("Wald intervals, estimate plus or minus %s standard errors on %d degrees of freedom. These are asymptotic: for a nonlinear model they are slightly optimistic compared with profile-likelihood intervals.",
r2(tq), df_res),
sprintf("The midpoint is estimated on the log10 scale and back-transformed, so its interval(%s to %s) is asymmetric around %s — which is the correct shape for a dose.",
fmt_num(ec50_lo), fmt_num(ec50_hi), fmt_num(ec50)),
if (n_zero > 0)
sprintf("%s zero-dose control row(s) were found, mean '%s' %s. They cannot sit on a log-dose axis, so they are excluded from the curve fit and drawn as their own series one decade below the lowest tested dose — that position is a drawing convention, not a measured dose.",
format(n_zero, big.mark = ","), resp_h, fmt_num(zero_mean))
else
sprintf("No zero-dose '%s' rows were present, so there is no untreated control to anchor a baseline on.", dose_h),
sprintf("The 4PL is compared against a three-parameter logistic(%s) and a log-linear model on AIC; %s has the lowest AIC here.",
status3, best_label),
lof_note,
"If nls does not converge from the self-start, the bounded port algorithm and a grid of perturbed starts are tried. If all fail, the analysis stops and reports the failure with a diagnosis. A linear or log-linear fit is never substituted for the curve.",
sprintf("This is a fitted description of how '%s' varies with observed '%s'. Unless the doses were assigned experimentally, the curve is associated with dose and does not by itself establish that changing '%s' causes the change in '%s'.",
resp_h, dose_h, dose_h, resp_h)),
stringsAsFactors = FALSE
)
metrics <- list(
`Observations Fitted` = n,
`Dose Levels` = n_levels,
`Midpoint Type` = ec_label,
`Midpoint Dose` = round(ec50, 5),
`Midpoint CI Low` = round(ec50_lo, 5),
`Midpoint CI High` = round(ec50_hi, 5),
`Hill Slope` = round(hill_signed, 3),
`Low Plateau` = round(lo, 3),
`High Plateau` = round(hi, 3),
`Residual Std Error` = round(rse, 4),
`Variance Explained` = if (is.finite(pseudo_r2)) round(pseudo_r2, 4) else NA_real_,
`Best Model By AIC` = switch(best_model, fourpl = "4PL",
threepl = "3PL", loglin = "log-linear"),
`Midpoint Extrapolated` = if (extrapolated) "yes" else "no"
)
extrap_clause <- if (extrapolated) {
sprintf(" The midpoint lies OUTSIDE the tested '%s' range of %s to %s, so this %s is an extrapolation beyond the doses actually observed and should be treated as a projection, not a measurement.",
dose_h, fmt_num(dose_min), fmt_num(dose_max), ec_label)
} else ""
group_clause <- if (n_groups >= 2) {
if (is.finite(group_test_p) && group_test_p < 0.05) {
sprintf(" Across '%s', the midpoints differ (%s): %s is a %s-fold shift, so the pooled curve averages over genuinely different curves and should not be read as any one group's response.",
group_h, group_note, ec_ratio_label, fmt_num(ec_ratio))
} else if (is.finite(group_test_p)) {
sprintf(" Across '%s', the midpoints are not distinguishable (%s), so one shared curve is a fair summary of all groups.",
group_h, group_note)
} else {
paste0(" ", group_note)
}
} else ""
json_output <- list(
answer = paste0(
"Four-parameter logistic dose-response fit of '", resp_h, "' on '",
dose_h, "' across ", format(n, big.mark = ","), " rows at ", n_levels,
" distinct positive dose levels: '", resp_h, "' ", dir_word, " from a low-dose plateau of ",
fmt_num(lo), " to a high-dose plateau of ", fmt_num(hi), ", with ",
ec_label, " = ", fmt_num(ec50), " (95% CI ", fmt_num(ec50_lo), " to ",
fmt_num(ec50_hi), ") and a Hill slope of ", r3(hill_signed), ".",
extrap_clause, group_clause,
" The curve leaves a residual standard error of ", fmt_num(rse),
" and ", best_label, " has the lowest AIC of the three models compared."
),
cards = lapply(
c("tldr", "overview", "preprocessing", "dose_response_curve",
"parameters", "model_comparison", "group_comparison",
"residual_diagnostics", "methods"),
function(cid) list(id = cid, metrics = metrics)
)
)
list(
initial_rows = initial_rows, final_rows = final_rows,
rows_removed = rows_removed,
dose_h = dose_h, resp_h = resp_h, group_h = group_h,
has_group_col = has_group_col,
n = n, n_levels = n_levels, n_zero = n_zero, zero_mean = zero_mean,
n_drop_na = n_drop_na, n_drop_nonpos = n_drop_nonpos,
lo = lo, hi = hi, le50 = le50, hill = h, hill_signed = hill_signed,
lo_ci = lo_ci, hi_ci = hi_ci, hill_ci = hill_ci_signed,
ec50 = ec50, ec50_lo = ec50_lo, ec50_hi = ec50_hi,
ec_label = ec_label, increasing = increasing, dir_word = dir_word,
span_obs = span_obs, extrapolated = extrapolated,
dose_min = dose_min, dose_max = dose_max, ctrl_x = ctrl_x,
plateau_low_seen = plateau_low_seen, plateau_high_seen = plateau_high_seen,
rse = rse, pseudo_r2 = pseudo_r2, df_res = df_res,
lof_p = lof_p, lof_note = lof_note, shapiro_p = shapiro_p,
sp_rho = sp_rho, sp_p = sp_p,
models_df = models_df, best_model = best_model, best_label = best_label,
status3 = status3, f_3pl_p = f_3pl_p, f_3pl_note = f_3pl_note,
curve_df = curve_df, params_df = params_df, residual_df = residual_df,
group_df = group_df, methods_df = methods_df,
n_groups = n_groups, group_test_p = group_test_p, group_note = group_note,
ec_ratio = ec_ratio, ec_ratio_label = ec_ratio_label,
group_excluded = group_excluded,
metrics = metrics, json_output = json_output
)
}