3
$\begingroup$

I'm trying to find out if my numeric predictors have a linear relation to the logit of my logistic regression. I tried to use the lrm fit in the rms package where I have used 3 knot cubic spline on all numeric predictors like this:

> fit <- lrm(y ~ rcs(x1,3)+rcs(x2,3)+.....)

There after I used anova on lrm fit. The main question is how do I use the result in anova(fit)?

My understanding is that the wald statistics are just the associated coefficients squared and dived by it's se. But what about the statistics for nonlinear terms here? are they the wald statistics for the coefficients for the squared predictors?

If none of the statistics are significant, can I conclude that there are no quadratic effect from my predictors?

$\endgroup$
7
$\begingroup$

These are all Wald tests, which assume that the sampling distribution of the vector of $\hat{\beta}$ has a multivariate normal distribution. Only in the special case where one is testing a single parameter does the Wald $\chi^2$ test equal the square of a Wald $z$-statistic; here $z = \frac{\hat{\beta_{j}}}{se}$ for a single coefficient $\beta_j$. The general Wald test is a "chunk test" involving multiple coefficients, and you can generalize this further by considering a general contrast with a null hypothesis of $H_{0}: C\beta = 0$. Some of the things that might be in the "chunks" are

  • nonlinear terms to get a test of linearity
  • nonlinear terms + linear term to get a test of flatness (association)
  • all terms involving one predictor, whether they are main effects or interaction effects, to get a general test of association such as whether age has an association with $Y$ for either sex group

Note that if a test is non-significant, it is not appropriate to remove the tested terms from the model as this causes bias, and especially makes confidence intervals too short and $p$-values too small.

The R rms package anova function makes it easy to see exactly which coefficients are being tested in any line of the ANOVA table. Scroll right to see this information on the far right of each table. For OLS we use $F$ tests instead of $\chi^2$. The model intercept corresponds to a subscript of $\beta$ of zero.

require(rms)
set.seed(123)
age <- rnorm(500, 50, 15)
treat <- factor(sample(c('a','b','c'), 500, TRUE))
bp  <- rnorm(500, 120, 10)
y   <- ifelse(treat=='a', (age-50)*.05, abs(age-50)*.08) + 3*(treat=='c') +
   pmax(bp, 100)*.09 + rnorm(500)

f   <- ols(y ~ treat*lsp(age,50) + rcs(bp,4))
Function(f)   # show algebraic form of fitted model.  Note rcs
              # is simplified so some redundant betas are added
function(treat = NA,age = NA,bp = NA) {-1.5357446+5.4522476*(treat=="b")+7.6742854*(treat=="c")+0.015671819*age+0.049487194*pmax(age-50,0)+0.095699259* bp-4.3486306e-05*pmax(bp-103.28133,0)^3+0.00020843892*pmax(bp-116.59859,0)^3-0.0002067844*pmax(bp-123.63285,0)^3+4.1831786e-05*pmax(bp-137.52664,0)^3+(treat=="b")*(-0.10304059*age+0.11755658*pmax(age-50,0))+(treat=="c")*(-0.084946042*age+0.085581901*pmax(age-50,0)) }

an <- anova(f); options(digits=3)
print(an, 'subscripts')

               Analysis of Variance          Response: y 

 Factor                                     d.f. Partial SS MS      F      P      Tested     
 treat  (Factor+Higher Order Factors)         6  1421.70    236.950 241.73 <.0001 1-2,8-11   
  All Interactions                            4    61.55     15.387  15.70 <.0001 8-11       
 age  (Factor+Higher Order Factors)           6   222.01     37.001  37.75 <.0001 3-4,8-11   
  All Interactions                            4    61.55     15.387  15.70 <.0001 8-11       
  Nonlinear (Factor+Higher Order Factors)     3   156.88     52.294  53.35 <.0001 4,10-11    
 bp                                           3   344.33    114.778 117.09 <.0001 5-7        
  Nonlinear                                   2     1.41      0.706   0.72 0.487  6-7        
 treat * age  (Factor+Higher Order Factors)   4    61.55     15.387  15.70 <.0001 8-11       
  Nonlinear                                   2    22.87     11.436  11.67 <.0001 10-11      
  Nonlinear Interaction : f(A,B) vs. AB       2    22.87     11.436  11.67 <.0001 10-11      
 TOTAL NONLINEAR                              5   157.75     31.550  32.19 <.0001 4,6-7,10-11
 TOTAL NONLINEAR + INTERACTION                7   194.53     27.790  28.35 <.0001 4,6-11     
 REGRESSION                                  11  1861.11    169.192 172.61 <.0001 1-11       
 ERROR                                      488   478.35      0.980                          

Subscripts correspond to:
 [1] treat=b        treat=c        age            age'           bp             bp'            bp''          
 [8] treat=b * age  treat=c * age  treat=b * age' treat=c * age'

print(an, 'dots')

                Analysis of Variance          Response: y 

 Factor                                     d.f. Partial SS MS      F      P      Tested     
 treat  (Factor+Higher Order Factors)         6  1421.70    236.950 241.73 <.0001 ..     ....
  All Interactions                            4    61.55     15.387  15.70 <.0001        ....
 age  (Factor+Higher Order Factors)           6   222.01     37.001  37.75 <.0001   ..   ....
  All Interactions                            4    61.55     15.387  15.70 <.0001        ....
  Nonlinear (Factor+Higher Order Factors)     3   156.88     52.294  53.35 <.0001    .     ..
 bp                                           3   344.33    114.778 117.09 <.0001     ...    
  Nonlinear                                   2     1.41      0.706   0.72 0.487       ..    
 treat * age  (Factor+Higher Order Factors)   4    61.55     15.387  15.70 <.0001        ....
  Nonlinear                                   2    22.87     11.436  11.67 <.0001          ..
  Nonlinear Interaction : f(A,B) vs. AB       2    22.87     11.436  11.67 <.0001          ..
 TOTAL NONLINEAR                              5   157.75     31.550  32.19 <.0001    . ..  ..
 TOTAL NONLINEAR + INTERACTION                7   194.53     27.790  28.35 <.0001    . ......
 REGRESSION                                  11  1861.11    169.192 172.61 <.0001 ...........
 ERROR                                      488   478.35      0.980                          
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.