# Continuous predictor variable has non-linear association with log odds of response

I'd like to know if a continuous predictor variable (sal) predicts a binary response (empty) in a logistic mixed effect model, but discovered there is a non-linear relationship between the predictor variable and the log odds of the response. My understanding is if this linearity assumption of logistic regression isn't met it means I can't use this predictor in the model, but is there an alternative?

# Station = location of data collection
# CYR = Calendar year
# Sal = Salinity
# Empty = Fish caught had empty stomach or not (empty=1)

# Salinity variable centered
logmod <- glmer(empty ~  center_sal + (1|Station) + (1|CYR), family=binomial(link = "logit"), data = c_neb,
control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))

prob <- predict(logmod, type = "response")
logit <- log(prob/(1-prob))

ggplot(c_neb, aes(center_sal, logit)) +
geom_point(size=0.5, alpha=0.5) +
geom_smooth(method = "loess") +
theme_bw()


The binary response is independent of time (so survival analysis wouldn't work). Would a GAMM with a binary response work ? Something like this?

gamm <- gamm(empty ~ s(sal) + s(Station, bs = "re") + s(CYR, bs = "re"), family="binomial", data=c_neb)


Data (random sample of 10 rows):

df <- c_neb[sample(nrow(c_neb), 10), ]

dput(df)

structure(list(CYR_Keyfield = c("C-2019-5-15-241", "C-2013-8-16-70",
"C-2015-6-20-122", "C-2013-7-14-70", "C-2016-11-11-280", "C-2009-7-11-67",
"C-2009-7-11-67", "C-2009-7-11-67", "C-2018-9-17-123", "C-2017-8-13-168"
), CYR = structure(c(10L, 5L, 6L, 5L, 7L, 1L, 1L, 1L, 9L, 8L), levels = c("2009",
"2010", "2011", "2012", "2013", "2015", "2016", "2017", "2018",
"2019", "2021"), class = "factor"), Month = c("5", "8", "6",
"7", "11", "7", "7", "7", "9", "8"), Day = c("15", "16", "20",
"14", "11", "11", "11", "11", "17", "13"), Station = structure(c(49L,
91L, 10L, 91L, 64L, 89L, 89L, 89L, 11L, 28L), levels = c("101",
"105", "106", "107", "111", "112", "117", "118", "119", "122",
"123", "124", "130", "133", "134", "135", "137", "143", "144",
"145", "146", "147", "156", "157", "158", "159", "167", "168",
"169", "171", "172", "173", "174", "175", "176", "20", "21",
"22", "224", "225", "226", "227", "229", "23", "237", "239",
"24", "240", "241", "253", "254", "255", "256", "257", "265",
"266", "267", "268", "269", "270", "271", "278", "279", "280",
"281", "282", "283", "284", "290", "291", "292", "294", "301",
"302", "303", "304", "312", "313", "315", "323", "40", "54",
"609", "618", "619", "621", "622", "65", "67", "68", "70", "71",
"73"), class = "factor"), Zone = c("Rankin", "West", "West",
"West", "Rankin", "West", "West", "West", "West", "West"), ID = c("20195241_3",
"2013870_24", "20156122_8", "2013770_15", "201611280_70", "2009767_42",
"2009767_40", "2009767_39", "20189123_18", "20178168_134"), Species = c("Cynoscion nebulosus",
"Cynoscion nebulosus", "Cynoscion nebulosus", "Cynoscion nebulosus",
"Cynoscion nebulosus", "Cynoscion nebulosus", "Cynoscion nebulosus",
"Cynoscion nebulosus", "Cynoscion nebulosus", "Cynoscion nebulosus"
), SL = c(36.6, 30.25, 41.36, 43.15, 68.95, 36.4, 47.88, 59.63,
36.06, 17.94), empty = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L), levels = c("0", "1"), class = "factor"), Latitude = c(25.075,
25.132, 25.025, 25.132, 25.12, 25.134, 25.134, 25.134, 25.012,
25.094), Longitude = c(-80.816, -80.941, -80.938, -80.941, -80.785,
-80.967, -80.967, -80.967, -80.926, -80.905), sal = c(38.43,
36.1, 46.7, 33.2, 32.06, 37.7, 37.7, 37.7, 36.59, 43), temp = c(30.41,
29.9, 31.3, 29.4, 24.3, 32.2999999999998, 32.2999999999998, 32.2999999999998,
32.5739999999998, 33.7999999999998), group = structure(c(1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L), levels = c("euhaline", "hyperhaline",
"polyhaline"), class = "factor")), row.names = c(NA, -10L), class = c("tbl_df",
"tbl", "data.frame"))

• Isn't this scatterplot backwards? What would modeling an explanatory variable like center_sal as a function of predicted response tell you?
– whuber
Dec 16, 2022 at 18:22
• Ah, yes, sorry! I'll fix it.
– Nate
Dec 16, 2022 at 18:23
• Your plot show some clustered structures. When modelling the processes going on here, then I would be more interested in capturing the essence of those clusters rather than trying to create a straightforward function/model that tries to blend all of those clusters together. Dec 19, 2022 at 13:45
• @SextusEmpiricus how would you begin to address clustered structures?
– Nate
Dec 20, 2022 at 13:57
• I would start making more plots including the other data as well. Dec 20, 2022 at 14:34

Using GAMMs might be a suitable choice. With GAMMs you can model the non-linear relationship using smooths (which are penalized splines) while accounting for the hierarchical structure of the data. One potential advantage of GAMMs over some alternative methods is that they estimate the amount of complexity that is needed from the data.

You can implement them using the bam() function from R's mgcv package. Alternatively, you can use the gamm() function, but while it has a few additional modelling capabilities (which I think you do not need) it is a bit slower.

That said, I suspect there are some dependencies in the data structure that you have not yet incorporated in the suggested GAMM code. Without more knowledge about the data and question(s) you want to answer, it is hard to say which random effects structure would be the most appropriate here. I suspect that with using only random intercepts for Station and CYR you have not yet captured how the (non-linear) effect of sal on the outcome might differ for each Station and each year.

In mgcv you can incorporate a non-linear random effect of sal by station like this: s(center_sal, Station, bs = "fs", m = 1)

Note that as such random smooths are unconstrained, meaning that they also model average differences of the outcome between levels of the random factor and therefore make random intercepts unnecessary.

Fully capturing the hierarchical structure of the data in the model is advisable especially if you want to obtain appropriate inferences (in your case about the presence of an effect of sal on the outcome). Additionally, using a complete random effects structure can improve the models (out of sample) predictions.

Putting these points together, the code for a potential model might look like this:

library(mgcv)

gamm <- bam(empty ~ s(center_sal) +                             # smooth for non-linear effect of center_sal
s(center_sal, Station, bs = "fs", m = 1) +  # random smooth for Station
s(center_sal, CYR, bs = "fs", m = 1),       # random smooth effect for Year
data = df,


After fitting the model you might want to check if the current model allows for sufficient complexity/wigglyness of the relationship. For this you can use the gam.check() function. It will help you decide whether you should increase the argument k of the s() function in the model (default value is 10, if too low try 20 and check again). While you are at it, you can also look at the plots that the function produces in order to check for potential heteroskedasticity (although I am not sure if the latter works for binomial models).

There are some which suggest making the model more parsimonious by removing unnecessary random effects. You can do that by looking at the p-value of the random smooth term in the model summary. A low p-value (e.g. smaller than .05) suggests that the random smooth might be helpful, a higher p-value suggests that it might not be necessary.

If running the model takes too long you can try the following: a) set discrete to TRUE to use a discretetized version of covariates, b) set nthreads to the number cores of your computer to run the model using multiple CPU-cores.

For answering your main question, I suggest looking at the p-value of the smooth term (that of s(center_sal)) in the model summary.

For visually inspecting the relationship you could either use the built-in plot() function of mgcv or the plot_smooth() function from the itsadug package. The latter allows for using simultaneous confidence intervals which might be useful, depending on the questions you want to answer.

For a general reference on GAMMs I suggest looking at the book of the author of mgcv:

• Wood, S. (2017) Generalized additive models: An introduction with R, 2nd. Ed.

Another interesting reference that provides a (short) introduction and offers some insights into which specification of random effects might be appropriate is the following:

• Soskuthy, M. (2021) - Evaluating generalised additive mixed modelling strategies for dynamic speech analysis

Finally I also found the articles and blog of G. Simpson to be a useful resource. To mention just one blog post: https://fromthebottomoftheheap.net/2021/02/02/random-effects-in-gams/

• An additive model is an interesting idea, but screening based on p-values distorts all downstream inferences unless considerable care is taken.
– Dave
Dec 19, 2022 at 13:45
• Absolutely, I understood the OPs question of deciding whether sal predicts empty to be the primary question and therefore not as screening followed by other inferential goals. Instead of looking at the p-value of the smooth one could also use a model comparison approach by looking at AIC values of models with and without sal as a predictor. However, that requires additional steps and potentially additional modeling decisions. How would you put "taking considerable care" into practice if the OPs intention is a screening? Dec 19, 2022 at 13:59
• How to be careful with stepwise regression sounds like a reasonable question to post!
– Dave
Dec 19, 2022 at 14:01
• Cross-validation approaches might come to mind (screen on part of the data, test on another). Dec 19, 2022 at 14:08
• Treating p-values close to the threshold differently than p-values far from the threshold might be another way to be cautious. Yet, I do not think it directly targets the screening plus downstream inference problem. But I guess treating p-values like that should be standard practice anyway. In GAMMs especially the p-values of smooths are "approximate" p-values anyway. For this the literature suggests to always specially "flag and report" p-values close to the threshold. Dec 19, 2022 at 14:11

Nonlinear basis functions are perfectly acceptable for generalized linear models like logistic regressions. You’re probably familiar with polynomial models, such as including a quadratic term in a linear regression. For modeling more bizarre curvature like you seem to have, splines could be your friend.

For reputable sources, Frank Harrell’s blog and Regression Modeling Strategies textbook talk about using nonlinear basis functions (chiefly splines) to model nonlinear relationships, and his rms R package implements splines through the rms::rcs function for so-called restricted cubic splines.

• I’m not sure about random effects in the rms package, but I don’t think you need to call rms::rcs in the package regression functions. I would expect the function to work within glmer.
– Dave
Dec 19, 2022 at 6:06
• Sorry, I'm confused. Does the smooth function s() mean that center_sal no longer needs to be linearly associated with the log odds (it can have a non-linear trend)? Only in a glmm is the linear association required?
– Nate
Dec 19, 2022 at 17:47

but discovered there is a non-linear relationship between the predictor variable and the log odds of the response

That plot only shows a non-linear relationship for the entire population. You might still have linear relationships within subgroups of the population. For instance you could have something like the following image

In the image this non-linear relationship for the entire population arises because other predictor variables (shown with colors in the image) correlate with the predictor variable in the plot.

• Did some plotting and the various "colors" are quite well mixed, no discernable (linear) pattern. Thank you though!
– Nate
Dec 20, 2022 at 15:25
• @Nate, in your graph you have a linear pattern by definition because your model to make it is linear (it's the other parameters that make the U-shaped relationship). It might be less visible because you have more groups and smaller groups. Anyway, maybe it does make sense to add a non-linear term in the model, but the current graph in the question does not show this very well. If you would compute another model with an additional quadratic term, then you can use an ANOVA or likelihood ratio test to see whether the extra non-linear term makes an improvement. Dec 20, 2022 at 15:27