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ályV 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.
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).
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.
nmle
, V rámci tejte knižnice je k dispozícii
funkcia Variogram()
, ktorá počíta výberový variogram, resp.
semi-variogram.
variogram()
dostupná v R knižnici gstat
(install.packages("gstat")
)
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}.
\]
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;
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;
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.
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.
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;
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.
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:
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;
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).