1

I'm having trouble fitting my model to a bacteria growth curve. I made a sample of my data below.

Code:

df <- data.frame(Bacteria = "EC",
                 Strains = c("AR1", "AR1", "AR1", "AR1", "AR1", "AR1", 
                             "AR2", "AR2",  "AR2", "AR2",  "AR2", "AR2",
                             "AR1", "AR1", "AR1", "AR1", "AR1", "AR1", 
                             "AR2", "AR2",  "AR2", "AR2",  "AR2", "AR2"),
                 Time_Points = c(0,4,6,8,12,16,0,4,6,8,12,16,0,4,6,8,12,16,0,4,6,8,12,16),
                 CFU = c(6.000000e+04, 7.100000e+06, 7.000000e+08, 2.166667e+09, 3.633333e+09, 2.666667e+09, 
                         1.333333e+04, 2.733333e+07, 9.266667e+08, 3.066667e+09, 6.200000e+09, 2.900000e+09,
                         0.000000e+00, 7.033333e+06, 3.766667e+08, 1.966667e+09, 1.366667e+09, 1.666667e+09, 
                         0.000000e+00, 4.866667e+06, 1.866667e+08, 7.800000e+08, 9.333333e+08, 1.733333e+09
                 ))
View(df)

#Make a subset of df of 1 strain

df.AR1 <- subset(df, Strains == "AR1")

View(df.AR1)

plot(df.AR1$Time_Points, df.AR1$CFU)

#make a function

k=4.510^9 N0=12 r=0.5 e_func <- function(t) k/(1 + ((k-N0)/N0)exp(-r*t))

time_ec <- (df.AR1$Time_Points) ECmodel.counts <- e_func(time_ec) + df.AR1$CFU

fit.ECmodel <- nls(ECmodel.counts~ k/(1 + ((k-N0)/N0)exp(-rTime_Points)), start=list(k= 4000000000, N0=1, r=0.05), data= df.AR1, trace=TRUE)

#Output: #3.415235e+19 (2.09e+00): par = (4e+09 1 0.05) #Error in nls(ECmodel.counts ~ k/(1 + ((k - N0)/N0) * exp(-r * Time_Points)), :

singular gradient

I've tried changing my starting parameters, but the same error occurs. I'm not sure where else to troubleshoot. It might be my "Time_Points," but I'm unsure.

Thank you!

Edit: singular gradient error solved

I've used the suggestions below, and I can fit the cure. However, I run into another problem with confidence intervals

Using code from Sal Mangiafico:

fit.ECmodel <- nls(ECmodel.counts~ k/(1 + ((k-N0)/N0)*exp(-r*Time_Points)), 
                   start=list(k= 2.332e+09, N0=3.769e+04, r=1.637e+00), 
                   data= df.AR1, trace=TRUE)

summary(fit.ECmodel)

plot(df.AR1$Time_Points, ECmodel.counts) lines(df.AR1$Time_Points, predict(fit.ECmodel))

confint(fit.ECmodel)

Output: > confint(fit.ECmodel) Waiting for profiling to be done... 3.437809e+18 (1.15e-01): par = (37685.5 1.637133) 3.433077e+18 (1.08e-01): par = (27873.03 1.681181) 3.428443e+18 (1.01e-01): par = (21384.36 1.721428) 3.423988e+18 (9.46e-02): par = (16878.15 1.758512) 3.419835e+18 (8.78e-02): par = (13627.03 1.792894) 3.419453e+18 (8.70e-02): par = (8793.148 1.85692) 3.414117e+18 (7.74e-02): par = (6436.856 1.908924) 3.408836e+18 (6.66e-02): par = (4981.086 1.953969) 3.407693e+18 (6.40e-02): par = (3023.672 2.032921) 3.398734e+18 (3.81e-02): par = (1743.85 2.135082) 3.393806e+18 (1.64e-03): par = (1563.866 2.172616) 3.393799e+18 (7.78e-04): par = (1456.776 2.184367) 3.393797e+18 (1.32e-06): par = (1459.651 2.184449) 2.176242e+22 (7.52e+01): par = (-34692.05 2.730645) Error in prof$getProfile() : singular gradient

Edit2: Confidence Interval Issues Solved with the below code library(nlstools); confint2(fit.ECmodel, level = 0.95)

KJack718
  • 13
  • 3
  • 1
    As a software solution, you might try the nlmrt package. It uses a different fitting algorithm that sometimes works better than nls in difficult cases. – Sal Mangiafico Mar 14 '23 at 15:26
  • 2
    This code has numerous errors that make it unexecutable! "view" should be "View" and it accesses df.AR1 before that data frame is created. As far as getting nls to run goes, standardize your data first. You unnecessarily create serious numerical problems by trying to relate data on a scale of $10^9$ to data on a scale of $1.$ Simply dividing ECmodel.counts by 1e9 solves the problem. – whuber Mar 14 '23 at 16:00
  • Apologies for the code error, but your solution to divide ECmodel.counts by 1e9 doesn't solve the problem; while I do get some results I still get a singular gradient error. – KJack718 Mar 14 '23 at 16:24
  • 1
    My system does not encounter that error. – whuber Mar 14 '23 at 17:12
  • It runs now; my starting values were too big. Thank you. – KJack718 Mar 14 '23 at 17:21
  • Are you just trying to get the confidence intervals for the parameters ? You might try library(nlstools); confint2(fit.ECmodel, level = 0.95) – Sal Mangiafico Mar 14 '23 at 18:59
  • For your plot, you're better off using many points for x which can be passed to predict(). For a simple example, try: library(rcompanion); plotPredy(data=df.AR1, x=Time_Points, y=ECmodel.counts, model=fit.ECmodel) – Sal Mangiafico Mar 14 '23 at 19:03
  • 1
    This is very helpful! Thank you for introducing me to these new packages and helping me solve my problems! – KJack718 Mar 14 '23 at 19:06

1 Answers1

3

It's basically a problem with the starting values, and the limitations of the fitting algorithm in nls().

The following works.

library(minpack.lm)

model = nlsLM(ECmodel.counts~ k/(1 + ((k-N0)/N0)exp(-rTime_Points)), start=list(k= 4000000000, N0=1, r=0.05), data= df.AR1, trace=TRUE)

summary(model)

Or, knowing the correct starting values will work with nls().

fit.ECmodel <- nls(ECmodel.counts~ k/(1 + ((k-N0)/N0)*exp(-r*Time_Points)), 
                   start=list(k= 2.332e+09, N0=3.769e+04, r=1.637e+00), 
                   data= df.AR1, trace=TRUE)

summary(fit.ECmodel)

Sal Mangiafico
  • 11,330
  • 2
  • 15
  • 35
  • How did you know what starting values to use? – KJack718 Mar 14 '23 at 17:21
  • In this case I got them from nlsLM() ! [Using your starting values]. ... But here's my general advice. You can set up a plot in Excel or LibreOffice Calc where you have cells with the parameters (starting values), and then have a column of x values and calculated y values based on your function and the cells with the parameters. Your plot will be essentially hot linked to these columns, so that you can change the parameter cells, and the plot will update. You can play with the parameter values until you get a reasonable starting fit. – Sal Mangiafico Mar 14 '23 at 18:45
  • Thank you for the advice! I will definitely use it. – KJack718 Mar 14 '23 at 18:55