Letný semester 2024-2025 | Cvičenie 8 | 14.04.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.
Predpokládajme základný (lineárny) regresný model s náhodnými
efektami, ktorý je vyjadrený prostredníctvom rovnice \[
\boldsymbol{Y}_i = \mathbb{X}_i\boldsymbol{\beta} +
\mathbb{Z}_i\boldsymbol{b}_i + \boldsymbol{\varepsilon}_i
\] kde \(\boldsymbol{Y}_i = (Y_{i1},
\dots, Y_{i n_i})^\top \in \mathbb{R}^{n_i}\) je vektor
opakovaných meraní vrámci \(i\)-teho
subjektu (pre \(i \in \{1, \dots,
N\}\)) a \(\boldsymbol{b}_i = (b_{i1},
\dots, b_{i q})^\top \in \mathbb{R}^q\) je vektor náhodných
(nepozorovaných) efektov v rámci \(i\)-teho subjektu, pre ktorý väčšinou
predpokládame, že \(\boldsymbol{b}_i \sim
N_q(\boldsymbol{0}, \mathbb{D})\) a náhodný vektor chýb \(\boldsymbol{\varepsilon}_{i}\) má opäť
mnohorozmerné normálne rozdelenie s nulovým vektorom stredných hodnôt a
určitou variančnou-kovariančnou štruktúrou. Jednotlivé vektory \(\boldsymbol{Y}_1, \dots, \boldsymbol{Y}_N\)
sú vzájomné nezávislé a chybový člen \(\boldsymbol{\varepsilon}_i\) sa často v
literatúre rozdeľuje na tzv. chyby merania a tzv. sériovú
koreláciu.
timeClas v
zásise nižšie).Nutná je opäť identifikácia subjektov. Pomocou parametru
type = sa volí variančná-kovariačná matica pre sériouvú
koreláciu opakovaných pozorovaní. V návode (help) procedúry pozri aj
parametre r a rcorr.
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;
time2 = time * time;
run;contrast statement a voľba (napr.)
chisc) covtest v
proc mixed statement)V nasledujúcom kóde je využitý približný Waldov test (resp. \(t\)-test a \(F\)-test) s Kenward-Roger aproximáciou počtu stupňov voľnosti.
proc mixed data = sm.data2 method = ml;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;
contrast 'Age vs. time' age*time2 1 / chisq;
run;
proc mixed data = sm.data2 method = ml;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;
contrast 'Age vs. time' age*time2 1,
age*time 1 / chisq;
run;
proc mixed data = sm.data2 method = ml;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;
contrast 'Age vs. time' age*time2 1,
age*time 1 / chisq;
contrast 'Quadratic time' time2 1 / chisq;
run;
proc mixed data = sm.data2 method = ml;
class id gender timeCls;
model EDSS = gender age time time2 / solution;
random intercept gender / type = vc subject = id;
repeated timeCls / type = simple subject = id;
contrast 'Overall time' time2 1,
time 1 / chisq;
run;Využijte metódu REML, prípadne alternatívne metódy aproximácie
stupňov voľnosti a porovnajte výsledky. (robustné odhady pomocou tzv.
sandwich odhadu pomocou parametru empirical v
proc statemente)
Test podmodelu založený na pomere vierohodnosti spočítanom cez REML
môže mať za následok aj zápornú hodnotu testovej štatistiky. Testy o
\(\boldsymbol{\beta} \in \mathbb{R}^p\)
parametri (t.j., štruktúra podmienenej strednej hodnoty) pomerom
vierhodnosti preto nie sú valídne v prípade, že je použitá metóda REML.
V pripade testov o paremtri \(\boldsymbol{\alpha} \in \mathbb{R}^q\) je
vhodnejšie vyžiť metódu REML, ktorá má za následok nestranné odhady
parametrov modelujúcich variačnu/kovariančnú štruktúru.
V následujúcej časti je uvažované á postupná inferencie o náhodných efektov. Niektoré testy ale nemaju logický význam…
proc mixed data = sm.data2 method = reml covtest;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept time time2 / type = un subject = id;
repeated timeCls / type = simple subject = id;
run;
proc mixed data = sm.data2 method = reml covtest;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept time time2 / type = vc subject = id;
repeated timeCls / type = simple subject = id;
run;
proc mixed data = sm.data2 method = reml covtest;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept time / type = vc subject = id;
repeated timeCls / type = simple subject = id;
run;
proc mixed data = sm.data2 method = reml covtest;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept time / type = un subject = id;
repeated timeCls / type = simple subject = id;
run;Test pomerom vierohodnosti:
/* Full model with random effect */
proc mixed data = sm.data2 method=ml noclprint noinfo covtest ic;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept time time2 / type = un subject = id;
repeated timeCls / type = simple subject = id;
ods output FitStatistics=full_fit;
run;
/* Reduced model without random effect */
proc mixed data = sm.data2 method=ml noclprint noinfo ic;
class id gender timeCls;
model EDSS = gender age time time2 age*time age*time2 / solution;
random intercept time / type = un subject = id;
repeated timeCls / type = simple subject = id;
ods output FitStatistics=reduced_fit;
run;
/* Extract -2 Log Likelihoods */
data full_ll;
set full_fit;
if Description = '-2 Res Log Likelihood' then call symputx('ll_full', Value);
run;
data reduced_ll;
set reduced_fit;
if Description = '-2 Res Log Likelihood' then call symputx('ll_reduced', Value);
run;
/* Compute LRT and p-value */
data lrt_result;
ll_diff = 1527.5 - 1430.0;
p_value = 1 - cdf('CHISQ', ll_diff, 1);
put "----------------------------------------";
put "Likelihood Ratio Test for Random Effect:";
put "LRT statistic = " ll_diff;
put "Degrees of freedom = " df;
put "P-value = " p_value;
put "----------------------------------------";
run;