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
______________________________________________
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.
