Letný semester 2023-2024 | Cvičenie 6 | 15.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ályZákladný lineárny regresný model s náhodnými efektami je vyjadrený pomocou rovnice
\[ \boldsymbol{Y} = \mathbb{X}\boldsymbol{\beta} + \mathbb{Z}\boldsymbol{b} + \boldsymbol{\varepsilon}, \] kde \(\boldsymbol{Y} = (Y_{11}, \dots, Y_{i m_1}, Y_{21}, \dots, Y_{nm_n})^\top \in \mathbb{R}^N\) predstavuje vektor závislej premennej \(Y\) nameranej jednak pre \(n \in \mathbb{N}\) nezávislých subjektov, ale zároveň aj pre \(m_i \in \mathbb{N}\) opakovaných pozorovaní vrámci každého subjektu \(i \in \{1, \dots, N\}\). Celkový počet pozorovaní (dĺžka vektoru \(\boldsymbol{Y}\)) je teda \(N = \sum_{i = 1}^n m_i\).
Regresná matica \(\mathbb{X}\) je
typu \(N \times p\) pre vektor
neznámych prametrov (pevných efektov) \(\boldsymbol{\beta} \in \mathbb{R}^p\) a
maticu \(\mathbb{Z} \in \mathbb{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é \(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, kde \[
\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 m_i})^\top\) je vektor opakovaných pozorovaní pre
daný subjekt \(i \in \{1, \dots, n\}\)
(ktorý môže byť obecne rôznej dĺžky pre rôzne subjekty \(i\)). Regresná matica je typu \(m-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) a regresná matica \(\mathbb{Z}_i\) je typu \(m_i \times q\).
Tzv. sériova korelácia úmožňuje modelovať individuálny profil daného subjektu ako určitý časovo závislý stochastický proces, ktorého konkrétnú realizáciu sledujeme v jednotlivých opakovaniach nameraných vrámci daného subjektu. Korelácia medzi jednotlivými dvoma časovými okamžikmi stochastického procesu 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 by by preto na jednotlivých pozíciach mala obsahovať prvky \(g(|t_{i j} - t_{i k}|)\), kde \(g(\cdot)\) je nejaká klesajúca 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/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;
data sm.data2;
set sm.data;
timeCls = time;
run;
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; 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;
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_{m_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)^{-m_i/2} |\mathbb{V}_i(\boldsymbol{\xi})|^{-1/2} exp\Big\{ -\frac{1}{2}(\boldsymbol{Y}_i - \mathbb{X}_i\boldsymbol{\beta})^\top \mathbb{V}_i^{-1}(\boldsymbol{\xi}) (\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{\xi}\) reprezentuje všetky neznáme parametre vrámci uvažovanej variančnej 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_{m}(\boldsymbol{\mu, \Sigma})\), s nejakou obecnou, pozitívne-definitnou variančnou-kovariančnou maticou \(\Sigma \in \mathbb{R}^{m \times m}\).
Ak by boli parametre \(\boldsymbol{\xi}\) 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{\xi})\). Vo väčšine aplikovaných
prípadoch je sú ale parametre \(\boldsymbol{\xi}\) a tým pádom aj
variančná-kovariančná štruktúra \(\mathbb{Z}_i\mathbb{D}\mathbb{Z}_i^\top +
\boldsymbol{\Sigma}_i\) neznáme a je nutné parametre \(\boldsymbol{\xi}\), 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čnj-kovariančnej matice) bez indukovania
vychlenia (bias) na z8klade toho, že namiesto skutočnej strednej hodnoty
(resp. skutočného vektoru stredných hodnot) pracujeme s príslušnými
empirickými protejškami (odhadmi) – 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})\). Restricted maximum likelihood
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á transofmácia do priestoru kolmého na
stĺpce \(\mathbb{X}\) je \((n - p)\) rozmerný linárny podpriestor v
\(\mathbb{R}^n\).
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_{m_i} (\mathbb{X}_i\boldsymbol{\beta}, \mathbb{V}_i). \] alebo celkový model pre všetkých \(n \in \mathbb{N}\) subjektov dohromady, teda \[ \boldsymbol{Y} \sim N_{N} (\mathbb{X}\boldsymbol{\beta}, \mathbb{V}(\boldsymbol{\xi})), \] kde \(N = \sum_{i = 1}^n m_i\) a \(\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úhrnne označíme ako \(\boldsymbol{\xi}\) (s príslušnou dimenziou).
Idea REML je opať odhadnúť neznáme parametre \(\boldsymbol{\xi}\) bez potreby prvotného
odhadovania neznámych parametrov v \(\boldsymbol{\beta} \in \mathbb{R}^p\).
Pôvodné data \(\boldsymbol{Y} \in
\mathbb{R}^N\) opať transformujeme ortogonálne vzȟladom na stĺpce
matice \(\mathbb{X}\) pomocou matice
\(\Delta\) \[
\Delta^\top \boldsymbol{Y} \sim N_{N - p}(\boldsymbol{0}, \Delta^\top
\mathbb{V}(\boldsymbol{\xi})\Delta)
\] a metódou maximálnej vierohodnosti pomocou transformovaných
dat \(\Delta^\top \boldsymbol{Y}\)
získame odhad neznámych parametrov \(\boldsymbol{\xi}\). Tento odhad sa nazýva
“restricted maximum likelihood” odhadnom a často sa v literatúre
označuje aj ako \(\widehat{\xi}_{REML}\).
Aj 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))Použijte data o pacientoch so sklerózou multiplex (prípadne môžete využiť aj iné data, ktoré sú longitudinálneho charakteru s aspoň troma opakovanými pozorovaniami aspoň u niektorých subjektov).
Riešenie si priprave vo svojom účte na SAS OnDemand na výuku cvičenia v pondelok, 22.04.2024.