Just like when running `anova()`

, the fact that you want to quantify which predictors are very important doesn’t mean that should lead to changes in the model. You are getting into the classic Grambsch-O’Brien problem discussed in RMS.

Dear all,

I’m developing a prediction model using different strategies (N = 122; DF = 10; continuous discrete response; using `lrm()`

):

- Full model unpenalized;
- Full model penalized (maximization of AIC);
- Full model penalized (penalty yielding effective DF of approx 122/15 = 8);
- Reduced pre-specified model (built using background knowledge);
- Reduced model (built with backwards elimitation of predictors using AIC);
- Reduced model from 5. estimated with shrinkage of coefficients using penalty from either 2. or 3.

I’m validating all these models with bootstrap and also producing calibration measures with bootstrap.

**Questions:**

Which performance indexes are most usefull/adequate for comparing these models? E.g., *Emax, Slope, Mean Absolute Error* etc

Which of these models is the most appropriate for making inferences about the relationships between predictor and response?

Thank you

That’s too many models, lowering the chance that AIC can select the best model with reasonable probability. I’d go with subject matter knowledge to the extent possible. Note that AIC is not valid in the stepwise variable selection context. I’d also make sparse principal components part of the game plan.

Good performance metrics into E90 and those described here.

Hi Professor, thank you for the comments/suggestions.

Actually, I was thinking of going with only model 4 in the paper (my firtst choice really). I would fit the other models and show their performance (calibrationwise) just for comparative purposes and to show the audience/readership how models like 1 and 5 are usually not appropriate for the task.

Please note that I’m not comparing these models against each other on the basis of their AIC. I used AIC only to estimate the penalty factor for the full penalized model (model 2) and to perform backwards variable selection (`fastbw(fit, rule='p', type='individual', sls=0.157`

) (model 5). Then I applied the penalty from model 2, or the penalty yielding an effective DF of 8 (model 3), to this reduced selected model*, resulting in model 6.

*as suggested by Steyerberg (page 259 of the 2nd edition of his book, *Clinical Prediciton Models*):

“For penalized estimates of the regression coefficients after selection, we can

apply the penalty factor that was identified as optimal for the full model, before

selection took place”

By *E90* do you mean the *0.9 Quantile of Absolute Error*, as reported by `calibrate()`

?

Thank you

Yes on the last question. And yes that makes sense if the penalty is estimated for the biggest model.

Thank you Professor. Much appreciated feedback.

Just to illustrate, here’s a plot with the E90 measures from each of these models.

The model yielding the lowest error is the model selected with `fastbw(fit, rule='p', type='individual', sls=0.157`

and then validated taking selection into account (4th model from left to right on the plot). This model outperforms the pre-specified reduced model (last model on the plot). Interestingly, if I take the selected 4th model and apply the penalty factor from the full penalized model, this results in worse peformance (5th model on the plot).

Is it safe to perform inference with this 4th (backwards step-down, validated) model?

By inference I mean interpreting *partial effects plots* and *contrasts*.

Minor suggestion: horizontal dot charts allow you to more easily read the results because the model labels will be in larger font next to the dots they represent.

I don’t know the answer to your question because

- Calibration, though super important, is only one dimension of accuracy
- You are comparing too many models to be able to avoid damaging model uncertainty

Ok,

Thank you for the great suggestion!

Then I believe that the safest/most appropriate approach here would be not to compare that many models, but simply fit a full penalized model. Perhaps compare (using AIC) 2 or 3 versions of this full model with varying degrees of complexity (in terms of interaction and nonlinear terms) and then select the best one.

**The question now becomes:** What degree of penalization to use?

The one that maximizes AIC or the one that yields an effective DF of approximately N/15, where N is the sample size. As we can see on the plot above, the latter aproach (which corresponds to less penalization) seems better (at least with regards to the E90 index).

Take a look at chapter 14 of RMS where I use effective AIC to choose a penalty. But beware that there is a sample size requirement to consider for either effective AIC or cross-validation to be reliable for selecting shrinkage parameters.

Thank you Professor.

What do you mean by “sample size requirement”?

I have 122 observations and an ordinal discrete response variable, with 36 distinct values. The `lrm()`

model has 10 df with 9 predictors.

When the sample size is small enough to need penalization it’s often too small to be able to guide you in selection of how much shrinkage to apply.

Interesting, and rather paradoxical.

I’ve come across this great paper of yours.

- Aren’t you also comparing too many models/strategies here?

I tried to replicate the approach of finding an optimal penalty factor through bootstrap, but couldn’t do it using the `rms`

functions*, so had to do it another way (using a `for`

loop and the `dplyr::sample_n`

function. The problem is that I need to find a way to ensure balanced resampling, so that there is at least one observation for each distinct value of Y in each bootstrap sample (something which is achieved with the argument `group`

when using `rms::validate`

or `rms::calibrate`

).

Anyway, I decided to try a different strategy (see plot attached):

To search for an optimum penalty factor by comparing the optimism-corrected (via bootstrap with `rms::validate`

) calibration slope yielded by different penalty factors, and choosing the penalty yielding the calibration slope that was closest to 1.

- Does that make sense?

*I don’t think there is a function to search for the penalty factor using bootstrap, correct?

Thank you again for your attention and invaluable commments.

Hi,

this post comes from this other one that was posted in the wrong place:

Professor Harrell’s answer:

You may have to combine infrequent levels. `Hmisc::combine.levels`

can help with that. Or use the `group=`

argument to `calibrate`

and `validate`

to do balanced bootstrap sampling. Specify something like `group=mydataset$FLUID2`

.

→ Thank you for your quick reply. I have another question with regards the use of `group=`

argument. Is it possible to use it for grouping more than one covariate at the same time?

I mean to do something similar to ‘group=c(mydataset$FLUID2,mydataset$TYPE)’.

note: I tested this way to use a vector to specify more than one covariate but it didn’t work.

Thank you

Right there is no built-in function for bootstrapping with penalization-finding.

You can use the `interact()`

function (or is it `interaction`

?) to paste together multiple factors to use in `group`

.

Thanks.

What about the strategy illustrated on the plot above? Any thoughts?

It’s useful but would be nice to use a comprehensive measure that rewarded both calibration and discrimination, such as Brier score, mean squared prediction error, mean absolute prediction error.

Hi guys,

I am little bit crazy trying to run a penalized cox ph regression model with functions of rms package. I’ve done it with glmnet package but I would like to analize the outputs with `rms`

and `Hmisc`

functions later.

Is it possible to do so with `cph()`

function? as well, how can be performed with cross-validation to found the optimum lambda which leads e.g. to a minimum mean cross-validated error.

Any help or do you have an example to follow?

Thank you

Hi professor Harrell,

I am now checking the cox Lasso regression on my full-original model. After optimizing the value of lambda with a cross-validation procedure, I get a penalized model with its coefficients (I used `penalized()`

package).

`fit <- penalized(formula_model, data = g_glmnet, lambda1=optL1fit$lambda, epsilon=1e-5, model = "cox")`

After doing that, I tried to run a new model with `cph()`

function, in order to get all the parameters in the format in which I have all the other models and to compare the statistics.

I want to force to have the same model with `cph()`

but I get weird results, as, for example LR X2 = 0.

This is how I am trying to do so:

`cph_lasso_fit <- cph(formula_model, data = g_glmnet, init=coefficients(fit,"all"),iter=0, tol=1e-13, x=TRUE,y=TRUE, surv=TRUE, ties = 'efron', id = id)`

I understand that I am doing something wrong or maybe it’s not possible what I am trying to do. Can you help me somehow in this issue?

Thank you so much in advance for your help

You can’t feed results of a penalized analysis into an unpenalized fit.