analiza Scorului de înclinație
context
când se estimează efectele tratamentului asupra unui rezultat binar în studiile observaționale, se întâmplă adesea ca tratamentele să nu fi fost repartizate aleatoriu subiecților. Dacă, de exemplu, pacienții mai bolnavi au fost deseori repartizați la tratament, în timp ce pacienții mai sănătoși nu au fost adesea tratați, o analiză simplă ar putea estima greșit gradul sau direcția unui efect de tratament.
o modalitate obișnuită de a încerca să se adapteze pentru potențiala părtinire datorată acestui tip de confuzie este prin utilizarea modelelor de regresie logistică multivariabile. O abordare alternativă este utilizarea analizei scorului de înclinație. În secțiunile următoare oferim un mic exemplu de set de date, apoi descriem și ilustrăm aceste metode alternative de analiză statistică. Ne concentrăm pe cel mai simplu exemplu în care pacienții sunt desemnați să primească fie tratament activ, fie control (adică 2 grupuri). La final menționăm pe scurt posibilele extensii la trei sau mai multe grupuri de tratament.
exemplu de date
următorul exemplu de set de date va fi utilizat pentru a ilustra conceptele de bază. Datele includ 400 de subiecți incluși într-un studiu retrospectiv de cohortă la bărbați cu vârste cuprinse între 40 și 70 de ani internați în spital cu suspiciune de infarct miocardic. Rezultatul interesului este mortalitatea de 30 de zile (deces=1). De interes este efectul posibil al administrării rapide a unui medicament mai nou de eliminare a cheagurilor (trt=1) față de o terapie standard (trt=0) asupra riscului de mortalitate. Covariatele relevante sunt un scor preexistent al factorului de risc (pe o scară de la 0 la 5, 5 fiind cel mai rău) și un scor de severitate a admiterii (pe o scară de la 0 la 10, 10 fiind cel mai rău). Iată datele pentru primii 12 subiecți:
age | male | risk | severity | trt | death |
---|---|---|---|---|---|
48 | 1 | 3 | 8 | 0 | 0 |
59 | 1 | 4 | 6 | 1 | 0 |
67 | 1 | 3 | 6 | 0 | 1 |
51 | 1 | 0 | 6 | 0 | 0 |
56 | 1 | 1 | 6 | 1 | 0 |
60 | 1 | 1 | 6 | 0 | 0 |
53 | 1 | 0 | 3 | 1 | 0 |
54 | 1 | 1 | 2 | 0 | 0 |
54 | 1 | 2 | 7 | 0 | 0 |
62 | 1 | 0 | 4 | 0 | 0 |
64 | 1 | 2 | 6 | 1 | 1 |
70 | 1 | 3 | 8 | 1 | 0 |
… | … | … | … | … | … |
datele pot fi, de asemenea, descărcate ca propensiune.csv sau numit direct în R folosind comanda:example <- read.csv("http://web.hku.hk/~bcowling/data/propensity.csv", header=TRUE)
un script r însoțitor pentru a rula toate analizele următoare poate fi găsit aici: propension.R.
analiza descriptivă
Un total de 192 (48%) pacienți au primit noul tratament (trt=1). Ratele de mortalitate de 30 de zile pentru pacienții tratați și netratați sunt rezumate în tabelul următor:
rezultat | trt=0 | trt=1 |
---|---|---|
a supraviețuit | 168 | 162 |
a murit | 40 | 30 |
rata mortalității pe 30 de zile | 19% | 16% |
o modalitate de a investiga efectul potențial al tratamentului este cu o estimare a diferenței de risc între cele două grupuri. Riscul relativ de mortalitate asociat tratamentului 1 este de 0, 375 / 0, 40, care este de 0, 81, sugerând un beneficiu ușor pentru terapia mai nouă.
o altă modalitate de a estima efectul tratamentului ar putea fi calcularea raportului de șanse, mai degrabă decât a riscului relativ. Raportul de Cote este (168×30) / (162×40), care este de 0,78, iar un interval de încredere de 95% poate fi calculat ca (0,46, 1,31).o a treia modalitate de a estima efectul tratamentului este de a analiza reducerea absolută a ratelor mortalității. Aici modificarea asociată tratamentului 1 este de -3,6% (de la 19,2% la 15,6%) și un interval de încredere de 95% este (-11.5%, 4,3%), adică o reducere de 12% sau o creștere de 4% a ratelor mortalității.
cu toate acestea, următoarele două parcele arată că subiecții cărora li s-a administrat tratamentul mai nou au fost puțin mai în vârstă decât cei cărora li s-a administrat terapia standard:
investigațiile ulterioare arată că, de asemenea, pare să existe diferențe între factorii de risc și starea actuală:
comparația formală a distribuțiilor acestor variabile explicative între cele două grupuri de tratament arată că diferențele de vârstă (testul t, ppp
regresie logistică
modelele multivariabile sunt adesea folosite pentru a evalua efectul tratamentului ajustarea pentru variabile explicative importante. Ajustarea pentru variabile explicative importante este necesară pentru a asigura comparabilitatea între grupurile de tratament și de control, iar dacă ajustarea nu este efectuată, atunci diferențele dintre grupuri pot duce la estimări părtinitoare ale efectului tratamentului.
tabelul de mai jos prezintă raportul de cote brute de tratament, apoi efectul ajustat pentru alte variabile explicative. Comparația criteriului de informare Akaike pentru fiecare model sugerează că scorul de risc și indicele de severitate nu îmbunătățesc semnificativ potrivirea, adică modelul 2 poate fi preferat față de modelul 3. Există o sugestie a unui beneficiu al tratamentului (deși nu este semnificativ statistic) și, de asemenea, o confuzie aparentă în funcție de vârstă, așa cum este suspectat de analizele descriptive de mai sus.
estimările dintr-un model care presupune efecte liniare ale covariatelor sunt foarte asemănătoare cu modelul 3 (rezultatele nu sunt prezentate).
Factor | n | Model 1 | Model 2 | Model 3 | |||
---|---|---|---|---|---|---|---|
95% iî | 95% IÎ | 95% IÎ | |||||
tratament 0 | 208 | 1,00 | – | 1,00 | 1,00 | – | |
tratament 1 | 192 | 0,78 | (0,46, 1,31) | 0.67 | (0.39, 1.15) | 0.62 | (0.35, 1.11) |
Age 40-49 | 95 | – | – | 1.00 | – | – | – |
Age 50-59 | 131 | – | – | 1.72 | (0.77, 3.82) | 1.26 | (0.52, 3.01) |
Age 60-70 | 175 | – | – | 2.62 | (1.23, 5.62) | 2.03 | (0.84, 4.95) |
Risk score 0 | 112 | – | – | – | – | 1.00 | – |
Risk score 1 | 103 | – | – | – | – | 3.06 | (1.34, 6.97) |
Risk score 2-3 | 132 | – | – | – | – | 1.33 | (0.54, 3.28) |
Risk score 4-5 | 53 | – | – | – | – | 2.64 | (0.95, 7.35) |
Severity index 0-3 | 108 | – | – | – | – | 1.00 | – |
Severity index 4 | 69 | – | – | – | – | 1.29 | (0.56, 2.96) |
Severity index 5 | 80 | – | – | – | – | 0.78 | (0.33, 1.87) |
Severity index 6 | 56 | – | – | – | – | 1.28 | (0.53, 3.08) |
Severity index 7-10 | 87 | – | – | – | – | 1.43 | (0.65, 3.16) |
AIC | 374 | 371 | 371 |
For completeness we could also use a non-linear regression model to check the shape of the effects of age, pre-existing risk and severity in the fully adjusted model. Scriptul r însoțitor conține codul pentru estimarea și trasarea funcțiilor spline corespunzătoare și nu le arătăm aici; observăm că efectele au fost destul de liniare.
propension score analysis
o abordare alternativă a analizei este de a încerca să imite condițiile unui studiu randomizat controlat (RCT). Într-un RCT, probabilitatea ca un participant să primească un anumit tratament este aceeași pentru toți participanții sau, în modele stratificate, depinde doar de variabilele explicative cunoscute ale unui pacient, cum ar fi vârsta, sexul etc. Cu alte cuvinte, vârsta, sexul pacientului (etc.) sunt suficiente informații pentru a ne spune probabilitatea pacientului de a primi tratamentul.
dacă într-un cadru de studiu observațional am avea toate informațiile disponibile profesioniștilor din domeniul sănătății care au atribuit tratament subiecților, ar trebui să putem recrea procesul decizional al acestora și să estimăm probabilitatea ca pacienții individuali să primească tratamentul. Această probabilitate este denumită scor de înclinație, iar în lucrarea lor seminală din 1983 Rosenbaum și Rubin au arătat că, atâta timp cât scorul de înclinație este o măsură adecvată a probabilității de a primi tratament, scorurile pot fi utilizate pentru a ajuta la estimarea efectelor cauzale ale tratamentului. Scorurile sunt utilizate pentru a echilibra variabilele de prognostic între grupurile tratate și netratate și există (cel puțin) patru modalități posibile de a face acest lucru:
- stratificați pacienții în grupuri (de exemplu, chintile) prin scorul de înclinație și comparați efectele tratamentului în fiecare strat.
- potriviți pacienții tratați și netratați și comparați perechile potrivite rezultate.
- ponderarea inversă a rezultatelor prin scorul de înclinație.
- ajustați pentru Scorul de propensiune într-un model de regresie logistică.
5.1 estimarea scorului de propensiune
scorul de propensiune este probabilitatea condiționată ca un subiect să fie tratat având în vedere variabilele explicative observate; intenția este ca această probabilitate unică să poată rezuma informațiile despre mecanismul de atribuire a tratamentului. Ar trebui apoi să putem obține estimări imparțiale ale efectelor tratamentului prin compararea subiecților care au avut probabilități similare de a primi un tratament (indiferent dacă l-au primit sau nu).
scorurile de propensiune sunt de obicei estimate folosind un model de regresie logistică multivariabilă.
în exemplul nostru, am montat un model de regresie logistică pentru a estima efectele vârstei, scorului de risc și indicelui de severitate asupra probabilității de a primi tratamentul 1, mai degrabă decât tratamentul 0. Constatăm că vârsta mai înaintată (p=0,05), scorul de risc mai mare (p=0,05) și indicele de severitate mai mare (p=0.01) sunt toate asociate cu o probabilitate mai mare de a primi tratament 1. Scorurile de propensiune variază de la 0,2 la 0,8 și comparăm distribuția scorurilor între cele două grupuri de tratament în figura de mai jos. Barele arată intervalul median și inter-quartile.
așa cum ar fi de așteptat, scorurile de înclinație (adică probabilitățile de a primi tratament) sunt în medie ușor mai mari în grupul de tratament. Putem vedea că există un grad bun de suprapunere, unde putem găsi indivizi în ambele grupuri de tratament pentru orice scor de înclinație între 0,2 și 0,8. Acest lucru este important, deoarece principiul esențial al analizei scorului de înclinație este că dacă găsim doi indivizi, unul în fiecare grup de tratament, ne putem imagina că acei doi indivizi au fost repartizați aleatoriu fiecărui grup în sensul că oricare dintre alocări este la fel de probabilă.
5.2 scorul propensiunii echilibrează grupurile?
în orice analiză a scorului de înclinație ar trebui să verificăm dacă propensityscore ne permite să echilibrăm distribuția variabilelor explicative. Există multe modalități de a verifica echilibrul ; de exemplu, am putea analiza distribuția unei variabile explicative în chintilele scorului de înclinație. În figura de mai jos se trasează intervalele de vârstă mediană și interquartilă în fiecare chintilă cu scor de înclinație:
fără ajustare (în general) există o discrepanță considerabilă. Cu toate acestea, în cadrul fiecărei quntile, distribuțiile sunt foarte strâns aliniate.
putem cuantifica diferențele inițiale prin calcularea statisticilor cu două eșantioane (adică un test t pentru diferențele de vârstă pe grupe de tratament). Acest lucru este echivalent cu găsirea statisticii t pentru tratament dintr-un model de regresie liniară (sau ANOVA) pentru vârstă față de grupul de tratament. Putem măsura în continuare diferențele după ajustarea scorului de propensiune, calculând Statisticile t pentru tratament dintr-un model de regresie liniară multivariabilă (sau ANOVA) pentru ajustarea vârstei pentru tratament, precum și ajustarea pentru chintilele scorului de propensiune. Statisticile t neajustate (cercuri umplute) și ajustate (cercuri deschise) sunt prezentate în figura de mai jos:
putem vedea că ajustarea scorului de înclinație elimină aproape toate diferențele inițiale de vârstă, scor de risc și indice de severitate între cele două grupuri de tratament.
5.3 ratele mortalității în cadrul scorului de propensitate chintile
am constatat anterior că covariatele sunt echilibrate în cadrul chintilelor scorului de propensitate. Rosenbaum și Rubin au arătat că efectul mediu al tratamentului în straturile scorului de înclinație este o estimare imparțială a efectului adevărat al tratamentului (cu condiția ca unele ipoteze să fie valabile). Am trasat ratele de mortalitate de 30 de zile (cu intervale de încredere de 95%) în funcție de grupul de tratament în fiecare chintilă cu scor de propensitate, mai jos:
ratele de mortalitate au fost, în general, mai mici în grupul căruia i s-a administrat tratamentul 1 (Albastru) decât în grupul căruia i s-a administrat tratamentul 0 (roșu), cu excepția T3 unde ratele au fost similare. Cu toate acestea, nu există dovezi puternice că efectele tratamentului variază în funcție de intervalul scorurilor de propensiune.
putem calcula diferența ratelor de mortalitate între grupurile de tratament în fiecare chintilă și putem obține efectul mediu al tratamentului ca medie ponderată între chintile. Figura de mai jos prezintă reducerea absolută a ratei mortalității pentru tratamentul 1 față de tratamentul 0 și media ponderată, cu intervale de încredere de 95% :
per total a existat o reducere absolută de 6% a ratei mortalității de 30 de zile pentru tratamentul 1 comparativ cu tratamentul 0, cu un interval de încredere destul de larg.
5.4 ratele de mortalitate între perechi potrivite de indivizi
o abordare alternativă este de a găsi perechi de subiecți, câte unul în fiecare grup de tratament, cu scoruri de înclinație foarte similare. Prin definirea scorului de înclinație, doi subiecți cu scoruri de înclinație similare ar trebui să fie, de asemenea, similari pe toate covariatele importante. Această procedură de potrivire este mai simplă din punct de vedere computațional decât potrivirea simultană pe toate covariatele importante.
folosind un algoritm de potrivire în datele exemplului, găsim 177 de perechi potrivite (adică 354 de indivizi) din cei 400 de subiecți originali. Am putut verifica dacă algoritmul de potrivire a atins echilibrul între grupuri prin compararea distribuțiilor covariatelor între cele două grupuri de tratament, între perechile potrivite. În subsetul corespunzător, au existat 23 de decese în grupul cu tratament 1 și 36 de decese în grupul cu tratament 0, ceea ce reprezintă o reducere absolută semnificativă statistic de 7,8% (interval de încredere 95%: -13,7%, -1,8%).
5.5 ponderarea inversă prin scorurile de înclinație
Rosenbaum descrie o utilizare alternativă a scorului de înclinație, ca factor de ponderare. Fără a intra în detalii despre derivare, el arată că rata mortalității așteptate dacă toți subiecții au fost alocați grupului de tratament 1 în loc de grupului 0 este egală cu E(YT/p), unde Y este variabila de rezultat, T este grupul de tratament și p este scorul de înclinație de a fi atribuit grupului de tratament 1. În mod similar, rata mortalității așteptate dacă toți indivizii sunt repartizați în grupul de tratament 0 este dată de E(Y(1-T)/(1-p)). Efectul cauzal mediu este apoi diferența dintre aceste două rate de mortalitate așteptate.
folosind scorurile de propensiune ca ponderi, am estimat că tratamentul 1 a fost asociat cu o reducere absolută de 6,5% (interval de încredere 95%: -13,9%, 1,8%) față de tratamentul 0.
5.6 ajustarea regresiei logistice pentru Scorul de propensiune
am estimat efectul tratamentului 1 vs tratamentul 0 într-un model de regresie logistică ajustarea pentru Scorul de propensiune (în chintile). Rata șanselor pentru tratamentul 1 a fost estimată la 0,65 (interval de încredere 95%: 0,37, 1,13). Am găsit un raport de cote estimat similar atunci când am adăugat variabilele explicative originale la model (adică. ajustat pentru Scorul de înclinație, vârstă, risc și severitate).
Rezumatul rezultatelor
ratele de mortalitate observate la 30 de zile au fost de 19% în grupul căruia i s-a administrat tratamentul 0 și de 16% în grupul căruia i s-a administrat tratamentul 1. O comparație a estimărilor din diferitele metode statistice este prezentată în tabelul de mai jos.
abordare | Diferență absolută | Rata cotelor | ||
---|---|---|---|---|
estimare | 95% IÎ | estimare | IÎ 95% | |
nicio ajustare | -3,6% | (-11.5%, 4.3%) | 0,78 | (0,46, 1,31) |
regresie logistică ajustarea pentru vârstă, scor de risc, și indicele de severitate | – | – | 0,62 | (0,35, 1.11) |
Stratifying by PS | -6.0% | (-25.8%, 13.7%) | – | – |
Matching by PS | -7.8% | (-13.7%, -1.8%) | 0.58 | (0.33, 1.04) |
Weighting by PS | -6.5% | (-13.9%, 1.8%) | 0.63 | (0.34, 1.11) |
Logistic regression adjusting for PS | – | – | 0.65 | (0.37, 1.13) |
în general, metodele scorului de propensiune dau rezultate similare modelului de regresie logistică. Aceasta este o constatare bine cunoscută din studiile empirice și de simulare anterioare .
notați ușoara discrepanță în semnificația statistică pentru metoda de potrivire, unde intervalul de încredere de 95% pentru raportul de Cote a fost calculat prin aproximarea standard și poate fi prea mare.
discuție
în secțiunile de mai sus, a fost descrisă și ilustrată utilizarea ajustării regresiei și a scorurilor de înclinație pentru analiza datelor observaționale. Este important de remarcat limitarea inevitabilă a datelor observaționale privind efectele tratamentului în comparație cu datele dintr-un studiu randomizat. Adică, metodele bazate pe ajustarea regresiei sau scorurile de înclinație din datele observaționale permit doar echilibrarea analizei peste covariatele cunoscute, în timp ce randomizarea echilibrează peste covariatele cunoscute și necunoscute.
atunci când se utilizează analiza scorului de înclinație, este vital să se verifice dacă factorii prognostici importanți sunt echilibrați de Scorul de înclinație – fără echilibru, teoria de bază eșuează. Cu toate acestea, dacă există un număr mare de predictori, este posibil să nu fie rezonabil să ne așteptăm la un echilibru perfect pentru fiecare, în același mod în care într-un RCT o comparație a factorilor de bază va găsi ocazional diferențe între grupuri din întâmplare.deoarece scorurile de propensiune trebuie să echilibreze distribuția variabilelor explicative între grupuri, uneori modelul va trebui să includă nu numai efecte principale, ci și termeni de interacțiune între variabilele explicative. Din fericire, modelul care este folosit pentru a estima scorurile de înclinație nu este de obicei centrul atenției și, prin urmare, nu trebuie să fie parsimonios – trebuie doar să permită echilibrul. Austin și colab. a efectuat un studiu amplu de simulare și a arătat că cele mai importante variabile care trebuie incluse într-un model de scor de înclinație (și pentru a asigura echilibrul) sunt acele variabile explicative asociate cu rezultatul interesului. Pe de altă parte, nu este esențial să se includă variabile care sunt asociate cu atribuirea tratamentului, dar care nu sunt asociate cu rezultatul.
o situație deosebit de potrivită pentru o abordare a scorului de înclinație este atunci când rezultatul interesului este rar, dar tratamentul este comun . În această situație, este posibil să nu existe prea multe date pentru a modela relația dintre rezultat și variabilele prognostice – o regulă comună este că trebuie observate 10 evenimente pentru fiecare variabilă prognostică (nivelul a) inclusă într – un model de regresie logistică multivariabilă) – în timp ce pot exista suficiente date pentru a construi un model bun pentru Scorul de înclinație. În acest caz, ajustarea folosind scorul de înclinație poate fi singura abordare viabilă a analizei.
un avantaj potențial al metodelor scorului de înclinație față de ajustarea regresiei este că poate fi mai ușor să se verifice dacă scorul de înclinație are variabile măsurate echilibrate între subiecții tratați și cei netratați, în timp ce este mai dificil să se judece dacă un model de regresie a fost specificat corect .
în cele din urmă, este important să rețineți că analiza scorului de înclinație estimează un efect de tratament diferit de ajustarea regresiei. Analiza scorului de înclinație estimează efectul marginal, în timp ce ajustarea regresiei estimează efectul condițional . Efectul marginal al tratamentului este interpretat la nivelul populației: cum s-ar schimba tratamentulnumărul total de rezultate observate în populație? Atunci când se utilizează un model de regresie logistică, efectul tratamentului condițional este modificarea șanselor rezultatului pentru un individ atunci când este expus la tratament în comparație cu primirea unui tratament, condiționată de variabilele explicative ale individului respectiv – adică efectul condițional este interpretat la nivel individual. Un exemplu numeric al acestui efect este dat în tabelul următor, unde o boală afectează 13.200 de indivizi. Majoritatea persoanelor sunt considerate ‘risc scăzut’, în timp ce un număr mic sunt’ risc ridicat’, cu rate de mortalitate de 5% și, respectiv, 25%, sub vechiul tratament. Un nou tratament extrem de eficient va reduce șansele de deces cu 80% (raportul cote condiționale este de 0,2), dar raportul cote la nivelul populației nu este de 0,2:
grup de risc | n | tratament vechi | tratament nou | rel. Risk | Odds Ratio | ||
---|---|---|---|---|---|---|---|
High risk | 1200 | 300 | (25%) | 75 | (6.25%) | 0.250 | 0.200 |
Low risk | 12000 | 600 | (5%) | 125 | (1.04%) | 0.208 | 0.200 |
Total | 13200 | 900 | (6.8%) | 200 | (1.52%) | 0.222 | 0.210 |
- Rosenbaum PR, Rubin DB. Rolul central al scorului de înclinație în studiile observaționale pentru efectele cauzale. Biometrika, 1983; 70:41-55. .
- Baser O. prea mult zgomot despre modelele scor propensiune? Compararea metodelor de potrivire a Scorului de înclinație. Valoare în sănătate, 2006;9(6):377-85..
- Rosenbaum PR. Ajustare directă bazată pe Model. Jurnalul Asociației Americane de Statistică, 1987; 82:387-94. .
- Shah BR, Laupacis A, Hux JE, Austin PC. Metodele scorului de înclinație au dat rezultate similare modelării tradiționale de regresie în studiile observaționale: o revizuire sistematică. Jurnalul de epidemiologie clinică, 2005;58(6):550-9..
- Austin PC, Grootendorst P, Anderson GM. O comparație a capacității diferitelor modele de scor de înclinație de a echilibra variabilele măsurate între subiecții tratați și netratați: un studiu Monte Carlo. Statistici în medicină, 2007; 26(4):734-53..
- Braitman LE, Rosenbaum PR. Rezultate Rare, tratamente comune: strategii analitice folosind scoruri de înclinație . Analele medicinei interne, 2002; 137:693-5. .
- Wang J, Donnan pt. Metode de scor de înclinație în studiile de siguranță a medicamentelor: practică, puncte forte și limitări. Farmacoepidemiologie și siguranța medicamentelor, 2001; 10(4):341-4. .
- Austin PC, Grootendorst P, Normand SL, Anderson GM. Condiționarea scorului de înclinație poate duce la o estimare părtinitoare a măsurilor comune ale efectului tratamentului: un studiu Monte Carlo. Statistici în medicină, 2007; 26(4):754-68..
Lectură suplimentară
- Austin PC. O evaluare critică a potrivirii scorurilor de înclinație în literatura medicală între 1996 și 2003. Statistici în medicină, 2008 (în presă)..
- D ‘ Agostino RB Jr.metode de scor de înclinație pentru reducerea părtinirii în comparația unui tratament cu un grup de control non-randomizat. Statisticeîn Medicină, 1998; 17(19):2265-81..
- Imbens GW. Rolul scorului de propensiune în estimarea funcțiilor doză-răspuns. Biometrika, 2000; 87(3):706-10..
- Rosenbaum PR, Rubin DB. Reducerea prejudecății în studiile observaționale folosind subclasificarea pe scorul de înclinație. Jurnalul Asociației Americane de Statistică, 1984; 79(387):516-24..
- Winkelmayer WC, Kurth T. Scoruri de înclinație: ajutor sau hype?Transplant De Dializă Nefrologică, 2004; 19:1671-3..
mulțumiri
datorită Eric Lau pentru a ajuta la dezvoltarea exemplul ilustrativ.
această lucrare este licențiată sub licența Creative Commons Attribution 3.0 Unported. Acesta a fost scris de Ben Cowling
această pagină a fost modificată ultima dată pe
Leave a Reply