How to run nlxb and wrapnls inside dplyr?

Question

I am trying to fit many nonlinear fits using wrapnls in parallel using dplyr and broom (and eventually mclapply), but I am getting a parsing evaluation error from nlxb:

Error in parse(text = joe) (from #11) : <text>:1:6: unexpected input
1: b1.10% <- 20

I get this error using both do and lapply approaches.

library(nlmrt)
library(dplyr)
library(purrr)
library(broom)

data_frame(x = seq(0, 200, 0.1),
           y = 1.2*exp(-(times - 10)^2/(2*4.2^2)) + 2.4*exp(-(times - 50)^2/(2*3.8^2)) + 5.3*exp(-(times - 80)^2/(2*5.1^2)) + rnorm(length(times), sd = 0.05)) %>%
  do({
    xl  <- quantile(.$x, 0.1, na.rm = TRUE)
    xm  <- quantile(.$x, 0.5, na.rm = TRUE)
    xh  <- quantile(.$x, 0.8, na.rm = TRUE)
    starts <- c(a1 = 5, a2 = 5, a3 = 5,
                b1 = xl, b2 = xm, b3 = xh,
                c1 = 5, c2 = 5, c3 = 5)
    fmla <- y ~ a1*exp(-(x - b1)^2/(2*c1^2)) + a2*exp(-(x - b2)^2/(2*c2^2)) + a3*exp(-(x - b3)^2/(2*c3^2))
    df <- data_frame(x = .$x, y = .$y)
    mod <- wrapnls(fmla, lower = 0, upper = 200, start = starts, data = df)
    tidy(mod)
  })

Is there any way around this?

Answer 1

The problem isn't with the do aspect, it's the code inside the do, so you can debug that part directly. The starts vector is getting the b# names concatenated with the quantiles:

names(starts)

## [1] "a1"     "a2"     "a3"     "b1.10%" "b2.50%" "b3.80%" "c1"     "c2"     "c3"  

Adding unname to the quantile calculation fixes the issue.

data_frame(x = seq(0, 200, 0.1),
           y = 1.2*exp(-(x - 10)^2/(2*4.2^2)) + 2.4*exp(-(x - 50)^2/(2*3.8^2)) + 5.3*exp(-(x - 80)^2/(2*5.1^2)) + rnorm(length(x), sd = 0.05)) %>%
  do({
    xl  <- quantile(.$x, 0.1, na.rm = TRUE) %>% unname()
    xm  <- quantile(.$x, 0.5, na.rm = TRUE) %>% unname()
    xh  <- quantile(.$x, 0.8, na.rm = TRUE) %>% unname()
    starts <- c(a1 = 5, a2 = 5, a3 = 5,
                b1 = xl, b2 = xm, b3 = xh,
                c1 = 5, c2 = 5, c3 = 5)
    fmla <- y ~ a1*exp(-(x - b1)^2/(2*c1^2)) + a2*exp(-(x - b2)^2/(2*c2^2)) + a3*exp(-(x - b3)^2/(2*c3^2))
    df <- data_frame(x = .$x, y = .$y)
    mod <- wrapnls(fmla, lower = 0, upper = 200, start = starts, data = df)
    tidy(mod)
  })

##   term  estimate   std.error  statistic p.value
## 1   a1  2.386492 0.007455097   320.1155       0
## 2   a2  5.296250 0.006437509   822.7174       0
## 3   a3  1.199384 0.007132559   168.1562       0
## 4   b1 49.997697 0.013702894  3648.6960       0
## 5   b2 80.004023 0.007150546 11188.5193       0
## 6   b3 10.077847 0.028644821   351.8209       0
## 7   c1  3.798829 0.013702940   277.2273       0
## 8   c2  5.094727 0.007150573   712.4921       0
## 9   c3  4.175235 0.028944448   144.2499       0