Nonliner models

Ioannis Kosmidis

Professor of Statistics

Department of Statistics, University of Warwick

ioannis.kosmidis@warwick.ac.uk

ikosmidis.com ikosmidis ikosmidis_

19 April 2024

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.

in Tukey (1986). Sunset Salvo. The American Statistician, 40 (1), p. 74.

\[ Y_i = x_i^\top \beta + \epsilon_i \,, \quad \epsilon_i \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0, \sigma^2) \]

- Parameters enter the model through a linear combination of coefficients and covariates.

\[ Y_i = \eta(x_i, \beta) + \epsilon_i \,, \quad \epsilon_i \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0, \sigma^2) \]

Parameters and covariates are allowed to have a nonlinear relationship.

The linear model results for \(\eta(x_i, \beta) = x_i^\top \beta\).

*Physical parameters*that have particular meaning in the subject-area where the model comes from. Estimating the value of physical parameters, then, contributes to scientific understanding.*Tuning parameters*that do not necessarily have any physical meaning. Their presence is justified as a simplification of a more complex underlying system. The aim when estimating them is to make the model represent reality as well as possible.

*Mechanistically*: prior scientific knowledge is incorporated into building a mathematical model for the mean response. That model can often be complex and \(\eta(x, \beta)\) may not be available in closed form.*Phenomenologically*(empirically): a function \(\eta(x, \beta)\) may be posited that appears to capture the non-linear nature of the mean response.

`?SMPracticals::calcium`

The uptake of calcium (in nmoles per mg) at set times (in minutes) by \(27\) cells in “hot” suspension.

Assume that the rate of growth is proportional to the calcium remaining, i.e. \[ \frac{d \eta}{d t} = (\beta_0 - \eta) / \beta_1 \tag{1}\]

Solving the differential equation (1) with initial condition \(\eta(0,\beta) = 0\), gives \[ \eta(t, \beta) = \beta_0 \left( 1 - \exp \left( - t / \beta_1 \right) \right) \]

\(\beta_0\) is the calcium uptake after infinite time.

\(\beta_1\) controls the growth rate of calcium uptake.

```
calc_lm1 <- lm(cal ~ time, data = calcium)
calc_lm2 <- lm(cal ~ time + I(time^2), data = calcium)
calc_nlm <- nls(cal~ beta0 * ( 1 - exp(-time/beta1)), data = calcium,
start = list(beta0 = 5, beta1 = 5))
newdata <- data.frame(time = seq(min(calcium$time), max(calcium$time), length.out = 100))
pred_lm1 <- predict(calc_lm1, newdata = newdata)
pred_lm2 <- predict(calc_lm2, newdata = newdata)
pred_nlm <- predict(calc_nlm, newdata = newdata)
plot(cal ~ time, data = calcium,
xlab = "Time (minutes)",
ylab = "Calcium uptake (nmoles/mg)",
bg = "#ff7518", pch = 21)
lines(newdata$time, pred_lm1, col = gray(0.8), lty = 1, lwd = 2)
lines(newdata$time, pred_lm2, col = gray(0.6), lty = 2, lwd = 2)
lines(newdata$time, pred_nlm, col = gray(0.4), lty = 3, lwd = 2)
legend("bottomright", legend = c("LM (linear)", "LM (quadratic)", "NLM"),
col = gray(c(0.8, 0.6, 0.4)), lty = 1:3, lwd = 2)
```

A comparison of the three models in terms of number of parameters, maximized log-likelihood value, and AIC and BIC returns

```
p loglik AIC BIC
LM (linear) 2 -28.701 63.403 67.290
LM (quadratic) 3 -20.955 49.910 55.093
NLM 2 -20.955 47.909 51.797
```

```
data("Theoph", package = "datasets")
plot(conc ~ Time, data = Theoph, type = "n",
ylab = "Concentration (mg/L)", xlab = "Time (hours)")
for (i in 1:30) {
dat_i <- subset(Theoph, Subject == i)
lines(conc ~ Time, data = dat_i, col = "grey")
}
points(conc ~ Time, data = Theoph,
bg = "#ff7518", pch = 21, col = "grey")
```

Theophylline is an anti-asthmatic drug. An experiment was performed on \(12\) individuals to investigate the way in which the drug leaves the body. The study of drug concentrations inside organisms is called *pharmacokinetics*.

An oral dose was given to each individual at time \(t = 0\), and the concentration of theophylline in the blood was then measured at 11 time points in the next 25 hours.

Compartmental models are a common class of model used in pharmacokinetics studies.

If the initial dosage is \(D\), then pharmacokinetic model with a first-order compartment function is \[ \eta(\beta, D,t) = \frac{D \beta_1 \beta_2}{\beta_3(\beta_2 - \beta_1)} \left( \exp \left( - \beta_1 t\right) - \exp \left( - \beta_2 t\right)\right)\ \tag{2}\] where

\(\beta_1 > 0\): the elimination rate which controls the rate at which the drug leaves the organism.

\(\beta_2 > 0\): the absorption rate which controls the rate at which the drug enters the blood.

\(\beta_3 > 0\): the clearance which controls the volume of blood for which a drug is completely removed per time unit.

Since all the parameters are positive, and their estimation will most probably require a gradient descent step (e.g. what some of the methods in `optim`

do), it is best to rewrite (2) in terms of \(\gamma_i = log(\beta_i)\). We can write \[
\eta'(\gamma, D, t) = \eta(\beta, D, t) = D \frac{\exp(-\exp(\gamma_1) t) - \exp(-\exp(\gamma_2) t)}{\exp(\gamma_3 - \gamma_1) - \exp(\gamma_3 - \gamma_2)}
\] where

\(\gamma_1 \in \Re\): the logarithm of the elimination rate.

\(\gamma_2 \in \Re\): the logarithm of the absorption rate.

\(\gamma_3 \in \Re\): the logarithm of the clearance.

We fit the model with predictor \(\eta'(\gamma, D_i,t_{ij})\) using nonlinear least-squares (`nls()`

in R).

```
Formula: conc ~ Dose * (exp(-exp(gamma1) * Time) - exp(-exp(gamma2) *
Time))/(exp(gamma3 - gamma1) - exp(gamma3 - gamma2))
Parameters:
Estimate Std. Error t value Pr(>|t|)
gamma1 0.39922 0.11754 3.397 0.000908 ***
gamma2 -2.52424 0.11035 -22.875 < 2e-16 ***
gamma3 -3.24826 0.07439 -43.663 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.459 on 129 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 3.709e-06
```

The estimates and estimated standard errors (using the delta method) for \(\beta_1\), \(\beta_2\) and \(\beta_3\) are

```
res <- residuals(pkm, type = "pearson")
ord <- order(ave(res, Theoph$Subject))
subj <- Theoph$Subject[ord]
subj <- factor(subj, levels = unique(subj), ordered = TRUE)
plot(res[ord] ~ subj,
xlab = "Subject (ordered by mean residual)",
ylab = "Pearson residual",
col = "#ff7518", pch = 21)
abline(h = 0, lty = 2)
```

Clear evidence of unexplained differences between individuals.

```
library("ggplot2")
st <- unique(Theoph[c("Subject", "Dose")])
pred_df <- as.list(rep(NA, nrow(st)))
for (i in seq.int(nrow(st))) {
pred_df[[i]] <- data.frame(Time = seq(0, 25, by = 0.2),
Dose = st$Dose[i],
Subject = st$Subject[i])
}
pred_df <- do.call("rbind", pred_df)
pred_df$conc <- predict(pkm, newdata = pred_df)
## Order according to mean residual
theoph <- within(Theoph, Subject <- factor(Subject, levels = unique(subj), ordered = TRUE))
fig_theoph <- ggplot(theoph) +
geom_point(aes(Time, conc), col = "#ff7518") +
geom_hline(aes(yintercept = Dose), col = "grey", lty = 3) +
facet_wrap(~ Subject, ncol = 3) +
labs(y = "Concentration (mg/L)", x = "Time (hours)") +
theme_bw()
fig_theoph +
geom_line(data = pred_df, aes(Time, conc), col = "grey")
```

Accounting for heterogeneity between individuals seems worthwhile.

\[ Y_{ij} = \eta(\beta + A b_i, x_{ij}) + \epsilon_{ij} \,, \quad \epsilon_{ij} \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0,\sigma^2)\,, \quad b_i \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0,\Sigma_b) \] where \(\Sigma_b\) is a \(q \times q\) covariance matrix and \(A\) is a \(p \times q\) matrix of zeros and ones, which determines which parameters are fixed and which are varying.

The linear mixed model is a special case of the nonlinear mixed model with \[ \eta(\beta + A b_i, x_{ij}) = x_{ij}^\top \left( \beta + A b_i\right) = x_{ij}^\top \beta + x_{ij}^\top A b_i = x_{ij}^\top \beta + z_{ij}^\top b_i \,. \]

A random intercept model results, if the first element of \(x_{ij}\) is \(1\) for all \(i\) and \(j\), \(q = 1\) and \(A = (1, 0, \ldots, 0)^\top\).

We fit a nonlinear mixed model that allows all the parameters to vary across individuals, i.e. \(A = I_3\) using `nmle()`

from the `nlme`

R package.

```
Nonlinear mixed-effects model fit by maximum likelihood
Model: fm
Data: Theoph
Log-likelihood: -173.32
Fixed: gamma1 + gamma2 + gamma3 ~ 1
gamma1 gamma2 gamma3
0.4514513 -2.4326850 -3.2144578
Random effects:
Formula: list(gamma1 ~ 1, gamma2 ~ 1, gamma3 ~ 1)
Level: Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
gamma1 0.6376932 gamma1 gamma2
gamma2 0.1310518 0.012
gamma3 0.2511873 -0.089 0.995
Residual 0.6818359
Number of Observations: 132
Number of Groups: 12
```

We fit a nonlinear mixed model that allows all the parameters to vary across individuals, i.e. \(A = I_3\) using `nmle()`

from the `nlme`

R package.

```
Nonlinear mixed-effects model fit by maximum likelihood
Model: fm
Data: Theoph
Log-likelihood: -173.32
Fixed: gamma1 + gamma2 + gamma3 ~ 1
gamma1 gamma2 gamma3
0.4514513 -2.4326850 -3.2144578
Random effects:
Formula: list(gamma1 ~ 1, gamma2 ~ 1, gamma3 ~ 1)
Level: Subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
gamma1 0.6376932 gamma1 gamma2
gamma2 0.1310518 0.012
gamma3 0.2511873 -0.089 0.995
Residual 0.6818359
Number of Observations: 132
Number of Groups: 12
```

Let’s consider the model with random effects with means \(\gamma_1\) and \(\gamma_3\), and just a population parameter for the logarithm of the absorption rate.

\[ Y_{ij} = \eta'\left( \begin{bmatrix} \gamma_1 + b_{i1} \\ \gamma_2 \\ \gamma_3 + b_{i3} \end{bmatrix}, D_i, t_{ij}\right) + \epsilon_{ij} \, \quad \epsilon_{ij} \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0,\sigma^2)\,, \quad (b_{i1}, b_{i3})^\top \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0,\Sigma_b) \] This corresponds to the general form pf the nonlinear mixed effects model with \[ A = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{bmatrix} \quad \text{and} \quad b_i = \begin{bmatrix} b_{i1} \\ b_{i3} \end{bmatrix} \]

A comparison to the model with all effects varying across individuals gives

```
conc_nlm <- pred_df$conc
conc_nlme_2 <- predict(pkmR_2, newdata = pred_df)
conc_nlme_3 <- predict(pkmR, newdata = pred_df)
pred_df_all <- pred_df[c("Subject", "Dose", "Time")]
pred_df_all <- rbind(
data.frame(pred_df_all, conc = conc_nlm, model = "NLM"),
data.frame(pred_df_all, conc = conc_nlme_2, model = "NLME(2)"),
data.frame(pred_df_all, conc = conc_nlme_3, model = "NLME(3)"))
fig_theoph +
geom_line(data = pred_df_all, aes(Time, conc, color = model))
```

The generalized nonlinear mixed effects model (GNLMM) assumes \[ Y_i \mid x_i, b_i \stackrel{\text{ind}}{\sim}\mathop{\mathrm{EF}}(\mu_i,\sigma^2)\,, \quad \begin{bmatrix} g(\mu_1)\\\vdots \\ g(\mu_n) \end{bmatrix} = \eta(\beta + A b_i, x_i) \,, \quad b_i \stackrel{\text{ind}}{\sim}{\mathop{\mathrm{N}}}(0,\Sigma_b) \] where \(\mathop{\mathrm{EF}}(\mu_i,\sigma^2)\) is an exponential family with mean \(\mu_i\) and variance \(\sigma^2 V(\mu_i) / m_i\).

Linear model, nonlinear model, linear mixed effects model, nonlinear mixed effects model, generalized linear model, and generalized nonlinear model.

As in GLMMs, the likelihood function may not have a closed form and needs approximation.

General-purpose optimizers may not converge to a global maximum of the likelihood.

Evaluating \(\eta(\beta, x)\) can be computationally expensive in some applications, like, for example, when \(\eta(\beta, x)\) is defined via a differential equation, which can only be solved numerically.

Materials | Link |
---|---|

Preliminary material | ikosmidis.com/files/APTS-SM-Preliminary |

Notes | ikosmidis.com/files/APTS-SM-Notes |

Slides: Introduction | ikosmidis.com/files/APTS-SM-Slides-intro |

Slides: Model selection | ikosmidis.com/files/APTS-SM-Slides-model-selection |

Slides: Beyond GLMs | ikosmidis.com/files/APTS-SM-Slides-beyond-glms |

Slides: Nonlinear models | ikosmidis.com/files/APTS-SM-Slides-nonlinear-models |

Slides: Latent variables | ikosmidis.com/files/APTS-SM-Slides-latent |

Lab 1 | ikosmidis.com/files/APTS-SM-Notes/lab1.html |

Lab 2 | ikosmidis.com/files/APTS-SM-Notes/lab2.html |

Ioannis Kosmidis - Statistical Modelling: Nonlinear models