Longitudinální a panelová data – NMST422

Letný semester 2024-2025 | Cvičenie 7 | 07.04.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ácie prostredníctvom emailu. 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. Symbol XXX v zdrojových kódoch 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


VII. Lineárne modely s náhodnými efektami (inferencia)

Predpokládajme lineárny regresný model pre spojitú závislú premennú \(Y\) (t.j., náhodnú veličinu, ktorá bola opakovane nameraná/realizovaná na \(N \in \mathbb{N}\) vzájomne nezávislých subjektoch. Navyše budeme predpokládať, že \(Y\) má podmienené normálne rozdelenie (t.j., základný lineárny regresný model s náhodnými efektami). Matematicky je tento fakt vyjadrený prostredníctvom zápisu \[ \boldsymbol{Y}_i = \mathbb{X}_i\boldsymbol{\beta} + \mathbb{Z}_i\boldsymbol{b}_i + \boldsymbol{\varepsilon}_i \] kde \(\boldsymbol{Y}_i = (Y_{i1}, \dots, Y_{i n_i})^\top \in \mathbb{R}^{n_i}\) je vektor opakovaných meraní vrámci \(i\)-teho subjektu (pre \(i \in \{1, \dots, N\}\)) a \(\boldsymbol{b}_i = (b_{i1}, \dots, b_{i q})^\top \in \mathbb{R}^q\) je vektor náhodných (nepozorovaných) efektov v rámci \(i\)-teho subjektu, pre ktorý väčšinou predpokládame, že \(\boldsymbol{b}_i \sim N_q(\boldsymbol{0}, \mathbb{D})\) a náhodný vektor chýb \(\boldsymbol{\varepsilon}_{i}\) má opäť mnohorozmerné normálne rozdelenie s nulovým vektorom stredných hodnôt a určitou variančnou-kovariančnou štruktúrou. Jednotlivé vektory \(\boldsymbol{Y}_1, \dots, \boldsymbol{Y}_N\) sú vzájomné nezávislé a chybový člen \(\boldsymbol{\varepsilon}_i\) sa často v literatúre rozdeľuje na tzv. chyby merania a tzv. sériovú koreláciu.

Poznámka

V praxi sa ale často stane, že predpoklad (mnohorozmerného) normálneho rozdelenia pre opakovane pozorovania – t.j., náhodné vektory \(\boldsymbol{Y}_i\), pre \(i =1, \dots, n\) je nerealistický a je nutné hľadať iný pravdepodobnostný model (napr. pretože sledovaná závislá premenná informuje výhradne len o úspechu/neúspechu liečby – binárna premenná – alebo sa všeobecne jedná o realizácie nejakej diskrétnej náhodnej veličiny). Ak je možné naviac postulovať (predpokladať) konkrétne rozdelenie pre závislú premennú (to znamená aj možnosť definovať celkovú vierohodnosť), tak je vhodné použíť tzv. zovšeobecnené lineárne modely s náhodnými efektami (generalized linear model with random effects). Jedná sa o rozšírenie triedy zovšeobecnených lineárnych regresných modelov (GLM) v podobnom zmysle, ako sú lineárne regresné modely s náhodnými efektami zovšeobecnením klasických lineárnych regresných modelov.

Ak nie je možné apriórne postulovať (predpokladať) nejaké vhodné pravdepodobnostné rozdelenie pre závislú premennú \(\boldsymbol{Y}\), ale je možné špecifikovať niektoré (podmienené) momenty náhodných vektorov \(\boldsymbol{Y}_i\), tak je vpohodné využiť iné odhadovacie metódy a postupy (napr. tzv. GEE metódy, ktoré budeme diskutovať neskôr).



Na rozdiel od klasických lineárnych regresných modelov, ktoré predpokládajú nezávislé pozorovania, je nutné pri modeloch s náhodnými efektami náležite zohľadniť korelačnú štruktúru v rámci opakovaných pozorovaní Pri analýze dat je preto nutné dbať na správnu špecifikáciu variančnej-kovariančnej (resp. korelačnej) štruktúry a správny popis jednotlivých zdrojov variability – t.j. variabilita medzi jednotlivými subjektami a variabilita opakovaných pozorovaní v rámci konkrétheho jedinca (v rámci konrétnych subjektov).


Podrobnejšie sa pozrieme na dve rôzne formulácie modelu s náhodnými efektami (LMM)

  • Marginálna formulácia modelu, kde predpokládame \(\boldsymbol{Y}_i \sim N(\mathbb{X}_i\boldsymbol{\beta}, \mathbb{Z}_i \mathbb{D}\mathbb{Z}_i^\top + \mathbb{R}_i)\)
  • Hierarchická formulácia modelu, kde predpokládame \(\boldsymbol{Y}_i | \boldsymbol{b}_i \sim N(\mathbb{X}_i\boldsymbol{\beta} + \mathbb{Z}_i \boldsymbol{b}_i, \mathbb{R}_i)\)



1. Teoretický a empirický variogram

V prípade opakovaných balancovaných pozorovaní je variančnú-kovariančnú štruktúru možné celkom dobre exploratívne analyzovať napr. pomocou výberovej korelačnej matice, prípadne vhodného scatterplot grafu (scatterplot grafov). V prípade, že celková dimenzia datoveho súboru nie je príliš veľka, celkom nápomocné sú tzv. pairwise scatterplots, ktoré postupne vykresľujú scatterplot pre všetky možné dvojice uvažovaných premenných. Avšak aj kompletný súbor všetkých možných vzájomných dvojíc scatterplotov neobsahuje komplexnú (mnohorozmernú) závislostnú štruktúru celkového datového súboru – vždy sa jedná len o určitú projekciu do dvojrozmerného priestoru—t.j., do 2D roviny grafu.

Pre niektoré uvažované premenné z datového súboru o pacientoch so sklerózou multiplex dostaneme v programe SAS prislušnú maticu scatterplotov následujúcim kódom:

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 sgscatter data=sm.data;
  title "Scatterplot Matrix (SM data)";
  matrix time age LEMsum timeBef EDSSini EDSS NEDA
         / group = gender diagonal = (histogram);
run;
title;

Analogický postup může byť použitý aj pre nebalancované data, ak je časová doména vhodne diskretizovaná, avšak vhodnejší postup je založený na analýze variogramu, resp. semi-variogramu.

Užitočné

  • V programe R je pre analýzu longitudinálnych dat a odhadovanie lineárneho regresného modelu s náhodnými efektami určená (napr.) knižnica nmle, V rámci tejto knižnice je k dispozícii funkcia Variogram(), ktorá počíta výberový variogram, resp. semi-variogram.
  • Alternatívnou možnosťou je napr. funkcia variogram() dostupná v R knižnici gstat (install.packages("gstat"))
  • Zaujímavý prehľad k variogramu a modelovaniu kovariancie obecně, je napr. dostupný na tejto stránke.



Pre pripomenutie zopakujme, že linárny regresný model s náhodnými efektami môže byť zapísaný (použitím vhodnej dekompozície nepozorovaných – latentných premenných – náhodných chýb) v tvare \[ \boldsymbol{Y}_i = \mathbb{X}_i\boldsymbol{\beta} + \mathbb{Z}_i\boldsymbol{b}_i + \boldsymbol{\varepsilon}_{1 i} + \boldsymbol{\varepsilon}_{2 i}, \] kde jednotlivé stochastické členy v zápise postupne predstavujú
  • \(\boldsymbol{b}_i\) – medzi-subjektovú variabilitu (between-subject variability);
  • \(\boldsymbol{\varepsilon}_{1 i}\) – tzv. chyby merania (measurement error);
  • \(\boldsymbol{\varepsilon}_{2 i}\) – tzv. sériovu koreláciu (serial correlation).



Ak budeme uvažovať model pre strednú hodnotu (t.j. pre časť \(\mathbb{X}_i\boldsymbol{\beta}\)) , tak chybovú (reziduálnu) časť můžeme následne odhadnúť pomocou \(N \in \mathbb{N}\) rezídui (vektorov rezídui) \[ \boldsymbol{r}_i = (r_{i 1}, \dots, r_{i n_i})^\top, \] kde \(r_{i j} = Y_{i j} - \mathbb{X}_i\widehat{\boldsymbol{\beta}}\), pre \(i = 1, \dots, N\). Následne lze predpokládať, že rezídua je možné charakterizovať pomocou modelu (resp. rozhladu) ako \[ \boldsymbol{r}_i = \mathbb{Z}_i\boldsymbol{b}_i + \boldsymbol{\varepsilon}_{1 i} + \boldsymbol{\varepsilon}_{2 i}. \]

Užitočné

  • Ak v regresnom modeli predpokládame pouze náhodný intercept, pre rozptyl (rezídua) to ma za následok predpoklad homoskedasticity (rozptyl je konštantný). Celkovú variančnú-kovariančnú maticu \(\boldsymbol{Y}_i\) resp. \(\boldsymbol{r}_i\) môžeme vyjadriť v tvare \[ Var(\boldsymbol{Y}_i) = Var(\boldsymbol{b}_i) = \nu^2 \mathbb{Z}_i\mathbb{Z}_i^\top + \sigma^2\mathbb{I}_{m_i} + \tau^2 \mathbb{H}_i, \] kde \(\mathbb{Z}_i\mathbb{Z}_i^\top\) je matica typu \(n_i \times n_i\), ktorá ma na všetkých pozíciach hodnotu jedna.


  • Korelácia medzi dvoma rezíduammi \(r_{ij}\) a \(r_{ik}\) (v rámci toho istého \(i\)-teho subjektu) je potom vyjadrená ako \[ \rho(|t_{ij} - t_{ik}|) = \frac{\nu^2 + \tau^2 g(|t_{ij} - t_{ik}|)}{\nu^2 + \sigma^2 + \tau^2}. \]


  • Pre \(j \neq k\) teda platí, že \[ \frac{1}{2}E(r_{ij} - r_{ik})^2 = \sigma^2 + \tau^2 (1 - g(|t_{ij} - t_{ik}|)) = v(u_{ijk}), \] kde posledný člen sa nazýva variogram (semi-variogram)—t.j., funkcia \(v(\cdot)\), ktorá závisí pouze na časových okamžikoch jednotlivých meraní, t.j. \(t_{i j}\) a \(t_{i k}\). Klesajúca funkcia \(g(\cdot)\) má za následok rastúci priebeh variogramu, pričom \(v(0) = \sigma^2\) a \(\lim_{u \to \infty} v(u) = \sigma^2 + \tau^2\).

  • Príslušný empirický odhad teoretického variogramu (semi-variogramu) môže byť priamočiaro využitý k exploratívnej analýze variančnej-kovariančnej štruktúry dat.



Samostatne

  • Pre ilustraciu budeme opäť uvažovať datový súbor pacientov so sklerózou multiplex, ale namiesto scatterplot matice využijeme variogram aplikovaný na rezídua. Pomocou SAS helpu pre štatistický program SAS sa podívajte na procedúru PROC VARIOGRAM.

  • Jednoduchá ukážka aplikácie variogramu v programe SAS (na základe odhadnutých rezídui z jednoduchého regresného modelu):

    /* Simple model for the mean */
proc glm data=sm.data;
model EDSS=time;
output out=out r=residual;
run;

/* Calculation of the variogram */
proc variogram data=out outpair=out;
coordinates xc=time yc=id;
compute robust novariogram;
var residual;
run;

data variogram;set out;
if y1=y2;vario=(v1-v2)**2/2; run;
data variance;set out;
if y1<y2; vario=(v1-v2)**2/2; run;

/* Calculation of the total variance */
proc means data=variance mean;
var vario;
run;

proc loess data=variogram;
ods output scoreresults=out;
model vario=distance;
score data=variogram;
run;

proc sort data = out; by distance; run; 

proc gplot data = out; 
plot vario*distance =1 p_vario*distance =2 / overlay haxis = axis1 vaxis = axis2 vref = 3.043 lvref=3;
symbol1 c=black v=dot h=1 mode=include; 
symbol2 c=red i=join w=2 mode=include; 
axis1 label=(h = 2 'Time lag') value =(h = 1.5) order = (0 to 8 by 1) minor=none; 
axis2 label=(h = 2 A = 90 'V(u)') value=(h=1.5) order = (0 to 4 by 0.2) minor = none; 
title h=3 'Semi-variogram';
run;



2. Odhadovanie modelu s náhodnými efektami v SAS

Základným nástrojom pre odhadovanie lineárnych regresných modelov s náhodnými efektami v programe SAS je procedúra PROC MIXED. Parametre sú odhadované pomocou metódy maximálnej vierohodnosti, pomocou metódy rezidálnej vierohodnosti (REML – default v programe SAS) a prípadne dalšími iterativnými numerickými algoritmami.



Pre správnu implementáciu modelu sú dôležité následújúce časti procedúry:

  1. CLASS Statement
    Definícia faktorových/diskrétných premenných. Pre program SAS je nutné definovať aj čas opakovaných pozorovaní (v rámci spojitého stochastického procesu \(W_I(t)\)) ako diskrétnú premennú–t.j., pozorovania v rámci diskreétnych časových okamžikoch \(t_{i1} < \dots < t_{i n_i}\)
  2. MODEL Statement
    Špecifikácia podmienenej strednej hodnoty \(\mathbb{X}_i \boldsymbol{\beta}\)–t.j. pevných efektov. Dodatočný parameter solution zabezpečí vypísanie odhadnutých pevných efektov–odhadnutej štruktúry podmienenej strednej hodnoty. Rôzne parametrizácie pre faktorové premenné je možne ovplyvniť voľbou parametra noint.
  3. RANDOM Statement
    Špecifikácia náhodných efektov a to vrátane náhodného interceptu, ktorý musí byť uvedený explicitne. Potrebná je aj špecifikácia subjektov, v rámci ktorých sa predpokladá nezávislosť. Parameter type = upresňuje variančnú-kovariančnú maticu \(\mathbb{D}\) pre vektor náhodných efektov. V návode (help) procedúry pozri aj parametre g, gcorr, v a vcorr.
  4. REPEATED Statement
    Špcifikácia poradia opakovaných pozorovaní v rámci konkrétneho subjektu. Opakované pozorovania sú špecifikované pomocou faktorovej/ordinálnej premennej (viď timeClas v zásise nižšie).Nutná je opäť identifikácia subjektov. Pomocou parametru type = sa volí variančná-kovariačná matica pre sériouvú koreláciu opakovaných pozorovaní. V návode (help) procedúry pozri aj parametre r a rcorr. 



data sm.data2;
set sm.data;
timeCls = time;
time2 = time * time;
run; 

proc mixed data = sm.data2 method = reml; 
class id gender timeCls;
model EDSS = gender age time time2 age*time  / solution;
random intercept  / type = vc subject = id;
repeated timeCls / type = simple subject = id;
run;



Samostatne

Všimnite si, že v prípade variogramu počítanom vyššie, boli využité rezidua z jednoduchého regresného modelu, kde EDSS záviselo pouze na čase (jednotlivé roky od začiatku liečby). Model, ktorý je prezentovaný vyššie (asi nie finálny model, ale určite model, ktorý je zo štatistického hľadiska významne lepší, než závislosť EDSS na čase), používa komplexnejšiu štruktúru pre modelovanie podmienenej strednej hodnoty.

  • Spočítajte príslušný variogram na základe rezídui z komplexného modelu a porovnajte oba variogramy.
  • Aké sú závery plynpce z variogramu na základe rezídui z komplexnejšieho modelu?
  • Pokúste sa model dopracovať a podmienenú strednú hodnotu modelovať čo najlepším modelom (v rámci frameworku lineárnych regresných modelov). Využijte rezídua modelu a opäť spočítajte príslušný variogram (semi-variogram).



3. Inferencia v modeloch s náhodnymi efektami

Otázka inferencie v lineárnom regresnom modeli s náhodnými efektami sa týká jednak štruktúry podmienenej strednej hodnoty (t.j., vektor nezámych parametrov \(\boldsymbol{\beta} \in \mathbb{R}^p\)), ale taktiež aj neznámych parametrov \(\boldsymbol{\alpha} \in \mathbb{R}^q\), ktoré parametrizuju podmienenú variančnú/kovariančnú štruktúru.

Niektoré (najčastejšie používané) základne inferenčné nástroje:

  • Inference pro parametre strednej hodnoty
    • približné Waldove testy
    • \(t\)-testy a \(F\)-testy s korekciou pre odhad \(\widehat{\boldsymbol{\alpha}}\)
    • testy pomerom vierohodnosti (nie pre REML)
    • robustné prístupy za podmienky správne špecifikovanej podmienenej strednej hodnoty


  • Inference pro parametre variančnej-kovariančnej štruktúry
    • približné Waldové testy
    • testy pomerom vierohodnosti



V nasledujúcom SAS kóde si všimnite úlohu jednotlivých parametrov (pomocou návodu zistite ich význam) –inferencia o parametroch \(\boldsymbol{\beta}\in \mathbb{R}^p\) (contrast statement a voľba (napr.) chisc) a tiež inferencia o parametroch \(\boldsymbol{\alpha} \in \mathbb{R}^q\) (voľba covtest v proc mixed statement).

V nasledujúcich častiach je postupne robená inferencia vzhľadomk k \(\boldsymbol{\beta}\in \mathbb{R}^p\) parametrom pomocou približného Waldovho testu, pomocou \(t\)-testu (\(F\)-testu resp.) s Kenward-Roger aproximáciou stupňov voľnosti, robustná inferencia, inferenia ohľadom náhodných efektov a individuaálne profily.

proc mixed data = sm.data2 method = ml; 
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;

contrast 'Age vs. time' age*time2 /  chisq;
run;


proc mixed data = sm.data2 method = ml; 
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;

contrast 'Age vs. time' age*time2 1,
                                    age*time 1 /  chisq;
run;


proc mixed data = sm.data2 method = reml empirical; 
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;

contrast 'Age vs. time' age*time 1, 
                                    age*time2 1 / chisq;
run;


proc mixed data = sm.data2 method = reml covtest; 
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;

contrast 'Age vs. time' age*time 1, 
                                    age*time2 1 / chisq;
run;


proc mixed data = sm.data2 method = reml; 
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id solution;
repeated timeCls / type = simple subject = id;
run;



Významnosť náhodných efektov je možné testovať aj pomocou testu pomerom vierohodnosti. Je nutné si ale uvedomiť, že vpodstate testujeme nulovú hypotézu v zmysle \[ H_0: \nu^2 = 0 \] teda príslušný náhodný efekt mú nulový rozptyl, pričom hodnota nula je na hranici parametrického priestoru pre rozptyl \(\nu^2 > 0\). Výsledné rozdelenie testovej štatistiky preto nemá \(\chi^2\) rozdelenie s jedným stupňom voľnosti (ako by tomu bolo v prípade testu oľadom hodnoty unitř parametrického prostoru), ale jedná sa o rovnoměrnou změs dvoch rozdelení – \(\chi^2\) rozdelenie s jedným stupňom voľnosti a \(\chi^2\) rozdelenie s nula stupňami voľnosti.

  /* Full model with random effect */
  proc mixed data = sm.data2 method=ml noclprint noinfo covtest;
    class id gender timeCls;
    model EDSS = gender age time time2 age*time age*time2;
    random intercept  / subject = id;
    repeated timeCls / type = simple subject = id;
    ods output FitStatistics=full_fit;
  run;

  /* Reduced model without random effect */
  proc mixed data = sm.data2 method=ml noclprint noinfo;
    class id gender timeCls;
    model EDSS = gender age time time2 age*time age*time2;
    repeated timeCls / type = simple subject = id;
    ods output FitStatistics=reduced_fit;
  run;

  /* Extract -2 Log Likelihoods */
  data full_ll;
    set full_fit;
    if Description = '-2 Res Log Likelihood' then call symputx('ll_full', Value);
  run;

  data reduced_ll;
    set reduced_fit;
    if Description = '-2 Res Log Likelihood' then call symputx('ll_reduced', Value);
  run;

  /* Compute LRT and p-value */
  data lrt_result;
    ll_diff = 3026.7 - 1765.9;
    df = 1;
      
    p0 = (ll_diff = 0);
    p1 = 1 - cdf('CHISQ', ll_diff, 1);

    /* Mixture p-value: 0.5 * p0 + 0.5 * p1 */
    p_value = 0.5 * p0 + 0.5 * p1;

    put "----------------------------------------";
    put "Likelihood Ratio Test for Random Effect:";
    put "LRT statistic = " ll_diff;
    put "Degrees of freedom = " df;
    put "P-value = " p_value;
    put "----------------------------------------";
  run;



Samostatne

Pomocou helpu k SAS procedure (napr. na tejto stránke) sa pokúste urobiť štatistickú inferenciu o neznámych parametroch (pre strednú hodnotu aj pre variančnú-kovariančnú štruktúru).

  • využijte približný Wald test aj presnejšie \(t\)-testy prípadne \(F\)-testy
  • využijte rôzne metódy pre aproxímáciu stupňov voľnosti v menovateli testovej štatistiky (\(t\)-testy a \(F\)-testy)
  • akým spôsobom je možné urobiť štatistickú inferenciu ohľadom variančnej-kovariančnej štruktúry?