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")
?lmerV 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.