Letný semester 2024-2025 | Cvičenie 6 | 31.03.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.
Základný lineárny regresný model s náhodnými efektami je (typicky) vyjadrený pomocou rovnice
\boldsymbol{Y} = \mathbb{X}\boldsymbol{\beta} + \mathbb{Z}\boldsymbol{b} + \boldsymbol{\varepsilon}, kde \boldsymbol{Y} = (Y_{11}, \dots, Y_{i n_1}, Y_{21}, \dots, Y_{Nn_N})^\top \in \mathbb{R}^\mathcal{N} predstavuje celkový vektor závislej premennej Y nameranej jednak pre N \in \mathbb{N} nezávislých subjektov, ale pre ka6d7 subjekt je postupne uvažovaných n_i \in \mathbb{N} opakovaných pozorovaní (pre i \in \{1, \dots, N\}). Celkový počet pozorovaní (dĺžka vektoru \boldsymbol{Y}) je teda \mathcal{N} = \sum_{i = 1}^N n_i.
Regresná matica \mathbb{X} je
typu '\mathcal{N} \times p pre
vektor neznámych prametrov (pevných efektov) \boldsymbol{\beta} \in \mathbb{R}^p a
maticu \mathbb{Z} \in
\mathbb{R}^{\mathcal{N} \times Nq}, ktorá prislúcha náhodnym
efektom \boldsymbol{b} =
(\boldsymbol{b}_1^\top, \dots, \boldsymbol{b}_N^\top)^\top \in
\mathbb{R}^{Nq}, kde \boldsymbol{b}_i = (b_{i 1}, \dots, b_{i q})^\top
\in \mathbb{R}^q reprezentuje tzv. ``subject-specific’’ náhodné
efekty pre každý subjekt, t.j. pre každé i
\in \{1, \dots, N\}. Všimnite si, že dimenzia (počet) náhodných
efektov je pre každý subjekt rovnaká, t.j. q
\in \mathbb{N}.
Lineárny regresný model s náhodnými efektami môžeme vyjadriť aj v
tzv. “subject-specific” formulácii, teda pomocou rovnice
\boldsymbol{Y}_i = \mathbb{X}_i\boldsymbol{\beta} +
\mathbb{Z}_i\boldsymbol{b}_i + \boldsymbol{\varepsilon}_i,
\boldsymbol{Y}_i = (Y_{i 1},
\dots, Y_{i n_i})^\top je vektor opakovaných pozorovaní pre
daný subjekt i \in \{1, \dots, N\}
(ktorý ale môže byť obecne rôznej dĺžky pre rôzne subjekty i). Regresná matica je typu n_i \times p (t.j. vektor neznámych
parametrov – pevných efektov \boldsymbol{\beta}\in \mathbb{R}^p má pre
všetky subjekty rovnakú dĺžku (resp. dimenziu) a jednotlivé jeho zložky
majú rovnakú interpretáciu). Regresná matica \mathbb{Z}_i je typu n_i \times q a podobne ako regresná
matica \mathbb{X}_i závisí na
indexe i \in \{1, \dots, N\} (t.j.,
``subject-specific’’ informácia).
Tzv. sériova korelácia úmožňuje modelovať individuálny profil daného subjektu ako určitý časovo závislý (nekonečne-rozmerný) stochastický proces, ktorý ale sledujeme len prostredníctvom jednotlivých (konečne mnoho) opakovaných pozorovaní nameraných v rámci daného subjektu. Korelácia medzi jednotlivými dvoma časovými okamžikmi stochastického procesu v rámci jedného subjekty by mala byť preto klesajúca funkcia vzhľadom k času, ktorý dané dva okamžiky oddeľuje/separuje. Príslušná korelačné matica sa preto často parametrizuje a jednotlivé prvky majú často tvar g(|t_{i j} - t_{i k}|), kde g(\cdot) je nejaká klesajúca (parametrická) funkcia (taková, že platí g(0) = 1) a t_{ij}, t_{i k} \in [0,T] sú okamžíky dvoch konkrétnych meraní uskutočnených pre subjekt i \{1, \dots, N\}.
V programe SAS slúži ako základný nástroj na fitovanie lineárnych
regresných modelov s náhodnými efektami procedúra
PROC MIXED
.
V programe R je možné využíť napr. funkciu
lmer()
z knižnice lme4.
install.packages("lme4")
library("lme4")
?lmer
V následnej časti využijeme opäť datový súbor s pacientami so sklerózou multiplex a pozrieme sa na porovnanie modelov odhadnutých pomocou programu R a pomocou programu SAS. Porovajte následujúce dva výstupy:
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;
data sm.data2;
set sm.data;
timeCls = time;
run;
proc mixed data = sm.data2 method = ml;
class gender(ref = "F") timeCls;
model EDSS = gender time*gender / s corrb;
repeated timeCls / subject = id;
random intercept / subject = id;
run;
sm <- read.csv(url("https://www2.karlin.mff.cuni.cz/~maciak/NMST422/sm_data2.csv"), header = T)
summary(lmer(EDSS ~ gender*time + (1|id), data = sm, REML = F))
Podobne porovnajte aj následujúce dva výstupy/modely:
proc mixed data = sm.data2 method = ml;
class gender timeCls;
model EDSS = gender time*gender / s corrb;
repeated timeCls / type = vc subject = id;
random intercept time / subject = id;
run;
options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))
summary(lmer(EDSS ~ gender*time + (1 + time || id), data = sm, REML = F))
lmer()
a k SAS
procedúre PROC MIXED
.
PROC MIXED
– konkrétne význam tzv. “statements” (MODEL,
RANDOM a REPEATED) v súvislosti s explicitnym matematickým zápisom
modelu vyššie.
Predchádzajúce modely (aj v programe R, aj v programe SAS) boli odhadnuté pomocou metódy maximálnej vierohodnosti. Nejedná sa ale o štandardnú metódu, ktoré sú defaultne používané. Ak budeme uvažovať marginálny model \boldsymbol{Y}_i \sim N_{n_i}\Big(\mathbb{X}_i\boldsymbol{\beta}, \mathbb{Z}_i\mathbb{D}\mathbb{Z}_i^\top + \boldsymbol{\Sigma}_i \Big). tak metóda maximálnej vierohodnosti za predpokladu normálneho modelu dáva vierohodnostnú rovnicu L(\mathcal{X}_N, \boldsymbol{\theta}) = \prod_{i = 1}^N \Big[ (2\pi)^{-n_i/2} |\mathbb{V}_i(\boldsymbol{\alpha})|^{-1/2} exp\Big\{ -\frac{1}{2}(\boldsymbol{Y}_i - \mathbb{X}_i\boldsymbol{\beta})^\top \mathbb{V}_i^{-1}(\boldsymbol{\alpha}) (\boldsymbol{Y}_i - \mathbb{X}_i\boldsymbol{\beta})\Big\}\Big], kde \mathcal{X}_N predstavuje namerané data \{(\boldsymbol{Y}_i, \mathbb{X}_i);~i = 1, \dots, N\}, pevné efekty sú reprezentované vektorom \boldsymbol{\beta} \in \mathbb{R}^p a vektor \boldsymbol{\alpha} \in \mathbb{R}^q reprezentuje všetky neznáme parametre v rámci uvažovanej variančno-kovariančnej štruktúry \mathbb{Z}_i\mathbb{D}\mathbb{Z}_i^\top + \boldsymbol{\Sigma}_i.
Je dobré si uvedomiť, že sa vpodstate jedná len o drobné zobecnenie problému, kde sledujeme mnohorozmerný náhodný výber \boldsymbol{Y}_1, \dots, \boldsymbol{Y}_N z mnohorozmerného normálneho rozdelenia N_{n}(\boldsymbol{\mu}, \Sigma), s nejakou obecnou, pozitívne-definitnou variančnou-kovariančnou maticou \Sigma \in \mathbb{R}^{n \times n} (kde b \in \mathbb{N} predstavuje počet opakovaných meraní v rámci každého uvažovaného subjektu).
Ak by boli parametre v \boldsymbol{\alpha} známe, tak odhad
neznámeho vaktoru \boldsymbol{\beta} \in
\mathbb{R}^p pomocou metódy maximálnej vierohodnosti by sa
redukoval na
\widehat{\boldsymbol{\beta}} = \Big(\sum_{i = 1}^n \mathbb{X}_i
\mathbb{W}_i \mathbb{X}_i\Big)^{-1} \cdot \Big( \sum_{i = 1}^n
\mathbb{X}_i\mathbb{W}_i\boldsymbol{Y}_i \Big),
kde pre jednoduchosť \mathbb{W}_i
= \mathbb{V}_i^{-1}(\boldsymbol{\alpha}). Vo väčšine
aplikovaných prípadoch je ale vektor \boldsymbol{\alpha} neznámy a tým pádom
je neznáma aj variančná-kovariančná štruktúra \mathbb{Z}_i\mathbb{D}\mathbb{Z}_i^\top +
\boldsymbol{\Sigma}_i. Preto je nutné nezáme parametre v \boldsymbol{\alpha}, ktoré
definujú/určujú \mathbb{Z}_i\mathbb{D}\mathbb{Z}_i^\top +
\boldsymbol{\Sigma}_i nejak odhadnúť. K tomuto účelu sa
najčastejšie používa niektorý z následujúcich postupov:
Restricted maximum likelihood je metóda pre odhadovanie
variančného parametru (resp. variančnej-kovariančnej matice) bez
indukovania vychýlenia (bias) na zaklade toho, že namiesto skutočnej
strednej hodnoty (resp. skutočného vektoru stredných hodnôt) pracujeme
len so stochastickým empirickým odhadom – t.j. do odhadovacej procedúry
sa vnáša dodatočná neurčitosť.
REML v linárnom regresnom modeli
Analogický problém sa objavuje aj v prípade (štandardného) lineárneho
regresného modelu, kde odhad rozptylu v tvare
\widehat{\sigma}^2 = (\boldsymbol{Y} -
\mathbb{X}\boldsymbol{\beta})^\top(\boldsymbol{Y} -
\mathbb{X}\boldsymbol{\beta})/N
je obecne vychýlený a nevychýlený odhad je až odhad definovaný
ako \frac{1}{N - p} (\boldsymbol{Y} -
\mathbb{X}\boldsymbol{\beta})^\top(\boldsymbol{Y} -
\mathbb{X}\boldsymbol{\beta}), kde p \in \mathbb{N} je počet nezámych
parametrov vo vektore \boldsymbol{\beta} \in
\mathbb{R}^p (resp. dimenzia/dĺžka vektoru \boldsymbol{\beta}). Restricted maximum
likelihood (REML) metóda vpodstate opäť hľadá vhodnú transformáciu
pôvodných dat \boldsymbol{Y}
pomocou vhodnej matice \Delta, tak
aby E[\Delta^\top \boldsymbol{Y}] =
\boldsymbol{0} (t.j., transformácia dat do priestoru kolmého na
priestor generovaný stĺpcami matice \mathbb{X}), čím sa vlastne vytratí
závislosť transformovaných dat \widetilde{\boldsymbol{Y}} = \Delta^\top
\boldsymbol{Y} na neznámom parametri (parametroch) strednej
hodnoty \mathbb{X}\boldsymbol{\beta}.
Transformované data majú opäť mnohorozmerné normálne rozdelenie a navyše
platí, že
\Delta^\top \boldsymbol{Y} \sim N_{N - p}(\boldsymbol{0},
\sigma^2\Delta^\top\Delta).
Keďže obecne predpokládame regulárnu maticu \mathbb{X} typu N \times p (ktorej stĺpce generujú
linárny podpriestor v \mathbb{R}^N,
ktorého dimenzia je p \in
\mathbb{N}), tak príslušná transfomácia do priestoru kolmého na
stĺpce matice \mathbb{X} je (N - p) rozmerný linárny podpriestor v
\mathbb{R}^N (ktorého projekčná
matica je napr. matica \mathbb{M} =
(\mathbb{I} -
\mathbb{X}(\mathbb{X}^\top\mathbb{X})^{-1}\mathbb{X}^\top)).
REML v linárnom regresnom modeli s náhodnými
efektami
V prípade lineárneho regresného modelu s náhodnými efektami môžeme formulovať jednak model pre jednotlivé subjekty \boldsymbol{Y}_i \sim N_{n_i} (\mathbb{X}_i\boldsymbol{\beta}, \mathbb{V}_i). alebo celkový model pre všetkých N \in \mathbb{N} subjektov dohromady (celkovo teda \mathcal{N} = \sum_{i = 1}^N n_i pozorovaní), teda \boldsymbol{Y} \sim N_{\mathcal{N}} (\mathbb{X}\boldsymbol{\beta}, \mathbb{V}(\boldsymbol{\alpha})), kde \mathbb{X} = (\mathbb{X}_1^\top, \dots, \mathbb{X}_N^\top)^\top. Variančné matice \mathbb{V}_i pre i = 1, \dots, N a \mathbb{V} závisia na neznámych parametroch, ktoré sú súhrnne označené ako \boldsymbol{\alpha} \in \mathbb{R}^q
Idea REML je odhadnúť neznáme parametre \boldsymbol{\alpha} bez potreby prvotného
odhadovania neznámych parametrov v \boldsymbol{\beta} \in \mathbb{R}^p.
Pôvodné data \boldsymbol{Y} \in
\mathbb{R}^\mathcal{N} je potrebné opať transformovať
ortogonálne vzȟladom na stĺpce matice \mathbb{X} pomocou vhodnej matice \Delta, tak aby
\Delta^\top \boldsymbol{Y} \sim N_{\mathcal{N} - p}(\boldsymbol{0},
\Delta^\top \mathbb{V}(\boldsymbol{\alpha})\Delta)
a násedne aplikovať metódou maximálnej vierohodnosti na
transformované data \Delta^\top
\boldsymbol{Y}, čím sa získa odhad neznámych parametrov \boldsymbol{\alpha}. Tento odhad sa
nazýva restricted maximum likelihood (REML) odhadnom a často sa
v literatúre označuje aj ako \widehat{\alpha}_{REML}.
V programe R aj v programe SAS odhadneme lineárne regresné modely pomocou metódy maximálnej vierohodnosti a následne pomocou restricted maximum likelihood. Jednotlivé výstupy porovnáme.
proc mixed data = sm.data2 method = ml;
class gender timeCls;
model EDSS = gender time*gender / s;
repeated timeCls / subject = id;
random intercept / subject = id;
run;
proc mixed data = sm.data2 method = reml;
class gender timeCls;
model EDSS = gender time*gender / s;
repeated timeCls / subject = id;
random intercept / subject = id;
run;
summary(lmer(EDSS ~ gender*time + (1|id), data = sm, REML = F))
summary(lmer(EDSS ~ gender*time + (1|id), data = sm))
Jedná z (asi mnohých) možností, ako v programe SAS získať variogram,
je využit procedúru PROC VARIOGRAM
– viď napr. tento
SAS
tutoriál.
V následujúcej časti SAS kódu sa pokúsíme variogram zostrojiť manuálne v postupnosti niekoľkých krokoch. Využijeme k tomu datový súbor s pacientami so sklerózou multiplex. Výstupom bude priemerný variogram spočítaný pre všetky uvažované subjekty.
proc sort data=sm.data2;
by id time;
run;
/* Create dataset with lagged differences */
data variogram;
set sm.data2;
by id;
/* Compute lags for EDSS */
lag0 = EDSS; /* Lag 0 (same time point) */
lag1 = lag1(EDSS); /* Lag 1 */
lag2 = lag2(EDSS); /* Lag 2 */
lag3 = lag3(EDSS); /* Lag 3 */
lag4 = lag4(EDSS); /* Lag 4 */
/* Compute squared differences */
if _n_ > 0 then do;
semi0 = (EDSS - lag0)**2 / 2; /* Lag 0 semivariance */
semi1 = (EDSS - lag1)**2 / 2; /* Lag 1 semivariance */
semi2 = (EDSS - lag2)**2 / 2; /* Lag 2 semivariance */
semi3 = (EDSS - lag3)**2 / 2; /* Lag 3 semivariance */
semi4 = (EDSS - lag4)**2 / 2; /* Lag 4 semivariance */
end;
/* Keep only relevant values */
keep SubjectID time semi0 semi1 semi2 semi3 semi4;
run;
proc print data = variogram;
run;
/* Compute mean semivariance for each lag */
proc means data=variogram mean;
var semi0 semi1 semi2 semi3 semi4;
output out=semivariances mean=semi0_mean semi1_mean semi2_mean semi3_mean semi4_mean;
run;
/* Reshape data for plotting */
data semivariogram;
set semivariances;
length Lag 8 Semivariance 8;
/* Convert to long format */
Lag = 0; Semivariance = semi0_mean; output;
Lag = 1; Semivariance = semi1_mean; output;
Lag = 2; Semivariance = semi2_mean; output;
Lag = 3; Semivariance = semi3_mean; output;
Lag = 4; Semivariance = semi4_mean; output;
keep Lag Semivariance;
run;
/* Plot the empirical variogram */
proc sgplot data=semivariogram;
series x=Lag y=Semivariance / markers lineattrs=(thickness=2);
scatter x=Lag y=Semivariance / markerattrs=(symbol=circlefilled size=10);
xaxis label="Time Lag";
yaxis label="Semivariance";
title "Empirical Variogram for EDSS (Pooled Across Subjects)";
run;
Alternatívna možnosť je podívať sa na individuálne variogrami jednotlivých subjektov. Tie v programe SAS vykreslíme do grafu napr. pomocou následujúceho kódu:
/* Reshape data for plotting */
data semivariogram;
set variogram;
length Lag 8 Semivariance 8;
/* Convert to long format */
Lag = 0; Semivariance = semi0; output;
Lag = 1; Semivariance = semi1; output;
Lag = 2; Semivariance = semi2; output;
Lag = 3; Semivariance = semi3; output;
Lag = 4; Semivariance = semi4; output;
keep id Lag Semivariance;
run;
/* Plot individual variograms */
proc sgplot data=semivariogram;
series x=Lag y=Semivariance / group=id markers lineattrs=(thickness=1);
xaxis label="Time Lag";
yaxis label="Semivariance";
title "Empirical Variogram for EDSS (Individual Profiles)";
run;
Použijte vhodné longitudinálne data podľa vlastného výberu (napr. datový súbor o pacientoch so sklerózou multiplex) a pomocou programu SAS urobťe následujúce:
Riešenie si priprave vo svojom účte na SAS OnDemand na výuku cvičenia v pondelok, 28.04.2025.