Longitudinální a panelová data – NMST422

Letný semester 2023-2024 | Cvičenie 5 | 08.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ály


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.., že počet opakovaných pozorovaní vrámci jedneho subjektu je pre rôzne subjekty rôzna). Z toho 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é matematické/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 nebalancovaných) dat (opakovaných meraní vrámci subjektov) 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}\) subjektov), matica \(\mathbb{X}\) je tzv. regresná matica modelu, vektor \(\boldsymbol{\beta} \in mathbb{R}^p\) predstavuje neznáme parametre, ktoré je potrebné odhadnúť a \(\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 m_1}, Y_{21}, \dots, Y_{nm_n})^\top \in \mathbb{R}^N\) vektor závislej premennej \(Y\) nameranej jednak pre \(n \in \mathbb{N}\) nezávislých subjektov, ale zároveň aj \(m_i \in \mathbb{N}\) opakovaných pozorovaní vrámci každého subjektu \(i \in \{1, \dots, N\}\). Celkový počet pozorovaní 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}\).

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. Dvoj-fázový (2-stage) regresný model pre longitudinálne data

Idea modelovania longitudinálnych dat pomocou dvoj-fázoveho postupu je založená na dvoch postupných krokoch:

  • Krok 1:
    Lineárny regresný model pre každý subjekt samostatne – model, modeluje vývoj vrámci opakovaných pozorovaní u každého subjektu;
  • Krok 2:
    Lineárny regresný model, ktorý vysvetlí variabilitu v odhadnutých regresných parametroch z individuálnych modelov z prvého kroku;

Uvažujúc značenie zavedené vyššie, v prvom kroku sa jedná o fitovanie \(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 m_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í vrámci subjektu) predpokládame napr. že platí \[ \boldsymbol{\varepsilon}_i = \left( \begin{array}{c} \varepsilon_{i 1}\\ \vdots\\ \varepsilon_{i m_i} \end{array} \right) \sim N_{m_i}(\boldsymbol{0}, \Sigma_i), \] kde \(\Sigma_i \in \mathbb{R}^{m_1 \times m_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_{m_i}(\boldsymbol{0}, \Sigma_i)\) popisuje tzv. within-subject variability v datach.



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č subjektvami identická.
  • Dimenzia \(q \in \mathbb{N}\) tym pádom dáva priame predpoklady na počet opakovaných pozorovaní vrá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}_{(m_i \times m_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 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 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, m_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_{m_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}^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}^{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}^{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 výhodami a nevýhodami jednotlivých formulácii.



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/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;

a pomcou procedúry proc mixed nodhadneme parametre príslušného lineárneho modelu s náhodnými efektami.

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;



Samostatne

    ´
  • Aký model je odhadnutý vyššie? 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. Č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.