CODEWITHSUNDEEP
rloopscurve-fittingnls
yaynikkiprograms
yaynikkiprograms

Reputation: 109

loop function with nls

I am struggling with looping nls functions. So here is an example data set for a single sample

dat<-read.table(text="time y
1 4.62
2 13.55
3 30.82
6 93.97
12 145.93
24 179.93", header = TRUE)
plot(data);lines(data)
model <- nls(y ~ Max * (1-exp(-k * (time - Lag))),data=dat,start=list(Max = 200, k = 0.1, Lag = 0.5))

but what if I want to apply model to multiple columns of samples? for example

dat<-read.table(text="time gluc starch solka
+ 1 6.32 7.51 1.95
+ 2 20.11 25.49 6.43
+ 3 36.03 47.53 10.39
+ 6 107.52 166.31 27.01
+ 12 259.28 305.19 113.72
+ 24 283.40 342.56 251.14
+ 48 297.55 353.66 314.22", header = TRUE)

How can I get R to solve for Max, k, and Lag for each sample (gluc, starch, solka)?

Upvotes: 0

Views: 618

Answers (2)

G. Grothendieck
G. Grothendieck

Reputation: 270075

In all the alternatives below we use these values:

long <- tidyr::pivot_longer(dat, -1, values_to = "y")
long$name <- factor(long$name)
st0 <- list(Max = 200, k = 0.1, Lag = 0.5)

1) nls grouped data Convert dat to long form and then use the grouped data feature of nls This solution is the most suitable among those presented here for testing whether certain parameters are common among the three names since it is easy to simply remove the subscript on a parameter if it is to be common among the names. The fitting itself does not use any packages but we show ggplot2 and lattice package graphics for plotting.

# get better starting values
model0 <- nls(y ~ Max * (1-exp(-k * (time - Lag))), long, start = st0)  
st <- with(as.list(coef(model0)), 
  list(Max = rep(Max, 3), k = rep(k, 3), Lag = rep(Lag, 3)))

model <- nls(y ~ Max[name] * (1-exp(-k[name] * (time - Lag[name]))), 
  long, start = st)
model

giving:

Nonlinear regression model
  model: y ~ Max[name] * (1 - exp(-k[name] * (time - Lag[name])))
   data: long
     Max1      Max2      Max3        k1        k2        k3      Lag1      Lag2 
306.48737 389.84657 361.82290   0.12214   0.03857   0.13747   1.38072   2.02205 
     Lag3 
  1.31770 
 residual sum-of-squares: 7167

Number of iterations to convergence: 8 
Achieved convergence tolerance: 9.186e-06

ggplot2 graphics could be done like this.

library(ggplot2)
fitdf <- transform(long, fit = fitted(model))
ggplot(fitdf, aes(x = time, y = y, color = name)) +
    geom_point() +
    geom_line(aes(y = fit))

A slightly different looking plot can be generated using lattice graphics which comes with R so the package does not have to be installed. The code is particularly compact.

library(lattice)
xyplot(fit + y ~ time | name, fitdf, type = c("l", "p"), auto.key = TRUE)

2) nlsList If you don't need to investigate common settings for parameters among the names then another possibility is to use nlsList in the nlme package (which comes with R so you don't have to install it). long and st0 are from above.

library(nlme)
fit <- nlsList(y ~ Max * (1-exp(-k * (time - Lag))) | name, long, start = st0)

giving an nlsList object whose 3 components are the three nls objects obtained by running nls for each name.

> fit
Call:
  Model: y ~ Max * (1 - exp(-k * (time - Lag))) | name 
   Data: long 

Coefficients:
            Max          k      Lag
gluc   306.4875 0.12214330 1.380713
solka  389.8449 0.03856544 2.022057
starch 361.8231 0.13747402 1.317698

Degrees of freedom: 21 total; 12 residual
Residual standard error: 24.43858

We can plot the data and fit:

levs <- levels(long$name)
col <- setNames(rainbow(length(levs)), levs) 
plot(y ~ time, long, col = col[name], pch = 20, cex = 1.5)
for(lv in levs) lines(fitted(fit[[lv]]) ~ time, dat, col = col[lv])
legend("bottomright", leg = levs, col = col, pch = 20, cex = 1.5)

screenshot

3) subset An approach which is similar to (2) is to perform three nls runs using subset= to select the data. This returns a named list of nls objects. st0 and long are from above. No packages are used.

fit <- Map(function(nm) nls(y ~ Max * (1-exp(-k * (time - Lag))), data = long, 
    start = st0, subset = name == nm), levels(long$name))

The graphics code in (2) also works here.

Upvotes: 1

Gregor Thomas
Gregor Thomas

Reputation: 146110

Build the formulas you want to use as strings:

outcomes = c("gluc", "starch", "solka")
my_formulas = paste(outcomes, "~ Max * (1-exp(-k * (time - Lag)))")
model_list = list()

for(i in seq_along(outcomes)) {
  model_list[[outcomes[i]]] = nls(
    as.formula(my_formulas[i], 
    data = dat,
    start = list(Max = 200, k = 0.1, Lag = 0.5)
  )
}

This will create a list of models, you can the access with, e.g., summary(model_list[[1]]) or summary(model_list[["solka"]])

Upvotes: 1

Related Questions