Longitudinální a panelová data – NMST422

Letný semester 2023-2024 | Cvičenie 7 | 22.04.2024

Prihlásenie k SAS OnDemand: https://www.sas.com/en_us/software/on-demand-for-academics.html

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


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

V doterajších častiach cvičenia nás zaujímal hlavne prípad, keď o závislej premennej \(Y\) (ktorá bola opakovane meraná na \(n \in \mathbb{N}\) vzájomne nezávislých subjektoch) môžeme predpokladať, že je spojitá a prípadne navyše aj podmienene normálne rozdelená (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 m_i})^\top \in \mathbb{R}^{m_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 vrámci \(i\)-teho subjektu. Väčšinou predpokládame, že \(\boldsymbol{b}_i \sim N_q(\boldsymbol{0}, \mathbb{D})\). Jednotlivé vektory \(\boldsymbol{Y}_1, \dots, \boldsymbol{Y}_n\) sú vzájomné nezávislé. Chybový člen \(\boldsymbol{\varepsilon}_i\) sa často v lieteratúre rozdeľuje na tzv. chyby merania a stzv. 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}\), tak je nutné využiť iné odhadovacie metódy a postupy (napr. tzv. GEE, 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í – čo platí aj pre lineárne regresné modely s náhodnými efektami (LMM) a samozrejme aj pre zovšeobecnené lineárne regresné modely s náhodnými efektami.

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 (vrámci konrétnych subjektov).



Podrobnejšie sa pozrieme na dve rôzne formulácie v súvislosti s lineárnymi regresnými modelmi 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 – 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/u63241636/sasuser.v94';
filename reffile '/home/u63241636/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";
  matrix time age LEMsum timeBef EDSSini EDSS NEDA
         / group=id;
run;
title;

Analogický postup může byť využitý aj pre nebalancované data, ak je časová doména vhodne diskretizovaná, avšak vhodnejší postup je založený na analýze tzv. 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 tejte 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útˇpomocou rezídui \[ \boldsymbol{r}_i = (r_{i 1}, \dots, r_{i m_i})^\top, \] kde \(r_{i j} = Y_{i j} = \mathbb{X}_i\widehat{\boldsymbol{\beta}}\). Následne lze predpokládať, že rezídua môžeme charakterizovať pomocou modelu (resp. rozhladu) \[ \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 + \sigma^2\mathbb{I}_{m_i} + \tau^2 \mathbb{H}_i. \]
  • 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) (funkcia \(v(\cdot)\)) a 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. Niektoré premenné zobrazíme pomocou jednoduchého scatterplot grafu.

    libname sm '/home/u63241636/sasuser.v94';
filename reffile '/home/u63241636/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";
  matrix time age LEMsum timeBef EDSSini EDSS NEDA
         / group=id;
run;
title;
  • Pomocou 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=red v=dot h=0.2 mode=include; 
symbol2 c=black 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 5 by 0.1) minor = none; 
title h=3 'Semi-variogram';
run;
  • The total variability is estimated as \(3.043\). Measurement error is rather negligible with respect to the overall variability. Serial correlation and random effects seem to be both substantial.



Samostatne

Consider some other longitudinal dataset and calculate the empirical version of the corresponding variogram (semi-variogram respectively). Compare different structures of the variograms and discuss the outcomes.



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 age*time age*time2 / solution;
random intercept time / 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 age*time age*time2 / solution;
random intercept time / type = vc subject = id;
repeated timeCls / type = simple subject = id;

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


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

contrast 'Gender vs. time' age*time 1, 
                           age*time2 1;
run;


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

contrast 'Gender 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 age*time age*time2 / solution;
random intercept time / type = vc subject = id;
repeated timeCls / type = simple subject = id;
run;


proc mixed data = sm.data2 method = reml; 
class id gender timeCls;
model EDSS = gender age age*time age*time2 / solution;
random intercept time / type = vc subject = id solution;
repeated timeCls / type = simple subject = id;
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?