Interpreting and troubleshooting nls in R with quadratic plateau model

Question

I am trying to run a quadratic plateau model on some proportion data where values are bound between 0 and 100. I would like some help troubleshooting some errors I have encountered, and correctly interpreting results as well as understanding the equation and how to write it out correctly. If anyone has experience with these models any help will be greatly appreciated as I have hit a wall.

Example data:

Days    Type    Area 
0   Abrasion    0
11  Abrasion    65.6513749
13  Abrasion    79.1887936
15  Abrasion    88.3947998
26  Abrasion    98.2726653
38  Abrasion    100
0   Abrasion    0
70  Abrasion    93.5047459
124 Abrasion    100
0   Abrasion    0
7   Abrasion    78.2666991
8   Abrasion    78.3624009
9   Abrasion    78.9448106
14  Abrasion    81.6443138
24  Abrasion    97.9969096
29  Abrasion    98.8788699
50  Abrasion    99.4708654
53  Abrasion    100
0   Laceration  0
8   Laceration  8.05965381
22  Laceration  67.1254163
83  Laceration  100
0   Laceration  0
8   Laceration  59.1650901
69  Laceration  96.1942307
74  Laceration  100
0   Laceration  0
49  Laceration  82.5396751
133 Laceration  100
0   Laceration  0
125 Laceration  100
0   Laceration  0
16  Laceration  48.5178133

X = Days Y = Area

I want to fit a quadratic plateau model to this data.

Code I am using:

###  Find reasonable initial values for parameters

fit.lm    = lm(Area ~ Days, data=healing)

a.ini     = fit.lm$coefficients[1]
b.ini     = fit.lm$coefficients[2]
clx.ini   = mean(healing$Area)


###  Define quadratic plateau function

quadplat = function(x, a, b, clx) {
  ifelse(x  < clx, a + b * x + (-0.5*b/clx) * x * x, 
         a + b * clx + (-0.5*b/clx) * clx * clx)}

###  Find best fit parameters


model = nls(Area ~ quadplat(Days, a, b, clx), 
            data = healing, 
            start = list(a   = a.ini, 
                         b   = b.ini, 
                         clx = clx.ini),
            trace = FALSE,
            nls.control(maxiter = 1000))

summary(model)

When I run this on some data it works fine but other times I get the following error:

Error in nls(Area ~ quadplat(Days, a, b, clx), data = healing,  : 
  singular gradient

I am unsure as to why I get this with some data and not with others. For example, when I run the Laceration subset the model runs fine. Model output:

Formula: Area ~ quadplat(Days, a, b, clx)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
a     1.2304     3.8509   0.320    0.753    
b     3.0869     0.5595   5.518 2.54e-05 ***
clx  62.7697    11.0592   5.676 1.80e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.86 on 19 degrees of freedom

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

I interpret this as the critical threshold where there is no statistical change in Y with an increase in X is 62.7697 days. Is that a correct interpretation?

Plot below:

To me, this plot looks good. However, when I run the same analysis with the abrasion subset I get the singular gradient error. Why might this be, is this because the data is not fitting well?

Please can someone with knowledge of nls help me by explaining exactly what this quadratic model is doing and why I might get an error. I don't want to 'black box' this analysis and I think I am missing key understanding. Additionally, if anyone out there is good at interpreting formulas, can you help me by writing up this code into a readable formula?

function(x, a, b, clx) {
  ifelse(x  < clx, a + b * x + (-0.5*b/clx) * x * x, 
         a + b * clx + (-0.5*b/clx) * clx * clx)}

Any information on this question is greatly appreciated or directions toward good resources on nls. I really need some help here and can attach my full dataset if needed.

Answer 1

We need better starting values. Fit a non-plateau model, model0, and use the parameters from that to fit all the data points giving model and then use a and b from that and a grid of values for clx (due to its problematic nature) giving model.Ab and model.La. (Note that it will not be able to produce fits from some of the grid's starting values resulting in error messages but nls2 will keep processing further starting values so those errors can be ignored.)

library(nls2)

# ensure data is sorted for plotting
o <- with(healing, order(Type, Days))
h <- healing[o, ]

# last argument specifies whether there is or is not a plateau
quadplat = function(x, a, b, clx, plat = TRUE) {
  if (plat) x <- pmin(x, clx)
  a + b * x + (-0.5*b/clx) * x * x
}

# fit no plateau model with all data
st <- c(a = 1, b = 1, clx = 1)
model0 <- nls(Area ~ quadplat(Days, a, b, clx, FALSE), h, start = st)

# fit all data model
model <- nls(Area ~ quadplat(Days, a, b, clx), h, start = coef(model0))
co <- coef(model)

We can now fit and plot the subset models using values computed above in the starting values.

if (exists("model.Ab")) rm(model.Ab)
model.Ab <- nls2(Area ~ quadplat(Days, a, b, clx), h, subset = h$Type == "Abrasion",
  start = data.frame(a = co[[1]], b = co[[2]], clx = 0:140))

if (exists("model.La")) rm(model.La)
model.La <- nls2(Area ~ quadplat(Days, a, b, clx), h, subset = h$Type == "Laceration",
  start = data.frame(a = co[[1]], b = co[[2]], clx = 0:140))
  
cols <- c(Abrasion = "red", Laceration = "blue")
plot(Area ~ Days, h, col = cols[Type], pch = 20, cex = 1.5)
lines(fitted(model.Ab) ~ Days, subset(h, Type == "Abrasion"), 
  col = cols["Abrasion"])
lines(fitted(model.La) ~ Days, subset(h, Type == "Laceration"), 
  col = cols["Laceration"])

(continued after graphics)

Alternate model

If it is ok to consider other models then this model has only two parameters, are easier to fit and despite having fewer parameters have lower residual sums of squares.

model.Ab2 <- nls(Area ~ a * (1 - exp(- b * Days)), h, 
   subset = Type == "Abrasion", start = c(a = 100, b = .1))
 
model.La2 <- nls(Area ~ a * (1 - exp(- b * Days)), h, 
   subset = Type == "Laceration", start = c(a = 100, b = .1))

# plot
cols <- c(Abrasion = "red", Laceration = "blue")
plot(Area ~ Days, h, col = cols[Type], pch = 20, cex = 1.5)
lines(fitted(model.Ab2) ~ Days, subset(h, Type == "Abrasion"), 
  col = cols["Abrasion"])
lines(fitted(model.La2) ~ Days, subset(h, Type == "Laceration"), 
  col = cols["Laceration"])

(continued after graphics)

One parameter model

If we fix a = 100 in the 2 parameter model of the last section we get a 1 parameter model which is not statistically distinguishable from the 2 parameter model. That is seen from the p value shown in the anovas which are greater than 0.05 indicating that we cannot reject the null hypothesis that the 1 and 2 parameter models describe the data equally well for each of the two subsets.

model.Ab3 <- nls(Area ~ 100 * (1 - exp(- b * Days)), h, 
   subset = Type == "Abrasion", start = c(b = .1))
 
model.La3 <- nls(Area ~ 100 * (1 - exp(- b * Days)), h, 
   subset = Type == "Laceration", start = c(b = .1))

anova(model.Ab3, model.Ab2)

anova(model.La3, model.La2)

Also note that the point at which it reaches y = 95, i.e. near plateau, is -log(1 - 95/100)/b (based on inverting the model equation). The numerator is approximately 3 so it reaches 95 at roughly 3/b.

Other

If m <- nls(...) then summary(m) will give standard errors of coefficients and other information.

Answer 2

Additionally, if anyone out there is good at interpreting formulas, can you help me by writing up this code into a readable formula?

function(x, a, b, clx) {
  ifelse(x  < clx, a + b * x + (-0.5*b/clx) * x * x, 
         a + b * clx + (-0.5*b/clx) * clx * clx)}

$$ f(x, a, b, x_{cl}) = \begin{cases} a + bx + (\frac{-0.5b}{x_{cl}}) \times x^2 , & \text{if}\ x < x_{cl} \\ a + bx_{cl} + (\frac{-0.5b}{x_{cl}}) \times {x_{cl}}^2 , & \text{otherwise} \end{cases} $$

which simplifies to:

$$ f(x, a, b, x_{cl}) = \begin{cases} a + bx \left( 1 - \frac{x}{2x_{cl}} \right) , & \text{if}\ x < x_{cl} \\ a + \frac{bx_{cl}}{2} , & \text{otherwise} \end{cases} $$

where I have substituted $x_{cl}$ for clx to make it more readable.