6
$\begingroup$

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:

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ (+1) (In geostatistics this is known as a "spherical variogram model.") Its second derivative is discontinuous and its parameter estimates, as given, will be very highly correlated: these mathematical properties make it less than desirable as a model. Indeed, it's not a good fit to your example data and using nls (least squares) to fit proportions is not statistically suitable. One solution, then, would be to reconsider both your model and the fitting procedure, rather than working to make this one run without errors. Perhaps the computational errors are telling you something important! $\endgroup$ – whuber Apr 8 at 14:19
5
$\begingroup$

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) screenshot

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) screenshot

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Wow, this is so helpful and you have helped me move forward with something I was really struggling with for a while. Thank you so much. Is is possible to generate standard errors and p values for the first models you suggested. I will also consider the second. Do you know how I might write out the quadratic equasion used in the first as a greek formula? Huge thank you again for your help here. $\endgroup$ – Freya Womersley Apr 11 at 12:37
  • $\begingroup$ Have fixed cut and paste error resulting in difference between code and description. Have also added One Parameter and Other sections which include some information that address some of the questions in your comment. $\endgroup$ – G. Grothendieck Apr 11 at 13:26
  • $\begingroup$ For reference, would these models be possible if my data was negatively decaying? $\endgroup$ – Freya Womersley Apr 12 at 17:51
  • $\begingroup$ It would depend on the actual data but you could try a * exp(-b * Days) or fix a to 100 if appropriate. For example, try curve(100 * exp(- .2 * x), 0, 50) $\endgroup$ – G. Grothendieck Apr 12 at 19:23
  • $\begingroup$ One last question. Can you help with the equation of this model? In the alternative one and parameter models you suggested, am I correct in the following equation: y = ae^-bx. In this case, what does a equate to as it is not Y when x = 0. $\endgroup$ – Freya Womersley Apr 22 at 14:01
4
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.