Longitudinální a panelová data – NMST422

Letný semester 2026 | Cvičenie 4 | 30.03.2026

Prihlásenie k SAS OnDemand: https://www.sas.com/en_us/software/on-demand-for-academics.html
Nutná je registrácia s vytvorenie účtu s vlastným identifikačným číslom a potvrdenie registrácie prostredníctvom (univerzitného) emailu zadaného pri registrácii. Identifikačné číslo užívateľa (vo forme uXXX, kde XXX je samotné číslo uživateľa) sa vyskytuje v jednotlivých SAS skriptoch uvedených nižšie (symbol XXX v skriptoch je potrebné 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.
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)

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

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 a navyše sú tieto pozorovania uskutočné vo vzájomne iných časových okamžikoch). Z tohto dôvodu nie je možné aplikovať mnohorozmerné štatistické postupy (napr. mnohorozmerný linárny regresný model – viď minulé cvičenie – ktorý predpokladá nezávislé a rovnako rozdelené vektory náhodných chýb medzi jednotlivými subjektami). 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 nutné pomocou dat 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{w} + \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 daného subjektu, \(i \in \{1, \dots, N\}\). Celkový počet pozorovaní je teda \(\mathrm{N} = \sum_{i = 1}^N n_i\).

Regresná matica \(\mathbb{X}\) je typu \(\mathrm{N} \times p\), pre vektor neznámych parametrov (pevných efektov) platí \(\boldsymbol{\beta} \in \mathbb{R}^p\) a matica \(\mathbb{Z}\) je typu \(\mathrm{N} \times Nr\) a prislúcha náhodnym efektom \(\boldsymbol{w} = (\boldsymbol{w}_1^\top, \dots, \boldsymbol{w}_N^\top)^\top \in \mathbb{R}^{Nr}\), kde \(\boldsymbol{w}_i = (w_{i 1}, \dots, w_{i r})^\top \in \mathbb{R}^r\) reprezentuje tzv. ``subject-specific’’ náhodné efekty (typicky sa predpokladá, že \(\boldsymbol{w}_i \sim N_r(\boldsymbol{0}, \mathbb{G})\)) 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. \(r \in \mathbb{N}\).

Pre pripomenutie, náhodné efekty v \(\mathbb{Z}\boldsymbol{b}\) vznikajú formálne, takzvanou mean-variance separáciou a rozlíšenim rôznych zdrojov variability. Konkrétne pre stochastický člen \(\varepsilon_{ij}\), ktorý predstavuje
náhodnú chybu \(i\)-teho subjektu v jeho \(j\)-tom pozorovaní, predpokládame, že platí

\[ \varepsilon_{ij} = \boldsymbol{z}_{ij}^\top \boldsymbol{w}_i + W_i(t_{ij}) + \omega_{ij}, \] kde \(\boldsymbol{z}_{i,j}\) predstavuje vektor vysvetľujúcich premenných pre subjekt \(i \in \{1, \dots, N\) v jeho \(j\)-tom meraní (väčšinou sa v praxi jedná o niektoré vybrané regresory z matice \(\mathbb{X}\)), \(\boldsymbol{\omega}_1, \dots, \boldsymbol{\omega}_N\) sú nezávislé náhodné vektory typicky z \(r\)-rozmerného normálneho rozdelenia \(N_r(\boldsymbol{0}, \mathbb{G})\), a \(W_i(t_{ij}) + \omega_{ij}\) predstavuje latentnú chybu (tzv. serial correlation \(+\) measurement error).

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 pri 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}^r\) 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 \(r \in \mathbb{N}\), ale tiež interpretácia jednotlivých zložiek je naprieč subjektami rovnaká.
Dimenzia \(r \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ň \(r \in \mathbb{N}\) opakovaných pozorovaní, aby bolo možné v modeli \(\boldsymbol{Y}_i = \mathbb{Z}_i\boldsymbol{\beta}_{i} + \boldsymbol{\varepsilon}_{i}\) odhadnúť neznámy parameter \(\boldsymbol{\beta}_i \in \mathbb{R}^r\).
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:142){
  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 142 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 čase time).

Následne sa môžeme graficky pozrieť na odhadnuté subject-specific parametre individuálných regresných modelov a prípadne pomocou funkcie lowess() (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], f = 1), col = 1, lwd = 2)
lines(lowess(BETA[BETA[,3] == 2, 1] ~ BETA[BETA[,3] == 2,4], f = 1), 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], f = 1), col = 1, lwd = 2)
lines(lowess(BETA[BETA[,3] == 2, 2] ~ BETA[BETA[,3] == 2,4], f = 1), col = 2, lwd = 2)
legend("topleft", legend = c("male", "female"), lwd = c(2,2), col = c(1,2))

V druhom kroku sú odhadnuté subject-specific parametre \(\widehat{\boldsymbol{\beta}_i}\) modelované pomocou druhého regresného modelu, ktorý vysvetľuje variabilitu medzi jednotlivými subjektami – t.j., between-subject variabilita.

Model lze matematicky formulovať ako \[ \boldsymbol{\beta}_i = \mathbb{K}_i\boldsymbol{\beta} + \boldsymbol{w}_i, \] pričom platí, že \(\boldsymbol{\beta}_i \in \mathbb{R}^r\), regresná matica \(\mathbb{K}_i \in \mathbb{R}^{r \times p}\) je opäť tzv. subject-specific (a je typu \(r \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{w}_i \sim N_r(\boldsymbol{0}, \mathbb{G})\) modelujú variabilitu medzi jednotlivými subjektami – t.j., tzv. between-subject variabilitu.

Z hľadiska lineárneho regresného modelu vysvetľujeme variabilitu v intercept parametroch BETA[,1] a v smerniciach BETA[,2] v závislosti na pohlaví BETA[,3] a prípadne veku BETA[,4].

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

## 
## Call:
## lm(formula = BETA[, 1] ~ BETA[, 3])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8827 -1.0124  0.1173  1.2483  3.0043 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.0315     0.4902   8.224 1.21e-13 ***
## BETA[, 3]    -0.1489     0.2779  -0.536    0.593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.511 on 140 degrees of freedom
## Multiple R-squared:  0.002046,   Adjusted R-squared:  -0.005082 
## F-statistic: 0.287 on 1 and 140 DF,  p-value: 0.593

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

## 
## Call:
## lm(formula = BETA[, 2] ~ BETA[, 3])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66969 -0.06969 -0.06969  0.08031  0.70513 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.09148    0.06763   1.353    0.178
## BETA[, 3]   -0.01089    0.03833  -0.284    0.777
## 
## Residual standard error: 0.2085 on 140 degrees of freedom
## Multiple R-squared:  0.0005766,  Adjusted R-squared:  -0.006562 
## F-statistic: 0.08077 on 1 and 140 DF,  p-value: 0.7767

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

## 
## Call:
## lm(formula = BETA[, 1] ~ BETA[, 4])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2669 -0.7942  0.2068  1.0221  2.9816 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.76644    0.51281   3.445 0.000755 ***
## BETA[, 4]    0.05683    0.01409   4.035 8.96e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.432 on 140 degrees of freedom
## Multiple R-squared:  0.1042, Adjusted R-squared:  0.09776 
## F-statistic: 16.28 on 1 and 140 DF,  p-value: 8.955e-05

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

## 
## Call:
## lm(formula = BETA[, 2] ~ BETA[, 4])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68090 -0.08836 -0.03874  0.06538  0.62027 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.103106   0.073105  -1.410   0.1606  
## BETA[, 4]    0.004973   0.002008   2.477   0.0145 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2041 on 140 degrees of freedom
## Multiple R-squared:  0.04197,    Adjusted R-squared:  0.03513 
## F-statistic: 6.134 on 1 and 140 DF,  p-value: 0.01445

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í v rámci konkrétneho subjektu \(i \in \{1, \dots, N\}\) sumarizovaný pomocou (``summary statistic’’) odhadnutého vektoru parametrov \(\widehat{\boldsymbol{\beta}_i} \in \mathbb{R}^r\) 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^{(i)});~j = 1, \dots, n_i\}\), kde \(\boldsymbol{z}_j^{(i)} \in \mathbb{R}^r\) je \(j\)-tý riadok matice \(\mathbb{Z}_i\). Matica \(\mathbb{Z}_i\) obsahuje len premenné, ktoré sa v rámci subjektu menia. Zaradenie tzv. cross-sekčných premenných, ktoré sú v rámci subjektu konštantné, by znamenalo lineárne závisle stlpce matice \(\mathbb{Z}_i\), keďže model obsahuje intercept parameter.

Model 2 Celkový regresný model pre subject-specific parametre \(\boldsymbol{\beta}_i\) a \(i = 1, \dots, N\) kde \[ \boldsymbol{\beta}_i = \mathbb{K}_i\boldsymbol{\beta} + \boldsymbol{w}_i \] resp. \(r\) samostatných regresných modelov (pre \(\ell \in \{1, \dots, r\}\)), kde \[ \beta_{i \ell} = \boldsymbol{k}_{\ell}^{(i)}\boldsymbol{\beta} + w_{i \ell}, \] kde vektor \(\boldsymbol{k}_\ell^{(i)}\) predstavuje \(\ell\)-tý riadok matice \(\mathbb{K}_i\). Matica \(\mathbb{K}_i \in \mathbb{R}^{r \times p}\) ma často jednoduchý tvar a zároveň \(r < p\) (t.j., počet náhodných efektov je menší, než počet pevných efektov). Pri vhodnej permutácii prvkov vektorov \(\boldsymbol{\beta}\) a \(\boldsymbol{\beta}_i\) a riadkov matice \(\mathbb{K}_i\) pre \(i = 1, \dots, N\) je možné maticu \(\mathbb{K}_i\) vyjadriť napr. v tvare \(\mathbb{K}_i = (\mathbb{I}_{r \times r}, \mathbb{0})\), kde \(\mathbb{I}_{r \times r}\) je jednotková matice typu \(r \times r\) a \(\mathbb{O} \in \mathbb{R}^{r \times (p - r)}\) je matica núl s príslušnými rozmermi (teda na prvých miestach sú náhodne efekty a potom následujú pevné efekty).

V praxi sa namiesto skutočných (neznámych) parametrov \(\beta_{i \ell}\) používajú 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, r\) postupných modelov.

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{w}_i\\ \end{array} \right\} \Longrightarrow \boldsymbol{Y}_i = \mathbb{Z}_i \mathbb{K}_i\boldsymbol{\beta} + \mathbb{Z}_i\boldsymbol{w}_i + \boldsymbol{\varepsilon}_i, \] čo je základná formulácia (definícia) lineárneho regresného modelu s náhodnými efektami \(\boldsymbol{w}_i \sim N_r(\boldsymbol{0}, \mathbb{G})\) a zároveň \(\boldsymbol{\varepsilon}_i \sim N_{n_i}(\boldsymbol{0}, \Sigma_i)\). Jedná sa ale o tzv. subject-specific model, teda pre konkrétny subjekt \(i \in \{1, \dots, N\}\). 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{w}_1, \dots, \boldsymbol{w}_N\).

Ak použijeme namiesto matice \(\mathbb{Z}_i\mathbb{K}_i\) regresnú maticu \(\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}^\mathrm{N}\), tak získame výsledný model v tvare \[ \boldsymbol{Y} = \mathbb{X}\boldsymbol{\beta} + \mathbb{Z}\boldsymbol{w} + \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}^{\mathrm{N} \times Nq}\) (prislúchajúca náhodným efektom \(\boldsymbol{w} = (\boldsymbol{w}_1^\top, \dots, \boldsymbol{w}_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

Ako by vyzerala matica \(\mathbb{K}_i\) v druhom kroku dvojsupňového modelu vyššie v prípade jednoduchého modelu s interceptom a smernicou (ako na prednáške)?
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{w}_i\) a rozdelenie náhodných efektov \(\boldsymbol{w}_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';  /* SAS library reference */
filename reffile '/home/uXXX/sasuser.v94/data/sm_data2.csv'; /* data file reference */

proc import datafile=reffile
    dbms=csv
    out=sm.data
    replace;
    getnames=yes;
run;
    
proc print datafile = sm.data; 
run;

a pomocou procedúry PROC MIXED nodhadneme parametre príslušného lineárneho regresného modelu (bez náhodných efektov a za predpokladu nezávislosti jednotlivých meraní a následne s explicitnou špecifikáciou štruktúry opakovaných pozorovaní – napr. AR(1) proces).

data sm.data2;
set sm.data;
timeCls = time;
run; 

/* independent observations */
proc mixed data = sm.data2 method = ml; 
class gender timeCls; 
model EDSS = gender time*gender / s; 
run; 

/* independent observations */
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;

Predchádzajúci model vychádza z tvaru \(Y_{ij} = \boldsymbol{X_{ij}}^\top\boldsymbol{\beta} + \varepsilon_{i,j}\), pričom predpokladá autokorelačnú štruktúru prvého rádu – AR(1) model pre náhodné chyby, teda \(cov(\varepsilon_{ij}, \varepsilon_{ik}) = \sigma^2 \rho^{|j - k|}\).

Následne lze špecifikovať maticu náhodných efektov \(\mathbb{Z}\) pomocu tzv. random statement v procedúre PROC MIXED, teda predpokládame model \(Y_{ij} = \boldsymbol{X}_{ij}^\top \boldsymbol{\beta} + w_{i} + \varepsilon_{ij}\), kde \(\varepsilon_{i,j} \sim N(0, \sigma^2)\) a \(w_i \sim N(n, \nu^2)\), pričom platí, že \(cor(Y_{ij}, Y_{ik}) = \frac{\nu^2}{\sigma^2 + \nu^2}\). Model je v SASe implementovaný pomocou random statement v procedure PROC MIXED:

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

V modeli vyššie bol použitý odhad pomocou REML (method = reml) pomocou ktorého získavame nestrané odhady variančnej štruktúry. Porovnajte výsledky predchádzajúceho modelu – špecialne získané odhady smerodatných chýb – s modelom fitovaným pomocou klasickej metódy maximálnej vierohodnosti (method = ml):

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

Porovnajte modely vyššie s následujúcim modelom ktorý kombinuje jednak autoregresnú štruktúru opakovaných pozorovaní a taktiež využíva náhodný intercep (prípadne využijte dodatočný parameter noint v model statement a modely porovnajte):

proc mixed data = sm.data2 method = reml; 
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;

Uvažovaný model je relatívne jednoduchý (čo sa týka špecifikácie podmienenej strednej hodnoty) a uvažované korelačné štruktúry pôsobia proti sebe: AR(1) v repeated statement predpokladá postupne klesajúcu lineárnu závislosť medzi \(Y_{ij}\) a \(Y_{ik}\) (dvoma opakovanými pozorovaniami). Na druhej strane, random statement predpokladá, že corelácia medzi ľubovolnými dvoma opakovanými pozorovaniami je konštantná.

Samostatne

Ako vyzerá príslušná regresná matica \(\mathbb{X}\) matica a matica náhodných efektov \(\mathbb{Z}\) v modeloch 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 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.