Letný semester 2024 | Cvičenie 3 | 11.03.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ályKorelované (resp. závislé) opakované (longitudinálne) pozorovania lze často vhodne analyzovať aj klasickými štatistickými metódami a postupmi, ktoré sú primárne určené pre analýzu nezávislých a rovnako rozdelených pozorovaní (náhodný výber – tzv. \(i.i.d.\) pozorovania).
Takúto ``jednoduchú’’ analýzu je ale nutné spracovať korektným postupom. Niektoré základné metódy v tomto smere a tiež ich konkrétnu aplikáciu na reálne data v programe SAS ilustrujeme nižšie.
V prvom rade je ale nutné správne pochopenie celkovej štruktúry
longitudinánlych dat. Ich formálna (matematická) reprezentácia môže byť
ale rôzna a nemusí byť jednoznačná a nie je ani ekvivalentná.
Už bolo spomenuté, že v terminológii longitudinálnych dat a ich
analýzy sa bežne používajú pojmy ako long-data, resp.
wide-data.
Z matematického hľadiska by bolo možne oba typy
reprezentovať nasledujúcim spôsobom.
Predpokládajme, že sledujeme opakovane (v čase) \(n \in \mathbb{N}\) vzájomne nezávislé subjekty a zaujíma nás konkrétna odozva – t.j. napr. náhodná veličina \(Y\). Jednotlivé (realizované) porozorvania pre každý subjekt môžeme reprezentovať ako vektor \(\boldsymbol{Y}_{i} = (Y_{i j}, \dots, Y_{i m_i})^\top\), kde \(m_i \in \mathbb{N}\) je počet opakovaných pozorovaní vrámci \(i\)-teho subjektu. Pre celkový počet pozorovaní \(N \in \mathbb{N}\) (počet sledovaných subjektov je značený ako \(n \in \mathbb{N}\)) samozrejme platí, že \[ N = \sum_{i = 1}^n m_i. \] Celkový vektor pozorovaní (t.j. long-data) je potom možné reprezentovať združene prostredníctvom výsledného vektoru \[ \boldsymbol{Y} = \Big( Y_{1 1}, \dots, Y_{1 m_1}, Y_{2 1}, \dots, Y_{n m_n}\Big)^\top = \Big(\boldsymbol{Y}_1~^\top, \dots, \boldsymbol{Y}_{n}~^\top\Big)^\top. \] Ak su naviac longitudinálne data balancované, resp. pre každý subjekt je odsledovaný rovnaký počet pozorovaní (niekedy sa navyše ešte predpokladá, že jednotlivé pozorovania sú uskutočnené v rovnakých časových okamžikoch), t.j. \(m _i = m \in \mathbb{N}\) pro \(\forall i \in \{1, \dots, n\}\) tak ``data’’ lze reprezentovat i v maticovom zápise (t.j. wide-data) ako \[ \mathbb{Y} = \left( \begin{array}{ccc} Y_{1 1} & \dots & Y_{1 m_1}\\ Y_{2 1} & \dots & Y_{2 m_2}\\ \ldots & \dots & \ldots\\ Y_{n 1} & \dots & Y_{n m_n} \end{array} \right) = \left( \begin{array}{c} \boldsymbol{Y}_1~^\top\\ \boldsymbol{Y}_2~^\top\\ \vdots \\ \boldsymbol{Y}_n~^\top\\ \end{array} \right). \] Jeden aj druhý zápis majú svoje zjavné výhody aj nevýhody (resp. obmedzenia). Poskúste sa výhody aj nevýhody oboch uvedených reprezentácii stručne, ale explicitne popísať. Prípadné výhody a nevýhody sú evidentné hlavnne pri fomulácii regresného modelu pre jednotlivé reprezentácie.
Regresné rozšírenie pre vektor \(\boldsymbol{Y}\)
Predchádzajúci model lze rozšíriť v zmysle regresného konceptu. Pre slcový vektor závislých pozorovaní \(\boldsymbol{Y} = \Big(\boldsymbol{Y}_1~^\top, \dots, \boldsymbol{Y}_{n}~^\top\Big)^\top\) a regresnú maticu \(\mathbb{X} \in \mathbb{R}^{N \times p}\) môžeme uvažovať regresný model v tvare \[ \boldsymbol{Y} = \mathbb{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \] kde \(\boldsymbol{\beta} \in \mathbb{R}^p\) je \(p\)-rozmerný (neznámy) vektor parametrov. Keďže je nutné správne zohľadniť variačnú a kovariačnú štruktúru pozorovaní medzi jednotlivými (nezávislými) subjektami a korelovanosť, resp. závislosť opakovaných meraní vrámci jedného subjektu, je potrebné pre vektor náhodných chýb predpokladať vhodné rozdelenie – napr. mnohorozmerné normálne rozdelenie \[ \boldsymbol{\varepsilon} = \left( \begin{array}{c} \varepsilon_{1 1}\\ \vdots \\ \varepsilon_{1 m_1}\\ \varepsilon_{2 1}\\ \vdots \\ \varepsilon_{n m_n} \end{array} \right) \sim N_{N} \left( \left( \begin{array}{c} 0\\ \vdots \\ 0 \end{array} \right), \left( \begin{array}{c} \Sigma_1 & 0 & \dots & 0\\ 0 & \Sigma_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \Sigma_n \end{array} \right) \right), \]
kde covariančna štruktúra (nulové prvky – resp. nulové matice
vhodných rozmerov) reprezentujú vzájomnú nezávislosť medzi jednotlivými
subjektvami a príslušné (pozitívne definitné) variačné kovariačné matice
\(\Sigma_{i} \in \mathbb{R}^{m_i \times
m_i}\) pre \(i \in \{1, \dots,
n\}\) sú obecně v tvare \[
\Sigma_{i} = \left(
\begin{array}{c}
\sigma_{1 1}^{(i)} & \sigma_{1 2}^{(i)} & \dots & \sigma_{1
m_i}^{(i)}\\
\sigma_{2 1}^{(i)} & \sigma_{2 2}^{(i)} & \dots & \sigma_{2
m_i}^{(i)}\\
\vdots & \vdots & \ddots & \vdots \\
\sigma_{m_i 1}^{(i)} & \sigma_{m_i 2}^{(i)} & \dots &
\sigma_{m_i m_i}^{(i)}\\
\end{array}
\right).
\]
Regresné rozšírenie pre maticu \(\mathbb{Y}\)
V prípade odozvy reprezentovanej v maticovom tvare je zápis modelu trochu odlišny. Analogický model, ako v predchádzajúcom prípade, môžeme formálne zapísať ako \[ \mathbb{Y} = \mathbb{X}\boldsymbol{B} + \boldsymbol{E}, \] kde tentokrát \(\mathbb{X} \in \mathbb{R}^{n \times p}\) a namiesto vektoru neznámych parametrov \(\boldsymbol{\beta} \in \mathbb{R}^p\) (ako v predchádzajúcom prípade) reprezentuje \(\boldsymbol{B} \in \mathbb{R}^{p \times m}\) maticu neznámych parametrov a \(\boldsymbol{E}\) je matica chybových členov, pričom platí (uvažujúc opäť normálny model), že \(\boldsymbol{E} = (\boldsymbol{\varepsilon}_{1}^\top, \dots, \boldsymbol{\varepsilon}_{n}^\top)^\top\), pričom náhodné vektory \(\boldsymbol{\varepsilon}_i = (\varepsilon_{i 1}, \dots, \varepsilon_{i m})^\top\) pre \(i = 1, \dots, n\) tvoria náhodný výber z mnohorozmerného normálneho rozdelenia s nulovým vektorom stredných hodnot a variačnou kovaričnou maticou \[ \Sigma = \left( \begin{array}{c} \sigma_{1 1} & \dots & \sigma_{1 m}\\ \vdots & \ddots & \vdots\\ \sigma_{m 1} & \dots & \sigma_{m m} \end{array} \right). \]
Najednoduchším príkladom štatistickej analýzy longitudinálnych dat
môže byť klasický jednovýberový párový \(t\)-test. Ten je čosto možné veľmi
efektívne aplikovať aj v prípade, že longitudinálne data, ktoré je
potrebné analyzovať, obsahujú vrámci jednotlivých subjektov výrazne
viacej, ako len dva opakované pozorovania. Keďže správna (korektná)
štatistická analýza longitudinálnych dat je pomerne zložitá, v praxi sa
celkom bežne používajú niektoré jednoduchšie postupy, ktorých cieľom je
určité zjednodušenie štruktúry dat a následná štatistická analýza
nezávislých (\(i.i.d\)) pozorovaní.
V následujúcej časti budeme pre ilustráciu uvažovať csv súbor) s datami o pacientoch so sklerózou multiplex. Zameriame sa hlavne na implementáciu párového \(t\)-testu v programe SAS. Ako závislú premennú budeme uvažovať veličinu EDSS (t.j., Expanded Disability Status Scale) .
Data načítame v SAS OnDemand a náslene upravíme do vhodnejšieho tvaru (pre aplikáciu pároveho \(t\)-testu) pomocou následujúcich príkazov:
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 print datafile = sm.data;
run;
proc sort data=sm.data; by id; run;
proc transpose data=sm.data out=sm.dataWide;
by id;
var EDSS;
run;
proc print datafile = sm.dataWide;
run;Pre ilustráciu budeme pomocou párového \(t\)-testu analyzovať zmenu EDSS v období piatich rokov po aplikácii novej liečby – Lemtrady. Implementácia párového \(t\)-testu v programe SAS je následovná:
proc ttest data=sm.dataWide alpha=0.05;
paired COL5*COL1;
run;Pre porovnanie, aplikujeme ešte klasický dvojvýberový \(t\)-test a vzájomne porovnáme výsledky. V
prvom rade je potrebné upraviť data vo forme dvoch (nezávislých)
náhodných výberov. K tomuto účelu z datového súboru
sm.dataWide vyberieme pouze relevantne premenné (t.j.
stĺpce COL1 a COL5):
data sm.data15;
set sm.dataWide;
keep COL1 COL5;
run;
proc print datafile = sm.data15;
run;Vytvorili sme nový datový súbor s názvom sm.data15.
Následne hodnoty z oboch stĺpcov združíme do jedného a príslušné riadky
označíme poradím merania:
data sm.twoSample;
length Sample $284.;
set sm.data15;
array numvars[*] _NUMERIC_;
do i = 1 to dim(numvars);
Sample = vname(numvars[i]);
EDSS = numvars[i];
output;
end;
keep Sample EDSS;
run;
proc print datafile = sm.twoSample;
run;Následne je už možné priamo aplikovať dvojvýberový \(t\)-test na datový súbor
sm.twoSample. Konkrétna implementácia v programe SAS je
následujúca:
proc ttest data=sm.twoSample alpha=0.05;
class Sample;
var EDSS;
run;Zakladné rozdiely medzi
long-data'' reprezentaciou (vektor odozvy $\boldsymbol{Y}$) awide-data’’
reprezentáciou (matica odozvy \(\mathbb{Y}\)) vrámci regresných modelov
formulovaných vyššie, ilustrujeme aj na reálnych datach. Pre tento účel
využijeme datový súbor cars.csv (data lze stiahnúť ako
csv súbor).
Následne data načítame do programu SAS a pomocou procedúry
print sa na data môžeme podívať.
libname car '/home/u63241636/sasuser.v94';
filename reffile '/home/u63241636/sasuser.v94/data/cars.csv';
proc import datafile=reffile
dbms=csv
out=car.data
replace;
getnames=yes;
run;
proc print datafile = car.data;
run;Budeme uvažovať jednoduchhý príklad, kde odozva bude dvojrozmerná a bude predstavovať dve ceny – výrobcom navrhnutú cenu (t.j., Manufacturer’s suggested retail price – MSRP) (premenná MSRP) a skutočnú cenu, za ktorú bolo auto zakúpené (premenná Invoice).
Obe uvedené premenné sú ale v datovom súbore reprezentované ako charakterové hodnoty. Pre účely regresného modelu je potrebné premenné vhodne upraviť, aby nové premenné predstavovali kvantitatívne hodoty (ceny).
To lze dosiahnúť napr. pomocou následujúceho kódu:
data car.data2;
set car.data;
MSRP2 = input(compress(translate(MSRP,"","$,")),best.);
Invoice2 = input(compress(translate(Invoice,"","$,")),best.);
run;
proc print datafile = car.data2;
run;Pomocou procedúry proc reg alebo procedúry
proc glm môžeme fitovať požadovaný regresný model:
proc reg data=car.data2;
model MSRP2 Invoice2 = EngineSize Horsepower;
run;
proc glm data=car.data2;
model MSRP2 Invoice2 = EngineSize Horsepower;
run;Modely vyššie porovnajte s následujúcimi modelmi:
proc reg data=car.data2;
model MSRP2 = EngineSize Horsepower;
run;
proc reg data=car.data2;
model Invoice2 = EngineSize Horsepower;
run;Pomocou nasledujúceho kódu môžeme zároveň otestovať napríklad nulovú hypotézu, že regresná priamka pre MSRP cenu a regresná priamka pre skutočnú cenu (Invoice) sú totožné.
proc reg data=car.data2;
model MSRP2 Invoice2 = EngineSize Horsepower;
mtest intercept, EngineSize, Horsepower, MSRP2-Invoice2;
run;Ako by analogicky štatistický test vyzeral a bol implementovaný v
programe SAS v prípade, že by sa pozorovania reprezetovali v tzv.
long-data type (t.j., pomocou vektoru \(\boldsymbol{Y}\) )?
Analogický a v určitom zmysle aj ekvivalentný model dostaneme ale aj pomocou tzv. long-data formatu, ked vhodne zostavíme regresný model pre vektor závislých pozorovaní \(\boldsymbol{Y}\).
Pre prehľadnosť vyberieme len relevantné premenné:
data car.dataLong;
set car.data2;
keep MSRP2 Invoice2 EngineSize Horsepower;
run;
proc print datafile = car.dataLong;
run;Data následne upravíme do vhodného formátu:
data car.twoSample;
length Sample $856.;
set car.dataLong;
array numvars[2] MSRP2 Invoice2;
do i = 1 to dim(numvars);
Type = i;
Price = numvars[i];
output;
end;
run;
proc print datafile = car.twoSample;
run;A pre účely implementácie potrebných interačných členov prípravíme interakčné premenné:
data car.twoSampleInt;
set car.twoSample;
Type2 = Type - 1;
EngineType = Type2 * EngineSize;
HorseType = Type2 * Horsepower;
run; Výsledný regresný model získame následujúcim spôsobom:
proc reg data=car.twoSampleInt;
model Price = Type2 EngineSize Horsepower EngineType HorseType;
run;Porovnajte odhadnuté parametre, príslušné smerodatné chyby a \(p\)-hodnoty príslušných štatistických testov v regresnom modeli, ktorý bol založený na wide-data type (na matici závislých pozorovaní \(\mathbb{Y}\)) a modeli, ktorý bol vybudovaný pre long-data typ (teda vektor \(\boldsymbol{Y}\) ).