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)
```
If you mean the *link function*, which translates from the data scale
to the linear predictor scale, you can get it via
b_model$family$linkfun
If you want the cycles at which the signal is equal to 0.5,
eta <- b_model$family$linkfun(0.5)
with(as.list(coef(b_model)), (eta-`(Intercept)`)/Cycles)
You might find the drc package useful as well.
cheers
Ben Bolker
On 4/15/25 08:57, Luigi Marongiu 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 at 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.
--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
* E-mail is sent at my convenience; I don't expect replies outside of
working hours.
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.luigi at 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 at 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.-------------- next part -------------- A non-text attachment was scrubbed... Name: signature.asc Type: application/pgp-signature Size: 603 bytes Desc: OpenPGP digital signature URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20250415/7ac3a87e/attachment.sig>