How do I know if my logistic model's numeric variables violate the linearity assumption on the log(target variable)?
I am using the code found at THIS LINK. They graphed on the same data but I, personally, only used my training data.
The variable age and pedigree is not linear and might need some transformations.
I'm ready to do so, but I'd first like to know why it is that they consider DiabetesPedigreeFunction & age as non-linear? What shape would you say qualifies the variable as linear?
You could use a smoothing spline approach (proposed by @AdamO) in a higher-dimensional context as well. The conditioning plots are fine if there are only three independent variables, but in a situation like the one you propose several different conditioning plots would be required and it would be difficult to derive much actionable intelligence from them. Instead, you could use a generalized additive model (GAM) to model flexible non-linearity on the link scale. Here's an example using the mgcv package in R.
First, load the data, specify the independent variables and make the formulas for the models. Note that wrapping a variable in s() changes the form from linear to (potentially) non-linear by using a smoothing spline to estimate the relationship with y.
library(tidyverse)
library(broom)
library(mgcv)
data("PimaIndiansDiabetes2", package = "mlbench")
PimaIndiansDiabetes2 <- na.omit(PimaIndiansDiabetes2)
# Fit the logistic regression model
ivs <- c("pregnant", "glucose", "pressure", "triceps", "insulin", "mass", "pedigree", "age")
forms <- sapply(ivs, function(x)
reformulate(ifelse(ivs == x, paste0("s(", ivs, ")"), ivs),
"diabetes"))
Next, we can estimate the null model where all effects are linear.
mod.null <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
insulin + mass + pedigree + age,
data = PimaIndiansDiabetes2,
family = binomial)
We could then get the AIC differences of each model to the null model.
smods <- lapply(forms, function(x)gam(x, data=PimaIndiansDiabetes2, family=binomial))
sort(sapply(smods, AIC) - AIC(mod.null))
#> age insulin mass pedigree triceps
#> -6.9494706510 -4.1703198157 -2.1917412788 -0.4979687138 0.0001220039
#> pressure glucose pregnant
#> 0.0001247168 0.0001520076 0.0001769288
Here, non-linear terms in age, insulin, mass and pedigree are improvements over the null model, though modeling non-linearity in pedigree doesn't seem to improve the model much. We could then start with the biggest improvement and then add in other non-linear trends as they are necessary. We cannot just automatically add non-linear trends for age, insulin and mass because they modeling unobserved non-linearity for one variable may change the way other variables relate to y. So, we start with the age variable.
m1 <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
insulin + mass + pedigree + s(age),
data = PimaIndiansDiabetes2,
family = binomial)
anova(mod.null, m1, test='Chisq')
#> Analysis of Deviance Table
#>
#> Model 1: diabetes ~ pregnant + glucose + pressure + triceps + insulin +
#> mass + pedigree + age
#> Model 2: diabetes ~ pregnant + glucose + pressure + triceps + insulin +
#> mass + pedigree + s(age)
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 383.00 344.02
#> 2 379.12 331.30 3.8846 12.725 0.01158 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We see from above that allowing a non-linear functional form for age produces a significantly better model. Moving on, we will test the linearity of the effect of insulin. One important caveat here is that we need to ensure that the only thing changing across the two models is the functional form for insulin. As such, we must hold fixed what would otherwise be a potentially changing non-linear relationship between age and y. To do this, we specify the degrees of freedom the relationship will use (k is the degrees of freedom + 1) and use fx=TRUE in the call to s(). This will ensure that the non-linearity in age will be of the same form across the two models.
summary(m1)
#>
#> Family: binomial
#> Link function: logit
#>
#> Formula:
#> diabetes ~ pregnant + glucose + pressure + triceps + insulin +
#> mass + pedigree + s(age)
#>
#> Parametric coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -8.3908109 1.2235167 -6.858 6.99e-12 ***
#> pregnant 0.0402072 0.0571396 0.704 0.4816
#> glucose 0.0375236 0.0058529 6.411 1.44e-10 ***
#> pressure -0.0042853 0.0119421 -0.359 0.7197
#> triceps 0.0099064 0.0173894 0.570 0.5689
#> insulin -0.0003452 0.0013406 -0.257 0.7968
#> mass 0.0627002 0.0274747 2.282 0.0225 *
#> pedigree 1.0854046 0.4313465 2.516 0.0119 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df Chi.sq p-value
#> s(age) 3.888 4.885 13.12 0.0215 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.375 Deviance explained = 33.5%
#> UBRE = -0.094205 Scale est. = 1 n = 392
m1f <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
insulin + mass + pedigree + s(age, k=5, fx=TRUE),
data = PimaIndiansDiabetes2,
family = binomial)
m2 <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
s(insulin) + mass + pedigree + s(age, k=5, fx=TRUE),
data = PimaIndiansDiabetes2,
family = binomial)
anova(m1f, m2, test='Chisq')
#> Analysis of Deviance Table
#>
#> Model 1: diabetes ~ pregnant + glucose + pressure + triceps + insulin +
#> mass + pedigree + s(age, k = 5, fx = TRUE)
#> Model 2: diabetes ~ pregnant + glucose + pressure + triceps + s(insulin) +
#> mass + pedigree + s(age, k = 5, fx = TRUE)
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 380.00 331.61
#> 2 372.22 313.97 7.7841 17.646 0.02141 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note here that the non-linearity in insulin is also significant - indicating that the model significantly improves when the functional form of the insulin variable is made more flexible. Following the same process as above, we can test for non-linearity in mass.
summary(m2)
#>
#> Family: binomial
#> Link function: logit
#>
#> Formula:
#> diabetes ~ pregnant + glucose + pressure + triceps + s(insulin) +
#> mass + pedigree + s(age, k = 5, fx = TRUE)
#>
#> Parametric coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -8.599060 1.381969 -6.222 4.90e-10 ***
#> pregnant 0.057864 0.059608 0.971 0.33168
#> glucose 0.034614 0.006207 5.577 2.45e-08 ***
#> pressure -0.001621 0.012620 -0.128 0.89782
#> triceps 0.007756 0.018136 0.428 0.66889
#> mass 0.065503 0.029002 2.259 0.02391 *
#> pedigree 1.229714 0.446129 2.756 0.00584 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df Chi.sq p-value
#> s(insulin) 8.1 8.784 13.13 0.106
#> s(age) 4.0 4.000 11.35 0.023 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.399 Deviance explained = 37%
#> UBRE = -0.10161 Scale est. = 1 n = 392
m2f <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
s(insulin, k=9, fx=TRUE) + mass + pedigree + s(age, k=5, fx=TRUE),
data = PimaIndiansDiabetes2,
family = binomial)
m3 <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
s(insulin, k=9, fx=TRUE) + s(mass) + pedigree + s(age, k=5, fx=TRUE),
data = PimaIndiansDiabetes2,
family = binomial)
anova(m2f, m3, test='Chisq')
#> Analysis of Deviance Table
#>
#> Model 1: diabetes ~ pregnant + glucose + pressure + triceps + s(insulin,
#> k = 9, fx = TRUE) + mass + pedigree + s(age, k = 5, fx = TRUE)
#> Model 2: diabetes ~ pregnant + glucose + pressure + triceps + s(insulin,
#> k = 9, fx = TRUE) + s(mass) + pedigree + s(age, k = 5, fx = TRUE)
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 373.00 315.00
#> 2 370.43 310.64 2.5729 4.3605 0.1735
Here, the difference is not statistically significant - indicating that the relationship is not significantly improved when non-linearity is allowed. Retaining the linear form (on the link scale) is appropriate. We can move on and do the same for pedigree.
m4 <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
s(insulin, k=9, fx=TRUE) + mass + s(pedigree) + s(age, k=5, fx=TRUE),
data = PimaIndiansDiabetes2,
family = binomial)
anova(m2f, m4, test='Chisq')
#> Analysis of Deviance Table
#>
#> Model 1: diabetes ~ pregnant + glucose + pressure + triceps + s(insulin,
#> k = 9, fx = TRUE) + mass + pedigree + s(age, k = 5, fx = TRUE)
#> Model 2: diabetes ~ pregnant + glucose + pressure + triceps + s(insulin,
#> k = 9, fx = TRUE) + mass + s(pedigree) + s(age, k = 5, fx = TRUE)
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 373 315
#> 2 373 315 0.00023655 0.00013726 0.001065 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
This is an interesting case because the result is significant, but as you can see the difference in residual degrees of freedom is about 0.0002 - indicating a model that is very nearly linear. Despite a significant result here, the linear relationship (on the link scale) should be retained. We can then estimate, summarize and plot the final model.
mod.final <- gam(diabetes ~ pregnant + glucose + pressure + triceps +
s(insulin) + mass + pedigree + s(age),
data = PimaIndiansDiabetes2,
family = binomial)
summary(mod.final)
#>
#> Family: binomial
#> Link function: logit
#>
#> Formula:
#> diabetes ~ pregnant + glucose + pressure + triceps + s(insulin) +
#> mass + pedigree + s(age)
#>
#> Parametric coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -8.7741665 1.3915436 -6.305 2.88e-10 ***
#> pregnant 0.0573195 0.0594034 0.965 0.33459
#> glucose 0.0357828 0.0061995 5.772 7.84e-09 ***
#> pressure -0.0007469 0.0126257 -0.059 0.95283
#> triceps 0.0069614 0.0180179 0.386 0.69923
#> mass 0.0659178 0.0290142 2.272 0.02309 *
#> pedigree 1.3001933 0.4485428 2.899 0.00375 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df Chi.sq p-value
#> s(insulin) 8.287 8.857 13.82 0.0921 .
#> s(age) 4.409 5.429 11.48 0.0554 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.406 Deviance explained = 37.4%
#> UBRE = -0.10374 Scale est. = 1 n = 392
par(mfrow=c(1,2))
plot(mod.final, select=1)
plot(mod.final, select=2)
Created on 2022-12-29 by the reprex package (v2.0.1)
Looking at the rug plots, you may want to re-evaluate whether these non-linear relationships are the potential result of some outliers (i.e., more regression diagnostics), but a procedure like this could provide more concrete advice about the appropriate functional forms of variables.
Glucose and Insulin are highly correlated. In my analysis, I dropped Insulin due to this information. Given that you've downloaded and inspected the data, do you see anything wrong with this step I took?
$\endgroup$
glucose and insulin are correlated at 0.58 and their VIFs are 1/67 and 1.56, respectively. I probably wouldn't remove a variable based on those numbers. In fact the highest VIF is for age at 2.13 - still probably not high enough to merit removing it. I would look at how influential the insulin values over 600 and age over 60 are.
$\endgroup$
Commented
Dec 30, 2022 at 16:37
Insulin due to the higher number of NA values, as well (48%), thinking that the lack of variance provided from the NA values impacts models on unseen data. Is it wise to still include it or does it come down to choice?
$\endgroup$
Insulin is a good idea or not. This is a bigger question you should ask (or search for other answers) - it probably shouldn't be answered in the comment of this question.
$\endgroup$
Commented
Dec 31, 2022 at 15:13
You're accidentally asking more than one question. The first is whether the trend relating a continuous exposure and a binary outcome. One possible way to inspect such a trend is to create a non-parametric smoothing function, such as a smoothing spline.
x <- seq(-3, 3, by=0.1)
y <- rbinom(length(x), 1, plogis(-2 + 1*x))
plot(x,y)
lines(smooth.spline(x=x, y=y), col='red')
lines(x, predict(glm(y ~ x, family=binomial), type='response'), col='blue')
You can also plot the GLM object itself, and pay particular attention to the plot of Pearson residuals versus fitted. Despite the infamous appearance, you can expect that the mean should be 0 and the variance should be constant across all predicted values. They give you a Loess smooth line to deal with this problem.
Departures from this could signify a misspecified logistic model.
The second question is whether the multivariate association holds, that is controlling for other predictors, U, V, W, does the X,Y relationship meet certain criteria (similarly, controlling for X,V,W does the U,Y relationship ... and so on...). For some reason, people seem to think that plotting the bivariate associations - like you've shown here - does something. Maybe, but dimensionality is complicated that way. A more powerful tool to inspect these relations is a coplot.
w <- abs(seq(-3, 3, by=0.1))
y2 <- rbinom(length(x), 1, plogis(-0.3 + 1 * x - 0.2*w))
mypanel <- function(x, y, ...) panel.smooth(x=x, y=y, span=1/2, ...)
coplot(y ~ x | w, panel=mypanel)
Notice what a terrible smoother LOESS is compared to a smoothing spline for binary response - this is because LOESS excludes outliers. Anyway, coplots work well for multivariate data because you expect the plots to be identical when there is no interaction between covariates. And, as above, you expect the smoother to show a nice logistic trend line relating the exposure to response.
