Régression logistique avec splines de régression dans R

12

J'ai développé un modèle de régression logistique basé sur des données rétrospectives d'une base de données nationale sur les traumatismes des traumatismes crâniens au Royaume-Uni. Le résultat clé est la mortalité à 30 jours (dénommée mesure de «survie»). D'autres mesures avec des preuves publiées d'effet significatif sur les résultats dans les études précédentes comprennent:

Year - Year of procedure = 1994-2013
Age - Age of patient = 16.0-101.5
ISS - Injury Severity Score = 0-75
Sex - Gender of patient = Male or Female
inctoCran - Time from head injury to craniotomy in minutes = 0-2880 (After 2880 minutes is defined as a separate diagnosis)

En utilisant ces modèles, étant donné la variable dépendante dichotomique, j'ai construit une régression logistique en utilisant lrm.

La méthode de sélection des variables du modèle était basée sur la littérature clinique existante modélisant le même diagnostic. Tous ont été modélisés avec un ajustement linéaire à l'exception de l'ISS qui a été modélisé traditionnellement par des polynômes fractionnaires. Aucune publication n'a identifié d'interactions significatives connues entre les variables ci-dessus.

Sur les conseils de Frank Harrell, j'ai procédé à l'utilisation de splines de régression pour modéliser l'ISS (cette approche présente des avantages dans les commentaires ci-dessous). Le modèle a donc été prédéfini comme suit:

rcs.ASDH<-lrm(formula = Survive ~ Age + GCS + rcs(ISS) +
    Year + inctoCran + oth, data = ASDH_Paper1.1, x=TRUE, y=TRUE)

Les résultats du modèle étaient:

> rcs.ASDH

Logistic Regression Model

lrm(formula = Survive ~ Age + GCS + rcs(ISS) + Year + inctoCran + 
    oth, data = ASDH_Paper1.1, x = TRUE, y = TRUE)

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
Obs          2135    LR chi2     342.48    R2       0.211    C       0.743    
 0            629    d.f.             8    g        1.195    Dxy     0.486    
 1           1506    Pr(> chi2) <0.0001    gr       3.303    gamma   0.487    
max |deriv| 5e-05                          gp       0.202    tau-a   0.202    
                                           Brier    0.176                     

          Coef     S.E.    Wald Z Pr(>|Z|)
Intercept -62.1040 18.8611 -3.29  0.0010  
Age        -0.0266  0.0030 -8.83  <0.0001 
GCS         0.1423  0.0135 10.56  <0.0001 
ISS        -0.2125  0.0393 -5.40  <0.0001 
ISS'        0.3706  0.1948  1.90  0.0572  
ISS''      -0.9544  0.7409 -1.29  0.1976  
Year        0.0339  0.0094  3.60  0.0003  
inctoCran   0.0003  0.0001  2.78  0.0054  
oth=1       0.3577  0.2009  1.78  0.0750  

J'ai ensuite utilisé la fonction calibrer dans le package rms afin d'évaluer la précision des prédictions du modèle. Les résultats suivants ont été obtenus:

plot(calibrate(rcs.ASDH, B=1000), main="rcs.ASDH")

Courbes d'étalonnage Bootstrap pénalisées pour sur-ajustement

Après avoir terminé la conception du modèle, j'ai créé le graphique suivant pour démontrer l'effet de l'année d'incident sur la survie, en basant les valeurs de la médiane sur des variables continues et le mode sur des variables catégorielles:

ASDH <- Predict(rcs.ASDH, Year=seq(1994,2013,by=1),Age=48.7,ISS=25,inctoCran=356,Other=0,GCS=8,Sex="Male",neuroYN=1,neuroFirst=1)
Probabilities <- data.frame(cbind(ASDH$yhat,exp(ASDH$yhat)/(1+exp(ASDH$yhat)),exp(ASDH$lower)/(1+exp(ASDH$lower)),exp(ASDH$upper)/(1+exp(ASDH$upper))))
names(Probabilities) <- c("yhat","p.yhat","p.lower","p.upper")
ASDH<-merge(ASDH,Probabilities,by="yhat")
plot(ASDH$Year,ASDH$p.yhat,xlab="Year",ylab="Probability of Survival",main="30 Day Outcome Following Craniotomy for Acute SDH by Year", ylim=range(c(ASDH$p.lower,ASDH$p.upper)),pch=19)
arrows(ASDH$Year,ASDH$p.lower,ASDH$Year,ASDH$p.upper,length=0.05,angle=90,code=3)

Le code ci-dessus a généré la sortie suivante:

Tendance de l'année avec inférieur et supérieur

Mes questions restantes sont les suivantes:

1. Interprétation des splines - Comment puis-je calculer la valeur de p pour les splines combinées pour la variable globale?

Fontaine Dan
la source
4
Bon travail. Pour afficher l'effet de l' année , je vous suggère de laisser les autres variables sont définies à des valeurs par défaut (médiane en mode continu, pour catégoriel) et varient d'une année sur l'axe x, par exemple plot(Predict(rcs.ASDH, Year)). Vous pouvez laisser d'autres variables varier, en créant différentes courbes, en faisant des choses comme plot(Predict(rcs.ASDH, Year, age=c(25, 35))).
Frank Harrell
1
Je ne sais pas pourquoi - mais je n'ai pas vu beaucoup d'exemples de courbes d'étalonnage à correction de biais dans la littérature. Semble être une bonne idée
Charles
1
Pour tester l'association globale avec plusieurs tests df, utilisez anova(rcs.ASDH).
Frank Harrell

Réponses:

8

Il est très difficile d'interpréter vos résultats lorsque vous n'avez pas prédéfini le modèle mais que vous vous êtes engagé dans des tests importants et des modifications de modèle en résultant. Et je ne recommande pas le test d'ajustement global que vous avez utilisé car il a une distribution dégénérée et non une distribution .χ2

Deux méthodes recommandées pour évaluer l'ajustement d'un modèle sont les suivantes:

  1. Courbe d'étalonnage lisse non paramétrique corrigée par le sur-ajustement Bootstrap (par exemple, * loess) pour vérifier la précision absolue des prédictions
  2. Examinez diverses généralisations du modèle, testant si une spécification de modèle plus flexible fonctionne mieux. Par exemple, effectuez un rapport de vraisemblance ou des tests de Wald sur des termes supplémentaires non linéaires ou d'interaction.χ2

Les splines de régression présentent certains avantages par rapport aux polynômes fractionnaires, notamment:

  1. Le prédicteur peut avoir des valeurs ; vous ne pouvez pas prendre de journaux de valeurs telles que requises par les FP0
  2. Vous n'avez pas à vous soucier de l'origine d'un prédicteur. Les FP supposent que zéro est une origine "magique" pour les prédicteurs, tandis que les splines de régression sont invariantes au décalage d'un prédicteur par une constante.

Pour plus d'informations sur les splines de régression et l'évaluation de la linéarité et de l'additivité, voir mes documents à http://biostat.mc.vanderbilt.edu/CourseBios330 ainsi que la fonction de rmspackage R. rcsPour les courbes d'étalonnage bootstrap pénalisées en cas de sur-ajustement, voir la rms calibratefonction.

Frank Harrell
la source
J'ai essayé de prédéfinir le modèle en créant un modèle complet avec toutes les variables cliniquement disponibles, y compris initialement et avec une pertinence connue pour le diagnostic et la variable de résultat. Il s'agissait de variables linéaires continues ou dichotomiques, à l'exception de l'ISS qui, selon des études antérieures, peut être modélisé avec des polynômes fractionnaires. La méthode que j'ai utilisée pour développer le modèle est, je crois, cohérente avec le "Multivariate Model-Building" de Willi Sauerbrei. Peut-être pourrais-je utiliser votre package rms en R afin d'évaluer la qualité globale de l'ajustement? Si oui, quelle formule recommanderiez-vous?
Dan Fountain
J'ai développé ma réponse pour répondre à certaines de ces questions.
Frank Harrell
Pourriez-vous recommander des packages pour l'exécution de 1 et 2 pour évaluer l'ajustement du modèle? Puis-je utiliser la fonction de validation dans rms? En termes de splines de régression, j'ai lu vos notes de cours et j'essaie actuellement de trouver une copie de votre livre (je l'achèterais si je le pouvais!) Seriez-vous en mesure de recommander un autre package / fonction R pour la construction de splines de régression pour les variables continues non linéaires dans ce modèle? Merci beaucoup pour toute votre aide jusqu'à présent.
Dan Fountain
Peut-être le paquet Terre basé sur les splines de régression adaptative multivariée de Friedman?
Dan Fountain
Une dernière question. La formule du score de gravité des blessures (ISS) est A ^ 2 + B ^ 2 + C ^ 2 où A, B et C sont différentes parties du corps avec un score de gravité de traumatisme indépendant de 1 à 5. Ainsi, le maximum est de 75 et le minimum de 1 dans cet ensemble de données. Compte tenu de cette formule, les polynômes fractionnaires seraient-ils la représentation la plus proche de la façon dont le score est réellement calculé par rapport aux splines de régression?
Dan Fountain