Check linearity in logistic regression

user_13321647 0
Viewed 322

In order to check linearity in logistic regression ->

Is independent1 and independent2variable linear related to the log-odds of depdendent?

I would like optimize this (working) calculations:

This is the code:

# Check Linearity  ---------------------------------------------------------

# quartiles of independent1
quantile(df$independent1, probs=c(0, 0.25, 0.5, 0.75, 1))


table(df$dependent[df$independent1<52])
table(df$dependent[df$independent1>=52 & df$independent1 < 60])
table(df$dependent[df$independent1>=60 & df$independent1 < 73])
table(df$dependent[df$independent1>=73 & df$independent1 < 91])

p1 <- mean(df$dependent[df$independent1<52])
p2 <- mean(df$dependent[df$independent1>=52 & df$independent1 < 60])
p3 <- mean(df$dependent[df$independent1>=60 & df$independent1 < 73])
p4 <- mean(df$dependent[df$independent1>=73 & df$independent1 < 91])
probs <- c(p1, p2, p3, p4)

# calculate the log-odds
logits <- log(probs/(1-probs))

# quartiles of independent1
q <- quantile(df$independent1, probs=seq(0,1,0.25))

# calculate median independent1 for each of the 4 groups
meds <- c( median(df$independent1[ df$independent1<q[2]]),
           median(df$independent1[ df$independent1>=q[2] & df$independent1<q[3]]),
           median(df$independent1[ df$independent1>=q[3] & df$independent1<q[4]]),
           median(df$independent1[ df$independent1>=q[4]])
)

plot(meds, logits, main="xxx",
     xlab = "independent1",
     ylab = "log-odds(dependent|independent1)", las=1)

For one variable this might be ok. But I have more independent variables. So how could I optimize this code (checking and plotting) for each independent variable (in this example independent1 and independent2)

My dataframe:

df <- structure(list(dependent = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), independent1 = c(84, 
49, 54, 75, 49, 70, 75, 42, 60, 72, 80, 73, 51, 61, 59, 78, 45, 
38, 78, 65, 91, 60, 39, 31, 42, 72, 41, 77, 73, 74, 39, 86, 71, 
55, 43, 75, 80, 75, 67, 74, 46, 70, 57, 66, 57, 72, 46, 52, 53, 
76, 57, 86, 67, 71, 57, 50, 76, 61, 41, 57, 62, 41, 64, 82, 53, 
75, 59, 38, 54, 56, 68, 63, 73, 26, 75, 76, 81, 46, 77, 53, 59, 
66, 51, 72, 80, 70, 39, 57, 62, 85, 84, 57, 73, 55, 70, 78, 66, 
69, 60, 51, 72, 68, 60, 62, 64, 44, 50, 59, 45, 81, 54, 68, 75, 
66, 54, 45, 52, 87, 44, 77, 49, 84, 68, 76, 82, 44, 58, 55, 69, 
33, 48, 62, 60, 76, 56, 73, 55, 58, 53, 53, 60, 52, 60, 41, 39, 
36, 38, 59, 54, 64), independent2 = c(23, 25, 34, 25, 31, 25, 
32, 19, 25, 28, 22, 18, 30, 26, 25, 25, 25, 19, 24, 27, 23, 28, 
39, 27, 30, 28, 22, 28, 25, 23, 18, 27, 27, 19, 25, 27, 26, 26, 
21, 26, 23, 28, 37, 32, 24, 32, 26, 23, 24, 27, 28, 25, 24, 22, 
34, 23, 35, 20, 29, 29, 21, 29, 25, 26, 23, 33, 25, 26, 29, 27, 
26, 28, 19, 22, 29, 22, 26, 35, 32, 29, 26, 23, 31, 30, 27, 28, 
23, 27, 34, 22, 24, 28, 21, 25, 18, 32, 21, 24, 31, 31, 24, 30, 
27, 23, 16, 26, 26, 19, 38, 21, 32, 34, 28, 19, 30, 24, 26, 24, 
40, 26, 15, 26, 28, 22, 25, 26, 31, 24, 26, 42, 26, 30, 28, 21, 
21, 19, 22, 20, 26, 31, 22, 25, 21, 20, 27, 27, 26, 29, 22, 24
)), row.names = c(NA, -150L), class = c("tbl_df", "tbl", "data.frame"
))
logistic-regression r
1 Answers

I'll demonstrate a somewhat different and decidedly more efficient method of splitting a variable that is to be used in a logistic regression model:

df$q41 <- with(df, cut(independent1, quantile(independent1), include = TRUE))
 # creates 4 level factor into roughly equally sized groups
 table(df$q41)
#--------------------
#[26,52] (52,60] (60,73] (73,91] 
#     39      37      39      35 

#Examine for "eyeball" trends in the log-odds of dependent 
fit1.q41 <- glm(dependent~q41+0, data=df, fam="binomial")
fit1.q41
#---------------------------
Call:  glm(formula = dependent ~ q41 + 0, family = "binomial", data = df)

Coefficients:
q41[26,52]  q41(52,60]  q41(60,73]  q41(73,91]  
    -3.638      -2.862      -2.918      -2.048  

Degrees of Freedom: 150 Total (i.e. Null);  146 Residual
Null Deviance:      207.9 
Residual Deviance: 65.52    AIC: 73.52

I chose to remove the intercept term because its presence prevented viewing the coefficient of the lowest group on the same scale as the upper 3. The coefficients are just the logits for the grouping I created. Compare:

> logits
[1] -3.555348 -2.740840 -2.970414 -2.169054
> coef(fit1.q41)
q41[26,52] q41(52,60] q41(60,73] q41(73,91] 
 -3.637586  -2.862201  -2.917771  -2.047693

I then tried to automate the process but ran into a bit of a problem because of the small numbers of events in one of the quartile groups, The ridiculously low coefficient for the lowest quartile in independent2 is from that lack of any events or "1"'s in that catergory. (An estimate log-odds of -19.566069 does rather point to a proportion of 0.)

 lapply( df[-1], function(x){cat(str(x)); IVq <- cut(x, quantile(x), include = TRUE); logits<-coef( summary(glm(df$dependent~IVq+0, fam="binomial"))); logits})
 num [1:150] 84 49 54 75 49 70 75 42 60 72 ...
 num [1:150] 23 25 34 25 31 25 32 19 25 28 ...
$independent1
            Estimate Std. Error   z value     Pr(>|z|)
IVq[26,52] -3.637586  1.0130639 -3.590678 3.298191e-04
IVq(52,60] -2.862201  0.7270292 -3.936845 8.256004e-05
IVq(60,73] -2.917771  0.7259663 -4.019155 5.840732e-05
IVq(73,91] -2.047693  0.5312796 -3.854266 1.160776e-04

$independent2
             Estimate   Std. Error     z value     Pr(>|z|)
IVq[15,23] -19.566069 1639.9716035 -0.01193074 9.904809e-01
IVq(23,26]  -3.091042    0.7229988 -4.27530783 1.908734e-05
IVq(26,28]  -2.397895    0.7385489 -3.24676555 1.167245e-03
IVq(28,42]  -1.856298    0.4808846 -3.86017349 1.133066e-04


> lapply( df[-1], function(x){ IVq <- cut(x, quantile(x), include = TRUE); table(IVq, df$dependent) })
$independent1
         
IVq        0  1
  [26,52] 38  1
  (52,60] 35  2
  (60,73] 37  2
  (73,91] 31  4

$independent2
         
IVq        0  1
  [15,23] 43  0
  (23,26] 44  2
  (26,28] 22  2
  (28,42] 32  5

At any rate I think I've demonstrated a more R-ish approach to calculating logits within quartiles. It also sets you up for a model comparison approach to examining departures from linearity as well as demonstrating possible pitfalls. If you had more events you might have considered looking at the change in deviance from the null model with the addition of a quartile factor on top of a simple linear model ... or even more powerfully using poly to create your comparison model.

In the past when working with datasets with adequate numbers of events I have chosen to split on the basis of quantiles calculated from an event==1 subset rather than letting the splits be based on the whole dataset.

user_18556770
Related