Letný semester 2022-2023 | Cvičenie 13 | 12.05.2023
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ályPROC GLIMMIX a PROC GEE v
Poissonovej regresiiJedná z ďalších konkrétnych volieb linkovacej funkcie \(g(\cdot)\) v obecnej GLM regresii je
logaritmus, teda \(g(x) = log(x)\), pre
\(x >0\), pričom táto voľba
linkovacej funkcie sa používa špeciálne pri modelovaní závislej
premennej, o ktorej predpokládame, že ma podmienene na vysvetľovaných
premenných Poissonovo rozdelenie s daným parametrom – jedná sa teda o
špcifickú verziu GLM modelu pre modelovanie tzv. Poissonových počtov
(t.j. ‘’Poisson counts data’’).
Už bolo spomínané, že v prípade GLM modelov pre korelované/opakované
(longitudinálne) data \(\{(Y_{ij},
\boldsymbol{X}_{ij});i = 1, \dots, n; j = 1, \dots, m_i\}\)
merané na \(n \in \mathbb{N}\)
nezávislých subjektov (pričom celkový počet pozorovaní je \(N = \sum_{i = 1}^n m_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 m_i})^\top
\] a variančnu-kovariačnu maticuu \(\mathcal{V}_i = Var \boldsymbol{Y}_i\). To
predstavuje tzv. stochastickú časť modelu a vedie to 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. ´
Korelovanosť opakovaných pozorovaní je zohľadnená buď prostredníctvom
náhodného efektu (resp. náhodných efektov), teda v rožírenej
špecifikácii podmienenej strednej hodnoty, kedy uvažujeme \[
\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),
\] pre náhodné efekty \(\boldsymbol{b}_i \sim N_{q}(\boldsymbol{0},
\mathbb{D})\), alebo v prostredníctvom špecifikácie štruktúry
variančnej-kovariančnej matice \(\mathcal{V}_i\), ktorú je ale na rozdiel
GLM modelov pre nekorelované pozorovania nutné odhadovať navyše a robí
sa to pomocou tzv. matice pracovných korelácii (resp. tzv.
``working correlation matrix’’).
V následujúcej ilustračnej časti budeme využívať mierne upravený
datový súbor pacientov so sklerózou multiplex, kde veličina EDSS
(Expanded DisabiLity Status Scale) je zaokrúhlená na celé čísla a pre
jednoduchosť budeme predpokládať, že táto veličina má Poissonovo
rozdelenie s konkrétnym parametrom, ktory závisí na vysvetľujúcich
premenných.
Podmienenú strednú hodnotu – parameter \(\lambda_{ij} > 0\) – teda modelujeme
pomocou linkovej funkcie \(g(x) =
log(x)\), teda dostávame \[
\lambda_{i j} = E\Big[Y_{i j} | \boldsymbol{X}_{ij}\Big]= \exp \{
\boldsymbol{X}_{ij}^\top\boldsymbol{\beta}\},
\] pre marginálny model (kde bude korelácia medzi pozorovaniami
ošetrená pomocou tzv. ‘’working correlation matice’’), resp. \[
\lambda_{i j} = E\Big[Y_{i j} | \boldsymbol{X}_{ij},
\boldsymbol{b}_i\Big]= \exp \{
\boldsymbol{X}_{ij}^\top\boldsymbol{\beta} +
\boldsymbol{Z}_{ij}\boldsymbol{b}_i\},
\] pre hierarchický (t.j., podmienený) model s náhodnými efektami
\(\boldsymbol{b}_i \sim N_q(\boldsymbol{0},
\mathbb{D})\), pričom pre \(Y_{ij}\) (t.j., \(j\)-te opakované pozorovanie v rámci \(i\)-tého subjektu) predpokládame platnosť
(marginálneho) modelu \[
Y_{ij} | \boldsymbol{X}_{ij} \sim Poiss(\lambda_{ij}).
\]
Ilustrácia na datach
Príslušný (upravený) datový súbor je k dispozícii ako csv súbor sm_data4.csv. Data načítajte do programu SAS (SAS OnDemand) a porovnajte následujúce modely:
libname sm '/home/u63241636/sasuser.v94';
filename reffile '/home/u63241636/sasuser.v94/data/sm_data4.csv';
proc import datafile=reffile
dbms=csv
out=sm.data
replace;
getnames=yes;
run;
proc print data=sm.data;
run;
/* hierarchicky (podmieneny) model s nahodnym interceptom*/
proc glimmix data=sm.data;
class id gender;
model EDSS = gender age EDSSini / dist=poisson solution;
random intercept / subject=id solution vcorr;
run;
/* marginalny model s pracovnou korelacnou maticou typu 'exch' */
proc genmod data=sm.data;
class id gender;
model EDSS = gender age EDSSini / d=poisson;
repeated subject = id / corrw mcorrb type=exch;
run;
/* marginalny model s pracovnou korelacnou maticou typu 'exch' */
proc gee data=sm.data;
class id gender;
model EDSS = gender age EDSSini / dist=poisson;
repeated subject=id / type=exch corrw mcorrb ;
run;V prípade GLM regresie pre Poissonové počty je v určitom zmysle (až na intercept parameter) interpretácia odhadnutých parametrov v hierarchickom modeli a v marginálnom modeli totožná. Rozdiel medzi modelmi sa prejaví pouze v interpretácii intercept parametru.
PROC GLIMMIX?
PROC MCMC (tzv. Monte Carlo Markov
Chain), ktorá Bayesovský model umožňuje odhadovaŤ (více podrobnosti
napr. v
tomto
článku).
PROC NLMIXEDPre 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=1 beta1=1 beta2=1;
lambda = exp(u + beta0 + beta1*age + beta2*EDSSini);
model EDSS ~ poisson(lambda);
random u ~ normal(0,10) subject=id;
run;Je dôležité porozumieť základným rozdielom medzi modelmi, kde
modelujeme nelineárne strednú hodnotu náhodnej veličiny \(Y_{ij}\) napr. model implementovaný pomocou
procedúry PROC NLMIXED vyššie) a strednú hodnotu
transformovanej náhodnej veličiny \(\log(Y_{ij})\) pomocou lineárneho
regresného modelu a klasickej procedúry PROC MIXED.