simple curve fitting in R-CodePudding

I am trying to find a fit for my data. But so far had no luck. Tried the logarithmic, different ones from the drc package .. but I am sure there must be a better one I just don't know the type. On a different note - I would be grateful for advice on how to go about curve hunting in general.

library(drc)    
df<-structure(list(x = c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
        20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 
        36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 
        52, 53, 54, 55), y = c(0.1066, -0.6204, -0.2028, 0.2621, 0.4083, 
        0.4497, 0.6343, 0.7762, 0.8809, 1.0029, 0.8089, 0.7845, 0.8009, 
        0.9319, 0.9414, 0.9505, 0.9323, 1.0321, 0.9381, 0.8975, 1.0929, 
        1.0236, 0.9589, 1.0644, 1.0411, 1.0763, 0.9679, 1.003, 1.142, 
        1.1049, 1.2868, 1.1569, 1.1952, 1.0802, 1.2125, 1.3765, 1.263, 
        1.2507, 1.2125, 1.2207, 1.2836, 1.3352, 1.1311, 1.2321, 1.4277, 
        1.1645), w = c(898, 20566, 3011, 1364, 1520, 2376, 1923, 1934, 
        1366, 1010, 380, 421, 283, 262, 227, 173, 118, 113, 95, 69, 123, 
        70, 80, 82, 68, 83, 76, 94, 101, 97, 115, 79, 98, 84, 92, 121, 
        97, 102, 93, 92, 101, 74, 124, 64, 52, 63)), row.names = c(NA, 
        -46L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(`47` = 47L), class = "omit"))
    
    fit <- drm(data = df,y ~ x,fct=LL.4(), weights = w)
    plot(fit)

CodePudding user response：

The basic idea is to understand how the selected function performs. Take a function you know (e.g. the logistic) and modify it. Or (even better) go to the literature and see which functions people use in your specific domain. Then create a user-defined model, play with it to understand the parameters, define good start values and then fit it.

Her a quick & dirty example of a user-defined function (with package growthrates). It can surely be made similarly with drc.

library("growthrates")

grow_userdefined <- function (time, parms) {
  with(as.list(parms), {
    y <- (K * y0)/(y0   (K - y0) * exp(-mumax * time))   shift 
    return(as.matrix(data.frame(time = time, y = y)))
  })
}

fit <- fit_growthmodel(FUN=grow_userdefined, 
  p = c(y0 = -1, K = 1, mumax = 0.1, shift = 1), 
  time = df$x, y = df$y)

plot(fit)
summary(fit)

It can of course be made better. As we have no exponential start at the onset, one can for example start with a simple saturation function instead of a logistic, e.g. something Monod-like. As said, the preferred way is to use a function related to the application domain.

modified logistic

CodePudding user response：

1) If we ignore the weights then y = a b * x c/x^2 seems to fit and is linear in the coefficients so is easy to fit.

fm <- lm(y ~ x   I(1 / x^2), df)
plot(y ~ x, df)
lines(fitted(fm) ~ x, df)

2) Another approach is to use two lines. We used the weights.

# starting values use lines fitted to 1st ten and last 10 points
fm_1 <- lm(y ~ x, df, subset = 1:10)
fm_2 <- lm(y ~ x, df, subset = seq(to = nrow(df), length = 10))
st <- list(a = coef(fm_1)[[1]], b = coef(fm_1)[[2]],
           c = coef(fm_2)[[1]], d = coef(fm_2)[[2]])

fm2 <- nls(y ~ pmin(a   b * x, c   d * x), df, start = st, weights = w)
plot(y ~ x, df)
lines(fitted(fm2) ~ x, df)