1

I am trying to fit a curve following a non-linear regression.

Here is my dataset:

stem_diameter <- c(15, 15, 16, 17, 19, 23, 23, 24, 24, 25, 25, 26, 27, 28, 29, 30, 30, 32, 32, 33, 34, 34, 35, 36, 36, 37, 38, 40, 41, 41, 42, 42, 46, 48, 48, 49, 51, 54, 55, 60)

total_biomass <- c(0.25, 0.65, 0.40, 0.40, 0.65, 0.60, 0.20, 0.60, 0.50, 0.35, 0.70, 0.65, 0.95, 0.60, 0.80, 0.90, 0.70, 1.15, 1.00, 1.70, 1.95, 1.15, 1.70, 1.25, 1.95, 1.20,
2.70, 2.00, 3.70, 2.35, 1.50, 2.50, 5.05, 4.10, 2.85, 2.15, 3.50, 4.80, 5.30, 2.95)

stem_data.df <- data.frame(stem_diameter, total_biomass)

I have tried fitting a linear model as below:

model1 <- lm(total_biomass ~ stem_diameter, data = stem_data.df)
summary(model1)
plot(total_biomass ~ stem_diameter)
lines(stem_diameter, predict(model1))

however I think a non-linear model is more appropriate for the data.

I then went on to do this:

m <- nls(total_biomass ~ a(stem_diameter^b), stem_data.df, start = list(a = -1.8, b = 0.1)) # power formula: y = a*x^b
summary(m)
plot(total_biomass ~ stem_diameter)
lines(stem_diameter, fitted(m), col = "green")

with a and b the coefficients from the above linear model, as I read elsewhere. The curve fitted on the plot however does not look how it should, leading me to think that there is a problem with the formula itself. The data itself appears to be fine, as I have been able to plot a curve in an excel scatter plot which had an R2 value of 0.79 and formula of y = 0.0009x^2.0598. Changing the values for a and b has no real effect on the curve either therefore I am not sure what the problem could be.

Can anyone explain what the problem may be?

Thanks in advance

G Jones
  • 43
  • 4

1 Answers1

2

I suggest you study this post: https://stats.stackexchange.com/a/255265/11849

Anyway,

fit1 <- lm(log(total_biomass) ~ log(stem_diameter), data = stem_data.df)

#you were missing a `*` and had bad starting values
m <- nls(total_biomass ~ a*(stem_diameter^b), stem_data.df, 
         start = list(a = exp(coef(fit1)[1]), b = coef(fit1)[2])) # power formula: y = a*x^b
summary(m)
plot(total_biomass ~ stem_diameter)
curve(predict(m, newdata = data.frame(stem_diameter = x)), add = TRUE)
Roland
  • 127,288
  • 10
  • 191
  • 288