Longitudinální a panelová data – NMST422

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

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).

• náhodný člen $\boldsymbol{b}_i$ predstavuje tzv. medzi-subjektovú variabilitu (between-subject variability)
• náhodný člen $\boldsymbol{\varepsilon}_{1 i}$ predstavuje tzv. chyby merania (measurement error)
• náhodný člen $\boldsymbol{\varepsilon}_{2 i}$ predstavuje tzv. sériovu koreláciu (serial correlation)

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\}$.

• exponenciálna sériová korelácia: $g(x) = exp\{-\phi x\}$, pre nejaký neznámy parameter $\phi > 0$
• tzv. gausovská sériová korelácia: $g(x) = exp\{- \phi x^2\}$, pre $\phi > 0$

1. Porovnanie odhadovaných modelov v SAS a R

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")
?lmer

V 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))

2. Restricted maximum likelihood – REML

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:

• metóda maximálnej vierohodnosti
• tzv. profilová vierohodnosť (profile likelihood)
• tzv. omezená maximálna vierohodnosť (restricted maximum likelihood)

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ť.

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))

3. Variogram v programe SAS

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;