Longitudinální a panelová data – NMST422

Letný semester 2025 | Cvičenie 5 | 24.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.

Doporučená literatúra a ďalšie užitočné materiály

Diggle, P. J., Heagerty, P., Liang, K. Y., & Zeger, S. (2002). Analysis of longitudinal data. Oxford university press.
Fitzmaurice, G. M., Laird, N. M., & Ware, J. H. (2012). Applied longitudinal analysis. John Wiley & Sons.
Hardin, J.W. and Hilbe, J.M. (2007). Generalized Linear Model and Extensions. StataPress.
Pinheiro, J., & Bates, D. (2006). Mixed-effects models in S and S-PLUS. Springer science & business media.
Härdle, K.H & Šimar, L. (2015.). Applied Multivariate Statistical Analysis, Springer-Verlag Berlin.
Jednoduchý (online) SAS tutorial (english)
Základný SAS OnDemand tutorial (english)
Stručný (Český) manuál uživatele SASu na stránke doc. Kulicha
Užitočné aj neužitočné príklady zdrojových kódov v SAS (english)

V. Lineárny regresný model s náhodnými efektmi

V praxi sa často štatistík stretáva s longitudinálnymi datami, ktoré nie sú balancované (tzv., počet opakovaných pozorovaní v rámci jedného subjektu je pre rôzne subjekty rôzna). Z tohto dôvodu nie je možné aplikovať mnohorozmerné štatistické postupy (napr. tie, ktoré boli explicitne zmienené na predchádzajúcich cvičeniach). Je nutné využiť iné stochastické/pravdepodobnostné modely a iné štatistické postupy, ktoré umožnia pracovať aj s nebalancovanými longitudinálnymi pozorovaniami.

Základným štatistickým (regresným) nástrojom pre analýzu longitudinálnych (nie nutne balancovaných) dat je tzv. (lineárny) regresný model s náhodnými efektmi. Jedná sa o rozšírenie klasického lineárneho regresného modelu, definovaného (maticovo) ako \[ \boldsymbol{Y} = \mathbb{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \] kde \(\boldsymbol{Y} = (Y_1, \dots, Y_N)^\top\) predstavuje vektor (nezávislých) pozorovaní závislej premennej/veličiny \(Y\) (celkovo pre \(N \in \mathbb{N}\) rôznych, t.j., nezávislých subjektov), matica \(\mathbb{X}\) je tzv. regresná (dizajnová )matica modelu a vektor \(\boldsymbol{\beta} \in mathbb{R}^p\) predstavuje neznáme parametre, ktoré je potrebné odhadnúť. Chybový člen \(\boldsymbol{\varepsilon} = (\varepsilon_1, \dots, \varepsilon_N)^\top\) predstavuje vektor nepozorovaných náhodných chýb a väčšinou sa predpokláda, že \[ \varepsilon_i \sim N(0, \sigma^2), \] pre nejaký nezámy parameter rozptylu \(\sigma^2 > 0\). Neznáme parametre \(\boldsymbol{\beta} = (\beta_1, \dots, \beta_p)^\top\) sa často označujú aj ako pevné efekty. U lineárneho regresného modelu s náhodnými efektmi navyše vystupujú tzv. náhodné efekty a celkový (marginálny) model je možné zapísať v (maticovom) tvare ako \[ \boldsymbol{Y} = \mathbb{X}\boldsymbol{\beta} + \mathbb{Z}\boldsymbol{b} + \boldsymbol{\varepsilon}, \] avšak v tejto formulácii predstavuje \(\boldsymbol{Y} = (Y_{11}, \dots, Y_{i n_1}, Y_{21}, \dots, Y_{Nn_N})^\top \in \mathbb{R}^{\sum n_i}(\) vektor závislej premennej \(Y\) nameranej jednak pre \(N \in \mathbb{N}\) nezávislých subjektov, ale zároveň aj \(n_i \in \mathbb{N}\) opakovaných pozorovaní v rámci každého subjektu \(i \in \{1, \dots, N\}\). Celkový počet pozorovaní je teda \(\mathcal{N} = \sum_{i = 1}^N n_i\).

Regresná matica \(\mathbb{X}\) je typu \(\mathcal{N} \times p\), pre vektor neznámych parametrov (pevných efektov) platí \(\boldsymbol{\beta} \in \mathbb{R}^p\) a matica \(\mathbb{Z}\) je typu \(\mathcal{N} \times Nq\) a 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}\).

Základný princíp lineárneho regresného modelu s náhodnými efektami môže byť dobre ilustrovaný pomocou tzv. dvoj-fázoveho regresného modelu.

1. Dvojfázový regresný model pre longitudinálne data

Idea modelovania longitudinálnych dat pomocou dvojfázoveho postupu je založená na dvoch samostatných (regresných) krokoch:

Krok 1:
Lineárny regresný model pre každý subjekt samostatne – model, ktorý popisuje vývoj v rámci opakovaných pozorovaní jedného konkrétneho subjektu (pričom sa niekedy predpokladá nezávislosť chýb jednodlivých opakovaných meraniach);

Krok 2:
Lineárny regresný model, ktorý popisuje variabilitu medzi jednotlivými subjektami pomocou regresného modelu pre odhadnuté parametre z individuálnych regresných modelov z prvého kroku;

Uvažujúc značenie zavedené vyššie, v prvom kroku sa jedná o \(N \in \mathbb{N}\) nezávislých regresných modelov (vzhľadom k nezávislosti jednotlivých subjektov), ktoré pre každý subjekt \(i \in \{1, \dots, N\}\) môžeme zapísať ako \[ \boldsymbol{Y}_i = (Y_{i 1}, \dots, Y_{i n_i})^\top = \mathbb{Z}_i\boldsymbol{\beta}_{i} + \boldsymbol{\varepsilon}_{i}, \] kde vektor neznámych parametrov \(\boldsymbol{\beta}_i \in \mathbb{R}^q\) je špecifický pre každý subjekt \(i \in \{1, \dots, N\}\) (teda \(\boldsymbol{\beta}_i\) sú obecně rôzne), \(\mathbb{Z}_i\) je príslušná regresná matica modelu a pre vektor chýb (vzhľadom ku korelovanosti/závislosti opakovaných pozorovaní v rámci subjektu) predpokládame napr. že platí \[ \boldsymbol{\varepsilon}_i = \left( \begin{array}{c} \varepsilon_{i 1}\\ \vdots\\ \varepsilon_{i n_i} \end{array} \right) \sim N_{n_i}(\boldsymbol{0}, \Sigma_i), \] kde \(\Sigma_i \in \mathbb{R}^{n_1 \times n_i}\) je pozitívne-definitná variačná-kovariančná matica (opäť obecně rôzna pre jednotlivé subjekty). Náhodný vektor \(\boldsymbol{\varepsilon}_i \sim N_{n_i}(\boldsymbol{0}, \Sigma_i)\) popisuje tzv. within-subject variability v datach (t.j., variabilitu v rámci jednotlivých subjektov).

Užitočné

Vektor neznámych parametrov je obecně odhadnutý rôzne pre rôzne subjekty, ale štruktúra parametru – t.j., príslušná dimenzia \(q \in \mathbb{N}\), ale tiež interpretácia jednotlivých zložiek je naprieč subjektami totožmá.
Dimenzia \(q \in \mathbb{N}\) tym pádom dáva priame predpoklady na počet opakovaných pozorovaní v rámci každého subjektu. Je dobré si uvedomiť, že týmto postupom potrebujeme mať pre každý subjekt aspoň \(q \in \mathbb{N}\) opakovaných pozorovaní, aby sme v modeli \(\boldsymbol{Y}_i = \mathbb{Z}_i\boldsymbol{\beta}_{i} + \boldsymbol{\varepsilon}_{i}\) dokázali odhadnúť neznámy parameter \(\boldsymbol{\beta}_i \in \mathbb{R}^q\) pre každý uvažovaný subjekt.
Predpoklad formulovaný pre variančnú kovariančnú maticu je pomerne bežný, avšak nie nevyhnutne nutný. Napr. pre určitý typ opakovaných pozorovaní je zmysluplné predpokladať, že \(\Sigma_i = \sigma_i^2 \mathbb{I}_{(n_i \times n_i)}\), resp. že dokonca platí \(\sigma_1^2 = \dots = \sigma_N^2\).
(zamyslite sa napr. nad závislostnou štruktúrou opakovaných meraní váhy/výšky u \(N \in \mathbb{N}\) subjektov, ktoré budeme opakovane merať každý deň, alebo každý rok).

Pre ilustráciu uvažujme datový súbor s opakovanými meraniami pacientov so sklerózou multiplex a pre každého pacienta samostatne uvažujme lineárny regresný model (v programe R) pre časovú závislost premennej EDSS. Z výsledných fitovaných regresných modelov nás ale zaujímajú hlavne odhadnuté neznáme (subject-specific) parametre. Nad rámec týchto parametrov zaznamenáme aj pohlavie každého pacienta (t.j., muž = 1 a žena = 2).

sm <- read.csv(url("https://www2.karlin.mff.cuni.cz/~maciak/NMST422/sm_data2.csv"), header = T)

BETA <- NULL
for (subject in 1:140){
  m <- lm(EDSS ~ time, data = sm[sm$id == subject,])
  if (sm$gender[sm$id == subject][1] == "M"){
    BETA <- rbind(BETA, c(m$coeff, 1, sm$age[sm$id == subject][1])) 
  } else {
    BETA <- rbind(BETA, c(m$coeff, 2, sm$age[sm$id == subject][1]))
  }
}

Odhadnuté regresné parametre pre všetkých 140 pacientov (každý z uvažovaných pacientov má k dispozícii aspoň dva opakované pozorovania a tiež platí, že \(\boldsymbol{\beta}_i \in \mathbb{R}^2\), pretože odhadujeme intercept a smernicu pre lineárnu závislosť EDSS' na časetime`).

Následne sa môžeme graficky pozrieť na odhadnuté `subject-specific'' parametre individuálných regresných modelov a prípadne pomocou funkcielowess()` (neparametrické vyhladzovanie dat) zohadní aj dodatočnú informáciu o pohlaví.

plot(BETA[,2] ~ BETA[,1], pch = 21, bg = BETA[,3], xlab = "Intercept", ylab = "Smernica")
lines(lowess(BETA[BETA[,3] == 1, 2] ~ BETA[BETA[,3] == 1,1]), col = 1, lwd = 2)
lines(lowess(BETA[BETA[,3] == 2, 2] ~ BETA[BETA[,3] == 2,1]), col = 2, lwd = 2)
legend("topleft", legend = c("male", "female"), lwd = c(2,2), col = c(1,2))

par(mfrow = c(1,2))
plot(BETA[,1] ~ BETA[,4], pch = 21, bg = BETA[,3], xlab = "Vek [roky]", ylab = "Intercept")
lines(lowess(BETA[BETA[,3] == 1, 1] ~ BETA[BETA[,3] == 1,4]), col = 1, lwd = 2)
lines(lowess(BETA[BETA[,3] == 2, 1] ~ BETA[BETA[,3] == 2,4]), col = 2, lwd = 2)
legend("topleft", legend = c("male", "female"), lwd = c(2,2), col = c(1,2))
plot(BETA[,2] ~ BETA[,4], pch = 21, bg = BETA[,3], xlab = "Vek [roky]", ylab = "Smernica")
lines(lowess(BETA[BETA[,3] == 1, 2] ~ BETA[BETA[,3] == 1,4]), col = 1, lwd = 2)
lines(lowess(BETA[BETA[,3] == 2, 2] ~ BETA[BETA[,3] == 2,4]), col = 2, lwd = 2)
legend("topleft", legend = c("male", "female"), lwd = c(2,2), col = c(1,2))

V druhom kroku sa odhadnuté subject-specific parametre \(\widehat{\boldsymbol{\beta}_i}\) (resp. variabilitu medzi subjektvami – t.j., between-subject variabilitu) vysvetliť pomocou druhého regresného modelu, ktorý je možné matematicky (teoretický model) formulovať ako \[ \boldsymbol{\beta}_i = \mathbb{K}_i\boldsymbol{\beta} + \boldsymbol{b}_i, \] pričom platí, že \(\boldsymbol{\beta}_i \in \mathbb{R}^q\), regresná matica \(\mathbb{K}_i \in \mathbb{R}^{q \times p}\) je opäť tzv. subject-specific (a je typu \(q \times p\)), vektor neznámych parametrov \(\boldsymbol{\beta} \in \mathbb{R}^p\) popisuje rozdiely medzi pacientmi (s analogickou interpretáciou, ako v štandardnom lineárnom regresnom modeli) a náhodné chyby \(\boldsymbol{b}_i \sim N_q(\boldsymbol{0}, \mathbb{D})\) modelujú variabilitu medzi jednotlivými subjektami – t.j., tzv. between-subject variabilitu.

Z hľadiska lineárneho regresného modelu nás zaujímajú následujúce regresné modely:

summary(lm(BETA[,1] ~ BETA[,3]))

## 
## Call:
## lm(formula = BETA[, 1] ~ BETA[, 3])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8676 -1.0338  0.1324  1.2662  3.0043 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.0014     0.5038   7.942 6.24e-13 ***
## BETA[, 3]    -0.1338     0.2842  -0.471    0.639    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.519 on 138 degrees of freedom
## Multiple R-squared:  0.001603,   Adjusted R-squared:  -0.005632 
## F-statistic: 0.2216 on 1 and 138 DF,  p-value: 0.6386

summary(lm(BETA[,2] ~ BETA[,3]))

## 
## Call:
## lm(formula = BETA[, 2] ~ BETA[, 3])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66969 -0.06969 -0.06214  0.08031  0.58031 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.054595   0.066645   0.819    0.414
## BETA[, 3]   0.007548   0.037593   0.201    0.841
## 
## Residual standard error: 0.2009 on 138 degrees of freedom
## Multiple R-squared:  0.000292,   Adjusted R-squared:  -0.006952 
## F-statistic: 0.04031 on 1 and 138 DF,  p-value: 0.8412

Samostatne

Interpretujte jednotlivé parametre (subject-specific parametre \(\boldsymbol{\beta}_i\), aj celkový parameter \(\boldsymbol{\beta}\)).
Aký je formálny záver vyplývajúci z výstupov lineárneho regresného modelu?
Ako by ste interpretovali celkový model pre závislosť EDSS na čase a na pohlaví pacienta?

2. Lineárny regresný model s náhodnými efektami

V predchádzajúcom dvoj-fázovom regresnom modelováni bol vektor opakovaných pozorování vrámci konkrétneho subjektu \(i \in \{1, \dots, N\}\) sumarizovaný pomocou (``summary statistic’’) odhadnutého vektoru parametrov \(\widehat{\boldsymbol{\beta}_i} \in \mathbb{R}^q\) a následne (v druhej fáze) jednotlivé odhadnuté parametre \(\widehat{\boldsymbol{\beta}_1}, \dots, \widehat{\boldsymbol{\beta}_n}\) boli sumarizované prostredníctvom odhadnutého vektoru parametrov \(\widehat{\boldsymbol{\beta}} \in \mathbb{R}^p\).

Uvažujeme teda dva lineárne regresné modely:

Model 1 (resp. \(N\) nezávislých modelov pre každý subjekt samostatne): \[ \boldsymbol{Y}_i = \mathbb{Z}_i\boldsymbol{\beta}_i + \boldsymbol{\varepsilon}_i, \] pre data v tvare \(\{(Y_{i j}, \boldsymbol{z}_j);~j = 1, \dots, n_i\}\), kde \(\boldsymbol{z}_j \in \mathbb{R}^q\).

Model 2 (resp. \(q\) samostatných regresných modelov): \[ \beta_{i \ell} = \mathbb{K}_i\boldsymbol{\beta}_\ell + b_{i \ell}, \] pričom namiesto skutočných (neznámych) parametrov \(\beta_{i \ell}\) používame empirické protějšky, t.j. odhady \(\widehat{\beta}_{i \ell}\) a teda príslušné data lze reprezentovať v tvare \(\{(\widehat{\beta}_{i \ell}, \boldsymbol{k}_\ell^{(i)});~i = 1, \dots, n\}\) pre \(\ell = 1, \dots, q\) postupných modelov. Vektor \(\boldsymbol{k}_\ell^{(i)}\) predstavuje \(\ell\)-tý riadok matice \(\mathbb{K}_i\).

Oba modely je možné uvažovať dohromady, resp. \[ \left. \begin{array}{c} \boldsymbol{Y}_i = \mathbb{Z}_i\boldsymbol{\beta}_i + \boldsymbol{\varepsilon}_i\\ \boldsymbol{\beta}_i = \mathbb{K}_i\boldsymbol{\beta} + \boldsymbol{b}_i\\ \end{array} \right\} \Longrightarrow \boldsymbol{Y}_i = \mathbb{Z}_i \mathbb{K}_i\boldsymbol{\beta} + \mathbb{Z}_i\boldsymbol{b}_i + \boldsymbol{\varepsilon}_i, \] čo je základná formulácia (definícia) lineárneho regresného modelu s náhodnými efektami \(\boldsymbol{b}_i \sim N_q(\boldsymbol{0}, \mathbb{D})\) a zároveň \(\boldsymbol{\varepsilon}_i \sim N_{n_i}(\boldsymbol{0}, \Sigma_i)\). Navyše sa štandardne predpokladá aj vzájomná nezávislosť medzi chybovými členmi, t.j. medzi náhodnými vektormi \(\boldsymbol{\varepsilon}_1, \dots, \boldsymbol{\varepsilon}_N, \boldsymbol{b}_1, \dots, \boldsymbol{b}_N\).

Ak označíme maticu \(\mathbb{Z}_i\mathbb{K}_i\) ako \(\mathbb{X}_i\) a združíme všetky subjekty \(i \in \{1, \dots, N\}\) do jedného modelu prostredníctvom vektoru závislých pozorovaní \(\boldsymbol{Y} = (\boldsymbol{Y}_1^\top, \dots, \boldsymbol{Y}_N^\top)^\top \in \mathbb{R}^\mathcal{N}\), tak získame výsledný model v tvare \[ \boldsymbol{Y} = \mathbb{X}\boldsymbol{\beta} + \mathbb{Z}\boldsymbol{b} + \boldsymbol{\varepsilon}, \] kde regresna matica \(\mathbb{X} \in \mathbb{R}^{\mathcal{N} \times p}\) (prislúchajúca pevným efektom) je definovaná ako \[ \mathbb{X} = (\mathbb{X}_1^\top, \dots, \mathbb{X}_N^\top)^\top \] a regresná matica \(\mathbb{Z} \in \mathbb{R}^{\mathcal{N} \times Nq}\) (prislúchajúca náhodným efektom \(\boldsymbol{b} = (\boldsymbol{b}_1^\top, \dots, \boldsymbol{b}_N^\top)\)), je definovaná ako \[ \mathbb{Z} = \left( \begin{array}{cccc} \mathbb{Z}_1 & \boldsymbol{0} & \dots & \boldsymbol{0}\\ \boldsymbol{0} & \mathbb{Z}_2 & \dots & \boldsymbol{0}\\ \vdots & \vdots & \ddots & \vdots\\ \boldsymbol{0} & \boldsymbol{0} & \dots & \mathbb{Z}_N \end{array} \right). \]

Samostatne

Overte dimenzie jednotlivých objektov (matíc a vektorov) v obecnej formulácii lineárneho regresného modelu s náhodnymi efektami.
Aká je interpretácia jednotlivých parametrov?
V súvislosti s lineárnym regresným modelom s náhodnými efektami sa v literatúre uvádzajú dve analogické, ale nie ekvivalentné formulácie: tzv. hierarchický model, ktorý špecifikuje podmienené rozdelenie \(\boldsymbol{Y}_i|\boldsymbol{b}_i\) a rozdelenie náhodných efektov \(\boldsymbol{b}_i\);
Druhou formuláciou je tzv. marginálny model, ktorý priamo špecifikuje rozdelenie náhodných vektorov \(\boldsymbol{Y}_i\). Zamyslite sa nad jednotlivými formuláciami a premyslite výhody a nevýhody jednotlivých zápisov.

Pre ilustráciu lineárneho modelu s náhodnými efektami využijeme opäť datový súbor s pacientami so sklerózou multiplex. Data načítame do programu SAS:

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;
    
proc print datafile = sm.data; 
run;

a pomcou procedúry ´proc mixed nodhadneme parametre príslušného lineárneho modelu (najprv bez náhodných efektov ale s explicitnou špecifikáciou štruktúry opakovaných pozorovaní – AR(1) proces).

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 / type = AR(1) subject = id;
run; 

proc mixed data = sm.data2 method = ml; 
class gender(ref = "F") timeCls;
model EDSS = gender time*gender / s;
repeated timeCls / type = AR(1) subject = id;
run;

Následne skusíme špecifikovať maticu náhodných efektov \(\mathbb{Z}\) pomocu tzv. ‘random statement’,

proc mixed data = sm.data2 method = reml; 
class gender(ref = "F") timeCls;
model EDSS = gender time*gender / s;
random intercept / subject = id v g cl solution;
run;

Porovnajte model vyššie s následujúcim modelom (a prípadne využijte parameter noint v “model statement”):

proc mixed data = sm.data2 method = ml; 
class gender(ref = "F") timeCls;
model EDSS = gender time*gender / s;
repeated timeCls / type = AR(1) subject = id;
random intercept / subject = id v g cl solution;
run;

Samostatne

Aké modely sú odhadnuté vyššie? Ako vyzerá príslušná regresná matica \(\mathbb{X}\) matica a matica náhodných efektov \(\mathbb{Z}\)? Aká je interpretácia jednotlivých parametrov?
Aká je variančna kovariačna štruktúra opakovaných pozorovaní? Explicitne zapíšte príslušnú variančnú-kovariačnú maticu.
Aké sú alternatvívne možnosti pre špecifikáciu neznámej variančnej-kovarianečnej matice?

Užitočné

Všimnite si, že hierarchický model priamo implikuje formuláciu marginálneho modelu. Opačná mplikácia ale zjavne neplatí. Prečo?

Samostatný úkol

Preštudujte si implementáciu SAS procedúry proc mixed napr. na tejto stránke. Čo je výstupom tejto funkcie a ako jednotlivé časti výstupu interpretovať?
Pomocou programu SAS sa pre Vami zvolené data (napr. datový súbor sm_data2.csv) pokúste implementovať dvoj-stupňový regresný model pre longitudinálne data.