## 1. The DNase data from 'nls', use all generic functions.
DNase1 <- subset(DNase, Run == 1)
set.seed(123)
DNase1$density <- sapply(DNase1$density, function(x) rnorm(1, x, 0.1 * x))
mod1 <- onls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),
data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))
print(mod1)
plot(mod1)
summary(mod1)
predict(mod1, newdata = data.frame(conc = 6))
logLik(mod1)
deviance(mod1)
formula(mod1)
weights(mod1)
df.residual(mod1)
fitted(mod1)
residuals(mod1)
vcov(mod1)
coef(mod1)
## 2a. Update model
DNase2 <- DNase1
DNase2$conc <- DNase2$conc * 2
mod2a <- update(mod1, data = DNase2)
print(mod2a)
## 2b. Example with a fixed parameter
## => Asym = 3.
mod2b <- onls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)),
data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1),
fixed = c(TRUE, FALSE, FALSE))
print(mod2b)
## 3. Multivariate example: matched curvature,
## low noise, decorrelated predictors
set.seed(123)
n <- 25
x1 <- runif(n, 1, 5)
x2 <- runif(n, 1, 5)
b1_true <- 5
b2_true <- 2
b3_true <- 1.5
z_true <- b1_true + b2_true * x1 + b3_true * x2^2
z <- z_true + rnorm(n, 0, 2)
x1 <- x1 + rnorm(n, 0, 0.5)
x2 <- x2 + rnorm(n, 0, 0.5)
DAT <- data.frame(x1 = x1, x2 = x2, z = z)
mod3 <- onls(z ~ b1 + b2 * x1 + b3 * x2^2, data = DAT,
start = list(b1 = 1, b2 = 1, b3 = 1), trace = TRUE)
print(mod3)
# \donttest{
## Reference tests comparing to pivotal literature
## 4. Example from odrpack_guide.pdf, 2.C.i, pages 39ff.
x <- c(0, 0, 5, 7, 7.5, 10, 16, 26, 30, 34, 34.5, 100)
y <- c(1265, 1263.6, 1258, 1254, 1253, 1249.8, 1237, 1218, 1220.6,
1213.8, 1215.5, 1212)
DAT <- data.frame(x, y)
mod4 <- onls(y ~ b1 + b2 * (exp(b3 * x) -1)^2, data = DAT,
start = list(b1 = 1500, b2 = -50, b3 = -0.1))
deviance_o(mod4) # 21.445 as on page 47
summary(mod4) # 1264.65481 (1.03492) / -54.01838 (1.583992) / -0.08785 (6.33222E-3) as on page 48
## 5. Example from Algorithm 676: ODRPACK, page 355 + 356.
x <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 85, 90, 95, 100, 105)
y <- c(4.14, 8.52, 16.31, 32.18, 64.62, 98.76, 151.13, 224.74, 341.35,
423.36, 522.78, 674.32, 782.04, 920.01)
DAT <- data.frame(x, y)
mod5 <- onls(y ~ b1 * 10^(b2 * x/(b3 + x)), data = DAT,
start = list(b1 = 1, b2 = 5, b3 = 100))
deviance_o(mod5) # 15.263 as on page 363
summary(mod5) # 4.4879 (0.56876) / 7.1882 (0.69504) / 221.8383 (37.2313) as on page 363
## 6. Example with bounds from simple_example.f90
## in https://www.netlib.org/toms/869.zip.
x <- c(0.982, 1.998, 4.978, 6.01)
y <- c(2.7, 7.4, 148.0, 403.0)
DAT <- data.frame(x, y)
mod6 <- onls(y ~ b1 * exp(b2 * x), data = DAT,
start = list(b1 = 2, b2 = 0.5),
lower = c(0, 0), upper = c(10, 0.9))
coef(mod6) # 1.4376 / 0.9 ## Different to reference 1.6334 / 0.9
deviance_o(mod6) # 0.1919 => lower RSS than original ODRPACK with 0.2674!
## 7. Example similar to Deming regression
## Comparison to XLstat
## https://help.xlstat.com/6650-run-deming-regression-compare-methods-excel
x <- c(9.8, 9.7, 10.7, 10.9, 12.4, 12.5, 12.8, 12.8, 12.9, 13.3,
13.4, 13.5, 13.7, 14.9, 15.2, 15.5)
y <- c(10.1, 11.4, 10.8, 11.3, 11.8, 12.1, 12.3, 13.6, 14.2, 14.4,
14.6, 15.3, 15.5, 15.8, 16.2, 16.5)
DAT <- data.frame(x, y)
mod7 <- onls(y ~ a + b * x, data = DAT, start = list(a = 2, b = 3))
print(mod7) ## -1.909 / 1.208 as on webpage
plot(mod7)
## 8. Linear multivariate model, using the closed-form Total Least Squares
## (TLS) solution from Golub & Van Loan (1980)
tls_fit <- function(X, y) {
X <- as.matrix(X)
n <- nrow(X); p <- ncol(X)
Xc <- scale(X, center = TRUE, scale = FALSE)
yc <- y - mean(y)
xbar <- colMeans(X); ybar <- mean(y)
Z <- cbind(Xc, yc)
SVD <- svd(Z)
v <- SVD$v[, p + 1L]
v_x <- v[1:p]; v_y <- v[p + 1L]
slope <- -v_x / v_y; intercept <- ybar - sum(slope * xbar)
list(intercept = intercept, slope = setNames(slope, colnames(X)),
singular_values = SVD$d)
}
set.seed(11)
n <- 40
x1_true <- runif(n, 0, 10)
x2_true <- runif(n, 0, 10)
b0_true <- 3; b1_true <- 1.5; b2_true <- -0.8
y_true <- b0_true + b1_true * x1_true + b2_true * x2_true
x1 <- x1_true + rnorm(n, 0, 0.5)
x2 <- x2_true + rnorm(n, 0, 0.5)
y <- y_true + rnorm(n, 0, 0.5)
DAT <- data.frame(x1 = x1, x2 = x2, y = y)
TLS <- tls_fit(DAT[, c("x1", "x2")], DAT$y)
mod8 <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT, start = list(b0 = 1, b1 = 1, b2 = 1))
TLS_vec <- c(b0 = TLS$intercept, b1 = TLS$slope[["x1"]], b2 = TLS$slope[["x2"]])
ONLS_vec <- coef(mod8)[c("b0", "b1", "b2")]
print(data.frame(TLS_closed_form = TLS_vec, onls = ONLS_vec,
abs_diff = abs(TLS_vec - ONLS_vec))) # all equal
## 9. Pearson (1901) / York (1966) -> "Pearson's data with York's weights"
# Intercept: 5.47991 (SE 0.29497)
# Slope: -0.48053 (SE 0.05799)
x <- c(0.0, 0.9, 1.8, 2.6, 3.3, 4.4, 5.2, 6.1, 6.5, 7.4)
y <- c(5.9, 5.4, 4.4, 4.6, 3.5, 3.7, 2.8, 2.8, 2.4, 1.5)
sd_x <- 1/sqrt(c(1000.0, 1000.0, 500.0, 800.0, 200.0, 80.0, 60.0, 20.0, 1.8, 1.0))
sd_y <- 1/sqrt(c(1.0, 1.8, 4.0, 8.0, 20.0, 20.0, 70.0, 70.0, 100.0, 500.0))
DAT <- data.frame(x = x, y = y)
mod9 <- onls(y ~ b0 + b1*x, data = DAT,
start = list(b0 = 5, b1 = -0.5),
sigma_x = sd_x, sigma_y = sd_y)
summary(mod9) # 5.47991 (0.29497) / -0.48053 (0.05799) as in paper
## 10. Daeron & Vermeesch (2024), Table 3 / Figure 2C toy example.
x <- c(9, 19, 31, 41)
y <- c(21, 31, 39, 49)
DAT <- data.frame(x = x, y = y)
mod10 <- onls(y ~ a + b * x, data = DAT, start = list(a = 10, b = 1), sigma_x = 1, sigma_y = 1)
summary(mod10) # 13.71 / 0.851 as in Table 3 of paper
## 11. Full predictor covariance (correlated predictor errors), compared to the closed-form
## generalized Total Least Squares (TLS) solution. Whitening the predictors with the Cholesky
## factor L of the precision matrix (t(L) %*% L = solve(Sigma)) and scaling y by 1/sigma_y
## turns the problem into plain TLS, which is solved by an SVD.
set.seed(123)
n <- 40
Sigma <- matrix(c(0.25, 0.15, 0.15, 0.16), 2) # correlation of predictor errors = 0.75
sigma_y <- 0.3
xt <- cbind(runif(n, 0, 10), runif(n, 0, 10))
E <- matrix(rnorm(2 * n), n) %*% chol(Sigma)
DAT <- data.frame(x1 = xt[, 1] + E[, 1], x2 = xt[, 2] + E[, 2],
y = 3 + 1.5 * xt[, 1] - 0.8 * xt[, 2] + rnorm(n, 0, sigma_y))
mod12 <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT,
start = list(b0 = 1, b1 = 1, b2 = 1), sigma_x = Sigma, sigma_y = sigma_y,
control = list(ftol = 1e-13, ptol = 1e-13))
L <- chol(solve(Sigma))
U <- as.matrix(DAT[, c("x1", "x2")]) %*% t(L) # whitened predictors
v <- DAT$y / sigma_y
Z <- cbind(scale(U, scale = FALSE), v - mean(v))
V <- svd(Z)$v[, 3]
w <- -V[1:2] / V[3]
gTLS <- c(sigma_y * (mean(v) - sum(w * colMeans(U))), sigma_y * drop(t(L) %*% w))
print(data.frame(gen_TLS = gTLS, onls = coef(mod12),
abs_diff = abs(gTLS - coef(mod12)))) # all equal
## 12. ODRPACK's separate weights WE (response) and WD (predictor).
## ODRPACK takes both as precisions (inverse variances), per observation. In onls(), WE is
## passed as 'weights' (with the default sigma_y = 1) and WD through sigma_x = 1/sqrt(WD)
## (for p > 1: an n x p matrix 1/sqrt(WD)). known_sigma = FALSE gives ODRPACK's scaling of
## the standard errors by the residual variance.
set.seed(123)
n <- 30
xt <- seq(0.5, 10, length.out = n)
WD <- runif(n, 0.5, 4) # predictor weights
WE <- runif(n, 0.5, 4) # response weights
x <- xt + rnorm(n, 0, 0.3/sqrt(WD))
y <- 2 * exp(-0.3 * xt) + 0.5 + rnorm(n, 0, 0.03/sqrt(WE))
DAT <- data.frame(x, y)
mod12 <- onls(y ~ a * exp(-b * x) + c, data = DAT, start = list(a = 1.5, b = 0.2, c = 0.3),
weights = WE, sigma_x = 1/sqrt(WD), known_sigma = FALSE)
summary(mod12)
# }
Run the code above in your browser using DataLab