# Custom link function needed for generalised linear model in R

I would like to use a generalised linear model to analyse data on the relationship between the size of a host and probability of parasitism in a wild population to determine the minimum host size at which parasitism is likely. The independent variable is continuous (host size) and the dependent variable is binomial (unparasitised/parasitised). A plot of the parasitism data resembles a logistic regression, except the maximum probability is not 1.

In other words, the asymptotic maximum probability of parasitism M in different populations is 0 < M < 1. It seems that I need to use a custom link function in R that looks something like family=binomial(link = M*logit).

Can someone suggest appropriate R code to do this?

My initial research suggested that this data set can be modelled by the four-parameter logistic model, but that is not the case.

#### Edit:

This is the outcome of nonlinear regression and a vanilla logistic regression:

Raw data look like this:

HS_data <- data.frame(
Size=c(4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14.5),
Total=c(1, 18, 50, 97, 80, 78, 73, 72, 32, 47, 49, 63, 65, 49, 27, 10, 6, 2, 2, 2),
Unparasitised=c(1, 18, 50, 97, 80, 78, 70, 71, 24, 38, 32, 42, 37, 29, 15, 7, 4, 1, 1, 1),
Parasitised=c(0, 0, 0, 0, 0, 0, 3, 1, 8, 9, 17, 21, 28, 20, 12, 3, 2, 1, 1, 1))


This is the non-linear regression model:

# Analysis using non-linear regression
# The analysis ignores the binary nature of the data  >
SizePref <- nls(Parasitised/Total ~ a/(1+exp(b*(Size+c))),
data = HS_data,
start=list(a=.45, b=-1.3, c=-11))

• I’m not sure why vanilla logistic regression would not for you; it never returns 0 or 1, but some M arbitrarily close to these limits. (Ex 0.0001 or 0.99, etc) May 31, 2023 at 5:33
• Unless you mean that for a given species, the max probability of parasitism is 0.7 for instance and you need to estimate the probability given that it’s less than 0.7. But I struggle to imagine a situation where’d you’d know this. May 31, 2023 at 5:35
• This is not a species specific situation. It refers to what occurs in a field population. So, if hosts are collected from the field, the overall % parasitism is never 100%. May 31, 2023 at 5:39
• @MikeKeller: It may be relevant to know if the formula proportion ~ M/(1 + exp(b * (Host_size + c))) has some scientific relevant background, or if it is picked because it normally fits the data well.
– Ute
Jun 2, 2023 at 8:03
• @MikeKeller, would it make sense to assume that the host population consists of a mix of resistent and susceptible individuals? The sigmoidal you fitted by nls looks indeed convincing :-) this would be a good justification for the second comment you got on this question
– Ute
Jun 2, 2023 at 10:26

## Custom link function for binomial linear models with upper bound

The standard model for success counts that depend some covariates are generalized linear models (glm) with logit or probit link. These will always give predicted success probabilities between 0 and 1, that can get arbitrarily close to 1.

In some applications, however, it makes sense to assume that the true success probabilities are bounded by some number $$M$$, as is the case in the present question, where individuals from some (animal?) population were classified into parasitized and not parasitized:
If we assume that there are two subpopulations, namely immune and susceptible individuals, then the probability of being parasitized should not be larger than the proportion of the susceptible subpopulation. We will call that proportion $$M$$, and let $$p_0(\boldsymbol{x})$$ be the success probability corresponding to covariates value $$\boldsymbol{x}$$ in the susceptible subpopulation ($$\boldsymbol{x}$$ can be a number or a vector $$(x_1,\dots,x_k)$$ if we have $$k$$ covariates).

Under the assumption that all individuals in the population have the same chance of being sampled, the probability of success for a randomly picked individual from the total population is then $$p(\boldsymbol{x})=M p_0(\boldsymbol{x})$$

Now consider a binomial glm with link function $$g_S$$ that describes the susceptible subpopulation. The glm implies that there exist coefficients $$\beta_0,\dots,\beta_k$$ such that $$p_0(\boldsymbol{x}) = g_S^{-1}(\eta),\quad \eta=\beta_0+\beta_1x_1 +\dots \beta_kx_k.$$ The corresponding link function $$g$$ for the total population must fulfil $$p(\boldsymbol{x}) = Mp_0(\boldsymbol{x}) = g^{-1}(\eta),\quad\text{thus}\quad g^{-1}(\eta)=Mg_S^{-1}(\eta).$$

This is why OP has asked for a custom link function of the form link = M*logit.

Note that the parameter $$M$$ is part of the model here. As a consequence, $$p$$-values or confidence interval calculated for the fitted model with standard methods are statistically valid. This is not exactly the case if the data are used to estimate $$M$$.

### Turning any link function into a custom link function with upper bound

Indeed, it is possible to define custom link functions in R to be used with glm (this has been used by @BenBolker to answer a question about Gamma glm on stackowerflow).
The required ingredients are the link function, its inverse, the derivative of the inverse and a function that defines the support of the inverse link function.

Writing $$\mu$$ for $$p(\boldsymbol{x})$$, we can obtain these functions from the original subpopulation link function $$g_S$$:

• the link function, $$g(\mu) = g_S(\mu / M)$$,
• the inverse link function, $$g^{-1}(\eta) = Mg_S^{-1}(\eta)$$,
• the derivative of the inverse link function $$(g^{-1})^\prime(\eta)= M(g_S^{-1})^\prime(\eta)$$,
• the support of $$g^{-1}$$ is the same as the support of $$g^{-1}_S$$.

The following R function Mscale_linkfunction takes the original link function $$g_S$$ (argument oli)and returns the link function $$g$$:

Mscale_linkfunction <- function(oli = "logit", M = 1){
if (class(oli) == "character")  oli <- make.link(oli)
structure(list(
linkfun = function(mu) oli$linkfun(mu/M), linkinv = function(eta) oli$linkinv(eta)*M,
mu.eta = function(eta) oli$mu.eta(eta)*M, valideta = oli$valideta,
name = paste0("M-scaled ",oli$name,", M = ",signif(M,3))), class = "link-glm" ) }  To use these constraint link functions in glm, you need provide a start value that ensures that the argument $$\mu/M$$ does not exceed 1. The wrapper code that finds $$M$$ from the data (in the last section) uses vanilla glm to find start values. #### Example Here is the analysis of the parasitism data with a probit link model with $$M=0.5$$: customlink <- Mscale_linkfunction("probit", 0.5) custom <- glm(cbind(Parasitised, Unparasitised) ~ Size, data = HS_data, family = binomial(link=customlink), start=c(0,0)) plot(Parasitised/Total ~ Size, data=HS_data, ylim = c(0,.8), xlab = "Size[mm]") lines(fitted(custom) ~ HS_data$Size, col = "blue")


### Using the data to determine the model parameter $$M$$

One option to determine the maximum probability from the data is to choose the best fitting model, according to AIC. The following code defines a function MaxProb_BinCLM that does this in a hands-on manner, and returns both the fitted glm and the value of $$M$$.

MaxProb_BinGLM <- function(formula, data, linkfoo, start = NULL){
glmM <- glm(formula, data = data,
family = binomial(linkfoo), start = start)
start <- glmM$coefficients M_AIC <- function(M) { Mlink <- Mscale_linkfunction(linkfoo, M) glmM <<- glm(formula, data = data, family = binomial(Mlink), start = start) start <<- glmM$$coefficients glmM$$aic } M <- optimize(M_AIC, c(0,1))$minimum
list(fitglm = glmM, M = M)
}


With a logit link, the function finds $$M=0.4036$$ for the present data.

parasites <- MaxProb_BinGLM(cbind(Parasitised, Unparasitised) ~ Size, data = HS_data, linkfoo = "logit")
parasites$M # gives 0.4036025 plot(Parasitised/Total ~ Size, data=HS_data, ylim = c(0,.8), xlab = "Size[mm]") lines(fitted(parasites$fitglm) ~ HS_data$Size, col = "blue")  The AIC for the fitted model should be corrected by adding 2, since now $$M$$ is a fitted model parameter. • +1. Just wanted to raise a (philosophical?) question on this. Do we ever want or need to mess with the functional form of the link function when it comes to addressing things like this? Shouldn't we always try to better fit on the linear-predictor side of the regression? In this approach, are we not choosing to prioritise the desired upper bound property over the actual fit/performance of the model? (Because we are still not addressing the non-linearity of the behaviour). Jun 21, 2023 at 0:52 • Yes, this approach needs to be complemented with a procedure to find the upper bound, if it is not a priori given. But apart from that, isnt glm in most casrs somewhat arbitrary, just like all modelling? – Ute Jun 21, 2023 at 9:06 • Does it make sense to impose a strict upper bound on the$p$but then to have observations above that bound? Isn't that contradictory? Jun 23, 2023 at 1:54 • @statsplease: It is not contradictory. The upper bound is for the whole (theoretical) population. The observations are only a random sample. Imagine you sit in a room with 10 people, and one has won in a lottery with a chance of 1 in 10000000. From the sample, you would estimate the probability to be 1 in 10. You can estimate the probability bound from the data, either by full maximum likelihood as suggested by Ben Bolker, or by looking for the best fit as in my approach. – Ute Jun 23, 2023 at 6:35 Here's another way to do it, as a generic maximum likelihood estimation problem: library(bbmle) mle2(Parasitised ~ dbinom(prob = plogis(maxpar)*plogis(a + b*Size), size = Total), data = HS_data, start = list(maxpar = 0, a = -5, b = 1))  • Advantages: flexible, doesn't require a nested estimation procedure (i.e. GLM within a loop to figure out the best maxpar estimate) • Disadvantages: needs starting values for all parameters. General nonlinear estimation is less efficient and less robust than iteratively reweighted least squares. Here's a slightly different formulation, using a linear sub-model for the logit-probability, that would be useful if you had a more complicated models with several predictors determining the probability (relative to maxpar). mle2(Parasitised ~ dbinom(prob = plogis(maxpar)*plogis(eta), size = Total), parameters = list(eta ~ Size), data = HS_data, start = list(maxpar = 0, eta = 0))  • Smart - MLE out of the box for complex models. I will try your handy package. – Ute Jun 23, 2023 at 7:20 • Ben Bolker's fully MLE approach gives virtually the same results as searching for minimal AIC, in the present simple case (coefficients and estimate for$M$deviate by less than 0.1 percent) - as it should be :-) – Ute Jun 23, 2023 at 17:44 I think there's a bit of a misunderstanding here on what a binomial GLM is. We can tweak your data, which might help your understanding. library(magrittr) x = data.frame( Size = c(4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14.5), Total = c(1, 18, 50, 97, 80, 78, 73, 72, 32, 47, 49, 63, 65, 49, 27, 10, 6, 2, 2, 2), Unparasitised = c(1, 18, 50, 97, 80, 78, 70, 71, 24, 38, 32, 42, 37, 29, 15, 7, 4, 1, 1, 1), Parasitised = c(0, 0, 0, 0, 0, 0, 3, 1, 8, 9, 17, 21, 28, 20, 12, 3, 2, 1, 1, 1)) x1 = x %>% dplyr::select(Size, Parasitised) %>% tidyr::uncount(weights = Parasitised) %>% dplyr::mutate(Y = 1) x2 = x %>% dplyr::select(Size, Unparasitised) %>% tidyr::uncount(weights = Unparasitised) %>% dplyr::mutate(Y = 0) y = rbind(x1, x2)  We now have the individual binary responses per row in the data. We can now pass to the model for fitting. I have used the convention that $$Y=1$$ indicates Parasitised. g = glm(Y ~ Size, family = binomial("logit"), data = y) summary(g) > Call: > glm(formula = Y ~ Size, family = binomial("logit"), data = rbind(x1, x2)) > > Deviance Residuals: > Min 1Q Median 3Q Max > -2.2432 -0.5598 -0.2687 -0.1842 2.4387 > > Coefficients: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -8.65640 0.70431 -12.29 <2e-16 *** > Size 0.76470 0.07106 10.76 <2e-16 *** > --- > Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > (Dispersion parameter for binomial family taken to be 1) > > Null deviance: 707.98 on 822 degrees of freedom > Residual deviance: 530.22 on 821 degrees of freedom > AIC: 534.22 > > Number of Fisher Scoring iterations: 6  We can see the fit results above. The fit indicates that Size increases the chances of being Parasitised. This is, of course, assuming this is the functional form you choose for the model. You have the freedom to impose whatever structure you want on the link-side of the regression (e.g. GAMs with splines). The true relationship between Size and being Parasitised could well be non-linear. This is how your fitted model looks. ggplot2::ggplot() + ggplot2::geom_point(data = y %>% dplyr::group_by(Size, Y) %>% dplyr::summarise(N = dplyr::n()), ggplot2::aes(Size, Y, size = N)) + ggplot2::geom_line(data = data.frame(Size = seq(0, 20, 0.01)) %>% dplyr::mutate(y = predict(g, newdata = ., type = "response")), ggplot2::aes(Size, y)) + ggplot2::labs(title = "Fitted logistic curve", y = latex2exp::TeX("$$E[Y|X]$$")) + ggplot2::scale_y_continuous(labels = scales::percent)  As you can see, the model does indeed predict probabilities across the full support of a Bernoulli random variable ($$0\leq p\leq 1$$). The point of fitting a binomial GLM is precisely to create a mapping between the probability parameter of the Bernoulli random variable to the regressor(s) you are interested in (here it is Size). For example, what your model here tells you is that, for a host of Size 15 (whatever the units are), your random variable is $$(\text{Parasitised}\,|\, \text{Size}=15)\sim \text{Bernoulli}\big(p = 1/(1 + \exp(-(-8.65640+15\times 0.76470))\big)$$ That is, $$(\text{Parasitised}\,|\, \text{Size}=15)\sim \text{Bernoulli}(p=0.943433)$$ So your model is telling you it believes that it is quite likely to be Parasitised. If this behaviour doesn't seem intuitive (based on theory outside of statistics) then there is probably something going on in the model. Have you ignored non-linearities where the model implicitly assumes larger Sizes lead to higher chances of Parasitised? Finally, it's just worth pointing out that the model does in general agree with your empirical observations of a 50% probability of being Parasitised across your empirical observations. If you look at the plot, the predicted probablity of being Parasitised is generally below 50% for most of the observations. Anything below a Size of 11.32 has a predicted probability less than or equal to 50%. I think what your question boils down to is the functional relationship between the regressor and the outcome is non-linear, but you haven't allow the GLM to account for this nature. The alternative approach using non-linear least-squares obviously is being allowed to account for this. If I just throw your data at a simple GAM (without any deep thought) I can already see the nature of the relationship is (surprise!) non-linear. g2 = mgcv::gam(Y ~ s(Size), family = binomial("logit"), data = y)  The blue line is the fit if we don't address the non-linear behaviour (irrespective of whether we manipulate the link function). If we replot our logistic fit with this new functional relationship for the link we get: I imagine you'll be happier now with this view of the fit. If you assess the fit of this new model you'll see if outperforms the old one in certain ways. (Included in the plot is the extrapolation of the fit beyond the spline. Here it's preserving the gradient beyond the support of the spline but you could equally choose any extrapolation). In general, if you've constructed a model that is sound theoretically: • Binary data with regressors - Binomial GLM (tick!) but something about your model behaviour doesn't seem intuitive: • then it's likely you're restricting it in some unfair way and it's being forced to work within those constraints. • This does not answer the original question (custom link function), anyway +1 for pointing out nonparametric alternatives :-) – Ute Jun 23, 2023 at 8:25 • Here is code to produce a simple plot comparing observed probabilities with the fitted curve: pagam <- mgcv::gam(cbind(Parasitised, Unparasitised) ~ s(Size), data = HS_data,family = binomial("logit")); plot(Parasitised/Total ~ Size, data=HS_data, ylim = c(0,.8), xlab = "Size[mm]"); lines(fitted(pagam) ~ HS_data$Size, col = "blue")  It looks good, and also here you see that the observed probabilities exceed the fitted values from the curve in the end.
– Ute
Jun 23, 2023 at 17:53
• Btw, as far as I interpret the question, the reason for asking was exactly the observation that a simple out-of-the-box link function does not do the job to explain the limited increase of parasitism probability with size sufficiently, and that this could be corrected in some way by multiplying the link function with the asymptotic upper bound (which usually is 1). And OP was right. They did not "ignore non-linearities" - a binomial glm with logit link will always give an S-shaped curve, it is just not visible in the plot ending in Size=14.
– Ute
Jun 23, 2023 at 18:10
• I agree the GAM does fit data. Regarding the last point. They did ignore the non-linear nature. I'm not saying a logit link function isn't non-linear (of course it is). The point is that on the linear-predictor side there is clearly non-linear behaviour at play which is being ignored. How is that debateable? I've updated the plot above to illustrate exactly that. The only thing that manipulating the link function achieves here is masking the poor fit. Jun 24, 2023 at 4:29
• I know my answer doesn't directly address the question. You'll find that a huge number of answers on this site are useful precisely because they point out why the approach in the question is misguided. Jun 24, 2023 at 4:31