Analyse af tid-til-event data
i store kohorte studier
analyseret med SAS/STAT
Jacob Simonsen,
Afdeling for Epidemiologisk Forskning
Statens Serum Institut
Register baseret forskning
Kendetegn:
• Baseret på landsdækkende registre, derfor store datasæt
• Udfald er oftest binære eller tid-til-event
• Prædefinerede og veldefinerede statistiske hypoteser
Centralpersonregistret
Psykiatriskregistret
Analyse data
Landspatientregistret
Lægemiddelregistret
Patalogiregistret
Sygesikringsregistret
Lægemiddelregistret
Cancerregistret
Tid til event modeller (overlevelses analyse)
I disse analyser gælder det om at finde sammenhænge mellem
forklarende faktorer og tid-til-event.
Tiderne kan være censoreret, dvs for nogle observationer kendes kun
en nedre grænse. Kan skyldes studietids ophør.
Og observationer kan være venstre trunkeret, dvs med forsinket entry.
Tid
Procedurer specialiseret til overlevelses analyse
Rate modeller: l(t)=f(Xb)
PROC PHREG
PROC ICPHREG
PROC SURVEYPHREG
Accelerated Failure time models: E(T)=f(Xb)
PROC LIFEREG
Ikke-parametriske analyser: Kaplan Meier kurver, Log rank test, etc
PROC LIFETEST
PROC ICLIFETEST
Percentil regression: fx Median(T)=Xb
PROC QUANTLIFE
Generelle procedure der også finder
anvendelse i overlevelses analyser
Generaliserede lineære modeller
PROC GENMOD
PROC GLIMMIX
PROC HPGENSELECT
Ikke lineære modeller
PROC HPNLMOD
PROC NLMIXED
Rate modeller
Her konstrueres modeller for raten (groft sagt:
1
rate=
)
𝑡𝑖𝑑−𝑡𝑖𝑙−𝑒𝑣𝑒𝑛𝑡
Mest populære rate modelle er
• Poisson regression (fuld parametrisk model):
l(t)=eX(t)b ,
X(t) er stykvise konstant
• Cox regression (semi parametrisk model):
l(t)=eX(t)b l0(t),
hvor l0(t) er en baggrundsrate
Mindre populært, men langt mere generel er
• Lin-Yings model:
l(t)=g(X(t)b) + h(Z(t)g) l0(t)
Poisson regression
Man kan nøjes med at regne på aggregerede persontider og events i kombinationer af forklarende faktorer:
Data på individ niveau
Aggregeret data
Alder
Køn
Eksponeret
0
20-24
Mand
Ja
7. Dec 85’
0
20-24
Mand
Nej
7. Dec 85’
7. Dec 86’
1
20-24
Kvinde
Ja
…
…
…
…
20-24
Kvinde
Nej
…
…
…
…
…
…
…
Person
ID
Køn
Eksponeret
Alder
Start
Slut
Event
1
Mand
Nej
20-24
1. Jan 78’
1. Jan 83’
1
Mand
Nej
24-25
1. Jan 83’
1
Man
Ja
24-25
…
…
…
…
…
…
Personår
Events
…
…
Med flere tidsafhænge faktorer er det et besværligt tælle-arbejde at danne det aggregerede data…
Heldigvis findes %stratify-makroen1 der udfører arbejdet.
1) Rostgaard K. Methods for stratification of person-time and events – a prerequisite for Poisson regression and SIR estimation. Epidemiologic Perspectives & Innovations : EP+I. 2008;5:7.
doi:10.1186/1742-5573-5-7.
Likelihood-funktionen får en form magen til Poisson fordelte data’s likelihood funktion.
Derfor kan procedurer til generaliserede lineære modeller benyttes.
PROC HPGENSELECT data=aggregeret;
CLASS alder koen eksponering;
MODEL events=alder koen eksponering/DIST=POISSON LINK=LOG offset=logpersontid;
RUN;
PROC HPGENSELECT er en multi-kerne procedure. Hurtig på store datasæt med mange forklarende variable.
Til gengæld har PROC GENMOD flere features end HPGENSELECT (fx REPEATED der tillader stokastisk afhængige
observationer).
Beregningstiden er proportional med antal kombinationer af forklarende variable.
Cox-regression
Inden analyse skal data opstilles på formen
Person
ID
Eksponeret
Start
(alder)
Slut
(alder)
event
1
Nej
54
56
0
1
Nej
56
57
0
1
Ja
57
57½
1
…
…
…
…
…
…
…
…
…
…
Bemærk, der er ikke en alder-klasse-variabel. Til gengæld siger start og slut
hvornår personen indgår i analysen.
Man kan selvfølgelig vælge en andent tidsakse en alder.
Data kan nu analyseres med PROC PHREG ...
Estimering i Cox-regression med PROC PHREG
PHREG anno 2002:
PROC PHREG data=mydata;
MODEL (entry exit)*event(0)=kon medicin;
RUN;
PHREG anno 2015:
Optimering af hastighed
Ved venstre trunkerede data
SAS/STAT 14.1 (9.4 M3)
PROC PHREG data=mydata fast;
CLASS kon medicin(REF=‘placebo’) site;
MODEL (entry exit)*event(0)=kon medicin/eventcode=1;
RANDOM site;
ASSESS PROPORTIONALHAZARDS;
WEIGHT myweight;
RUN;
Multiple event typer.
SAS/STAT 13.1
(SAS 9.4 M1)Class statement
(fra version 9.2)
Tillader stokastisk afhængige data
SAS/STAT 9.3
Goodness of fit
SAS 9.2
Vægte – smart ved aggregering
SAS/STAT 13.2 (9.4 M2)
Aggregering i Cox-regression
I Cox regression er det tilstrækkeligt at vide hvor mange der er i risiko og antal døde på hvert riskset
riskset
Tid
Der findes SAS-makroer til dannelse af det aggregerede datasæt.
Eksponeret
Riskset
Tid
(målt i alder)
dummytid
Antal
Nej
54
1
1
Nej
54
2
9
Ja
54
2
10
Nej
57
2
9
Ja
57
1
1
Ja
57
2
9
%coxaggregate(data=,entry=,exit=,covariate=..)
PROC PHREG data=aggregated;
CLASS eksponeret;
MODEL dummytime*dummytime(2)= eksponeret;
WEIGHT antal;
STRATA riskset;
RUN;
Cox regression med tidsafhængie faktorer
Uden aggregering:
Her tillades eksponerings-effekteen at ændre værdi ved fx tid=10. Eksponerings effekten kan ”splittes” således:
PROC PHREG data=ikke_aggregeret_data;
tidlig=(eksponeret=”yes”)*(tid<=10);
sen=(eksponeret=”yes”)*(tid>10);
MODEL tid*event(0)=tidlig sen;
RUN;
Dette give desværre en beregningstid der er kvadratisk med antal individer.
λ
? aggregering
Med
0.8
Aggregeringen
gør det muligt at definere en tidsmarkør i et datastep. Herefter modelleres tidsafhængigheden som en interaktion:
%coxaggregate(data=,entry=,exit=,covariate=..)
0.7
Data aggregated;
0.6
set aggregated;
if (riskset_tid<=10) then tid2=”tidlig”;
0.5
else tid2=”sen”;
Run;
0.4
PROC
PHREG data=aggregated;
0.3
CLASS eksponeret tid2/param=glm;
MODEL dummytime*dummytime(2)= eksponeret*tid2;
0.2
WEIGHT antal;
STRATA riskset;
0.1
hazardratio eksponeret/at(tid2=all);
0.0
RUN;
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Dette reducere beregningstiden til at
være lineær afhængig af antal individer.
Eksponeret
Ueksponeret
17
18
19
20
Tid
Beregningstid i Cox-regression
n=antal individer.
Venstre-trunkering
Tidsafhængige faktorer
Beregningstid
uden aggregering
Beregningstid
Med aggregering
Nej
Nej
O(n)
O(n)
Nej
Ja
O(n2)
O(n)
Ja
Nej
O(n2)
eller O(n) med FAST option
O(n)
Ja
Ja
O(n2)
O(n)
Lin Yings Ikke lineære model
Lad os forestille os at raten formen
l(t)=g(Xb) + h(Xa) l0(t)
•
•
•
•
b og a er parametre, der skal estimeres.
X og Z er vektorer med forklarende variable
l0(t) er en baggrundsrate (parametrisk eller ikke-parametrisk).
g og h er selvvalgte link-funtioner.
En ekstrem fleksibel model!
Modellen tillader additive komponenter, g(Xb), som udtrykker rate-differencer
og multiplikative komponenter, h(Za), som udtrykker rate-ratioer.
Estimering med parametrisk baggrundsrate
Antag modellen
l(t)=g(Xb) + h(Za) l0(t)
Med fx g(x)=x og h(x)=ex, og at baggrundsraten fx har formen:
𝑘 𝑡 𝑘−1
l0(t)=𝑟 𝑟
.
Estimering foretages ved maksimere log-likelihood-funktionen.
Hver observation har et bidrag på formen
li(b,g,s)=1event log(li(t))som er et rent parametrisk udtryk.
𝑡𝑒𝑥𝑖𝑡
l (s)
𝑡𝑒𝑛𝑡𝑟𝑦 i
𝑑𝑠,
Bidraget til likelihood-funktionen kan programmeres ind i HPNLMOD:
PROC HPNLMOD DATA=simulation ;
lambda=beta1*w1+exp(gamma1*x1)*(k/r)*(t/r)**(k-1)
integrale=t*beta1*w1+exp(gamma1*x1)*(t/r)**k
likelihood=event*log(lambda)-integrale;
BOUNDS k r>0;
MODEL t ~ GENERAL (likelihood);
RUN;
Chancerne for konvergens øges hvis man har et fornuftigt bud på startværdier
parameters beta1 0.25 gamma1 0.5
k 1 r 1;
Har man tilfældige effekter i modellen kan man benytte sig af PROC NLMIXED.
Estimering med ikke-parametrisk baggrundsrate
li(t)=g(Xib) + h(Zia) l0(t)
Nu antages ikke nogen parametrisk form på l0(t).
Modellen er derfor mere fleksibel på trods af færre parametre.
Parametrene kan estimeres ved at løse et sæt estimerings ligninger
(Se Lin-Yings artikel for den præcise form1):
𝐹 b, g =
𝑡
0
𝑓 a, b, 𝑑𝑎𝑡𝑎 𝑠 𝑑𝑀𝑖 (𝑠) = 𝟎 .
Generelt kan dette ikke opstilles som et maksimerings problem.
- PROC HPNLMOD eller nogen anden procedure virker derfor ikke hér.
Heldigvis kan alt løses i SAS/BASE (kræver dog lidt arbejde )
1Lin,
D. Y.; Ying, Zhiliang. Semiparametric Analysis of General Additive-Multiplicative Hazard Models for Counting Processes. Ann. Statist. 23 (1995), no. 5, 1712--1734. doi:10.1214/aos/1176324320.
http://projecteuclid.org/euclid.aos/1176324320.
Numerisk løsning til estimerings ligningen
𝐹 b, g =
𝑡
0
𝑓 a, b, 𝑑𝑎𝑡𝑎 𝑠 𝑑𝑀𝑖 (𝑠) = 𝟎
Kan løses med Newton Raphsons algoritme
−1
𝛿𝐹
𝜽𝑖 = 𝜽𝑖−1 −
(𝜽𝑖−1 )
𝐹 𝜽𝑖−1 ,
𝛿𝜽
Hvor ”𝜽” er hele parameter vektoren 𝜽 = a b .
Dvs, vælg passende startværdi 𝜽0 og iterér indtil der er konvergens.
Sådan kan en Newton-Raphson algoritme implementeres i et datastep
Proc sort data=mydata;
by t;
Run;
Parametre og afledte erklæres i 2x2 arrays;
Data estimater;
/***initialiser parametre****/
array parameters{&dimension.,1} _temporary_;
array afledte{&dimension.,&dimension.} _temporary_;
do until (konvergens=1);
***sæt integrand og afledte til 0****;
do i=1 to nobs;
set mydata point=i;
****opdater hjælpevariable
end;
****Udregn F(q) og dF(q)/dq
****Newton Rapson opdatering af parametre;
*** if (parameterændring<delta) then konvergens=1;***;
end;
keep parametre;
Run;
Newton-Raphson algoritme
Datasættet løbes igennem,
Hvorved integralet og dens afledte bliver løst
−1
𝛿𝐹
𝜃𝑖 = 𝜃𝑖−1 −
𝐹 𝜃
𝛿𝜃
(Mere om matrix-operationer på næste slide)
Algoritmen stoppes når
Der er opnået konvergens
Matrix funktioner
Matrix-funktioner kan ikke kaldes direkte fra et datastep,
men de kan gøres (permanent) tilgængelige via PROC FCMP således..
option cmplib=function.func;
libname function 'd:\sasdata\sasfunctions';
proc fcmp outlib=function.func.matrix;
subroutine invers(m[*,*],inv[*,*]);
outargs m,inv;
call inv(m,inv);
endsub;
run;
Herefter funktionen kaldes fra datastep:
array A{4,4} _temporary_;
array B{4,4} _temporary_;
call invers(A,B);
Tilsvarende med andre matrix funktioner (addér, multiplicér osv.).
PROC FCMP kan også bruges til at definere de to link-funktioner…
proc fcmp outlib=work.function.gh;
function g(x);
y=x;
return (y);
endsub;
function h(x);
y=exp(x);
return (y);
endsub;
run;
På samme måde skal 1. og 2. afledte af g og h defineres.
Hvis estimerings-datasteppet pakkes ind i en makro,
%macro lin_ying(data,model…);
%mend;
så skal brugeren kun definere link-funktioner g og h (og afledte) med
PROC FCMP,
og derefter kalde makroen
%lin_ying(data=simulation,…)
Eksempel
100.000 simulerede levetider, med en rate på formen
l(t)= xb + 𝑒 za l0(t),
hvor baggrundsraten er
𝑘 𝑡 𝑘−1
l0(t)=𝑟 𝑟
,
Link-funktionerne er g(x)=x, og h(x)=ex
z og x er binære (0 eller 1),
b=0.25
a=0.25
(dvs en Weibull form)
Eksempel – fortsat
Sammenligning mellem makro og PROC HPNLMOD
Begge giver gode estimater hvis modellen er specificeret korrekt
Sande
værdi
PROC HPNLMOD
Lin-Ying makro
Parametrisk model
Semi-parametrisk model
b
0.25
0.2513 (0.0019 std err)
0.2515 (0.0022)
a
0.25
0.2541 (0.0107)
0.254 (0.0125)
Beregningstid (real)
5.0 sec
7:08.92
Beregningstid (CPU)
1:01.98
7:09:21
• PROC HPNLMOD er (meget) hurtigere end makroen. Lin-Ying makroen har dog lineær beregningstid.
• HPNLMOD kræver en fuld-parametrisk model. Dvs flere model-antagelser, hvilket kan medføre bias.
• Makroen er mindre følsom over for valg af start værdier, da den semi-parametriske baggrundsrate sikrer pænt fit fra start.
-Hvorimod de potentielt mange parametre i en fuld-parametrisk model gør valg af startværdier svær/umulig.
Slut