Regresja logistyczna z splajnami regresji w R.

12

Opracowuję model regresji logistycznej oparty na danych retrospektywnych z krajowej bazy danych dotyczących urazów głowy w Wielkiej Brytanii. Kluczowym rezultatem jest 30-dniowa śmiertelność (oznaczona jako miara „przetrwania”). Inne miary z opublikowanymi dowodami znaczącego wpływu na wyniki poprzednich badań obejmują:

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)

Używając tych modeli, biorąc pod uwagę dychotomiczną zmienną zależną, zbudowałem regresję logistyczną za pomocą lrm.

Metodę selekcji zmiennych modelowych oparto na istniejącej literaturze klinicznej modelującej tę samą diagnozę. Wszystkie zostały modelowane z dopasowaniem liniowym, z wyjątkiem ISS, który został tradycyjnie modelowany za pomocą wielomianów ułamkowych. Żadna publikacja nie zidentyfikowała znanych znaczących interakcji między powyższymi zmiennymi.

Zgodnie z radą Franka Harrella zastosowałem splajny regresji do modelu ISS (zalety tego podejścia podkreślono w komentarzach poniżej). Model został więc wstępnie określony w następujący sposób:

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

Wyniki modelu to:

> 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  

Następnie użyłem funkcji kalibracji w pakiecie rms w celu oceny dokładności prognoz z modelu. Następujące wyniki zostały osiągnięte:

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

Krzywe kalibracji bootstrap karane za przeregulowanie

Po zakończeniu projektowania modelu stworzyłem następujący wykres, aby zademonstrować wpływ roku zdarzenia na przeżycie, opierając wartości mediany w zmiennych ciągłych i tryb w zmiennych kategorialnych:

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)

Powyższy kod zaowocował następującymi danymi wyjściowymi:

Trend roku z dolną i górną

Moje pozostałe pytania to:

1. Interpretacja splajnu - Jak obliczyć wartość p dla splajnów połączonych dla zmiennej ogólnej?

Fontanna Dan
źródło
4
Dobra robota. Dla wyświetlenia efektu Roku sugeruję, aby pozwolić innym ustawiać wartości domyślne (mediana dla ciągłego, tryb dla jakościowego) i zmieniać Rok na osi x, np plot(Predict(rcs.ASDH, Year)). Możesz zmieniać inne zmienne, tworząc różne krzywe, wykonując takie czynności plot(Predict(rcs.ASDH, Year, age=c(25, 35))).
Frank Harrell,
1
Nie wiem dlaczego - ale w literaturze nie widziałem wielu przykładów krzywych kalibracyjnych skorygowanych o błąd systematyczny. Wydaje się, że to dobry pomysł
Charles
1
Aby przetestować ogólne powiązanie z wieloma testami DF, użyj anova(rcs.ASDH).
Frank Harrell,

Odpowiedzi:

8

Bardzo trudno jest zinterpretować wyniki, jeśli nie określono wcześniej modelu, ale przeprowadzono znaczące testy i wynikające z niego modyfikacje modelu. I nie polecam zastosowanego globalnego testu dopasowania, ponieważ ma on zdegenerowany rozkład, a nie .χ2

Dwa zalecane sposoby oceny dopasowania modelu to:

  1. Łagodna, nieparametryczna (np. * Lessowa) krzywa kalibracji bootstrapu korygująca przepełnienie w celu sprawdzenia absolutnej dokładności prognoz
  2. Zbadaj różne uogólnienia modelu, sprawdzając, czy bardziej elastyczna specyfikacja modelu działa lepiej. Na przykład wykonaj testy prawdopodobieństwa lub testy Wald dodatkowych warunków nieliniowych lub interakcji.χ2

Istnieją pewne zalety splajnów regresji w porównaniu z wielomianami ułamkowymi, w tym:

  1. Predyktor może mieć wartości ; nie możesz brać dzienników takich wartości, jakie są wymagane przez FP0
  2. Nie musisz się martwić o pochodzenie predyktora. FP zakładają, że zero jest „magicznym” źródłem dla predyktorów, podczas gdy splajny regresji są niezmienne dla przesunięcia predyktora o stałą.

Więcej informacji na temat wypustami regresji i liniowość i oceny addytywności zobaczyć moje materiały informacyjne na http://biostat.mc.vanderbilt.edu/CourseBios330 jak również R rmspakiet rcsfunkcji. Aby zapoznać się z krzywymi kalibracji bootstrap karanymi za przeregulowanie, patrz rms calibratefunkcja.

Frank Harrell
źródło
Podjęłam próbę wstępnego określenia modelu przy tworzeniu pełnego modelu ze wszystkimi zmiennymi dostępnymi klinicznie, w tym początkowo o znanym znaczeniu dla zmiennej diagnozy i wyniku. Wszystkie były liniowymi zmiennymi ciągłymi lub dychotomicznymi, z wyjątkiem ISS, które zidentyfikowane w poprzednich badaniach można modelować za pomocą ułamkowych wielomianów. Metoda, którą zastosowałem do opracowania modelu, który moim zdaniem jest zgodny z „Wielowymiarowym budowaniem modelu” Willi Sauerbrei. Być może mógłbym użyć twojego pakietu rms w R, aby ocenić ogólną poprawność dopasowania? Jeśli tak, jaką formułę poleciłbyś?
Dan Fountain
Rozszerzyłem swoją odpowiedź, aby rozwiązać niektóre z tych kwestii.
Frank Harrell,
Czy możesz polecić pakiety do wykonania 1 i 2, aby ocenić dopasowanie modelu? Czy mogę użyć funkcji sprawdzania poprawności w RMS? Jeśli chodzi o splajny regresji, czytam notatki z wykładów i obecnie próbuję znaleźć kopię twojej książki (kupiłbym ją, gdybym mógł!) Czy byłbyś w stanie polecić kolejny pakiet / funkcję R do budowy splajnów regresji dla nieliniowych zmiennych ciągłych w tym modelu? Wielkie dzięki za dotychczasową pomoc.
Dan Fountain
Być może pakiet ziemski oparty na wielowymiarowych wielowypustowych regresjach adaptacyjnych Friedmana?
Dan Fountain
Ostatnie pytanie. Wzór na wskaźnik ciężkości urazu (ISS) to A ^ 2 + B ^ 2 + C ^ 2, gdzie A, B i C to różne części ciała z niezależnym wskaźnikiem nasilenia urazu 1-5. Zatem maksimum wynosi 75, a minimum 1 w tym zbiorze danych. Biorąc pod uwagę ten wzór, czy ułamkowe wielomiany byłyby dokładniejszym odwzorowaniem sposobu obliczania wyniku w porównaniu do splajnów regresji?
Fontanna Dan