Longitudinální a panelová data – NMST422

Letný semester 2024-2025 | Cvičenie 11 | 12.05.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


XI. GLM modely s náhodnými efektami (GLMM)

Už bolo spomínané, že v prípade GLM modelov (t.j., marginálnych modelov) pre korelované/opakované (longitudinálne) data \(\{(Y_{ij}, \boldsymbol{X}_{ij});i = 1, \dots, N; j = 1, \dots, n_i\}\) merané na \(N \in \mathbb{N}\) nezávislých subjektov (pričom celkový počet pozorovaní je \(\mathcal{N} = \sum_{i = 1}^N n_i\)) špecifikujeme analogické podmienky ako v prípade klasických GLM modelov pre nezávislé pozorovania, t.j, formu podmienenej strednej hodnoty \[ \mu_{i j} = E\Big[Y_{i j} | \boldsymbol{X}_{ij}\Big]= g^{-1}(\boldsymbol{X}_{ij}^\top\boldsymbol{\beta}); \] resp. vyjadrené vektorovo/maticovo pre jednotlivé subjekty ako \[ E\Big[\boldsymbol{Y}_{i}| \mathbb{X}_{i}\Big] = \boldsymbol{\mu}_i = (\mu_{i1}, \dots, \mu_{i n_i})^\top \] s variančnou-kovariačnou maticou \(\mathcal{V}_i = Var \boldsymbol{Y}_i\), ktorá zohľadňuje korelácie opakovaných pozorovaní. Odhad neznámych parametrov vedie na riešenie nelineárnych rovníc, ktoré špecifikujú konkrétny vzťah medzi prvými dvoma teoretickými momentami a rovnice \[ \sum_{i = 1}^N \frac{\partial \boldsymbol{\mu}_{i}^\top}{\partial \boldsymbol{\beta}} \mathcal{V}_i^{-1}(\boldsymbol{Y}_{i} - \boldsymbol{\mu}_i) = \boldsymbol{0}. \] sa štandardne riešia pomocou vhodných iteračných algoritmov – napr. Newton-Raphson algorithmus. ´

Pridanie náhodných efektov do GLM modelu (tzv. GLMM modely) ma za následok rozšírenie vyššie uvedenej špecifikácie podmienenej strednej hodnoty, kedy uvažujeme pre podmienenú strednú hodnotu tvar \[ \mu_{i j} = E\Big[Y_{i j} | \boldsymbol{X}_{ij}, \boldsymbol{b}_i\Big]= g^{-1}(\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + \boldsymbol{Z}_{ij}\boldsymbol{b}_i), \] kde \(\boldsymbol{b}_i \sim N_{r}(\boldsymbol{0}, \mathbb{D})\) je vektor náhodných efektov a \(\mathbb{Z}_i = (\boldsymbol{X}_{i1}, \dots, \boldsymbol{X}_{i n_i})^\top \in \mathbb{R}^{n_i \times r}\) je príslušná matica dizajnú pre pevné efekty (a pre konkrétny subjekt \(i \in \{1, \dots, N\}\)).

V závislosti na konkrétnom predpokladanom rozdelení (z rodiny rozdelení exponenciálneho typu) platí príslušný vzťah medzi podmieneným prvým momentom, teda \(\mu_{ij}\), a podmieneným druhým momentom \(v_{ij} = Var[Y_{ij}|\boldsymbol{X}_{ij}, \boldsymbol{b}_i] = \phi v(\mu_{ij})\), kde \(v(\cdot)\) je variančná funkcia a \(\phi > 0\) je tzv. overdisperzný parameter.

Výsledné odhady sa štandardne získavajú metódou maximálnej vierohodnosti, avšak s využitím rôzncyh aproximačných a iteratívnych postupov (z dôvodu značnej zložitosti výslednej združenej vierohodnosti).



Samostatne



1. Proc GLIMMIX (SAS)

GLM modely s náhodnými efektami (tzv. hierarchický GLM model) je možné v programe SAS odhadovať pomocou procedúry PROC GLIMMIX. Marginálna interpretácia modelu ale nie je totožná s interpretáciou marginálnych GLM modelov (viď nižšie).

Pre jednoduchú ilustráciu využijeme opäť data s pacientami so sklerózou multiplex a ako závislú premennú budeme uvažovť ukazateľ aktivity nemoci – NEDA (No Evidence of Disease Activity).

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 data=sm.data; 
run;

Porovnajte následujúce dva modely

proc genmod data=sm.data DESCENDING;
class id gender;
model NEDA = gender age time EDSS / dist=binomial;
run;

proc glimmix data=sm.data;
class id gender;
model NEDA = gender age time EDSS / dist=binomial solution;
run;

V oboch prípadoch sa jedna o modely bez náhodných efektov a za predpokladu nezávislých pozorovaní. Samotný odhadovací algoritmus je ale v oboch prípadoch iný, ale výsledné modely sú totožné. Odhadnuté parametre sú získane pomocou metódy maximálnej vierohodnosti iteratívnym riešením skórových rovníc. Modely ale nezohľadŇujú opakované pozorovania.

Pre zohľadnenie opakovaných/korelovaných pozorovaní (v prípade oboch procedúr), je možné použíť napr. následujúci kód:

proc genmod data=sm.data DESCENDING;
    class id gender;
    model NEDA = gender age time EDSS / dist=binomial;
    repeated subject=id / type=exch corrw;
run;

proc glimmix data=sm.data;
    class id gender;
    model NEDA = gender age time EDSS / dist=binomial solution;
    random intercept / subject=id vcorr;
run;

V oboch modeloch vyššie sa predpokladá rovnaká korelácia medzi ľubovoľnými dvoma opakovanými pozorovaniami v rámci jedného subjektu. V prípade PROC GENMOD je tento fakt reprezentovaná pomocou parametru type = exch a v prípade PROC GLIMMIX je toto zohľadnené zahrnutím náhodného interceptu. Výsledná odhadnutá korelačná matica ale tomuto požiadavku v prípade PROC GLIMMIX nezodpovedá… Prečo je tomu tak?

Pre pripomenutie, \[ Cor (Y_{ij}, Y_{ik}) = \frac{Cov(Y_{ij}, Y_{ik})}{\sqrt{Var Y_{ij} Var Y_{ik}}}, \] pričom oba rozptyly v menovateľi navyše závisia na podmienených stredných hodnotách \(\mu_{ij} = E[Y_{ij}|\boldsymbol{X}_{ij}, b_i]\) a \(\mu_{ik} = E[Y_{ik}|\boldsymbol{X}_{ik}, b_i]\) prostredníctvom vzťahu \(Var Y_{ij} = \phi v(\mu_{ij})\), kde \(v(\cdot)\) je variančná funkcia a \(\phi > 0\) je tzv. overdisperzný parameter.



Procedúra PROC GLIMMIX ale umožňuje pridanie viacerých náhodných efektov v random statement – viď model nižšie:

proc glimmix data=sm.data method =laplace;
class id gender;
model NEDA = gender age time EDSS / dist=binomial solution;
random intercept time / type = vc subject=id solution g;
run;
run;



Samostatne

  • GLM model s náhodným interceptom vyššie modifikujte a využijte komplexnejšiu štruktúru náhodných efektov.
  • Použijte rôzne nastavenia pre voľbu pracovnej korelačnej matice (parameter type = ...).
  • Porovnajte jednotlivé modely a pokúste sa interpretovať odhadnuté parametre.



3. SAS procedúra PROC NLMIXED

Pre odhadovanie nelineárnych regresných modelov (s náhodnými efektami, t.j. pre korelované/opakované pozorovania) je v programe SAS k dispozícii procedúra PROC NLMIXED – podrobná dokumentácia je k dispozícii na stránke
https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/statug/statug_nlmixed_toc.htm>



Principiálne sa ale jedná o iný model, keď modelujeme podmienenú strednú hodnotu náhodnej veličiny \(Y_{ij}\) ako \[ E\Big[Y_{i j} | \boldsymbol{X}_{ij}\Big]= \exp \{\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} \} \] v prípade marginálneho modelu, do ktorého lze zapracovať aj náhodný efekt \[ E\Big[Y_{i j} | \boldsymbol{X}_{ij}, u_i\Big]= \exp \{\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + u_i \}. \]

Príslušná implementácia v programe SAS by vyzerala následovne:

proc nlmixed data=sm.data;
parms beta0=0 beta1=0 beta2=0 beta3 = 0;
 lambda = exp(u + beta0 + beta1*age + beta2*time + beta3*EDSSini);
 model EDSS ~ poisson(lambda);
 random u ~ normal(0,1) subject=id;
run;

proc nlmixed data=sm.data;
parms beta0=0 beta1=0 beta2=0 beta3 = 0;
 lambda = exp(u + beta0 + beta1*age + beta2*time + beta3*EDSSini);
 model EDSS ~ poisson(lambda);
 random u ~ normal(0,10) subject=id;
run;



Je dôležité porozumieť základným rozdielom medzi jednotlivými modelmi, kde modelujeme nelineárne strednú hodnotu náhodnej veličiny \(Y_{ij}\) – napr. model implementovaný pomocou procedúry PROC NLMIXED v tvare \[ Y_{ij} = exp(\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + b_i) + \varepsilon_{ij} \] prípadne modelujeme strednú hodnotu transformovanej náhodnej veličiny \(\log(Y_{ij})\) pomocou lineárneho regresného modelu (a klasickej procedúry PROC MIXED) v tvare \[ \log Y_{ij} = \boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + b_i + \varepsilon_{ij}, \] alebo pomocou GLM modelu a vhodnej linkovacej funkcie modelujeme strednú hodnotu pomocou lineárneho prediktoru \[ log \big( E[Y_{ij} | \boldsymbol{X}_{ij}, b_i] \big)= \boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + b_i. \]



Samostatne

  • V čom spočíva zakladný rozdiel medzi modelmi uvedenými/popisanými vyššie?
  • Pokúste sa v programe SAS implementovať uvedené modely a jednotlivé modely následne porovnajte.
  • Aká je výsledná interpretácia odhadnutých parametrov v jednotlivých modeloch?