Longitudinální a panelová data – NMST422

Letný semester 2025 | Cvičenie 2 | 03.03.2025

Prihlásenie k SAS OnDemand: https://www.sas.com/en_us/software/on-demand-for-academics.html
Nutná je registrácia s vytvorením vlastného účtu s jedinečným identifikačným číslom a potvrdenie registrácia prostredníctvom emailu zadaného pri registrácii. Identifikačné číslo užívateľa (vo forme uXXX, kde XXX je samotné číslo uživateľa) sa objavuje v niektorých následujúcich SAS skriptoch. Príslušný symbol XXX je potrebné vždy nahradiť príslušným identifikačným číslom užívateľa.

Doporučená literatúra a ďalšie užitočné materiály

Diggle, P. J., Heagerty, P., Liang, K. Y., & Zeger, S. (2002). Analysis of longitudinal data. Oxford university press.
Fitzmaurice, G. M., Laird, N. M., & Ware, J. H. (2012). Applied longitudinal analysis. John Wiley & Sons.
Hardin, J.W. and Hilbe, J.M. (2007). Generalized Linear Model and Extensions. StataPress.
Pinheiro, J., & Bates, D. (2006). Mixed-effects models in S and S-PLUS. Springer science & business media.
Jednoduchý (online) SAS tutorial (english)
Základný SAS OnDemand tutorial (english)
Stručný (Český) manuál uživatele SASu na stránke doc. Kulicha
Užitočné aj neužitočné príklady zdrojových kódov v SAS (english)

II. Exploratívna analýza longitudinálnych dat

Podobne ako v prípade štandardných (lineárnych) regresných modelov (resp. akejkoľvek štatistickej analýzy – t.j. aplikácie konkrétneho stochastického/pravdepodobnostného modelu), je štatistická analýza dat rozdelená do dvoch samostatných, ale vzájomne prepojených časti – exploratívnej analýzy a konfirmačnej analýzy (exploratory vs. confirmatory part). Exploratívna analýza využíva rôzne popisné charakteristiky (t.j., empirické odhady teoretických charakteristík) a nástroje pre vizualizáciu dat. Konfirmačná analýza (i.e., štatistická inferencia) hľadá pomocou vhodných pravdepodobnostných metód odpoveď na konkrétnu (expertnú) otázku (t.j., testuje konkrétnu nulovú hypotézu).

Hlavnou výhodou longitudinálnych dat (a modelov určených na ich analýzu) je to, že na rozdiel od (viacmenej jednoduchého/priamočiarého) cross-sekčného porovnania dvoch (rôznych a disjunktných) subpopulácii umožňujú odhadovať a porovnávať (testovať) aj (výrazne komplexnejšiu a netriviálnu) zmenu v rámci konkrétnych subjektov a to vzhľadom k (časovo/priestorovo) sa meniacím sa hodnotám nezávislých premenných popisujúcih daný subjekt. Klasická regresná analýza cross-sekčných dat umožňuje len porovnávať dve podskupiny, ktoré sa vzájomne líšia vzhľadom k hodnotám (uvažovaných) nezávislých premenných. Longitudinálne data a modely umožňujú porovnanie jednak statických (medzi-populačných) cross-sekčných rozdielov a tiež dynamických zmien v rámci konkrétnuch subjektov. Nejedná sa ale o kauzálnu analýzu.

Pre \(N \in \mathbb{N}\) nezavislých subjektov, z ktorých každý subjekt \(i \in \{1, \dots, N\}\) sledujeme opakovane \(n_i \in \mathbb{N}\) krát, môžeme data reprezentovať v tvare \(\{(Y_{ij}, \boldsymbol{X}_{ij}^\top)^\top;~i = 1, \dots, N; j = 1, \dots, n_{i}\}\). V prípade reprezentácie v zmysle \[ \{(\boldsymbol{Y}_{i}^\top, \boldsymbol{X}_{i 1}^\top, \dots, \boldsymbol{X}_{i n_i}^\top)^\top;~i = 1, \dots, N\}, \] kde \(\boldsymbol{Y}_i = (Y_{i1}, \dots, Y_{i n_i})^\top\), môžeme dokonca (pri niektoých dodatočných predpokladoch) hovoriť o náhodnom výbere (nezávisl é a rovnako (združene) rozdelené náhodné vektory). Je nutné si ale uvedomiť, že rozmer (dĺžka) náhodných vektorov \((\boldsymbol{Y}_{i}^\top, \boldsymbol{X}_{i 1}^\top, \dots, \boldsymbol{X}_{i n_i}^\top)^\top\) je pre rôzne \(i \in \{1, \dots, N\}\) rôzna.

Obyčajný linárny regresný (tzv. “cross-sectional”) model môžeme formulovať ako \[ Y_{i 1} = \boldsymbol{x}_{i 1}^\top \boldsymbol{\beta}_{CS} + \varepsilon_{i 1}, \qquad i = 1, \dots, N, \] pre vektor neznámych parametrov \(\boldsymbol{\beta}_{CS} \in \mathbb{R}^p\). Príslušný ohad \(\widehat{\boldsymbol{\beta}}_{CS}\) lze po jednotlivých zložkách interpretovať ako odhadovaný očakávaný rozdiel (i.e., odhad rozdielu v teoretických stredných hodnotách) medzi dvoma podskupinami, ak je rozdiel príslušného regresoru (jeden konkrétny element v \(\boldsymbol{x}_{i 1}\)) jednotkový a ostatné hodnoty regresorov (zvyšných elementov \(\boldsymbol{x}_{i 1}\)) sú rovnaké.

Na druhej strane, jednoduchý model pre longitudinálne data môžeme zapísať aj v tvare \[ Y_{i j} = \boldsymbol{x}_{i 1}^\top \boldsymbol{\beta}_{CS} + (\boldsymbol{x}_{i j} - \boldsymbol{x}_{i 1})^\top \boldsymbol{\beta}_{L} + \varepsilon_{i j}. \qquad i = 1, \dots, N; j = 1, \dots, n_i, \] kde celkový počet pozorovaní je \(\sum_{i = 1}^N n_{i}\). Model ale okrem odhadu pre vektorový parameter \(\boldsymbol{\beta}_{CS}\) (ktorý má ekvivalentnú interpretáciu ako v prípade obýčajného regresného modelu vyššie) umožňuje odhadnúť aj vektorový parameter \(\boldsymbol{\beta}_L\). Interpretácia tohto parametru je ale už specifická pre daný subjekt a vyjadruje zmenu v rámci daného subjektu \(i \in \{1, \dots, N\}\). Príslušný odhad \(\widehat{\boldsymbol{\beta}}_{L}\) môžeme po jednotlivých zložkách interpretovať ako očakávanú zmenu v rámci daného subjektu v čase, ak sa príslušný regresor zmení v čase o jednotku.

Užitočné

Niekoľko jednoduchých doporučení, ako vhodne vizualizovať longitudinálne data:

Obrázok/graf by mal zobrazovať čo najviac dôležitej a hlavne relevantnej informácie a to výrazne nad rámec jednoduchých popisných charakteristík triviálne spočítaných spočítaných z dat;
Obrázok/graf by mal vizualizovať podstatnú informáciu a súvislosti v datach relevantne s expertnou otázkou – t.j., nulovou hypotézou;
Väčšinou je užitočné zobraziť jednak klasické porovnanie v zmysle štandardného lineárneho regresného modelu (cross-sectional), ale tiež longitudinálne profily (resp. zmeny vrámci sledovaných subjektov);
Obrázok/graf by mal byť pomerne jednoduchý, transparentný a ľahko pochopiteľný (obrázok je vždy v texte doplnený o popisok, ktorý dostatočne komplexne daný obrazok popisuje);

Samostatne

Uvedený longitudinálny model vyššie je samozrejeme možné formulovať aj alternatívnymi spôsobmi. Parameter \(\boldsymbol{\beta}_{CS}\) v podstate porovnáva/modeluje disjunktné subpopulácie v čase prvých (t.j., tzv. baselinových pozorovaní) – t.j., pri hodnotách regresorov \(\boldsymbol{x}_{i1}\), kde \(i = 1, \dots, N\). Aké by boli alternatívne formulácie ekvivalentného modelu?

1. Longitudinálne profily v SAS

Pre ilustráciu použijeme datový súbor pig_weights.csv, ktorý udáva váhu u 48 ošípaných, ktoré boli postupne merané po dobu 9 týždňov
(tzv. follow-up period, teda obdobie, počas ktorého boli subjekty opakovane sledované, je 9 týždňov) – avšak nejedná sa nutné o stejných 9 týždňov. Do programu SAS data načítame následujúcim príkazom:

libname weight '/home/uXXX/sasuser.v94';
filename reffile '/home/uXXX/sasuser.v94/data/pig_weights.csv';

proc import datafile=reffile
    dbms=csv
    out=weight.data
    replace;
    getnames=yes;
run;
    
proc print datafile = weight.data; 
run;

Samostatne sa podívajte na niektoré základné popisné charakteristiky spočítané z dat. Pre ilustráciu niektorých popisných charakteristík poslúži napr. aj následujúci graf. Pokúste sa ale vytvoriť aj niekoľko ďalších grafov a obrázkov samostatne.

proc sgscatter data=weight.data;
   matrix weight week /
      ellipse=(type=mean)
      diagonal=(histogram kernel);
run;

V nasledujúcej časti sa špecificky zameriame na vizualizáciu datového súboru vhľadom k jeho časovo-zavislej štruktúre a tzv. subject-specific longitudinálnych profilov.

Porovnajte následujúce grafické výstupy:

proc sgscatter  data=weight.data;
   plot weight * week ;
run;

proc sgscatter  data=weight.data;
   plot weight * week  / group = id;
run;

proc sgplot data=weight.data;
    series y=weight x=week /group=id;
run;

Podívajte sa na základnú syntax procedrúry PROC sgscatter napr. na tejto stránke a analogicky tiež na syntax procedúry PROC sgplot napr. na tejto stránke.

Samostatne

Aké rôzne informácie lze vyčítať z jednotlivých grafov?
Aké informácie, ktoré by nás v súvislosti s datami mohli zaujímať, naopak z grafov nedokážeme získať?
Ktorý graf je pre vizualizáciu longitudinálneho charakteru dat najviac vhodný?

2. Transformácia dat a štandardizované profily

Jednoduché longitudinálne profily nie vždy dobre poslúžia k tomu, aby sme získali dôležité informácie o (časovom) vývoji konkrétneho jednotlivca vzhľadom k celkovej skupine. Špeciálne v prípadoch, keď je uvažovných subjektov príliš mnoho, je niekedy užitočné zvážiť aj graf s tzv. štandardizovanými longitudinálnymi profilmi.

Jedna z možnosti, ako štandardizované longitudinálne profily získať v programe SAS, je v niekoľkých krokoch uvedená v následujúcich bodoch:

K vytvoreniu štandardizovaných profilov potrebujeme v prvom rade vhodné popisné charakteristiky spočítané z dat. Konkrétne výberové priemery a smerodatné odchýlky v jednotlivých (časových) meraniach. Vytvorenie popisných charakteristík (agregované data vzhľadom k jednotlivým týždňom follow-up obdobia) získame napr. následujúcim spôsobom:
```
proc means data=weight.data mean std;
  class week;
  var weight;
  output out=summary mean = MeanValue std = StdValue;
run;
```
```
proc print datafile = summary; 
run;
```
V ďalšom kroku potrebujeme priradiť jednotlivé agregované (empirické) charakteristiky k jednotlivým pozorovaniam (t.j., vhodne prepojiť dva datove súbory). V programe SAS je pre tento krok kľúčové správne usporiadanie jednolivých pozorovaní – v tomto konkrétnom prípade ide o usporiadanie oboch datových súborov podľa premennej, ktorá identifikuje jednotlivé týždne (premenná week):
```
proc sort data=weight.data; by week; run;
proc sort data=summary; by week; run;
```
Zlúčenie oboch datových súborov dohromady a opätovne preusporiadanie podľa jednotlivých subjektov (ako v pôvodných datach):
```
data weight.dataStd;
  merge summary weight.data ;
  by week;
run;

proc sort data=weight.dataStd; by id; run;
```
V aktuálnom datovom súbore už v podstate priamočiaro a jednoducho (pomocou lineárnej kombinácie už existujúcich premenných) vytvoríme novú premennú, ktorá bude obsahovať štandardizované profily. Na záver štandardizované profily vykreslíme:
```
data weight.dataStd;
    set weight.dataStd; 
    weightStd = (weight - MeanValue)/StdValue;
run;
```
```
proc print datafile = weight.dataStd; 
run;
```
```
proc sgplot data=weight.dataStd;;
    series y=weightStd x=week /group=id;
run;
```
Následujúci obrázok uvádza jednak (v ľavom paneli) klasické individuálne profily (teda subject-specific profiles) a tiež (v pravom paneli) štandardizované individuálne profily (teda standardized individual profiles).

Užitočné

Dôležitou časťou exploratívnej analýzy je postupne overovanie podmienok a predpokladov, ktoré sú postulované/predpokládané daným pravdepodobnostným modelom, ktorý na analýzu dát využívame. Niektoré technické predpoklady ale nie je možné overiť priamo z podkladového datového súboru.
Graf štandardizovaných rezíduí je možné napríklad použíť k nahliadnutiu na rezídua, prípadne na overenie niektorých predpokladov súvisiacich s náhodnými chybami.
```
ods graphics on;

proc loess data=weight.dataStd;
   model weightStd=week;
run;
```

Samostatne

Modifikujte predchádzajúci zdrojový kod tak, aby ste získali opäť štandardizované longitudinálne profily, ale štandardizácia bude nie časovo-špecifická, ale subjektovo-špecifická.
Aké základné rozdiely medzi oboma postupmi je možné sledovať? Aký je rozdiel medzi grafom, ktorý obsahuje základné (klasické) longitudinálne profily a grafom, ktorý obsahuje štandardizované profily?
Pokúste sa vymyslieť a navrhnúť iný spôsob vizualizácie longitudinálnych dat, ktorý bude zameraný na nejakú inú charakteristiku časovo závislých pozorovaní.
Pokúste sa vytvoriť novú kategorickú premennú – napr. veľký a malý jedinci – a následne využiť túto premennú v grafe.

3. Cross-sectional vs. longitudinálna štruktúra dat

V tejto časti použijeme o niečo komplexnejší datový súbor o pacientoch so sklerózou multiplex (roztroušená skleróza). Podkladový datový súbor je možné stiahnúť ako csv súbor (následne je potrebné súbor uploadovať do SAS OnDemand a pomocou nasledujúceho SAS kódu načítať).

libname sm '/home/uXXX/sasuser.v94';
filename reffile '/home/uXXX/sasuser.v94/data/sm_data2.csv';

proc import datafile=reffile
    dbms=csv
    out=sm.data
    replace;
    getnames=yes;
run;
    
proc print datafile = sm.data; 
run;

V praxi sa niektoré premenné rozlíšujú do dvoch skupín: tzv. subject-specific variables a tzv. time-dependent variables. V prvom prípade sa jedná o premenné, ktoré sú v čase nemenné (napr. pohlavie pacienta, vek pacienta v čase podstúpenia liečby, zdravotný stav pacienta v čase podstúpenia liečby, alebo aj napr. výskyt/nevýskyt komplikácii počas sledovania). V tom druhom prípade sa jedná o premenné, ktoré sa v rámci konkrétneho jedinca v priebehu opakovaných meraní môžu meniť (napr. hodnota EDSS, alebo opäť výskyt komplikácii zaznamenaný samostatne pri každej kontrole).

Cross-sekčná závislosť

Ide o zisťovanie a vyšetrovanie vzťahov (závislostnej štruktúry) v tzv. subject-specific premenných – porovnávanie vrámci nezávislých subjektov (podskupín, sub-populácii). V prípade datového súboru s pacientami so sklerózou multiplex – napríklad EDSS (t.j., Expanded Disability Status Scale) v závislosti na pohlaví pacienta:

proc loess data = sm.data;
   model EDSS = age / clm smooth=0.1 0.5 residual;;
   by gender;
run;

Longitudinálna závislosť

Jedná sa o zisťovanie vplyvu/zmeny v rámci opakovaných pozorovaní – korelované pozorovania vrámci daného subjektu. Opäť v prípade datového súboru s pacientami so sklerózou multiplex – napríklad EDSS (a jeho vývoj) v závislosti na čase:

proc loess data = sm.data;
   model EDSS = time / clm smooth=0.1 0.5 residual;
run;

Nemusí sa ale nutne jednať pouze o čas ako takový – t.j., premennú time. Analogicky móžeme zisťovať napr. zmenu v EDSS v závislosti na zmene hodnoty v premennej relapse:

proc sort data=sm.data; by relapse; run;

proc boxplot data=sm.data;
   plot EDSS*relapse;
run;

Samostatne

Podrobne sa podívajte na cross-sekčnú a longitudinálnu štruktúru datového súboru sm_data2.csv;
Pokúste sa vytvoriť (v programe SAS) niekoľko vhodných obrázkov, z ktorých by bola zrejmá závislosť EDSS na tzv. “subjec-specific” premenných a tiež závislosť na zmene v čase – tzv. “time-dependent” premenných;
Použijte vlastný (iný) datový súbor a sa pokúste vypracovať pre data analogickú exploratívnu analýzu zohľadňujúc oba typy premenných;

4. Variančná-kovariančná štruktúra

Na záver sa opäť vrátime k datovému súboru o váhach ošípaných, ktorý sme navyše doplnili o štandardizované profily. Vytvorené štandardizované longitudinálne profily použijeme k tomu, aby sme nahliadli do časovo-závislej štruktúry.

proc print datafile = weight.dataStd; 
run;

Existuje samozrejme množstvo rôznych spôsobov, ako empiricky vyšetrovať a zisťovať prípadnú závislosť medzi pozorovaniami. Predpokládame, že jednotlivé súbjekty sú vzájomne nezávisle (malo by to vyplývať z povahy samotného experimentu). Opakované pozorovania v rámci konkretného jedinca ale nezávislé samozrejme nie sú. Pre následné modelovanie (t.j., “confirmatory analysis”) je ale dôležité závislostnú štruktúru opakovaných pozorovaní pochopiť a následne správne v modeli ošetriť.

V prvom kroku vhodne pripravíme data. Pri analýze longitudinálnych dat sa často používa terminológia, ktorá odkazuje na tzv. long-data štruktúru dat, prípadne tzv. wide-data typ dat. Data uložené v datovom súbore weight.dataStd sú reprezentované v tzv. long-data type. Pomocou PROC transpose ale data môžeme reprezentovať vo wide-data formáte:

proc sort data=weight.dataStd; by id; run;

proc transpose data=weight.dataStd out=weight.dataWide;
    by id;
    var weightStd;
run;

Wide-data formát už môžeme priamo využiť pre spočítanie variančenej kovariačnej (resp. korelačnej) matice

proc corr data=weight.dataWide;
var COL1 COL2 COL3 COL4 COL5 COL6 COL7 COL8 COL9;
run;

prípadne tzv. matice scatterplotov, ktorá časovo-závislostnú štruktúru reprezentuje vizuálne:

proc sgscatter data=weight.dataWide;
  title "Scatterplot Matrix -- Repeated observations";
  matrix COL1 COL2 COL3 COL4 COL5 COL6 COL7 COL8 COL9;
run;
title;

Predchádzajúci výstup lze také porovnať s následujúcim výstupom:

proc sgscatter data=weight.dataWide; 
matrix COL1 COL2 COL3 COL4 COL5 COL6 COL7 COL8 COL9 / diagonal=(histogram kernel);
run;

Samostatne

Aké konrétne závery lze získať z pohľadu na odhadnutú korelačnú maticu, pripadne na graf príslušných scatterplotov?
Ako presne chápať a interpretovať grafy v matici scatterplotov?
Aké podstatné následky maju zistenia v empirickej časti analýzy vzhľadom k následnej konfirmačnej analýze?
Premyslite, akým iným spôsobom by bolo možné vizualizovať závislosť/nezávislosť medzi dvoma konkrétnymi premennými.

Samostatný úkol (a)

Pripomeňme, že pre ľubovolnú náhodnú veličinu \(X \sim F\) z nejakého (absolútne) spojitého rozdelenia s distribučnou funkciou \(F\) platí, že \[ P[X \leq x] = F(x), \qquad \forall x \in \mathbb{R}. \] Zároveň ale dostávame, že platí tiež \[ P[F(X) \leq x] = P[X \leq F^{-1}(x)] = F\Big(F^{-1}(x)\Big) = x, \qquad \textrm{pre} \quad x \in [0,1], \] kde \(F^{-1}(u)\) pre \(u \in [0,1]\) je príslučná kvantilová funkcia náhodnej veličiny \(X\) (t.j., inverzná funkcia k distribučnej funkcii \(F\)). Inými slovami, transformovaná náhodná veličina \(F(X)\) má rovnomerné (spojité) rozdelenie na intervale \([0,1]\).

Ako by sa vyššie uvedená (teoretická) vlastnosť dala využiť v (empirickej) exploratívnej analýze pre posúdenie závislosti/nezávislosti nejakých náhodných veličín?
Pomocou programu SAS sa pokúste daný postup implementovať a aplikovať na príklade dvoch konkrétnych (spojitých) premenných, u ktorých chcete posúdiť ich vzájomnú závislosť/nezávislosť.

Samostatný úkol (b)

Na webovej stránke http://faculty.washington.edu/heagerty/Books/AnalysisLongitudinal/datasets.html sú umiestnené niektoré data z knihy Diggle, P. J., Heagerty, P., Liang, K. Y., & Zeger, S. (Analysis of longitudinal data, 2002). Vyberte si jeden datový súbor a vypracujte základnú exploratívnu analýzu pre longitudinálny typ dat. Dôležité zistenia, vhodné grafy s príslušnými popisnými charakteristikami, prípadne závery so stručnou interpretáciou vo forme PDF dokumentu zašlite emialom na maciak[AT]karlin.mff.cuni.cz do dalšieho cvičenia, t.j., v pondelok, 10.03.2025 (do 14:00).