[R] What is the most useful way to detect nonlinearity in lo

[Frank E Harrell Jr]

(Ted Harding) wrote:
> On 05-Dec-04 Patrick Foley wrote:
> 
>>It is easy to spot response nonlinearity in normal linear
>>models using plot(something.lm).
>>However plot(something.glm) produces artifactual peculiarities
>>since the diagnostic residuals are constrained  by the fact
>>that y can only take values 0 or 1.
>>What do R users find most useful in checking the linearity
>>assumption of logistic regression (i.e. log-odds =a+bx)?
>>
>>Patrick Foley
>>patfoley at csus.edu
> 
> 
> The "most useful way to detect nonlinearity in logistic
> regression" is:
> 
> a) have an awful lot of data
> b) have the x (covariate) values judiciously placed.
> 
> Don't be optimistic about this prohlem. The amount of
> information, especially about non-linearity, in the binary
> responses is often a lot less than people intuitively expect.
> 
> This is an area where R can be especially useful for
> self-education by exploring possibilities and simulation.
> 
> For example, define the function (for quadratic nonlinearity):
> 
>   testlin2<-function(a,b,N){
>     x<-c(-1.0,-0.5,0.0,0.5,1.0)
>     lp<-a*x+b*x^2; p<-exp(lp)/(1+exp(lp))
>     n<-N*c(1,1,1,1,1)
>     r<-c(rbinom(1,n[1],p[1]),rbinom(1,n[2],p[2]),
>          rbinom(1,n[3],p[3]),rbinom(1,n[4],p[4]),
>          rbinom(1,n[5],p[5])
>         )
>     resp<-cbind(r,n-r)
>     X<-cbind(x,x^2);colnames(X)<-c("x","x2")
>     summary(glm(formula = resp ~ X - 1,
>             family = binomial),correlation=TRUE)
>   }
> 
> This places N observations at each of (-1.0,0.5,0.0.5,1.0),
> generates the N binary responses with probability p(x)
> where log(p/(1-p)) = a*x + b*x^2, fits a logistic regression
> forcing the "intercept" term to be 0 (so that you're not
> diluting the info by estimating a parameter you know to be 0),
> and returns the summary(glm(...)) from which the p-values
> can be extracted:
> 
> The p-value for x^2 is testlin2(a,b,N)$coefficients[2,4]}
> 
> You can run this function as a one-off for various values of
> a, b, N to get a feel for what happens. You can run a simulation
> on the lines of
> 
>   pvals<-numeric(1000);
>   for(i in (1:1000)){
>     pvals[i]<-testlin2(1,0.1,500)$coefficients[2,4]
>   }
> 
> so that you can test how often you get a "significant" result.
> 
> For example, adopting the ritual "sigificant == P<0.05,
> power = 80%", you can see a histogram of the p-values
> over the conventional "significance breaks" with
> 
> hist(pvals,breaks=c(0,0.01,0.03,0.1,0.5,0.9,0.95,0.99,1),freq=TRUE)
> 
> and you can see your probability of getting a "significant" result
> as e.g. sum(pvals < 0.05)/1000
> 
> I found that, with testlin2(1,0.1,N), i.e. a = 1.0, b = 0.1
> corresponding to log(p/(1-p)) = x + 0.1*x^2 (a possibly
> interesting degree of nonlinearity), I had to go up to N=2000
> before I was getting more than 80% of the p-values < 0.05.
> That corresponds to 2000 observations at each value of x, or
> 10,000 observations in all.
> 
> Compare this with a similar test for non-linearity with
> normally-distributed responses [exercise for the reader].
> 
> You can write functions similar to testlin2 for higher-order
> nonlinearlities, e.g. testlin3 for a*x + b*x^3, testlin23 for
> a*x + b*x^2 + c*x^3, etc., (the modifications required are
> obvious) and see how you get on. As I say, don't be optimistic!
> 
> In particular, run testlin3 a few times and see the sort of
> mess that can come out -- in particular gruesome correlations,
> which is why "correlation=TRUE" is set in the call to
> summary(glm(...),correlation=TRUE).
> 
> Best wishes,
> Ted.

library(Design)  # also requires Hmisc
f <- lrm(sick ~ sex + rcs(age,4) + rcs(blood.pressure,5))
# Restricted cubic spline in age with 4 knots, blood.pressure with 5
anova(f)  # automatic tests of linearity of all predictors
latex(f)  # see algebraic form of fit
summary(f)# get odds ratios for meaningful changes in predictors

But beware of using the tests of linearity.  If non-significant results 
cause you to reduce an effect to linear, confidence levels and type I 
errors are no longer preserved.  I use tests of linearity mainly to 
demonstrate that effects are more often nonlinear than linear, given 
sufficient sample size.  I.e., good analysts are needed.  I usually 
leave non-significant nonlinearities in the model.

-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University