1

Sorry for this, as I am new to R. I noticed that rms::calibrate with backward generated factors in final model and their coefficients are quite different from glm with step() function under the same criteria or the lrm() with the final factors.

Are coefficients of the final factors in "Approximate Estimates after Deleting Factor" different from the coefficients by these final factors in a lrm function?

I have already noticed to use the tyep = "individual" for the fast backward to behave like step. But the coefficient I got are still different.

Here I defined a full model with some variables and did the calibration with backward stepwise by p-value:

mod1.cal <- lrm(Death~ Age + Gender + F + IH + SH + M + A + C + BS + PSY, data = data, x = T, y = T)
set.seed(2020)
cal1 <- rms::calibrate(mod1.cal, method = 'boot', B = 100, bw = T, rule = 
'p', sls = 0.05, type = 'individual', data = data)

The results were:

        Backwards Step-down - Original Model

        Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
PSY     0.18   1    0.6716 0.18     1    0.6716 -1.82
F       2.84   1    0.0921 3.02     2    0.2212 -0.98

Approximate Estimates after Deleting Factors

             Coef     S.E.  Wald Z         P
Intercept    -7.74397 0.572982 -13.515 0.000e+00
Age           0.06096 0.006425   9.489 0.000e+00
Gender=Male   0.86812 0.202052   4.297 1.735e-05
IH=Present    0.52113 0.193214   2.697 6.994e-03
SH=Present    0.70469 0.196477   3.587 3.350e-04
M=Presnet     0.88153 0.286214   3.080 2.070e-03
A=Present     0.60711 0.213466   2.844 4.454e-03
C=Present    -1.02122 0.324155  -3.150 1.630e-03
BS=Abnormal   1.79088 0.291809   6.137 8.401e-10

Factors in Final Model

[1] Age    Gender IH    SH    M    A   C   BS

Then I applied the factos in final model again with lrm for a model:

mod1.final <- lrm(Death~ Age + Gender + IH + SH + M + A + C + BS, data = data, x = T, y = T)

And the results were:

Logistic Regression Model

lrm(formula = Death ~ Age + Gender + IH + SH + M + A + C + BS, 
data = data, x = T, y = T)

                   Model Likelihood     Discrimination    Rank Discrim.    
                      Ratio Test           Indexes           Indexes       
Obs          1395    LR chi2     326.85    R2       0.365    C       0.850    
  Alive       1184    d.f.             8    g        1.964    Dxy     0.700    
  Death        211    Pr(> chi2) <0.0001    gr       7.126    gamma   0.700    
max |deriv| 1e-10                          gp       0.178    tau-a   0.180    
                                        Brier    0.095                     

          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept    -7.7719 0.5702 -13.63 <0.0001 
Age           0.0612 0.0064   9.57 <0.0001 
Gender=Male   0.8708 0.2015   4.32 <0.0001 
IH=Present    0.5217 0.1926   2.71 0.0068  
SH=Present    0.7047 0.1968   3.58 0.0003  
M=Present     0.8856 0.2874   3.08 0.0021  
A=Present     0.6071 0.2131   2.85 0.0044  
C=Present    -1.0251 0.3252  -3.15 0.0016  
BS=Abnormal   1.7998 0.2920   6.16 <0.0001

And I checked the coef of the mod1.final:

coef(mod1.final)
Intercept     Age          Gender=Male  IH=Present   SH=Present   M=Present 
A=Present     C=Present    BS=Abnormal 
-7.7719073    0.0612293    0.8707690    0.5216842    0.7046757    0.8855560    
 0.6070599   -1.0250605    1.7997968

Can someone help on this.

Thanks

Add up the results from validation for questions:

          index.orig training    test optimism index.corrected   n
 Dxy           0.7889   0.8004  0.7752   0.0252          0.7637 100
 R2            0.4818   0.5008  0.4645   0.0363          0.4455 100
 Intercept     0.0000   0.0000 -0.0986   0.0986         -0.0986 100
 Slope         1.0000   1.0000  0.9105   0.0895          0.9105 100
 Emax          0.0000   0.0000  0.0387   0.0387          0.0387 100
 D             0.3220   0.3361  0.3084   0.0277          0.2943 100
 U            -0.0014  -0.0014  0.0015  -0.0030          0.0015 100
 Q             0.3235   0.3376  0.3069   0.0307          0.2928 100
 B             0.0802   0.0768  0.0824  -0.0056          0.0858 100
 g             2.4779   2.6190  2.3844   0.2346          2.2433 100
 gp            0.2016   0.2041  0.1984   0.0056          0.1959 100
Stephan Kolassa
  • 95,027
  • 13
  • 197
  • 357
Orson
  • 37
  • 5

1 Answers1

2

The answer to your question is contained in the title. Backward step-down as implemented by the fastbw() function used by calibrate() in the rms package "prints approximate parameter estimates for the model after deleting variables," as the help page says. If you were doing ordinary least squares the coefficients themselves would be exact, but with logistic regression neither the coefficients nor the standard errors reported are exact.

You will notice that the "approximate" coefficients are quite close to those for the lrm() model, differing only in their third significant figures, so this is not a big practical problem.

Be sure to pay attention to the validate() and calibrate() results themselves, as they provide important information about the likely generalizability of your model when applied to new data. My guess is that the "optimism" in the coefficient estimates (provided by validate()) will be much greater than these very minor differences you found in the regression coefficients.

Finally, as "Death" is your outcome variable do consider whether you should be performing survival analysis that considers time to death, not just whether or not death has occurred.

EdM
  • 57,766
  • 7
  • 66
  • 187
  • Thank you so much. Hope to clearify me following doubts: 1) about the model for final presentation. I should still present the coefficients of the logistic model rather than the coefficients from the calibrate() or validate ()? 2) I added the results from the validate(). does the Slope 0.9105 refers to the "optimism" you mentioned in the coefficient estimate? And if I want to shrink my coefficients, can I use this slope to times with the coefficients of the factors in logistic regression directly? or I should calculate with the coefficients from the calibrate() or validate() results? – Orson Mar 16 '20 at 03:36
  • 1
    @Orson if your main interest is prediction and you are willing to use shrinkage then you should consider ridge regression (with shrinkage built in) including all predictors instead of backward selection. [Harrell says](http://hbiostat.org/doc/rms.pdf) in Section 4.11: "Shrinkage is far superior to variable selection." Including all predictors avoids the omitted-variable bias that is [intrinsic in logistic regression](https://stats.stackexchange.com/q/113766/28500) if an omitted predictor (like your `F`) might have a relationship with outcome even if it isn't "significant" in your data set. – EdM Mar 16 '20 at 14:47
  • @Orson for reporting it won't matter much in your case as both point estimates and standard errors are very close for the "approximate" and final models. I'd report `mod1.final` while making it clear that the model was developed by backward elimination and reporting the slope optimism of 8.95%. Calling `plot(cal1)` will give you a calibration curve that might be more reliable than just shrinking coefficients directly to correct for optimism. – EdM Mar 16 '20 at 14:56
  • Thank you so much for clearing these doubts. Really helpful:) You are right, I am also using ridge and lasso here for comparisons. I am using lasso because I still want to make a model simpler for some clinical interpretation and application. – Orson Mar 16 '20 at 15:46
  • Hope to ask a bit more about the λ in lasso. For another model by lasso, I splitted my data here by time with the first 3 years as training data and the last year for a temporal validation. I found when using λ(1se), val.prob gave me a slope of 1.47 with intercept of 0.87 while λ(min) gave me much better results with a slope of 1.03 with intecept of 0.15. I saw lots of people use 1se over min, as it usually gives a much simpler model. I suppose I should go with min here? – Orson Mar 16 '20 at 15:49
  • And a minor quesiton. When using cv.glmnet() for binary logistic regression (family="binomial"), should the type.measure = "class" or "deviance"? As I found both the 1se and min λ are quite different under these two settings. I tried both, found the val.prob on the last year gave slightly different slopes. And the type.measure = "class" seems to show slightly better results in slope and intercept close to the ideal line on the final year temporal validation. – Orson Mar 16 '20 at 16:19
  • @Orson see [this page](https://stats.stackexchange.com/q/138569/28500) for min versus 1se for $\lambda$. See [this page](https://stats.stackexchange.com/q/255747/28500) for why a continuous measure like deviance is preferred to misclassification error; I think that "class" assumes a probability cutoff of 0.5, which isn't always appropriate depending on your tradeoff between false-positives and false-negatives. If you still have questions after looking at these links, please post as a new question to make it easier for later visitors to the site to find what they are looking for. – EdM Mar 16 '20 at 17:51
  • Get it. Thank you:) – Orson Mar 17 '20 at 00:33