Hello Luigi, Great follow-up — Looks like you’re on the right track using nls() for nonlinear regression. You're fitting a logistic-like sigmoidal model (as in Rutledge’s paper), I think both errors you’re encountering are signs of bad starting values or a maybe model formulation issue.
Breakdown of Your Model Attempt Your model formula: > Fluo ~ (MaxFluo / (1 + exp(-(Cycle - (HalfFluoCycle / Slope)))) + 1) This expression has perhaps two issues: 1. Parameter placement: In the paper, the term inside exp() should look like -(Cycle - HalfFluoCycle)/Slope — not Cycle - (HalfFluoCycle / Slope). 2. +1 additive constant: You called it backgroundSignal, but it’s hard-coded as +1. That constant should be a parameter (e.g., Bkg), otherwise the model fit has no way to adjust it. Alternative Model Form: Here’s a likely corrected version of the model based on Rutledge’s sigmoidal fit: > mod1 <- nls(Fluo ~ MaxFluo / (1 + exp(-(Cycle - HalfFluoCycle)/Slope)) > + > Bkg, > data = df, > start = list(MaxFluo = max(df$Fluo), > HalfFluoCycle = 20, > Slope = 2, > Bkg = min(df$Fluo))) Notes on nls() Convergence nls() is very sensitive to starting values. Always inspect your plot and roughly estimate parameters: MaxFluo: Use max(Fluo) Bkg: Use min(Fluo) HalfFluoCycle: Approximate where the curve inflects Slope: Try values like 1, 2, 5 depending on steepness You can also add trace=TRUE to watch convergence progress: > mod1 <- nls(..., trace=TRUE) About SSmicmen >From what I've gathered - this is meant for fitting the Michaelis-Menten model >(enzyme kinetics), not sigmoidal growth. That's why mod2 fails — wrong >self-starting model function. Visualizing the Fit After fitting: > plot(Fluo ~ Cycle, data = df) > lines(df$Cycle, predict(mod1), col = "red", lwd = 2) You might also want to compute the inverse (e.g., solving for Cycle given a fluorescence value) — that would involve solving the nonlinear equation manually or using uniroot(). All the best, gregg On Wednesday, April 23rd, 2025 at 8:23 AM, Luigi Marongiu <marongiu.lu...@gmail.com> wrote: > > > Further on this, I have found a formula from a paper from Rutledge > (DOI: 10.1093/nar/gnh177), which I rendered as > `MaxFluo / (1+ exp(-(Cycle-(HalfFluoCycle/Slope)))) + backgroundSignal` > I then see that one can use `nls` to fit non-linear regression, so I tried: > `df <- data.frame(Cycle=1:35, Fluo=c(15908.54, 16246.92, 16722.43, 17104.29, > 16794.93, 17031.44, 17299.05, 17185.49, 17362.28, 17127.43, 17368.96, > 17513.19, 17593.95, 17626.37, 18308.29, 18768.12, 19955.26, 22493.85, > 27088.12, 36539.44, 53694.18, 84663.18, 138835.64, 223331.89, 336409.69, > 457729.88, 561233.12, 637536.31, 688985.88, 723158.56, 746575.62, 766274.75, > 776087.75, 785144.81, 791573.81) ) plot(Fluo~Cycle, df) mod1 = > nls(Fluo~(MaxFluo / (1+ exp(-(Cycle-(HalfFluoCycle/Slope)))) + 1), data=df, > start=list(MaxFluo=max(df$Fluo), HalfFluoCycle=15, Slope=0.1)) mod2 = > nls(Fluo~SSmicmen(Fluo, Cycle), data=df)` > but I got errors in both cases: > ``` > > > mod1 > > Error in nlsModel(formula, mf, start, wts, scaleOffset = scOff, > nDcentral = nDcntr) : > singular gradient matrix at initial parameter estimates > > > mod2 > > Error in qr.qty(QR.rhs, .swts * ddot(attr(rhs, "gradient"), lin)) : > NA/NaN/Inf in foreign function call (arg 5) > ``` > How can I properly set this regression model? > Thank you > > > On Wed, Apr 16, 2025 at 7:08 AM Luigi Marongiu marongiu.lu...@gmail.com wrote: > > > Thank you. This topic is more complicated than anticipated. Best regards, > > Luigi > > > > On Tue, Apr 15, 2025 at 11:09 PM Andrew Robinson a...@unimelb.edu.au wrote: > > > > > A statistical (off-topic!) point to consider: when the GLM was fitted, > > > you conditioned on x and let y be the random variable. Therefore the > > > model supports predictions of y conditional on x. You’re now seeking to > > > make predictions of x conditional on y. That’s not the same thing, even > > > in OLS. > > > > > > It might not matter for your application but it’s probably worth thinking > > > about. Simulation experiments might shed some light on that. > > > > > > Cheers, Andrew > > > > > > -- > > > Andrew Robinson > > > Chief Executive Officer, CEBRA and Professor of Biosecurity, > > > School/s of BioSciences and Mathematics & Statistics > > > University of Melbourne, VIC 3010 Australia > > > Tel: (+61) 0403 138 955 > > > Email: a...@unimelb.edu.au > > > Website: https://researchers.ms.unimelb.edu.au/~apro@unimelb/ > > > > > > I acknowledge the Traditional Owners of the land I inhabit, and pay my > > > respects to their Elders. > > > > > > On 16 Apr 2025 at 1:01 AM +1000, Gregg Powell via R-help > > > r-help@r-project.org, wrote: > > > > > > Take a look at this Luigi... > > > > > > # The model is: logit(p) = β₀ + β₁*Cycles > > > # Where p is the probability (or normalized value in your case) > > > > > > # The inverse function would be: Cycles = (logit⁻¹(p) - β₀)/β₁ > > > # Where logit⁻¹ is the inverse logit function (also called the expit > > > >function) > > > > > > # Extract coefficients from your model > > > intercept <- coef(b_model)[1] > > > slope <- coef(b_model)[2] > > > > > > # Define the inverse function > > > inverse_predict <- function(p) { > > > # p is the probability or normalized value you want to find the >cycles > > > for > > > # Inverse logit: log(p/(1-p)) which is the logit function > > > logit_p <- log(p/(1-p)) > > > > > > # Solve for Cycles: (logit(p) - intercept)/slope > > > cycles <- (logit_p - intercept)/slope > > > > > > return(cycles) > > > } > > > > > > # Example: What cycle would give a normalized value of 0.5? > > > inverse_predict(0.5) > > > > > > This function takes a probability (normalized value between 0 and 1) and > > > returns the cycle value that would produce this probability according to > > > your model. > > > Also: > > > This works because GLM with binomial family uses the logit link function > > > by default > > > The inverse function will return values outside your original data range > > > if needed > > > This won't work for p=0 or p=1 exactly (you'd get -Inf or Inf), so you > > > might want to add checks > > > > > > All the best, > > > Gregg > > > > > > On Tuesday, April 15th, 2025 at 5:57 AM, Luigi Marongiu > > > marongiu.lu...@gmail.com wrote: > > > > > > I have fitted a glm model to some data; how can I find the inverse > > > function of this model? Since I don't know the y=f(x) implemented by > > > glm (this works under the hood), I can't define a f⁻¹(y). > > > Is there an R function that can find the inverse of a glm model? > > > Thank you. > > > > > > The working example is: > > > `V = c(120.64, 66.14, 34.87, 27.11, 8.87, -5.8, 4.52, -7.16, -17.39, > > > -14.29, -20.26, -14.99, -21.05, -20.64, -8.03, -21.56, -1.28, 15.01, > > > 75.26, 191.76, 455.09, 985.96, 1825.59, 2908.08, 3993.18, 5059.94, > > > 6071.93, 6986.32, 7796.01, 8502.25, 9111.46, 9638.01, 10077.19, 10452.02, > > > 10751.81, 11017.49, 11240.37, 11427.47, 11570.07, 11684.96, 11781.77, > > > 11863.35, 11927.44, 11980.81, 12021.88, 12058.35, 12100.63, 12133.57, > > > 12148.89, 12137.09) df = data.frame(Cycles = 1:35, Values = V[1:35]) M = > > > max(df$Values) df$Norm = df$Values/M df$Norm[df$Norm<0] = 0 b_model = > > > glm(Norm ~ Cycles, data=df, family=binomial) plot(Norm ~ Cycles, df, > > > main="Normalized view", xlab=expression(bold("Time")), > > > ylab=expression(bold("Signal (normalized)")), type="p", col="cyan") > > > lines(b_model$fitted.values ~ df$Cycles, col="blue", lwd=2)` > > > > > > ______________________________________________ > > > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > > > https://stat.ethz.ch/mailman/listinfo/r-help > > > PLEASE do read the posting guide > > > https://www.R-project.org/posting-guide.html > > > and provide commented, minimal, self-contained, reproducible code. > > > > -- > > Best regards, > > Luigi > > > > > -- > Best regards, > Luigi
signature.asc
Description: OpenPGP digital signature
______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide https://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
