Longitudinální a panelová data – NMST422

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?

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.