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
© Copyright 2024