Letný semester 2025-2026 | Cvičenie 8 | 04.05.2026

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

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


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

Už bolo spomínané, že v prípade marginálnych GLM 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 subjektoch (napr. GEE) , 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 – teda štruktúru 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 = (g^{-1}(\boldsymbol{X}_{i1}^\top\boldsymbol{\beta}), \dots, g^{-1}(\boldsymbol{X}_{i n_{i}}\boldsymbol{\beta}))^\top \] a variančnu-kovariačnú maticu \(\mathcal{V}_i = Var \boldsymbol{Y}_i\), ktorá je implikovaná väzbou medzi prvým a druhým momentom v obecnom rozdelení exponenciálneho typu a zohľadňuje korelačnú štruktúru opakovaných pozorovaní.

Odhad neznámych parametrov \(\boldsymbol{\beta} \in \mathbb{R}^p\) vedie na riešenie nelineárnych rovníc \[ \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}, \] ktoré berú do úvahy prvé dva teoretické momenty a štandardne sa riešia pomocou vhodných iteračných metód – napr. Newton-Raphson algorithmus. ´

Pridanie náhodných efektov do GLM modelu (tzv. GLMM modely) má za následok rozšírenie vyššie uvedenej špecifikácie podmienenej strednej hodnoty, kedy uvažujeme pre podmienenú strednú hodnotu (hierarchický) model \[ \mu_{i j} = E\Big[Y_{i j} | \boldsymbol{X}_{ij}, \boldsymbol{w}_i\Big]= g^{-1}(\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + \boldsymbol{Z}_{ij}\boldsymbol{w}_i), \] kde \(\boldsymbol{w}_i \sim N_{r}(\boldsymbol{0}, \mathbb{G})\) je \(r\)-rozmerný 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 dizajnu 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, ktorý lze obecne vyjádriť v tvare \[ v_{ij} = Var[Y_{ij}|\boldsymbol{X}_{ij}, \boldsymbol{w}_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ôznych aproximačných a iteratívnych postupov (z dôvodu nelinearity linkovej funkcie \(g(\cdot)\), \(r\)-rozmerného integrálu vzhľadom k náhodným efektom a teda značnej výpočetnej zložitosti pri maximalizácii 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 odhadovaných napr. pomocou procedúry PROC GENMOD, prípadne pomocou zovšeobecnených odhadovacích rovníc (GEE) – teda procedúry PROC GEE.

Pre jednoduchú ilustráciu použijeme data s pacientami so sklerózou multiplex a ako závislú premennú budeme uvažovať ukazateľ aktivity nemoci – t.j., binárnu veličinu NEDA (No Evidence of Disease Activity). Pre jednoduchosť je podmienená stredná hodnota modelovaná len pomocou informácie o pohlaví pacienta (muž/žena), veku pacienta (v rokocho), čase od nástupu na liečbu (premenná time) a hodnoty EDSS (ktorá bola v predchádzajúcich cvičeniach používaná ako závislá premenná).

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 method=mmpl;
class id gender;
model NEDA = gender age time EDSS / dist=binomial solution;
run;



V oboch prípadoch sa jedná 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ú ekvivalentné.

Odhadnuté parametre v prípade PROC GLIMMIX sú získané pomocou metódy maximálnej vierohodnosti a odhady v rámci PROC GENMOD sú získané 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 modelov), je nutné prislušné procedúry upraviť a doplniť požadovanú korelačnú štruktúru.

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á.

V prípade regresného modelu (s náhodným interceptom) pre normálne rozdelené data (t.j. identita ako linková funkcia) by korelácia medzi ľubovolnými dvoma opakovanými pozorovaniami bola skutočne rovnaká. V prípade GLM modelu s logit linkovou funkciou to ale nie je pravda – resp. na príslušnej logit škále je korelácia konštatná, ale na škále pozorovaných dat, teda \(Y_{ij}\), je korelácia už obecne rôzn a navýše závisí na podmienených stredných hodnotách (t.j., konštantný aditívny posun v rámci náhodného interceptu na logit škále má rôzny efekt v závislosti na tom, na akej hodnote sa v zmysle odhadovanej pravdepodobnosti na logit krivke nachádzame).

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 menovateii navyše závisia na podmienených stredných hodnotách \(\mu_{ij} = E[Y_{ij}|\boldsymbol{X}_{ij}, w_i]\) a \(\mu_{ik} = E[Y_{ik}|\boldsymbol{X}_{ik}, w_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.

Dodatočný parameter type = ... nemá teda tú istú funkcionalitu, ako v prípade marginálnych modelov, alebo klasických lineárnych modelov s náhodnými efektami. Parameter určuje korelačnú štruktúru náhodných efektov, ale neimplikuje korelačnú štruktúru samotných opakovaných pozorovaní. Tá je určená samotným predpokládaným rozdelením.



Parameter type = ... má teda význam iba v prípade viac-rozmerných náhodných efektov, čo procedúra PROC GLIMMIX aumožňuje v rámci random statement:

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;

Je dobré si ale všimnúť, že algorithmus založený priamo na reziduálnej vierohodnosti (defaultné nastavenie v SAS je method = RSPL) nekonverguje k výsledným odhadom. Na druhej strane, method = laplace, ktorá využíva podmienenú (hierarchickú) vierohodnosť neumožňuje zobraziť korelačnú maticu pre opakované pozorovania (t.j., využit parameter vcorr), pretože tá vychádza z marginálnej štruktúry modelu.

Metóda Laplace je vo všeobecnosti lepšia pre tento prípad – binárne data.



Samostatne

  • Pokúste sa model vyššie rozšíriť a uvažovať komplexnejšiu štruktúru pevných, či 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 (náhodný intercept) \[ E\Big[Y_{i j} | \boldsymbol{X}_{ij}, w_i\Big]= \exp \{\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + w_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} + w_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} + w_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}, w_i] \big)= \boldsymbol{X}_{ij}^\top\boldsymbol{\beta} + w_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?