Download Report

Kursus i basal statistik:
Ventetidsanalyse
(“Overlevelsesanalyse”)
Birthe Lykke Thomsen
Det Nationale Forskningscenter for Arbejdsmiljø
Ventetidsanalyse
Baggrund
Åreknuder i spiserøret er en alvorlig følgesygdom, der kan optræde
hos patienter med levercirrose. Hvis en af åreknuderne brister, er
patienten i alvorlig risiko for at dø af den indvendige blødning.
Forskningsspørgsmål: Vi ønsker at undersøge, om
skleroterapi kan bruges til at forebygge/udskyde
forekomsten af reblødning hos patienter med
blødende åreknuder i spiserøret p.g.a. levercirrose.
2
Ventetidsanalyse
Baggrund
Data:
Randomiseret studie af patienter, der har blødt
én gang fra åreknuder i spiserøret. Respons er,
hvornår patienten rebløder.
(EVASP, 1984)
3
Ventetidsanalyse
Baggrund
Særlige egenskaber ved denne slags data
• Respons: tid til en given hændelse – her reblødning; men det
kunne også være død, sygdom, helbredelse, tilbagefald,
graviditet, . . .
– gør det særligt interessant at lave simple modeller for
sandsynligheden for, at noget sker på et givet tidspunkt
forudsat, at personen på det tidspunkt er under risiko for, at
det sker
– disse betingede sandsynligheder kaldes rater eller hazards
– typisk kigger man efter multiplikative eﬀekter på raten – det
kaldes Cox’s proportional hazards model
– hvis (og kun hvis) raten er konstant, kan den estimeres ved
antal hændelser/samlet risikotid
(standard formel fra epidemiologi), og eﬀekter kan så
estimeres ved Poisson regression
4
Ventetidsanalyse
Baggrund
Særlige egenskaber ved denne slags data fortsat
• Forsinket indgang og censurering: Vi har kun en større eller
mindre bid af hvert individs liv med i vores studie
– forsinket indgang: nogle personer er ikke med fra start
– censurering: for nogle af individerne ved vi kun, at det
endnu ikke var sket, da de gik ud af studiet
• og vi må kun bruge den tid, hvor vi ville tælle
hændelsen med, når/hvis den skete
5
Ventetidsanalyse
Baggrund
Særlige egenskaber ved denne slags data fortsat
• Oftest ingen anelse om fordelingen af tid til hændelsen
– den er i praksis stort set aldrig normalt fordelt
(for censurerede, normalt fordelte data, f.eks. outcomes med
detektionsgrænse, bruges PROC LIFEREG; se s. 65–83 i disse noter til
selvstudium)
– det kan nogle gange være rimeligt at tro, at raten er konstant
– en forudsætning for mange epidemiologiske vurderinger
• men det er altid en forudsætning, at censurering er uafhængig af
outcome givet de kendte, forklarende variable – censurering må
ikke være prædiktivt for, om hændelsen snart vil ske
6
Ventetidsanalyse
Baggrund
Særlige metoder kan håndtere denne slags data
Fordi raten er i fokus, ser man på hvert tidspunkt, et ad gangen, og
opsummerer derefter i en
• Kaplan-Meier kurve
• kumuleret (summeret) rate
• Cox regression
Når man bruger disse metoder, så er estimaterne for raten
og Cox regressionen korrekte, uanset om der er forsinket
indgang og censureringer eller ej! Forudsat, selvfølgelig, at
censureringen er uafhængig af outcome...
Når der er forsinket indgang, og hændelsen sker for alle, der under risiko på
det givne (tidlige) tidspunkt, bliver Kaplan-Meier estimatet lig 0, og
dermed forkert for resten af tidsperioden (forklaring følger).
7
Ventetidsanalyse
Kaplan-Meier
Reblødning i de to grupper: Kaplan-Meier og kumuleret rate
uden henholdsvis med punktvise konfidensgrænser
8
Ventetidsanalyse
Kaplan-Meier
Fortolkning af Kaplan-Meier og kumuleret rate
Hvis hændelsen er “død af hvilken som helst årsag”, så viser
Kaplan-Meier kurven overlevelsessandsynligheden som funktion af
tiden t, dvs. hvor stor en andel, man kan regne med vil overleve
mindst til tid t.
Hvis hændelsen er hvad som helst andet, og død har medført
censurering, så giver Kaplan-Meier kurven ingen fornuftig mening!
Den kumulerede rate har ingen direkte, intuitiv fortolkning. Den er
ikke lig sandsynligheden for, at hændelsen sker inden det givne
tidspunkt – men det er en tilnærmelse, når den kumulerede rate er
lille. Hældningen på den kumulerede rate er lig med selve raten.
Rater og kumulerede rater er korrekte, også når død medfører
censurering!
9
Ventetidsanalyse
Kaplan-Meier
10
Reblødning i de to grupper – standard plots
Trappekurve: et trin ned hver gang,
der sker en hændelse.
S(T ) ≈ exp(−R(T ))
R(T ) ≈ − ln(S(T ))
(Stykkevis) konstant rate giver
(stykkevis) lineær kumuleret rate R(t).
For tidsintervaller med lineær R(t) kan
selve raten estimeres som forskellen
mellem startværdi og slutværdi for
R(t) divideret med længden af
tidsintervallet.
Høj rate=stejl stigning,
lav rate=flad kurve.
Ventetidsanalyse
Kaplan-Meier
Beregning af Kaplan-Meier og kumuleret rate
(alt under bølgelinien er skjult for os, og kun de tynde, grønne personer er med fra start)
11
Ventetidsanalyse
Kaplan-Meier
Beregning af Kaplan-Meier og kumuleret rate
På en given dag t observerer vi for hver gruppe g (kaldes
sædvanligvis “stratum”, flertal “strata”, i ventetidsanalyse)
1. ng (t) individer totalt
2. mg (t) individer, der begynder at rebløde
hvilket giver den daglige reblødningsrate
mg (t)
rg (t) =
ng (t)
Kaplan-Meier estimatet Sg (T ) på dag T for gruppe g beregnes ved
at gange leddene 1 − rg (t) sammen for alle dage t op til og med dag T .
Nelson-Aalen estimatet for den kumulerede reblødningsrate Rg (T )
på dag T for gruppe g beregnes ved at lægge de daglige
reblødningsrater rg (t) sammen for alle dage t op til og med dag T .
12
Ventetidsanalyse
Kaplan-Meier
13
Estimation af Kaplan-Meier kurver v.h.a. PHREG
Respons:
“Det tidspunkt, hændelsen sker”;
men det sker ikke for alle, så det er
nødvendigt at bruge en kombination
af 2 responsvariable:
“tid” og “hvad skete der”.
Datasættet SKL indeholder (bl.a.)
DAY: tid for udgang
BLD:
1 for reblødning,
0 for censurering
SKLERO:
1 for skleroterapigruppen,
0 for kontrolgruppen
NB: Alle værdierne for censurering skal skrives
i parentesen efter status-variablen! (her BLD)
ODS GRAPHICS ON;
PROC PHREG DATA=skl PLOTS(OVERLAY=ROW CL)=S;
MODEL Day*Bld (0) = ; /*NB ingen X-var*/
STRATA Sklero;
/*Hvis ODS GRAPHICS ikke duer, så lav datasæt:*/
BASELINE OUT=km SURVIVAL=KMcurves
LOWER=LowerBound UPPER=UpperBound
/ METHOD=PL CLTYPE=LOGLOG;
RUN;
SYMBOL1 C=BLACK V=NONE I=STEPLJ MODE=INCLUDE;
SYMBOL2 C=RED
V=NONE I=STEPLJ MODE=INCLUDE;
PROC GPLOT DATA=km;
PLOT kmcurves*day=sklero;
RUN;
Ventetidsanalyse
Kaplan-Meier
SAS’s ODS GRAPHICS figur
14
Ventetidsanalyse
Kaplan-Meier
15
Ekstra SAS kode for at få konfidensgrænser med i egen figur
SAS kode for nødvendige
datamodifikationer:
DATA KM_m_CL; SET km;
IF sklero=0 THEN DO;
type=1; curve=kmcurves; OUTPUT;
type=2; curve=lowerbound; OUTPUT;
type=3; curve=upperbound; OUTPUT;
END;
IF sklero=1 THEN DO;
type=4; curve=kmcurves; OUTPUT;
type=5; curve=lowerbound; OUTPUT;
type=6; curve=upperbound; OUTPUT;
END;
RUN;
Eksempel på kald af GPLOT:
AXIS1 LABEL=(H=2) VALUE=(H=1.75);
AXIS2 LABEL=(H=2 A=90 "S(t) = Kaplan-Meier")
VALUE=(H=1.75);
SYMBOL1
SYMBOL2
SYMBOL3
SYMBOL4
R=1
R=2
R=1
R=2
C=BLACK
C=BLACK
C=RED
C=RED
V=NONE
V=NONE
V=NONE
V=NONE
I=STEPLJ
I=STEPLJ
I=STEPLJ
I=STEPLJ
L=2 W=5 MODE=INCLUDE;
L=4 MODE=INCLUDE;
L=1 W=5 MODE=INCLUDE;
L=41 MODE=INCLUDE;
PROC GPLOT DATA=KM_m_CL;
PLOT curve*day=type
/ HAXIS=AXIS1 VAXIS=AXIS2 NOLEGEND;
RUN;
(SAS-kode for plot af kumuleret rate kan findes bagest i disse noter)
Ventetidsanalyse
Log rank
16
Log rank test (svarer til en-sidet variansanalyse)
Grupper kan sammenlignes v.h.a. et log rank test.
Princip (for 2 grupper):
Vi forudsætter (“nul-hypotesen”), at der ingen forskel er på de 2
grupper, og betinger for hvert ’døds’tidspunkt med
• det observerede totale antal ’døde’ m(ti ) (= m1 (ti ) + m2 (ti ))
• antallet af personer under risiko på det givne tidspunkt i hver af
de 2 grupper, n1 (ti ) og n2 (ti ) (n1 (ti ) + n2 (ti ) = n(ti ))
For gruppe 1 kan vi så for hvert ’døds’tidspunkt ti beregne
• det forventede antal ’døde’ E1 (ti ) = m(ti ) ·
• variansen på antallet af ’døde’
n2 (ti ) m(ti )·(n(ti )−m(ti ))
1 (ti )
V (ti ) = nn(t
·
n(ti ) ·
(n(ti )−1)
i)
n1 (ti )
n(ti )
Ventetidsanalyse
Log rank
Log rank test fortsat
De forventede antal ’døde’ E1 (ti ), henholdsvis varianserne V (ti ), adderes
for alle ’døds’tidspunkter ti til E1 , henholdsvis V . Desuden tælles det
totale antal observerede ’døde’ M1 i gruppe 1.
Log rank teststørelsen
(M1 − E1 )2
Log Rank Chi Square =
,
V
er χ2 -fordelt med 1 frihedsgrad.
Tilnærmelse, der også kan bruges for mere end 2 grupper, her G grupper:
G
∑
(Mg − Eg )2
Log Rank Chi Square ≈
,
E
g
g=1
er χ2 -fordelt med G − 1 frihedsgrader (bemærk, at alle grupper bidrager
til summen).
17
Ventetidsanalyse
Styrkeberegning
Styrkeberegning for ventetidsdata
Summen af varianserne, V , brugt ved beregningen af log rank
teststørrelsen, bruges også ved styrkeberegninger:
Hvis man ønsker at kunne påvise en rate ratio på RR for aktiv versus
kontrol med et tosidet signifikansniveau på α og en styrke på 1−β, så skal
( U +U )2
β
α/2
man inkludere tilstrækkeligt mange personer til, at V ≥
.
ln(RR)
(Ux =(1−x)-fraktilen i en normal fordeling, eksempelvis U2.5% =1.96,
U5% =1.645, U10% =1.282, U20% =0.842)
Hvis randomiseringsratioen (aktiv:kontrol) er a:1, og hændelsen er sjælden,
a
er V ≈ (a+1)
2 · M , hvor M =total antal hændelser, dvs. der skal inkluderes
tilstrækkeligt mange personer til, at
(
)
Uα/2 + Uβ 2
(a + 1)2
·
M≥
a
ln(RR)
18
Ventetidsanalyse
Log rank
Log rank test ved hjælp af PROC PHREG
Log rank teststørrelsen kan beregnes som score-teststørrelsen i
PROC PHREG med brug af option TIES=DISCRETE:
PROC PHREG DATA=skl;
MODEL day*bld(0) = sklero / TIES=DISCRETE;
RUN;
Hvis der er mere end 2 grupper, er det nødvendigt at bruge et CLASS statement,
eksempelvis
PROC PHREG DATA=skl;
CLASS sklero / PARAM=GLM;
MODEL day*bld(0) = sklero / TIES=DISCRETE;
RUN;
19
Ventetidsanalyse
Log rank
20
Dele af output fra PROC PHREG
The PHREG Procedure
Model Information
Data Set
WORK.SKL
Dependent Variable
Day
Censoring Variable
Bld
Censoring Value(s)
0
Ties Handling
DISCRETE
Summary of the Number of Event and Censored Values
Percent
Total
Event
Censored
Censored
187
91
96
51.34
Model Fit Statistics
Without
With
Criterion
Covariates
Covariates
-2 LOG L
738.406
737.488
:
:
:
Testing Global Null Hypothesis: BETA=0
Test
Chi-Square DF
Pr > ChiSq
Likelihood Ratio
0.9175
1
0.3381
Score
0.9174
1
0.3382 ← Log rank teststørrelse=0.9174, p=0.3382
Wald
0.9124
1
0.3395
Ventetidsanalyse
Log rank
Output fra PROC PHREG fortsat
Analysis of Maximum Likelihood Estimates
Parameter Standard
Hazard
Variable DF Estimate
Error Chi-Square Pr>ChiSq Ratio
Sklero
1 -0.20261
0.21212
0.9124
0.3395 0.817
Log rank testet er dårligt til at detektere forskelle, der
ændres markant med tiden, som eksempelvis en dårligere
korttidsprognose samtidigt med en bedre langtidsprognose!
21
Ventetidsanalyse
Opsummering 1. del
Opsummering 1. del
Ventetidsdata karakteristika:
• “et tidspunkt ad gangen”: vi laver modeller for raten
• censurering og forsinket indgang
• ingen fordelingsantagelser
Typisk grafik:
• Kaplan-Meier (pas på!)
• kumuleret rate – tjek for stykkevis linearitet
Numerisk test (’p-værdi’) for sammenligning af grupper (det,
der svarer til ensidet variansanalyse):
• log rank test – dårligt til at detektere tidsvarierende eﬀekter
22
Ventetidsanalyse
Cox model
23
Proportionale rater
Kvantificering af behandlingseﬀekten:
r(t; sklero) = r(t; medicinsk) · B
Eﬀekt af ascites (væske i bughulen):
r(t; ascites) = r(t; ingen ascites) · A
Kombineret:
r(t; sklero, ascites) = r(t; medicinsk, ascites) · B
= r(t; medicinsk, ingen ascites) · A · B
= r(t; medicinsk, ingen ascites) · exp(a + b)
med a=ln(A) og b=ln(B).
0 ∼
∼
Sættes X1 = {1
ingen ascites
0
og
X
=
{
2
ascites
1
∼
∼
medicinsk
skleroterapi
så er
r(t; sklero, ascites) = r0 (t) · exp(aX1 + bX2 )
Ventetidsanalyse
Cox model
Cox’s regression model
Denne model kaldes Cox’s regression model:
r(t; X1 , X2 , . . . , Xk ) = r0 (t) · exp(b1 X1 + b2 X2 + . . . + bk Xk )
Summen i eksponenten, b1 X1 + b2 X2 + . . . + bk Xk , kaldes nogle
gange det “prognostiske indeks”.
Hvis vi transformerer med den naturlige logaritme og skriver b0 (t) i
stedet for ln(r0 (t)), fås noget, der ligner de sædvanlige
regressionsmodeller:
ln(r(t; X1 , X2 , . . . , Xk )) = b0 (t) + b1 X1 + b2 X2 + . . . + bk Xk
Bemærk, at logaritmen til den underliggende rate r0 (t) svarer til den
sædvanlige intercept parameter (der altså afhænger af tiden her)!
24
Ventetidsanalyse
Cox model
En positiv værdi for en regressionsparameter bj betyder, at høje
værdier af den tilsvarende kovariat Xj hører sammen med høj rate:
For uønskede hændelser betyder det, at høje værdier forværrer
prognosen, så vær meget forsigtig med at bruge ordene
“positiv/negativ” om eﬀekter, brug eksempelvis “gavnlig/skadelig”.
25
Ventetidsanalyse
Cox model
26
PROC PHREG DATA=skl;
MODEL day*bld(0) = ascites bilirub sklero / RL;
RUN;
--------------------------------------------------------Summary of the Number of Event and Censored Values
Total
177
Event
87
Censored
90
Percent
Censored
50.85
Variable DF
Ascites
1
Bilirub
1
Sklero
1
Parameter Standard
Estimate
Error Chi-Square Pr>ChiSq
0.18072 0.22721
0.6326
0.4264
0.00476 0.00112
18.1500
<.0001
-0.21924 0.21801
1.0113
0.3146
Variable
Ascites
Bilirub
Sklero
Hazard
Ratio
1.198
1.005
0.803
95% Hazard
Confidence
0.768
1.003
0.524
Ratio
Limits
1.870
1.007
1.231
Ventetidsanalyse
Specielle aspekter ved Cox modellen
Forklarende variable i ventetidsanalyser
En variabel kan indgå i Cox regressionsanalyser på to essentielt
forskellige måder:
1. som en regressionsvariabel, dvs. i det prognostiske index, hvad
enten det er som
(a) en lineær variabel, eller som
(b) en kategorisk (CLASS) variabel
2. som en stratifikationsvariabel for den underliggende rate
De to muligheder, regressionsvariabel vs. stratifikationsvariabel, er
helt forskellige, og valget har ikke alene konsekvenser for både
modellering og fortolkning af den aktuelle variabel, men påvirker
også estimationen af eﬀekten af de øvrige variable!
27
Ventetidsanalyse
Specielle aspekter ved Cox modellen
Estimation i en Cox model
For den dag tj , hvor patient j reblødte, beregner vi sandsynligheden
for, at det var netop patient j, der reblødte, betinget med, at der
“skulle” ske en reblødning den dag i det stratum, som patient j
tilhører:
exp(b · Xj (tj ))
∑
.
i med i ℜj exp(b · Xi (tj ))
Her betegner ℜj alle de patienter (i’erne), som var under risiko for
reblødning i det samme stratum som j på netop den dag tj ,
hvor j reblødte.
Disse bidrag ganges sammen for alle reblødningstidspunkter, og b
estimeres som den værdi, b̂, som maksimerer dette totale produkt,
der kaldes “Cox’s partielle likelihood”.
28
Ventetidsanalyse
Specielle aspekter ved Cox modellen
29
Forklarende variable i ventetidsanalyser: Regressionsvariable
Regressionsvariable (både lineære og kategoriske variable) indgår
udelukkende via det prognostiske index, dvs. i exp(bX), og ikke i
den underliggende rate r0 (t). Med andre ord, raterne antages at være
proportionale for forskellige værdier af variablen, X:
r(t; X = 0) = r0 (t)
og
r(t; X = x + 1) = r0 (t) · exp( b·(x + 1) ) = r(t; X = x) · exp(b)
eller
r(t; X = referencegruppe) = r0 (t)
og
r(t; X = sammenligningsgruppe) = r0 (t) · exp(b)
Konsekvens:
1. eﬀekten kan beskrive med et enkelt tal (HR=exp(b)),
2. men dette tal er kun meningsfuldt, hvis antagelsen om
proportionale rater holder (tilnærmelsesvis).
Ventetidsanalyse
Specielle aspekter ved Cox modellen
30
Forklarende variable i ventetidsanalyser: Stratifikationsvariable
For stratifikationsvariable afhænger den underliggende rate af
hvilken værdi variablen har, dvs. forskellen på individer med X = 1
og individer med X = 2 ændres med tiden:
r(t; X = 1) = r1 (t)
og
r(t; X = 2) = r2 (t)
Konsekvens:
1. vi kan ikke beskrive eﬀekten med et enkelt tal,
2. stratifikationsvariable er nødt til at være kategoriske variable
med et begrænset antal forskellige værdier.
Ventetidsanalyse
Specielle aspekter ved Cox modellen
Stratifikation versus Interaktion
Stratifikation er ikke det samme som interaktion! Eﬀekten af
stratifikationsvariablen ændres med tiden, men eﬀekten afhænger ikke
af de øvrige variable, og deres eﬀekter antages at være identiske i de
forskellige strata — i modsætning til den epidemiologiske brug af
udtrykket “stratificerede analyser”, som betyder helt separate
analyser for hver værdi af stratifikationsvariablen!!
Interaktion betyder, at eﬀekten af en variabel, for eksempel bilirubin,
afhænger af værdien af en anden variabel, for eksempel behandling,
men eﬀekterne ændres ikke med tiden.
(Stratifikation og interaktion kan sagtens kombineres.)
31
Ventetidsanalyse
Transformation af regressionsvariable
Behov for transformation af regressionsvariable
− Grafiske vurderinger af sammenhængen er mere kompliceret at
lave end for normalt fordelte outcomes, fordi vi ikke kan tegne
direkte
+ Principielle kriterier og simple numeriske vurderinger er som i
enhver anden statistisk model
32
Ventetidsanalyse
Transformation af regressionsvariable
Kriterier for transformation af regressionsvariable
Kriterier for valg af parametrisering/transformation
• Biologisk/medicinsk begrundelse (bedst, men sjældent muligt).
– Raten vokser eksponentielt med utransformerede kovariater,
mens en logaritmetransformation af en kovariat betyder, at
raten vokser med en fast faktor eller procentdel, hver gang
kovariaten vokser med f.eks. 10%.
• Transformationer brugt af andre (sammenlignelighed).
• Den “bedst mulige” transformation for de aktuelle data — pas på:
signifikansen vil blive overvurderet, og man kan ikke gå ud fra, at
den samme transformation er den mest optimale i andre data —
men det kan være et fornuftigt kriterium for konfounderne.
33
Ventetidsanalyse
Transformation af regressionsvariable
Kriterier for transformation af regressionsvariable
Kriteria for valg af parametrisering/transformation fortsat
• Se på fordelingen af den forklarende variabel: Nogle få ekstreme
værdier af den forklarende variabel kan have urimeligt stor
indflydelse på konklusionerne, medmindre variablen transformeres
nogle få ekstremt høje → loga (x),
nogle få ekstremt lave → exp(x/c) (sjældent brugt).
Det bør altid tjekkes, om den valgte transformation
giver en rimeligt lineær sammenhæng.
Ved at vælge XX = log(x)/ log(1.1) som kovariat, får man, at exp(b̂)
(Hazard Ratio) direkte estimerer den faktor, som raten skal ganges med
for en 10% forskel i kovariaten.
34
Ventetidsanalyse
Transformation af regressionsvariable
Fordelingen af serum bilirubin
PROC UNIVARIATE DATA=skl; *- Tjek, om fordelingen er skæv -*;
VAR bilirub;
HISTOGRAM / HEIGHT=5;
RUN;
35
Ventetidsanalyse
Transformation af regressionsvariable
Fordelingen af log(serum bilirubin)
DATA skl; SET skl; log2Bilirub=LOG2(bilirub); RUN;
PROC UNIVARIATE DATA=skl;
VAR log2Bilirub;
HISTOGRAM / HEIGHT=5;
RUN;
36
Ventetidsanalyse
Transformation af regressionsvariable
Test af behov for transformation af regressionsvariable
Simpelt, numerisk test: Inkluder både en utransformeret og en
transformeret version af variablen samtidigt – f.eks. X og log2 (X) –
for at se, om der er et klart svar på hvilken, der er den bedste
prediktor. Kræver en fornuftigt valg af transformation –
logaritmetransformation er ofte et oplagt valg, fordi det svarer til
stabil eﬀekt af relative forskelle, hvilket ofte er det intuitivt naturlige
valg.
37
Ventetidsanalyse
Transformation af regressionsvariable
Test af behov for transformation af en variabel
Inklusion af X og log(X) samtidigt (lige meget hvilken logaritme), giver
flere svar på en gang (gælder alle statistiske modeller, også almindelig
regression og logistisk regression):
1. Hvis log(X) er insignifikant, kan sammenhængen antages at være lineær
2. Hvis X er insignifikant, kan sammenhængen beskrives som konstant eﬀekt af
relative (%-vise) forskelle
3. Hvis begge er insignifikante, kan man selv vælge ud fra hvad, der er nemmest
at formidle,
4. Hvis begge er signifikante, må man undersøge det nærmere. Tjek eﬀekterne
af den lavest mulige og den højest mulige fordobling (eller 10-dobling eller
10% forøgning eller en anden relevant procentvis forskel).
NB: Det enkelte parameterestimat i denne model kan ikke fortolkes:
“Ændring i outcome, når X fordobles/10-dobles/ændres 10% samtidigt
med, at X holdes fast”, giver ikke mening!
38
Ventetidsanalyse
Transformation af regressionsvariable
Test af behov for transformation af serum bilirubin
Laveste observerede bilirubin er 4, højeste er 545 (=2 · 272.5).
Inklusion af serum bilirubin utransformeret såvel som
log2 -transformeret i SAS, med estimater for eﬀekterne af den lavest
mulige og den højest mulige fordobling i de aktuelle data:
PROC PHREG DATA=skl;
MODEL day*bld(0) = Sklero Bilirub log2Bilirub;
ESTIMATE "Mindst dbl" Bilirub 4
Log2Bilirub 1 / EXP ;
ESTIMATE "Størst dbl" Bilirub 272.5 Log2Bilirub 1 / EXP ;
Begge: TEST Bilirub=log2Bilirub=0;
RUN;
Ovenstående TEST-statement beregner Wald’s test for, at begge
regressionsparametre er lig med 0, så man kan udelade begge variable på
en gang, altså at serum bilirubin intet betyder for reblødningsraten.
39
Ventetidsanalyse
Transformation af regressionsvariable
40
Udpluk af output
Variable
DF
Sklero
1
Bilirub
1
log2Bilirub 1
Parameter
Estimate
-0.18290
-0.0001959
0.48004
Standard
Error Chi-Square
0.21596
0.7172
0.00231
0.0072
0.18152
6.9939
Pr>ChiSq
0.3971
0.9325
0.0082
Linear Hypotheses Testing Results
Wald
Label
Chi-Square
DF
Pr>ChiSq
Begge
22.8760
2
<.0001
Estimate
Standard
Label
Estimate
Error z Value
Mindst dbl 0.4793 0.1738
2.76
Størst dbl 0.4267 0.4873
0.88
Pr>|z| Exponentiated
0.0058
1.6149
0.3812
1.5321
Ventetidsanalyse
Transformation af regressionsvariable
Plot af sammenhængen med bilirubin som en lineær spline
Lineær spline giver både grafik og et numerisk test (dvs. p-værdi)
(Eksempel til selvstudium med test og estimater af hældningerne i de
enkelte intervaller – velegnet til tabel i artikler – kan findes s. 84-93 i
disse noter, og generaliserbar “automatisk” SAS-kode findes i det
samlede SAS program bagerst i noterne)
41
Ventetidsanalyse
Transformation af regressionsvariable
Estimation med log2 (serum bilirubin)
PROC PHREG DATA=skl;
MODEL day*bld(0) = sklero log2bilirub / RL;
RUN;
--------------------------------------------------Parameter Standard
Variable
DF Estimate
Error Chi-Square Pr>ChiSq
Sklero
1 -0.18373
0.21575
0.7252
0.3944
log2Bilirub 1
0.46716
0.09706
23.1656
<.0001
Analysis of Maximum Likelihood Estimates
Hazard
95% Hazard Ratio
Variable
Ratio
Confidence Limits
Sklero
0.832
0.545
1.270
log2Bilirub 1.595
1.319
1.930
Fortolkning af eﬀekten af serum bilirubin: En dobbelt så stor serum
bilirubin værdi ved baseline svarer til en næsten 60% højere
reblødningsrate for personer i samme behandlingsgruppe.
42
Ventetidsanalyse
Opsummering 2. del
Opsummering 2. del
Standard valg af regressionsmodel:
• Cox’s proportional hazards regressionsmodel
Forskellige måder variable kan indgå på:
• Regressionsvariable indgår via det prognostiske index
– lineære variable
– klassevariable
– interaktion
• Stratifikationsvariable
– helt forskellige rater, men samme eﬀekt af alle andre variable
Transformation af regressionsvariable:
• Kriterier for valg af transformation er helt ’som sædvanligt’
• Simpelt numerisk test (dvs. en p-værdi) er helt ’som sædvanligt’
• Grafisk test med tilhørende p-værdi er mere besværligt (lineær
spline)
43
Ventetidsanalyse
Særlig modelkontrol for Cox modellen
Modelkontrol specielt for Cox modellen
Cox modellen forudsætter proportionale rater,
R(t; X) = R0 (t) exp(bX), og dermed
ln(R(t; X)) = ln(R0 (t)) + bX
Grafisk tjek af proportionale rater: Stratificer for hver variabel, en ad
gangen, og plot
ln(Rstratum (t))
= ln( − ln(Sstratum (t)) )
som funktion af tiden. Kurverne skal være nogenlunde parallelle.
For kontinuerte variable skal man dels vælge nogle passende
opdelingspunkter, dels huske at den kontinuerte variabel stadig skal
indgå i det prognostiske index (se hvordan i programeksemplet allerbagest i
disse noter)
44
Ventetidsanalyse
Særlig modelkontrol for Cox modellen
Plot af log(kumuleret rate) for grafisk tjek af proportionalitet
Lineær tidsakse
Logaritmisk tidsakse
45
Ventetidsanalyse
Særlig modelkontrol for Cox modellen
46
SAS kode til grafisk tjek af proportionale rater
PROC PHREG DATA=skl;
MODEL Day*bld(0) = log2bilirub;
STRATA sklero;
BASELINE OUT=check LOGLOGS=LCumRate / METHOD=CH;
RUN;
SYMBOL1 C=BLACK V=NONE I=STEPLJ L=1 MODE=INCLUDE;
SYMBOL2 C=RED
V=NONE I=STEPLJ L=1 MODE=INCLUDE;
AXIS2 LABEL=(A=90);
PROC GPLOT DATA=check; *- Først et plot med almindelig lineær tidsakse -*;
PLOT LCumRate*day=sklero / VAXIS=AXIS2 NOLEGEND;
RUN;
AXIS1 LOGBASE=10 /* Log-skaleret akse med hak ved 1,2,..,10,20,..,100,200,.. */
INTERVAL=PARTIAL ; /* Aksen slutter ved den største observerede værdi */
PROC GPLOT DATA=check; *- Så et plot med logaritmisk tidsakse -*;
WHERE 0<Day; *- Bemærk, at Day=0 skal fjernes pga. logaritmen -*;
PLOT LCumRate*Day=Sklero / HAXIS=AXIS1 VAXIS=AXIS2 NOLEGEND;
GOPTIONS DEVICE=; *- SKAL med for SAS Enterprise, ellers helt unødvendig! -*;
RUN;
Ventetidsanalyse
Særlig modelkontrol for Cox modellen
47
Numerisk test af proportionalitet ved brug af tidsafhængige variable
Vælg passende opdeling (her 14 og 105 dage), hvor rate ratioen tillades at
være forskellig for de forskellige tidsintervaller. Dette gøres v.h.a.
tidsafhængige dummy variable (her Skl_0_14, Skl_14_105 og SklFra105).
Proportionalitetsantagelsen testes ved at teste, om de tilsvarende
parameterestimater kan antages at være ens:
PROC PHREG DATA=skl;
MODEL day*bld(0) = log2bilirub Skl_0_14 Skl_14_105 SklFra105;
IF
day<=14 THEN Skl_0_14=sklero;
ELSE Skl_0_14=0;
IF 14<day<=105 THEN Skl_14_105=sklero; ELSE Skl_14_105=0;
IF 105<day
THEN SklFra105=sklero; ELSE SklFra105=0;
TestProp: TEST
RUN;
Skl_0_14
=
Skl_14_105
=
SKLfra105;
For hvert reblødningstidspunkt beregnes variablene inden i PHREG for hver enkelt af de
patienter, der er under risiko på netop det tidspunkt (i’erne i nævneren på s. 29). I
beregningerne for hver enkelt patient i bruger SAS de værdier, som variablene (her
sklero) har for den enkelte patient i, undtagen for tids-responsvariablen (her day),
som er lig det aktuelle reblødningstidspunkt tj (reblødningstidspunktet for patient j
på s. 29)
Ventetidsanalyse
Særlig modelkontrol for Cox modellen
48
Del af output fra PROC PHREG
Parameter Standard
Chi- Pr>
Hazard 95% Hazard Ratio
Variable
DF Estimate
Error Square ChiSq Ratio Confidence Limits
log2Bilirub 1 0.45629 0.09640 22.4022 <.0001 1.578
1.306
1.906
Skl_0_14
1 0.15759 0.31013 0.2582 0.6114 1.171
0.637
2.150
Skl_14_105 1 -0.46164 0.37626 1.5053 0.2199 0.630
0.301
1.318
SklFra105
1 -0.68164 0.57135 1.4234 0.2329 0.506
0.165
1.550
Linear Hypotheses Testing Results
Wald
Label
Chi-Square
DF
Pr > ChiSq
TestProp
2.5117
2
0.2848
Mønster som forventet biologisk, men ingen statistisk signifikans
Ventetidsanalyse
Tidsskalaer
Virkelighedstjek!
Tid fra stop af seneste blødning har stor betydning for
sandsynligheden for at begynde at rebløde her og nu
• dels er det ikke muligt at rebløde, mens den foregående blødning
stadig er i gang,
• dels mindskes sandsynligheden for at rebløde markant med, hvor
længe det er siden, man blødte sidst.
Vi har altså to komplikationer:
1. Patienterne er ikke under risiko fra tid 0 for at blive observeret
med start på en reblødning i vores studie.
2. ’Tid’ er ikke bare tid – der er flere tidsskalaer at vælge imellem.
49
Ventetidsanalyse
Tidsskalaer
Tidsskalaer
Eksempler på tidsskalaer
• alder
• kalendertid
• tid siden sygdommen begyndte
• tid fra en eller anden hændelse af stor betydning for raten (her
tid fra stop af foregående blødning)
• tid fra randomisering (ofte et ikke-optimalt valg)
Den eneste forskel for det enkelte individ er definitionen af tid=0,
men det kan gøre en væsentlig forskel, fordi det kan have stor
betydning for hvem, der ellers var “under risiko” i studiet, da
hændelsen skete for et givet individ.
50
Ventetidsanalyse
Tidsskalaer
51
Forskellige tidsskalaer
Kalendertid:
Tid siden indgang i studiet:
Ventetidsanalyse
Tidsskalaer
52
Eksempler på forskellige tidsskalaer ved reblødning
Tid siden randomisering
B
Tid siden indgang i studiet:
B
A
50 år gammel
50 år gammel
1. blødnings start
1. blødnings start
1. blødning slutter
randomiseret
50 år gammel
A
50 år gammel
1. blødnings start
1. blødnings start
1. blødning slutter
randomiseret
1. blødning slutter
randomiseret
1. blødning slutter
2. blødnings start
2. blødnings start
randomiseret
2. blødnings start
2. blødnings start
Ventetidsanalyse
Tidsskalaer
Tidsskalaens betydning for Cox modellen
Hidtil har vi brugt tid siden randomisering som tid i Cox modellen –
men reblødningsraten afhænger langt mere af, hvor længe siden det
er, patienten holdt op med at bløde.
• Fordelen ved Cox modellen er, at den tillader en vilkårlig
sammenhæng mellem den underliggende rate og den valgte
tidsskala.
• Ratioen mellem raterne for to vilkårlige patienter i samme
stratum på et hvilket som helst tidspunkt afhænger kun af
regressionsvariablene, altså ikke af tidspunktet.
• Egenskaber ved en relevant tidsskala: Der skal være en god grund
til at tro, at raten ændres med tid siden tid 0 på samme måde for
alle individer i samme stratum.
53
Ventetidsanalyse
Tidsskalaer
Valg af tidsskala
Vælg den tidsskala, hvor raten varierer mest (eller mest ujævnt) med tiden
inden for det tidsinterval, hvor man har observationer!
Andre tidsskalaer kan indgå via regressionsvariable i Cox modellen og/eller
via stratifikation. Hvis betydningen af en alternativ tidsskala ikke følger
mønstret “en måned mere betyder altid det samme”, så er det nødvendigt
enten at stratificere eller at bruge tidsafhængige regressionsvariable (se det
numeriske test for proportionalitet s. 47 for et eksempel på SAS
programmering med tidsafhængige variable).
Brug aldrig en tidsskala, hvor man skal konstruere pseudotidspunkter for nogle af personerne. Det er et helt sikkert bevis
på et forkert valg af tidsskala!
54
Ventetidsanalyse
Tidsskalaer
Forsinket indgang
Årsag: Der skal ske et eller andet bestemt for individerne, før det ville
tælle med i vores analyse, hvis de oplevede vores “outcome” hændelse.
Dette sker på forskellige tidspunkter for forskellige individer, og det er ikke
altid sket før tid 0 i den valgte tidsskala.
Eksempler:
• for nogle patienter sker randomiseringen efter tid 0 i den valgte
tidsskala
• i analyser med tid siden randomisering kan patienterne ikke rebløde,
før de er holdt op med den første blødning
• nogen af de forklarende variable kræver en speciel undersøgelse, og
nogen af patienterne er nødt til at vente på denne undersøgelse
• for at blive inkluderet skal personerne være i live og “raske” ved
studiestart
55
Ventetidsanalyse
Tidsskalaer
56
Forsinket indgang i SAS: Option ENTRY=
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = log2bilirub sklero
/ ENTRY=t_entry RL;
RUN;
--------------------------------------------Model Information
Data Set
WORK.SKL
Entry Time Variable
t_entry
Dependent Variable
tNotBld
Censoring Variable
Bld
Censoring Value(s)
0
Ties Handling
BRESLOW
Percent
Total
Event
Censored
Censored
149
86
63
42.28
:
Analysis of Maximum Likelihood Estimates
Parameter Standard
Hazard 95% Hazard
Variable
DF Estimate
Error Chi-Square Pr>ChiSq
Ratio Confidence
log2Bilirub 1
0.43431 0.09580
20.5534 <.0001
1.544
1.280
Sklero
1 -0.16470 0.21682
0.5770
.4475
0.848
0.555
Ratio
Limits
1.863
1.297
Ventetidsanalyse
Konkurrerende afgangsårsager
Virkelighedstjek 2!
Når vi overhovedet interesserer os for at mindske reblødningsraten, er
det fordi reblødning er en af de væsentligste dødsårsager blandt
patienter, der mindst en gang har blødt fra åreknuder i spiserøret.
Men patienterne kan jo stadig dø, inden de når at rebløde!
57
Ventetidsanalyse
Konkurrerende afgangsårsager
Konkurrerende afgangsårsager
Individerne kan udgå af flere årsager, her reblødning eller død.
Konsekvenser:
• Teknisk: Andre afgangsårsager end den, der er i fokus,
behandles som censureringer
• Fortolkning: Raten og de estimerede eﬀekter har den samme
fortolkning som før –
MEN
• Kaplan-Meier kurven kan ikke fortolkes som sandsynligheden for
at undgå den hændelse, der er i fokus
Hvis flere slags hændelser er af interesse (som her), så analyseres hver
enkelt slags hændelse separat med de andre slags hændelser
behandlet som censureringer.
58
Ventetidsanalyse
Konkurrerende afgangsårsager
59
Separate analyser af de to typer af hændelser
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0)=sklero log2bilirub ascites / ENTRY=t_entry RL;
RUN;
Parameter Standard
Hazard 95% Hazard Ratio
Variable
DF Estimate Error Chi-Sq. Pr>ChiSq Ratio Confidence Limits
Sklero
1 -0.19124 0.22021 0.7542 0.3851
0.826
0.536
1.272
log2Bilirub 1 0.42240 0.09677 19.0542 <.0001
1.526
1.262
1.844
Ascites
1 0.15762 0.22776 0.4789 0.4889
1.171
0.749
1.829
------------------------------------------------------------------------PROC PHREG DATA=skl;
MODEL tNotBld*dead(0)=sklero log2bilirub ascites / ENTRY=t_entry RL;
RUN;
Parameter Standard
Hazard 95% Hazard Ratio
Variable
DF Estimate Error Chi-Sq. Pr>ChiSq Ratio Confidence Limits
Sklero
1 0.17358 0.35173 0.2435 0.6217
1.190
0.597
2.370
log2Bilirub 1 0.50353 0.14482 12.0890 0.0005
1.655
1.246
2.198
Ascites
1 0.93763 0.38166 6.0354 0.0140
2.554
1.209
5.396
Ventetidsanalyse
Konkurrerende afgangsårsager
60
Sandsynligheden for at være i live uden reblødning
Hvis de forskellige hændelser kombineres, fås en vurdering af eﬀekten
af de forklarende variable på tid til den første af disse hændelser sker:
DATA skl; SET skl;
status = bld + 2*dead;
RUN;
PROC PHREG DATA=skl;
MODEL tNotBld*status(0) = sklero log2bilirub ascites / ENTRY=t_entry RL;
RUN;
-------------------------------------------------------------------Parameter Standard
Hazard 95% Hazard Ratio
Variable
DF Estimate
Error Chi-Sq. Pr>ChiSq Ratio Confidence Limits
Sklero
1 -0.08715 0.18555 0.2206 0.6386
0.917
0.637 1.319
log2Bilirub 1 0.44557 0.08044 30.6819 <.0001
1.561
1.334 1.828
Ascites
1 0.37333 0.19303 3.7405 0.0531
1.453
0.995 2.121
Ventetidsanalyse
Konkurrerende afgangsårsager
Konkurrerende afgangsårsager
Sandsynligheden for, at reblødningen sker efter præcis t dage, er lig
med
S(t; X) · r(t; X)
hvor S(t; X) er sandsynligheden for at være i live uden reblødning
indtil (men ikke med) dag t. Sandsynligheden for at have observeret
en reblødning inden tid T fås ved at summere disse for alle t ≤ T
Konsekvens:
• Faktorer, som ikke påvirker raten for en bestemt slags hændelse,
kan have en stærk eﬀekt på sandsynligheden for at opleve
hændelsen via en eﬀekt på raten for en konkurrerende hændelse
og dermed en eﬀekt på sandsynligheden for at være under risiko.
• Eﬀekten på sandsynligheden for at opleve hændelsen kan i
ekstreme tilfælde være modsat rettet eﬀekten på selve raten.
61
Ventetidsanalyse
Rate ratio, odds ratio og risk ratio
Ventetidsanalyse vs. analyse af sandsynlighed for hændelsen
• Meningsfuld definition af “sandsynlighed for hændelsen” kræver,
at man vælger et fast tidsrum (glemmes desværre ofte, med deraf
følgende problematisk fortolkning og konklusion)
• Ventetidsanalyse (rate ratio) har større styrke end analyser af
sandsynligheden for, at hændelsen er sket inden det bestemte,
faste tidspunkt – uanset om man laver logistiske (odds ratio) eller
log-lineære (risk ratio) regressionsanalyser
• Analyser af “sandsynlighed for hændelsen” er direkte forkerte, når
der er censureringer – det gælder både logistisk regression (odds
ratio) og log-lineær regression (risk ratio). Censureringer må
hverken ses som “ingen hændelse” eller helt udelades af analysen!
Kun ventetidsanalyse er korrekt!
62
Ventetidsanalyse
Rate ratio, odds ratio og risk ratio
Sammenligning af rate ratio, odds ratio og risk ratio
Alle tre kan desværre blive kaldt det upræcise “relativ risiko” i artikler!
Simple modeller: Hvis sandsynlighederne er meningsfuldt defineret
(fast tidsrum og ingen censurering), viser rate (“hazard”) ratio (HR)
altid stærkere eﬀekt end odds ratio (OR), som igen viser stærkere
eﬀekt end risk ratio (RR). Matematisk relation:
Enten er
• 1 < RR < OR < HR
eller
• 1 = RR = OR = HR
eller også er
• 1 > RR > OR > HR
Multiple regressionsmodeller: Modellerne passer ikke sammen!
Eﬀekterne kan ikke være multiplikative på alle 3 skalaer, så (mindst)
2 ud af 3 er forkerte! Hvis hændelsen er meget sjælden, ligner de dog
hinanden
63
Ventetidsanalyse
Opsummering 3. del
Opsummering 3. del
Modelkontrol af Cox model:
• proportionale rater
− grafisk vurdering ved plot af log-transformeret kumuleret rate
− numerisk test med tidsafhængige variable
Tidsskala:
• Valg af tidsskala
• Forsinket indgang
Konkurrerende afgangsårsager:
• separate analyser med de andre slags hændelser som censurering
• eﬀekten på sandsynligheden kan være meget anderledes end eﬀekten
på raten
Ventetidsanalyser vs. logistisk eller log-lineær regression
• Ventetidsanalyser har størst styrke
• Hvis der er censureringer, så er logistisk regression og log-lineær
regression direkte forkert
64
Censurerede, normalt fordelte data
Censurerede, normalt fordelte data
Eksempel med SAS programbidder
– til selvstudium –
65
Censurerede, normalt fordelte data (data med detektionsgrænse)
Data med detektionsgrænse
For lave værdier er der “venstrecensurering”, dvs. man ved bare, at
værdien er under en given grænse:
Baggrundsstøj eller begrænsning i måleudstyrets/assayets følsomhed
har medført, at man ikke kan skelne mellem meget lave værdier.
Man tør godt gøre antagelser om fordelingen af de sande målinger,
eksempelvis at residualerne er normalt fordelt (eventuelt efter
logaritmetransformation af målingerne).
66
Censurerede, normalt fordelte data (data med detektionsgrænse)
Eksempel på venstrecensurerede data
Målinger af NO2 indendørs og udendørs (Raaschou-Nielsen et al., 1997)
Vi har 85 sæt af samhørende mål for NO2
1. udenfor gadedøren
2. i soveværelset
med en detektionsgrænse på 0.75.
Vi ønsker at undersøge, hvor stor indflydelse udendørsniveauet har på
indendørsniveauet.
67
Censurerede, normalt fordelte data (data med detektionsgrænse)
Eksempel på venstrecensurerede data
Samhørende mål for NO2 inde og ude
(under detektionsgrænsen er plottet som lig med detektionsgrænsen)
68
Censurerede, normalt fordelte data (data med detektionsgrænse)
Estimation af sammenhæng
Hvad med bare at udelade de “ukendte”? Pas på: Det er selektion
baseret på responsvariablen!
DATA no2; SET no2;
IF UPCASE(UnderDetekt)="JA" THEN inde=.;
ude_25 = ude - 2.5; * Centrering af variabel ;
RUN;
PROC REG DATA=no2;
Udensmaa: MODEL inde = ude_25;
RUN;
69
Censurerede, normalt fordelte data (data med detektionsgrænse)
The REG Procedure
Model: Udensmaa
Dependent Variable: INDE
Analysis of Variance
Sum of
Mean
Source
DF Squares Square F Value Pr>F
Model
1 9.18865 9.18865 107.71 <.0001
Error
58 4.94781 0.08531
Corrected Total 59 14.13645
Root MSE
0.29207
Dependent Mean 1.52430
Coeff Var
19.16120
Variable
Intercept
ude_25
R-Square 0.6500
Adj R-Sq 0.6440
Parameter Estimates
Parameter Standard
DF Estimate
Error t Value Pr>|t|
1
1
1.60065
0.60009
0.03842
0.05782
41.66 <.0001
10.38 <.0001
70
Censurerede, normalt fordelte data (data med detektionsgrænse)
Fit uden data under detektionsgrænsen
Samhørende mål for NO2 inde og ude
Duer ikke på grund af bias.
Vi mangler alle de laveste indendørsværdier, så linien ligger for højt
for de lave udendørsværdier, og spredningen undervurderes!
71
Censurerede, normalt fordelte data (data med detektionsgrænse)
Estimation af sammenhæng
Kan vi ikke bare sætte alle observationer under detektionsgrænsen lig
med detektionsgrænsen (ligesom i tegningen)?
DATA no2; SET no2;
IF UPCASE(UnderDetekt)="JA" THEN inde=0.75;
RUN;
PROC REG DATA=no2;
Naiv: MODEL inde = ude_25;
RUN;
72
Censurerede, normalt fordelte data (data med detektionsgrænse)
73
The REG Procedure
Model: Naiv
Dependent Variable: INDE
Analysis of Variance
Source
DF
Sum of
Squares
Mean
Square
F Value Pr>F
Model
1 18.15521 18.15521 229.66 <.0001
Error
83 6.56128 0.07905
Corrected Total 84 24.71649
Root MSE
0.28116
Dependent Mean 1.29656
Coeff Var
21.68511
R-Square 0.7345
Adj R-Sq 0.7313
Parameter Estimates
Parameter Standard
Variable DF Estimate
Error t Value Pr>|t|
Intercept 1 1.55732
0.03502 44.48 <.0001
ude_25
1 0.64260
0.04240 15.15 <.0001
Censurerede, normalt fordelte data (data med detektionsgrænse)
Fit for naiv model
Samhørende mål for NO2 inde og ude
Duer heller ikke på grund af bias.
Vi overvurderer alle de laveste indendørsværdier, så linien ligger stadig
for højt for de lave udendørsværdier, og spredningen undervurderes!
74
Censurerede, normalt fordelte data (data med detektionsgrænse)
Estimation af sammenhæng
Alternativ (lidt bedre, men stadig tvivlsomt): Observationer under
detektionsgrænsen sættes lig et gæt på gennemsnitsværdien for
observationer under detektionsgrænsen ( 23 ·detektionsgrænsen svarer til, at
et histogram over de sande værdier ville ligne en trekant fra 0 op til
detektionsgrænsen):
DATA no2; SET no2;
IF UPCASE(UnderDetekt)="JA" THEN justeret = (2/3)*0.75;
ELSE justeret = inde;
RUN;
PROC REG DATA=no2;
Adhoc: MODEL justeret = ude_25;
RUN;
(Hvis man tror på, at alle værdier mellem 0 og detektionsgrænsen er lige sandsynlige,
så er 12 ·detektionsgrænsen et relevant gæt. Hvis man tror, at histogrammet ville følge
en parabel, er 34 ·detektionsgrænsen et relevant gæt.)
75
Censurerede, normalt fordelte data (data med detektionsgrænse)
76
The REG Procedure
Model: Adhoc
Dependent Variable: justeret
Analysis of Variance
Source
DF
Sum of
Squares
Mean
Square F Value Pr>F
Model
1 23.92209 23.92209
Error
83 8.72937 0.10517
Corrected Total 84 32.65146
Root MSE
0.32430
Dependent Mean 1.22303
Coeff Var
26.51639
227.45 <.0001
R-Square 0.7326
Adj R-Sq 0.7294
Parameter Estimates
Parameter Standard
Variable DF Estimate
Error t Value Pr>|t|
Intercept 1
1.52235 0.04039
37.69 <.0001
ude_25
1
0.73763 0.04891
15.08 <.0001
Censurerede, normalt fordelte data (data med detektionsgrænse)
Fit med gæt på gennemsnitsværdi
Samhørende mål for NO2 inde og ude
Vi ved ikke, om linien ligger for højt (for højt gæt på gennemsnit)
eller for lavt (for lavt gæt på gennemsnit) for de lave
udendørsværdier, men spredningen undervurderes stadig.
77
Censurerede, normalt fordelte data (data med detektionsgrænse)
Korrekt estimation af sammenhæng
Optimal udnyttelse af data opnås ved at inkludere præcis den viden,
vi har:
De censurerede data er mindre end detektionsgrænsen.
Princip:
• For hver af observationerne under detektionsgrænsen kan det
udregnes, hvad sandsynligheden er for, at observationen ligger under
detektionsgrænsen for givne værdier af alle modelparametrene, både
regressionsparametrene og variansparameteren. Denne sandsynlighed
indgår i estimationen.
• De “almindelige” observationer indgår på sædvanlig vis.
Alle modelparametrene estimeres som de værdier, der alt i alt passer bedst
med det observerede, inklusive at observationerne under detektionsgrænsen
faktisk lå under denne grænse.
78
Censurerede, normalt fordelte data (data med detektionsgrænse)
I SAS skal der være to variable, der viser de faktiske nedre henholdsvis
øvre grænser for responsvariablens værdi.
• For observerede værdier er nedre og øvre grænse lig observationen.
• For observationer under detektionsgrænsen er den nedre grænse
ukendt, hvilket angives som et . (missing) i SAS, og den øvre grænse er
lig med detektionsgrænsen.
DATA no2; SET no2;
IF UPCASE(UnderDetekt)="JA"
THEN DO;
nedre = .; oevre = 0.75;
END;
ELSE DO;
nedre = inde; oevre = inde;
END;
RUN;
PROC LIFEREG DATA=no2;
MODEL (nedre, oevre) = ude_25 / DIST=NORMAL NOLOG;
RUN;
79
Censurerede, normalt fordelte data (data med detektionsgrænse)
The LIFEREG Procedure
Model Information
Data Set
Dependent Variable
Dependent Variable
Number of Observations
Noncensored Values
Right Censored Values
Left Censored Values
Interval Censored Values
Name of Distribution
Log Likelihood
Algorithm converged.
WORK.NO2
nedre
oevre
85
60
0
25
0
Normal
-35.88065877
Effect
ude_25
80
Type III Analysis of Effects
Wald
DF
Chi-Square
Pr > ChiSq
1
177.8626
<.0001
Analysis of Parameter Estimates
Standard 95% Confidence
Parameter DF Estimate Error
Limits
Intercept 1 1.5203 0.0431
1.4359 1.6047
ude_25
1 0.7845 0.0588
0.6692 0.8997
Scale
1 0.3403 0.0320
0.2830 0.4092
Parameter
Intercept
ude_25
Scale
Chi-Square Pr > ChiSq
1245.07
<.0001
177.86
<.0001
113.08
<.0001
Censurerede, normalt fordelte data (data med detektionsgrænse)
Optimalt fit
Samhørende mål for NO2 inde og ude
81
Censurerede, normalt fordelte data (data med detektionsgrænse)
Estimation af standard deviation
scale = maximum likelihood estimat for standard deviationen (SD)
(=residualspredning=prediktionsspredning). For at få noget, der er
sammenligneligt med det sædvanlige estimat (“Root MSE” i SAS
output), skal der justeres:
√
n
SD* = scale ·
n−k−1
(n = antal observationer, k = antal kovariater).
√
Her fås SD*= 0.340 · 85
83 = 0.344.
82
Censurerede, normalt fordelte data (data med detektionsgrænse)
83
Sammenligning af resultater
Oversigt over resultaterne af de 4 analyser:
Model
Udensmaa
Naiv
Adhoc
Optimal
Ude_25
Parameter Standard
Estimate
Error
0.600
0.058
0.643
0.042
0.738
0.049
0.785
0.059
SD ("Root MSE")
0.292
0.281
0.324
0.344
Lineære splines
84
SAS kode til lineære splines
– til selvstudium –
Lineære splines
85
Konstruktion af en lineær spline
Lineære splines
86
Lineær spline, princip:
Vi konstruerer nogle variable, der tilsammen summerer til serum
bilirubin-variablen eventuelt fratrukket et passende tal. Hver variabel
fanger den del af variationen i serum bilirubin-værdierne, der ligger i et
præ-defineret interval.
Af hensyn til den statistiske styrke kan det anbefales at lægge intervallerne,
så der er ca. lige mange hændelser i hvert interval.
Her er valgt 4 variable, der summerer til bilirubin minus medianen blandt
reblødere. Interval-endepunkterne/knudepunkterne (engelsk “knots”) er
placeret ved kvartilerne blandt rebløderne: 26, 47 og 73
• b_u26 fanger variationen i bilirubin for værdier under 26
• b_26_47 fanger variationen i bilirubin for værdier mellem 26 og 47
• b_47_73 fanger variationen i bilirubin for værdier mellem 47 og 73
• b_o73 fanger variationen i bilirubin for værdier over 73
Lineære splines
87
Beregning af observerede percentiler blandt reblødere
PROC UNIVARIATE DATA=skl PCTLDEF=3; WHERE bld=1; VAR bilirub; RUN;
(PCTLDEF=3 giver observerede percentiler)
Quantile
100% Max
99%
95%
90%
75% Q3
50% Median
25% Q1
10%
5%
1%
0% Min
Estimate
495
495
177
146
73
47
26
15
12
4
4
Mulig kode for de nødvendige ekstra variable i SAS
Lineære splines
DATA skl; SET skl;
IF NOT MISSING(bilirub) AND bilirub<=26 THEN DO; *- de laveste bilirubinværdier -*;
b_u26 = bilirub - 26; *- beregner hvor meget under 26 -*;
b_26_47 = 26 - 47; *- starter med bundværdien 26-47 = -21 -*;
b_47_73 = 0;
b_o73 = 0;
END;
IF 26<bilirub<=47 THEN DO; *- bilirubin mellem laveste kvartil og median -*;
b_u26 = 0; *- stopper på topværdien 0 -*;
b_26_47 = bilirub - 47;
b_47_73 = 0;
b_o73 = 0;
END;
IF 47<bilirub<=73 THEN DO; *- bilirubin mellem median og højeste kvartil -*;
b_u26 = 0;
b_26_47 = 0;
b_47_73 = bilirub - 47;
b_o73 = 0;
END;
IF 73<bilirub THEN DO; *- de højeste bilirubinværdier -*;
b_u26 = 0;
b_26_47 = 0;
b_47_73 = 73 - 47; *- stopper på topværdien 73-47 = 26 -*;
b_o73 = bilirub - 73;
END;
RUN;
88
Lineære splines
89
Kort “smart” SAS-kode, der gør det samme:
DATA skl;
SET skl;
IF NOT MISSING(bilirub) THEN DO;
b_u26 = MIN(bilirub-26, 0);
b_26_47 = MIN( MAX(26-47,bilirub-47), 0);
b_47_73 = MAX( MIN(73-47,bilirub-47), 0);
b_o73 = MAX(bilirub-73, 0);
END;
RUN;
MIN finder den mindste værdi og MAX finder den største værdi af de
variable eller tal, der står i parentesen
Lineære splines
90
Estimation og test af lineær spline
PROC PHREG DATA=skl;
MODEL day*bld(0) =
b_u26 b_26_47 b_47_73 b_o73
sklero
/ RL;
TestLine: TEST b_u26 = b_26_47 = b_47_73 = b_o73;
RUN;
Lineære splines
91
Estimation og test af lineær spline fortsat
Variable DF
b_u26
1
b_26_47
1
b_47_73
1
b_o73
1
Sklero
1
Parameter
Estimate
-0.01390
0.03250
0.03483
0.00162
-0.13197
Standard
Error
0.03276
0.02161
0.01412
0.00161
0.21954
Chi-Sq. Pr>ChiSq
0.1801
0.6712
2.2618
0.1326
6.0838
0.0136
1.0189
0.3128
0.3613
0.5478
Hazard
Ratio
0.986
1.033
1.035
1.002
0.876
95% Hazard Ratio
Confidence Limits
0.925
1.052
0.990
1.078
1.007
1.065
0.998
1.005
0.570
1.348
Linear Hypotheses Testing Results
Label
Wald Chi-Square DF
Pr > ChiSq
TestLine
15.7811
3
0.0013
“Parameter Estimate” er hældningen for ln(rate ratio) inden for
hvert af intervallerne. “Hazard Ratio” er derfor et mål for den
interval-specifikke dosis-respons relation.
Lineære splines
Plot af lineær spline: Datasæt, der skal plottes
5 og 95 percentilerne for bilirubin blandt reblødere er 12 og 177
DATA Plot;
DO bili=12, 26, 47, 73, 177; * Knækpunkter og endepunkter;
* Variabelværdier svarende til knæk- og endepunkter beregnes ;
b_u26=MIN(bili-26,0);
b_26_47=MIN( MAX(26-47,bili-47), 0);
b_47_73=MAX( MIN(73-47,bili-47), 0);
b_o73=MAX(bili-73,0);
* Parameterestimaterne fra outputtet ganges på de tilsvarende variable ;
* for at udregne den prædikterede værdi af det Prognostiske Index
;
* (eksponenten i Cox modellen)
;
pi = -0.01390*b_u26 +0.03250*b_26_47
+0.03483*b_47_73 +0.00162*b_o73;
rr = EXP(pi); * PI omregnes til Rate Ratio ;
OUTPUT; * Gem værdierne af variablene på SAS-datasættet Plot ;
END;
RUN;
92
Lineære splines
93
Plot af lineær spline
Eksempel på et kald af PROC GPLOT med SYMBOL og AXIS statements:
SYMBOL1 C=BLACK V=CIRCLE I=JOIN L=1 MODE=INCLUDE;
AXIS1 LABEL=(F="Andalus" "Bilirubin") ;
AXIS2 LABEL=(F="Andalus" A=90 "Rate ratio")
LOGBASE=10 INTERVAL=PARTIAL MAJOR=(H=1) MINOR=(H=1);
PROC GPLOT DATA=plot;
PLOT rr*bili / HAXIS=AXIS1 VAXIS=AXIS2;
GOPTIONS DEVICE=; *- SKAL med for SAS Enterprise -*;
RUN;
Eksempel på program til ventetidsanalyse
Samlet SAS program eksempel
– til selvstudium –
94
Eksempel på program til ventetidsanalyse
95
*-----------------------------------------------------------------------*;
*- AXIS-statements og LEGEND-statement til brug i resten af programmet -*;
*-----------------------------------------------------------------------*;
GOPTIONS FTEXT="Times New Roman";
/* bestemmer skrifttypen i GPLOT */
*--- Til almindelig, lineært skaleret X-akse ---*;
AXIS1 LABEL=(H=2.5)
/* bestemmer skriftstørrelsen af aksens tekst */
VALUE=(H=2);
/* skriftstørrelsen af aksens talværdier */
*--- Til logaritmisk skaleret X-akse ---*;
AXIS11 LABEL=(H=2.5) VALUE=(H=2)
LOGBASE=10
/* log-skaleret med hak ved 1,2,..,10,20,..,100,200,.. */
INTERVAL=PARTIAL
/* aksen slutter ved den største værdi i data */
MAJOR=(H=1)
/* størrelsen på hakkene ved tal */
MINOR=(H=1);
/* størrelsen på hakkene mellem tallene */
*--- Til almindelig, lineært skaleret Y-akse ---*;
AXIS2 LABEL=(A=90 /* skriver teksten lodret, dvs. langs Y-aksen */ H=2.5)
VALUE=(H=2);
*--- Til logaritmisk skaleret Y-akse ---*;
AXIS22 LABEL=(A=90 H=2.5) VALUE=(H=2)
LOGBASE=10 INTERVAL=PARTIAL MAJOR=(H=1) MINOR=(H=1);
*--- LEGEND-statement, så legend er nemmere kan læses ---*;
LEGEND1 LABEL=(H=1.5) VALUE=(H=1.5 JUSTIFY=LEFT);
Eksempel på program til ventetidsanalyse
*------------------------------------------------------*;
*- Kaplan-Meier kurver i 2 versioner:
-*;
*- ODS GRAPHICS og en, hvor man selv styrer akser mv. -*;
*------------------------------------------------------*;
ODS GRAPHICS ON;
PROC PHREG DATA=skl PLOTS(OVERLAY=ROW CL)=S;
MODEL tNotBld*Bld(0) =
/
ENTRY=t_entry /* Skal med, hvis forsinket indgang ved tid "t_entry" */
;
STRATA Sklero;
BASELINE OUT=km
/* Lav datasæt KM til egne plots
SURVIVAL=KMcurves
/* Kaplan-Meier estimat
LOWER=LowerBound UPPER=UpperBound
/* Punktvise konfidensgrænser
/ METHOD=PL
/* Estimeres ved Kaplan & Meier’s Product Limit metode
CLTYPE=LOGLOG;
/* Holder konfidensgrænserne inden for 0-1
RUN;
96
*/
*/
*/
*/
*/
Eksempel på program til ventetidsanalyse
*- Helt simpelt Kaplan-Meier uden konfidensgrænser, -*;
*- men hvor man selv kan styre akser mm. til artikel -*;
SYMBOL1 C=BLACK V=NONE I=STEPLJ MODE=INCLUDE;
SYMBOL2 C=RED
V=NONE I=STEPLJ MODE=INCLUDE;
*- AXIS og LEGEND definitionerne står i starten af programmet -*;
PROC GPLOT DATA=km;
PLOT kmcurves*tNotBld=sklero
/
HAXIS=AXIS1 VAXIS=AXIS2 LEGEND=LEGEND1
;
RUN; QUIT;
97
Eksempel på program til ventetidsanalyse
*- Gør data klar til Kaplan-Meier MED konfidensgrænser, -*;
*- og hvor man selv kan styre akser mm. til artikel
-*;
DATA KM_m_CL; SET km;
LABEL Curve = "S(t) = Kaplan-Meier";
IF sklero=0 THEN DO;
type=1; curve=kmcurves; OUTPUT;
type=2; curve=lowerbound; OUTPUT;
type=3; curve=upperbound; OUTPUT;
END;
IF sklero=1 THEN DO;
type=4; curve=kmcurves; OUTPUT;
type=5; curve=lowerbound; OUTPUT;
type=6; curve=upperbound; OUTPUT;
END;
RUN;
98
Eksempel på program til ventetidsanalyse
99
*- Plot af figuren: -*;
SYMBOL1
SYMBOL2
SYMBOL3
SYMBOL4
R=1
R=2
R=1
R=2
C=BLACK
C=BLACK
C=RED
C=RED
V=NONE
V=NONE
V=NONE
V=NONE
I=STEPLJ
I=STEPLJ
I=STEPLJ
I=STEPLJ
L=2 W=5 MODE=INCLUDE;
L=4 MODE=INCLUDE;
L=1 W=5 MODE=INCLUDE;
L=41 MODE=INCLUDE;
*- AXIS definitionerne står i starten af programmet -*;
PROC GPLOT DATA=KM_m_CL;
PLOT curve*tNotBld=type
/ HAXIS=AXIS1 VAXIS=AXIS2
NOLEGEND;
/* Fravælger symbolforklaring */
RUN; QUIT;
Eksempel på program til ventetidsanalyse
*------------------------------------*;
*- Kumuleret rate, igen 2 versioner -*;
*------------------------------------*;
ODS GRAPHICS ON;
PROC PHREG DATA=skl PLOTS(OVERLAY=ROW CL)=CUMHAZ;
MODEL tNotBld*Bld(0) =
/
ENTRY=t_entry /* Skal med, hvis der er forsinket indgang */
;
STRATA Sklero;
BASELINE
OUT=CumulatedHazard
LOGSURV=LogSurvival
LOWER=LowerBound UPPER=UpperBound
/ METHOD=CH /* Beregner den kumulerede rate som sum af dagsraterne */
CLTYPE=LOGLOG;
RUN;
100
Eksempel på program til ventetidsanalyse
101
DATA CumulatedHazard;
SET CumulatedHazard;
KumRate = -LogSurvival; *- jf. noterne s. 10 -*;
LowerKumRate = -log(UpperBound); *- Konfidensgrænserne er for survival, -*;
UpperKumRate = -log(LowerBound); *- så de skal også transformeres.
-*;
RUN;
DATA CH_m_CL;
SET CumulatedHazard;
LABEL Curve = "Nelson-Aalen kum. rate";
IF sklero=0 THEN DO;
type=1; curve=KumRate; OUTPUT;
type=2; curve=LowerKumRate; OUTPUT;
type=3; curve=UpperKumRate; OUTPUT;
END;
IF sklero=1 THEN DO;
type=4; curve=KumRate; OUTPUT;
type=5; curve=LowerKumRate; OUTPUT;
type=6; curve=UpperKumRate; OUTPUT;
END;
RUN;
Eksempel på program til ventetidsanalyse
SYMBOL1
SYMBOL2
SYMBOL3
SYMBOL4
R=1
R=2
R=1
R=2
C=BLACK
C=BLACK
C=RED
C=RED
V=NONE
V=NONE
V=NONE
V=NONE
102
I=STEPLJ
I=STEPLJ
I=STEPLJ
I=STEPLJ
L=2 W=5 MODE=INCLUDE;
L=4 MODE=INCLUDE;
L=1 W=5 MODE=INCLUDE;
L=41 MODE=INCLUDE;
*- AXIS definitionerne står i starten af programmet -*;
PROC GPLOT DATA=CH_m_CL;
PLOT curve*tNotBld=type
RUN; QUIT;
/ HAXIS=AXIS1 VAXIS=AXIS2 NOLEGEND;
*-----------------*;
*- Log-rank test -*;
*-----------------*;
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = sklero
/ TIES=DISCRETE /* For at beregne log-rank HELT korrekt */
ENTRY=t_entry /* Skal med, hvis der er forsinket indgang */
RL;
/* Konfidensgrænser på Hazard Ratio */
RUN;
Eksempel på program til ventetidsanalyse
103
*--------------------------*;
*- Cox’s regressionsmodel -*;
*--------------------------*;
DATA skl;
SET skl;
log2Bilirub=LOG2(bilirub);
RUN;
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = log2bilirub sklero
/ RL ENTRY=t_entry; *- ENTRY= ved forsinket indgang -*;
RUN;
/*
/*
/*
/*
/*
/*
/*
For små datasæt bør man bruge options
/ RL=PL TIES=EXACT ;
(og selfølgelig
ENTRY=...
hvis relevant).
RL=PL giver mere korrekte konfidensgrænser for Hazard Ratio for
små datasæt. Det gør næsten ingen forskel her.
TIES=EXACT beregner likelihoodfunktionen helt korrekt, når der er
flere hændelser på samme tid -- men det kan tage MEGET lang tid,
så brug det kun på små datasæt!
*/
*/
*/
*/
*/
*/
*/
Eksempel på program til ventetidsanalyse
*---------------------------------------------------*;
*- Er transformation af serum bilirubin nødvendig? -*;
*---------------------------------------------------*;
*- Find min og max for bilirubin i datasættet -*;
PROC MEANS MIN MAX DATA=skl;
VAR bilirub;
RUN;
*- Simpelt tjek: Inklusion af både utransformeret -*;
*- og log-transformeret bilirubin på en gang.
-*;
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = bilirub log2bilirub sklero
/ RL
ENTRY=t_entry ; *- ENTRY= ved forsinket indgang -*;
/*--- Kan man helt udelade bilirubin? ---*/
TestBoth: TEST Bilirub=Log2Bilirub=0;
/*--- Er der stor forskel på betydningen af en fordobling ---*/
/*--- i den lave ende og en fordobling i den høje ende?
---*/
ESTIMATE "Mindst dbl" Bilirub 4
Log2Bilirub 1 / EXP ;
ESTIMATE "Størst dbl" Bilirub 272.5 Log2Bilirub 1 / EXP ;
RUN;
104
Eksempel på program til ventetidsanalyse
105
*---------------------------------------------------------------*;
*- Vurdering v.h.a. lineær spline (koden kan bruges til enhver -*;
*- anden variabel -- man skal bare estatte ’bilirub’ med det
-*;
*- andet variabelnavn på denne slide OG den næste slide).
-*;
*---------------------------------------------------------------*;
******-
Find relevante percentiler blandt personer med events, fordi styrken
til at estimere hældningen i et interval afhænger af hvor mange
events, der er blandt personer med bilirubin i det interval.
Percentilerne gemmes direkte på datasættet PCTLDATA som
xpctl5, xpctl25, xpctl50, xpctl75, xpctl95, xpctl33_3 og xpctl66_7
(de sidste to skal bruges ved grafisk tjek af proportionalitet senere)
-*;
-*;
-*;
-*;
-*;
-*;
PROC UNIVARIATE DATA=skl PCTLDEF=3 /* PCTLDEF=3 giver observerede percentiler */
NOPRINT; *- NOPRINT: uden output, kun nyt datasæt -*;
WHERE bld=1; *- skal kun med, fordi det er en ventetidsanalyse med BLD som event -*;
VAR bilirub; *- BILIRUB kunne erstattes af en hvilken som helst variabel -*;
OUTPUT OUT=pctldata
PCTLPRE=Xpctl
PCTLPTS=5 25 50 75 95
33.3 66.7;
RUN;
*- Se det nye datasæt: -*;
PROC PRINT DATA=pctldata; RUN;
Eksempel på program til ventetidsanalyse
106
*- Beregne BILIRUB spline hjælpevariable, der skal bruges i PHREG. -*;
*- Husk: BILIRUB kunne erstattes af en hvilken som helst variabel! -*;
DATA skl;
SET skl;
IF _N_=1 THEN SET pctldata; *- lægger percentilerne ind i SKL datasættet -*;
*- Beregning af hjælpevariablene med den "smarte" kode -*;
IF NOT MISSING(bilirub) THEN DO;
b_u25pctl=MIN(bilirub-Xpctl25,0);
b_25_50pctl=MIN(MAX(Xpctl25-Xpctl50,bilirub-Xpctl50),0);
b_50_75pctl=MAX(MIN(Xpctl75-Xpctl50,bilirub-Xpctl50),0);
b_o75pctl=MAX(bilirub-Xpctl75,0);
END;
RUN;
Eksempel på program til ventetidsanalyse
*- Estimation og test -*;
PROC PHREG DATA=skl
OUTEST=parms /* gemmer parameterestimaterne på datasættet PARMS */
/* i variable med de samme navne som de tilsvarende */
/* variable i analysen (f.eks. b_u25pctl)
*/
;
MODEL tNotBld*bld(0) =
b_u25pctl
b_25_50pctl
b_50_75pctl
b_o75pctl
sklero
/
RL ENTRY=t_entry; *- ENTRY= ved forsinket indgang -*;
TestLine: TEST b_u25pctl = b_25_50pctl = b_50_75pctl = b_o75pctl;
RUN;
107
Eksempel på program til ventetidsanalyse
108
*- Beregning af punkter til plot af lineær spline -*;
DATA Plot;
MERGE pctldata parms; *- Begge data indeholder kun 1 obs (=1 linje) -*;
LABEL
bili = "Bilirubin" /* Hvis man har erstattet BILIRUB med en anden variabel, skal denne
label rettes. Resten kunne principelt bruges uden ændringer */
rr = "Rate ratio"
;
*- Man behøver kun beregne støttepunkter i enderne og ved knækkene -*;
DO bili=xpctl5, xpctl25, xpctl50, xpctl75, xpctl95;
*- De nye dummy-variable skal hedde noget lidt anderledes, fordi de -*;
*- "gamle" navne jo "er brugt" af PHREG til parameterestimaterne.
-*;
b_u25=MIN(bili-Xpctl25,0);
b_25_50=MIN(MAX(Xpctl25-Xpctl50,bili-Xpctl50),0);
b_50_75=MAX(MIN(Xpctl75-Xpctl50,bili-Xpctl50),0);
b_o75=MAX(bili-Xpctl75,0);
*- Beregning af det prognostiske index PI -*;
pi = b_u25pctl*b_u25 + b_25_50pctl*b_25_50 + b_50_75pctl*b_50_75 + b_o75pctl*b_o75;
*- Beregning af rate ratio (hazard ratio): -*;
rr = EXP(pi);
*- Læg observationen ud på Plot datasættet: -*;
OUTPUT;
END;
RUN;
Eksempel på program til ventetidsanalyse
*- Tegning af lineær spline -*;
SYMBOL1 C=BLACK V=CIRCLE I=JOIN L=1 MODE=INCLUDE;
*- AXIS definitionerne står i starten af programmet -*;
PROC GPLOT DATA=plot;
PLOT rr*bili / HAXIS=AXIS1 VAXIS=AXIS22;
GOPTIONS DEVICE=; *- SKAL med for SAS Enterprise, -*;
*- er ellers helt unødvendig!
-*;
RUN; QUIT;
109
Eksempel på program til ventetidsanalyse
110
*-------------------------------------------------------------------------*;
*- Test af vekselvirkning og estimation af effekterne ved vekselvirkning -*;
*-------------------------------------------------------------------------*;
PROC PHREG DATA=skl;
CLASS Sklero (REF="0") / PARAM=GLM;
MODEL tNotBld*bld(0) = sklero sklero*log2bilirub log2bilirub
/ ENTRY=t_entry ; *- Skal med, hvis der er forsinket indgang -*;
HAZARDRATIO Sklero / CL=PL; *- Effekt af Sklero ved mean af log2bilirub -*;
HAZARDRATIO log2Bilirub / CL=PL;
*- Effekterne af log2(bilirubin) -*;
RUN;
Eksempel på program til ventetidsanalyse
*------------------------------------------------------*;
*- Grafisk tjek for proportionalitet for behandlingen -*;
*------------------------------------------------------*;
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = log2bilirub / RL ENTRY=t_entry;
STRATA sklero;
BASELINE OUT=check LOGLOGS=LCumRate / METHOD=CH;
RUN;
SYMBOL1 C=BLACK V=NONE I=STEPLJ L=1 MODE=INCLUDE;
SYMBOL2 C=RED
V=NONE I=STEPLJ L=1 MODE=INCLUDE
R=1; *- SAS kan huske R=2 fra KM og CH plottene -*;
*- AXIS definitionerne står i starten af programmet -*;
*- Almindelig lineær tidsakse -*;
PROC GPLOT DATA=check;
PLOT LCumRate*tNotBld=sklero / HAXIS=AXIS1 VAXIS=AXIS2 NOLEGEND;
RUN; QUIT;
*- Logaritmisk tidsakse -*;
PROC GPLOT DATA=check;
WHERE 0<tNotBld; *- Bemærk, at tNotBld=0 skal fjernes pga. logaritmen -*;
PLOT LCumRate*tNotBld=Sklero / HAXIS=AXIS11 VAXIS=AXIS2 NOLEGEND;
GOPTIONS DEVICE=; *- SKAL med for SAS Enterprise, ellers helt unødvendig! -*;
RUN; QUIT;
111
Eksempel på program til ventetidsanalyse
112
*---------------------------------------------------*;
*- Grafisk tjek for proportionalitet for bilirubin -*;
*---------------------------------------------------*;
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = log2bilirub sklero / RL ENTRY=t_entry;
STRATA bilirub (36, 63); *- opdeling v 36 og 63 (33.3 og 66.7 pctl blandt reblødere) -*;
BASELINE OUT=check2 LOGLOGS=LCumRate / METHOD=CH;
RUN;
SYMBOL1
SYMBOL2
SYMBOL3
*- AXIS
C=BLACK V=NONE
C=RED
V=NONE
C=BLUE V=NONE
definitionerne
I=STEPLJ L=1 MODE=INCLUDE;
I=STEPLJ L=2 MODE=INCLUDE;
I=STEPLJ L=41 MODE=INCLUDE;
står i starten af programmet -*;
*- Almindelig lineær tidsakse -*;
PROC GPLOT DATA=check2;
PLOT LCumRate*tNotBld=bilirub / HAXIS=AXIS1 VAXIS=AXIS2 LEGEND=LEGEND1;
RUN; QUIT;
*- Logaritmisk tidsakse -*;
PROC GPLOT DATA=check2;
WHERE 0<tNotBld; *- Bemærk, at tNotBld=0 skal fjernes pga. logaritmen -*;
PLOT LCumRate*tNotBld=Sklero / HAXIS=AXIS11 VAXIS=AXIS2 LEGEND=LEGEND1;
GOPTIONS DEVICE=; *- SKAL med for SAS Enterprise, ellers helt unødvendig! -*;
RUN; QUIT;
Eksempel på program til ventetidsanalyse
*------------------------------------------------------*;
*- Numerisk test af proportionalitet for behandlingen -*;
*------------------------------------------------------*;
*- De kumulerede rater kan approksimeres nogenlunde ved rette linier i
*- intervallerne 0-14 dage, 14-105 dage og over 105 dage. Hvis raterne
*- er (nogenlunde) konstante i de intervaller, er de automatisk også
*- (nogenlunde) proportionale inden for hvert af disse intervaller!
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = log2bilirub Skl_0_14 Skl_14_105 SklFra105
/
ENTRY=t_entry /* Skal med, hvis der er forsinket indgang */
RL
;
IF
tNotBld<=14 THEN Skl_0_14 =
sklero; ELSE Skl_0_14 =
0;
IF 14<tNotBld<=105 THEN Skl_14_105 = sklero; ELSE Skl_14_105 = 0;
IF 105<tNotBld
THEN SklFra105 = sklero; ELSE SklFra105 = 0;
TestProp: TEST Skl_0_14 = Skl_14_105 = SklFra105;
RUN;
113
-*;
-*;
-*;
-*;
Eksempel på program til ventetidsanalyse
114
*---------------------------------------------------------*;
*- Numerisk test af proportionalitet for serum bilirubin -*;
*---------------------------------------------------------*;
*- I stedet for at kigge efter stykkevis linearitet kan man vælge at -*;
*- opdele i tidsintervaller med cirka lige mange events. Fordelen ved -*;
*- denne tilgang er, at man får cirka lige meget statistisk styrke i -*;
*- hvert interval (jf. styrkeberegningen for ventetidsanalyser).
-*;
*- Tertiler for event-tidspunkter -*;
PROC UNIVARIATE DATA=skl PCTLDEF=3 NOPRINT; *- NOPRINT: uden output -*;
WHERE bld=1;
VAR tNotBld;
OUTPUT
OUT=TidsTertil
PCTLPRE=Tid
PCTLPTS=33.3 66.7;
RUN; *- Kig på WORK.TidsTertil og vælg tider tæt på tertilerne, her 7 og 49 -*;
PROC PHREG DATA=skl;
MODEL tNotBld*bld(0) = sklero Bili_0_7 Bili_7_49 BiliFra49
/ RL
ENTRY=t_entry; *- ENTRY= hvis der er forsinket indgang -*;
IF
tNotBld<=7 THEN Bili_0_7 = log2bilirub; ELSE Bili_0_7 = 0;
IF 7<tNotBld<=49 THEN Bili_7_49 = log2bilirub; ELSE Bili_7_49 = 0;
IF 49<tNotBld
THEN BiliFra49 = log2bilirub; ELSE BiliFra49 = 0;
TestProp: TEST Bili_0_7 = Bili_7_49 = BiliFra49;
RUN;
Eksempel på program til ventetidsanalyse
115
*- Alternativt: Læg tidstertilerne (variablene Tid33_3 og Tid66_7) ind på -*;
*- analysedatasættet, og lad SAS tage sig af resten "helt automatisk".
-*;
DATA skl_m_TidsTertil;
SET skl;
IF _N_=1 THEN SET TidsTertil;
RUN;
PROC PHREG DATA=skl_m_TidsTertil;
MODEL tNotBld*bld(0) = sklero BiliTidligt BiliMellem BiliSent
/ RL
ENTRY=t_entry; *- ENTRY= hvis der er forsinket indgang -*;
IF
tNotBld<=Tid33_3 THEN BiliTidligt= log2bilirub; ELSE BiliTidligt= 0;
IF Tid33_3<tNotBld<=Tid66_7 THEN BiliMellem= log2bilirub; ELSE BiliMellem= 0;
IF Tid66_7<tNotBld
THEN BiliSent=
log2bilirub; ELSE BiliSent=
0;
TestProp: TEST BiliTidligt = BiliMellem = BiliSent;
RUN;