Nonlinear regression - Modeling with R series

Non linear regression

Options
- Linearize by transforming the data
- Fit polynomials or complex splines
- Fit nonlinear functions

Libraries

library(deSolve)
library(nlme)
library(plotly)
#library(bbmle)
library(reshape2)

Using nonlinear least squares

Simulate data

set.seed(123)
x <- seq(0, 50, 1)
y <- ((runif(1, 10, 20) * x) / (runif(1, 0, 10) + x)) + rnorm(51, 0, 1)
plot(x, y,
     main = "Simulated data",
     xlab = "Independent variable",
     ylab = "Dependent variable",
     col = 9,
     las = 1)

Using Michaelis-Menten equation (enzymatic reaction rate)

\[v=\frac{V_{\text{max}}\left[S\right]}{K_M + \left[S\right]}\tag{1}\]

Start values initialized to $1$

m <- nls(y ~ a * x / (b + x), start = list(a = 1, b = 1))

cor(y, predict(m))

## [1] 0.9420521

plot(x, y,
     main = "Simulated data",
     xlab = "Independent variable",
     ylab = "Dependent variable",
     col = 9,
     las = 1)
lines(x, predict(m),
      lty = 2,
      col = "red",
      lwd = 3)

The right starting values are required
- Start by eyeballing the plot and using your knowledge of equations

set.seed(123)
y <- runif(1, 5, 15) * exp(-runif(1, 0.01, 0.05) * x) + rnorm(51, 0, 0.5)
plot(x, y,
     main = "Simulated data",
     xlab = "Independent variable",
     ylab = "Dependent variable",
     col = 9,
     las = 1)

Use these values:
- a is $y$-value at $x=0$
- b is the decay rate

Our equation is $ Y=Ae^{-BX}$ and $b_{\text{initial}}=-2\frac{\ln{2}}{a_{\text{initial}}}$

a_start <- 8
b_start <- -2 * log(2)/a_start

m <- nls(y ~ a * exp(-b * x), start = list(a = a_start,
                                           b = b_start))
cor(y ,predict(m))

## [1] 0.9758937

plot(x, y,
     main = "Simulated data",
     xlab = "Independent variable",
     ylab = "Dependent variable",
     col = 9,
     las = 1)
lines(x, predict(m),
      lty = 2,
      col = "red",
      lwd = 3)

Using SelfStarting function

Using logistic population growth model given by equation (2)

\[\frac{dN}{dt}=RN\left(1-\frac{N}{K}\right)\tag{2}\]

Simulate data

log_growth <- function(time, st, params){
  with(as.list(c(st, params)),{
    dN <- R * N * (1 - N/K)
    return(list(c(dN)))
  })
}

params <- c(R = 0.2, K = 1000)  # Parameters for logistic growth
N_ini <- c(N = 1)  # Initial values
times <- seq(0, 50, 1)  # Time steps to evaluate the ODE

out <- ode(N_ini, times, log_growth, params)  # ODE

set.seed(123)
N_obs <- out[, 2] + rnorm(51, 0, 50)  # Add random noise

N_obs <- ifelse(N_obs < 1, 1, N_obs) # Restrict smallest value of N_obs to be 1

plot(times, N_obs,
     main = "Simulated data",
     xlab = "Time",
     ylab = "Dependent variable",
     col = 9,
     las = 1)

R has a built-in function to estimate the starting values of the parameters of a logistic equation SSLogis

\[N_t = \frac{\alpha}{1 + e^{\frac{x_{\text{mid}}-t}{\text{scale}}}}\tag{3}\]

Use getInitial which calculates some initial guesses for the parameters given the data
Pass a self-starting model to this function

ss <- getInitial(N_obs ~ SSlogis(times,
                                 alpha,
                                 xmid,
                                 scale),
                 data = data.frame(N_obs = N_obs,  # Pass data set as a data frame
                                   times = times))

ss

##      alpha       xmid      scale 
## 985.747083  34.169811   4.810827

SSLogis uses different parametrization than our problem, so some algebra is required

K_start <- ss["alpha"]
R_start <- 1 / ss["scale"]
N0_start <- ss["alpha"] / (exp(ss["xmid"] / ss["scale"]) + 1)

Formula $N_{t}=\frac{KN_{0}e^{Rt}}{K+N_{0}e^{Rt}-1}$

log_formula <- formula(N_obs ~ K * N0 * exp(R * times) / (K + N0 * (exp(R * times) - 1)))  # Logistic growth

m<-nls(log_formula,
       start = list(K = K_start,
                    R = R_start,
                    N0 = N0_start))
summary(m)

## 
## Formula: N_obs ~ K * N0 * exp(R * times)/(K + N0 * (exp(R * times) - 1))
## 
## Parameters:
##           Estimate Std. Error t value Pr(>|t|)    
## K.alpha  985.74708   27.58937  35.729   <2e-16 ***
## R.scale    0.20786    0.01366  15.216   <2e-16 ***
## N0.alpha   0.81049    0.34976   2.317   0.0248 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 42.73 on 48 degrees of freedom
## 
## Number of iterations to convergence: 0 
## Achieved convergence tolerance: 1.145e-06

cor(N_obs, predict(m))

## [1] 0.993136

plot(times, N_obs,
     main = "Simulated data",
     xlab = "Time",
     ylab = "Dependent variable",
     col = 9,
     las = 1)
lines(times, predict(m),
      col = "red",
      lty = 2,
      lwd = 3)

Effect of predictors on non-linear regression

Use logistic growth function with parameters $r$ being the growth rate, $K$ being the equilibrium population size, and $n_0$ being the initial population size
Three experiments (three conditions)

Ks <- c(100, 200, 150) 
n0 <- c(5, 5, 6)
r <- c(0.15, 0.2, 0.15)
time <- 1:50

logF <- function(time, K, n0, r){
  d <- K * n0 * exp(r*time) / (K + n0 * (exp(r*time) - 1))
  return(d)
}

set.seed(123)
dat <- data.frame(Group = character(),
                  Time = numeric(),
                  Count = numeric())

for(i in 1:3){
  Ab <- logF(time = time, K = Ks[i], n0 = n0[i], r = r[i])
  tmp <- data.frame(Group = paste0("G", i), Time = time,
                    Count = Ab + rnorm(time, 0, 5))
  dat <-rbind(dat, tmp)
}
head(dat)

##   Group Time     Count
## 1    G1    1  2.960164
## 2    G1    2  5.482371
## 3    G1    3 15.418436
## 4    G1    4  9.103423
## 5    G1    5 10.671537
## 6    G1    6 20.036875

plot_ly(data = dat,
        x = ~Time,
        y = ~Count,
        color = ~Group,
        type = "scatter",
        mode = "markers") %>% 
  layout(title = "Growth over time for three groups",
         yaxis = list(zeroline = F))

lF <- formula(Count ~ K * n0 * exp(r * Time) / (K + n0 * (exp(r * Time) - 1)) | Group)

m <- nlsList(lF, data = dat, start = list(K = 200, n0 = 5, r = 0.15))

dat$Pred <- predict(m)
head(dat)

##   Group Time     Count      Pred
## 1    G1    1  2.960164  6.320132
## 2    G1    2  5.482371  7.218905
## 3    G1    3 15.418436  8.234450
## 4    G1    4  9.103423  9.378669
## 5    G1    5 10.671537 10.663731
## 6    G1    6 20.036875 12.101777