Letný semester 2025-2026 | Cvičenie 6 | 20.04.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


VI. Lineárne modely s náhodnými efektami (inferencia)

Predpokládajme lineárny regresný model pre spojitú závislú premennú \(Y \in \mathbb{R}\) (t.j., náhodnú veličinu opakovane meranú/sledovanú pre \(N \in \mathbb{N}\) vzájomne nezávislých subjektov). Navyše budeme predpokladať, že \(Y\) má(hierarchické) podmienené normálne rozdelenie. Matematicky je tento fakt vyjadrený prostredníctvom zápisu \[ \boldsymbol{Y}_i | \mathbb{X}_i, \boldsymbol{w}_i \sim N_{n_i}\Big( \mathbb{X}_i\boldsymbol{\beta}_i + \mathbb{Z}_i \boldsymbol{w}_i, \mathbb{R}_i\Big) \] kde navyše platí (hierarchia modelu), že \(\boldsymbol{w}_i \sim N_{r}(\boldsymbol{0}, \mathbb{G})\) (pričom náhodné vektory \(\boldsymbol{w}_1, \dots, \boldsymbol{w}_N\) sú vzájomne nezávislé). Samotný (lineárny) regresný model (s náhodnými efektami) môže byť vyjadrený pomocou rovnice \[ \boldsymbol{Y}_i = \mathbb{X}_i\boldsymbol{\beta} + \mathbb{Z}_i\boldsymbol{w}_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í v rámci \(i\)-teho subjektu (pre \(i \in \{1, \dots, N\}\)), \(\boldsymbol{w}_i = (w_{i1}, \dots, w_{i r})^\top \in \mathbb{R}^r\) je vektor náhodných efektov a pre náhodný vektor chýb \(\boldsymbol{\varepsilon}_{i}\) sa predpokladá mnohorozmerné normálne rozdelenie s nulovým vektorom stredných hodnôt a určitou variančnou-kovariančnou štruktúrou – pozitívne definitná variančná-kovariančná matica \(\mathbb{R}_i \in \mathbb{R}^{n_i \times n_i}\). 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.

Poznámka

V praxi sa ale stáva, že predpoklad podmieneného (mnohorozmerného) normálneho rozdelenia pre opakovane pozorovania – t.j., náhodné vektory \(\boldsymbol{Y}_i\), pre \(i =1, \dots, n\) je nerealistický a je nutné hľadať iný pravdepodobnostný model (napr. pretože sledovaná závislá premenná informuje výhradne len o úspechu/neúspechu liečby – binárna premenná – alebo sa všeobecne jedná o realizácie nejakej diskrétnej náhodnej veličiny). Ak je možné naviac postulovať (predpokladať) konkrétne rozdelenie pre závislú premennú (to znamená aj možnosť definovať celkovú vierohodnosť), tak je vhodné použíť tzv. zovšeobecnené lineárne modely s náhodnými efektami (generalized linear model with random effects – GLMM). Jedná sa o rozšírenie triedy zovšeobecnených lineárnych regresných modelov (GLM) v podobnom zmysle, ako sú lineárne regresné modely s náhodnými efektami zovšeobecnením klasických lineárnych regresných modelov.

Ak nie je možné apriórne postulovať (predpokladať) nejaké vhodné pravdepodobnostné rozdelenie pre závislú premennú \(\boldsymbol{Y}_i \in \mathbb{R}^{n_i}\), ale je možné špecifikovať niektoré (podmienené) momenty náhodných vektorov \(\boldsymbol{Y}_i\), tak je vhodné využiť iné odhadovacie metódy a postupy (napr. tzv. GEE metódy, ktoré budeme diskutovať neskôr).

Na rozdiel od klasických lineárnych regresných modelov, ktoré predpokládajú nezávislé pozorovania, je nutné pri modeloch s náhodnými efektami náležite zohľadniť korelačnú štruktúru v rámci opakovaných pozorovaní. Pri analýze dat je preto nutné dbať na správnu špecifikáciu variančnej-kovariančnej (resp. korelačnej) štruktúry a správny popis jednotlivých zdrojov variability – t.j. variabilita medzi jednotlivými subjektami a variabilita opakovaných pozorovaní v rámci konkrétnych subjektov.


Podrobnejšie sa pozrieme na dve rôzne formulácie modelu s náhodnými efektami (LMM)

  • Marginálna formulácia modelu, ktorá vychádza z predpokladu \(\boldsymbol{Y}_i | \mathbb{X}_i \sim N_{n_i}(\mathbb{X}_i\boldsymbol{\beta}, \mathbb{Z}_i \mathbb{G}\mathbb{Z}_i^\top + \mathbb{R}_i)\)
  • Hierarchická formulácia modelu, postavené na predpoklade \(\boldsymbol{Y}_i | \mathbb{X}_i, \boldsymbol{w}_i \sim N_{n_i}(\mathbb{X}_i\boldsymbol{\beta} + \mathbb{Z}_i \boldsymbol{w}_i, \mathbb{R}_i)\)



1. Teoretický a empirický variogram

V prípade opakovaných balancovaných pozorovaní je variančnú-kovariančnú štruktúru možné celkom dobre exploratívne analyzovať napr. pomocou výberovej korelačnej matice (resp. výberovej autokorelačnej funkcie), prípadne vhodného scatterplot grafu (scatterplot grafov).

V prípade, že celková dimenzia datoveho súboru nie je príliš veľka, celkom nápomocné sú tzv. pairwise scatterplots, ktoré postupne vykresľujú scatterplot pre všetky možné dvojice uvažovaných premenných. Avšak aj kompletný súbor všetkých možných vzájomných dvojíc scatterplotov neobsahuje komplexnú (mnohorozmernú) závislostnú štruktúru celkového datového súboru – vždy sa jedná len o určitú projekciu do dvojrozmerného priestoru—t.j., do 2D roviny grafu—a určitú stratu celkovej informácie.

Pre niektoré uvažované premenné z datového súboru o pacientoch so sklerózou multiplex získame príslušnú maticu scatterplotov v programe SMS pomocou následujúceho kódu:

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 sgscatter data=sm.data;
  title "Scatterplot Matrix (SM data)";
  matrix time age LEMsum timeBef EDSSini EDSS NEDA
         / group = gender diagonal = (histogram);
run;
title;

Ak je potrebné zohľadniť viac, ako len jednú referečnú skupinu, je nutné jednotlivé kategorické premenné sumarizovať pomocou jednej premnnej a následne použíť analogicky SAS kód. Napríklad, pre zobrazenie podobného scatterplotu bodov, ale rozlíšujúc obe pohlavia aj 0/1 záznam o aktivite nemoci (NEDA), je možné použiť následujúci kód:

data sm.data2;
  set sm.data;
  group2 = catx("_", gender, NEDA);
run;

proc sgscatter data=sm.data2;
  title "Scatterplot Matrix (SM data)";
  matrix time age LEMsum timeBef EDSSini EDSS NEDA
         / group = group2 diagonal = (histogram) 
           markerattrs=(symbol=circlefilled)
           transparenc = 0.6;
run;
title;

Analogický postup může byť použitý aj pre nebalancované data, ak je časová doména vhodne diskretizovaná, avšak vhodnejší postup je založený na analýze variogramu, resp. semi-variogramu.

Užitočné

  • V programe R je pre analýzu longitudinálnych dat a odhadovanie lineárneho regresného modelu s náhodnými efektami určená (napr.) knižnica nmle, V rámci tejto knižnice je k dispozícii funkcia Variogram(), ktorá počíta výberový variogram, resp. semi-variogram.
  • Alternatívnou možnosťou je napr. funkcia variogram() dostupná v R knižnici gstat (install.packages("gstat"))
  • Zaujímavý prehľad k variogramu a modelovaniu kovariancie obecně, je napr. dostupný na tejto stránke.



Linárny regresný model s náhodnými efektami môže byť po dekompozícii nepozorovaných, latentných náhodných chýb, zapísany aj v podrobnejšom tvare ako v \[ \boldsymbol{Y}_i = \mathbb{X}_i\boldsymbol{\beta} + \mathbb{Z}_i\boldsymbol{w}_i + \boldsymbol{\varepsilon}_{1 i} + \boldsymbol{\varepsilon}_{2 i}, \] kde jednotlivé stochastické členy v zápise postupne predstavujú
  • \(\boldsymbol{w}_i\) – medzi-subjektovú variabilitu (between-subject variability);
  • \(\boldsymbol{\varepsilon}_{1 i}\) – tzv. sériovu koreláciu (serial correlation);
  • \(\boldsymbol{\varepsilon}_{2 i}\) – tzv. chyby merania (measurement error).



Ak budeme uvažovať model pre strednú hodnotu (t.j. pre časť \(\mathbb{X}_i\boldsymbol{\beta}\)) , tak chybovú (reziduálnu) časť můžeme následne odhadnúť pomocou \(N \in \mathbb{N}\) rezídui (vektorov rezídui) \[ \boldsymbol{r}_i = (r_{i 1}, \dots, r_{i n_i})^\top, \] kde \(r_{i j} = Y_{i j} - \mathbb{X}_i\widehat{\boldsymbol{\beta}}\), pre \(i = 1, \dots, N\). Následne lze predpokládať, že rezídua je možné charakterizovať pomocou modelu (resp. rozhladu) ako \[ \boldsymbol{r}_i = \mathbb{Z}_i\boldsymbol{w}_i + \boldsymbol{\varepsilon}_{1 i} + \boldsymbol{\varepsilon}_{2 i}. \]

Užitočné

  • Ak v regresnom modeli predpokládame len náhodný intercept, pre rozptyl (rezídua) to má za následok predpoklad homoskedasticity (rozptyl je konštantný). Celkovú variančnú-kovariančnú maticu \(\boldsymbol{Y}_i\), resp. \(\boldsymbol{r}_i\), môžeme vyjadriť v tvare \[ Var(\boldsymbol{Y}_i) = Var(\boldsymbol{r}_i) = \nu^2 \mathbb{Z}_i\mathbb{Z}_i^\top + \sigma^2\mathbb{I}_{m_i} + \tau^2 \mathbb{H}_i, \] kde \(\mathbb{Z}_i\mathbb{Z}_i^\top\) je matica typu \(n_i \times n_i\), ktorá má na všetkých pozíciach hodnotu jedna.


  • Korelácia medzi dvoma rôznymi rezíduami \(r_{ij}\) a \(r_{ik}\) (v rámci toho istého \(i\)-teho subjektu) je potom vyjadrená ako \[ \rho(|t_{ij} - t_{ik}|) = \frac{\nu^2 + \tau^2 g(|t_{ij} - t_{ik}|)}{\nu^2 + \sigma^2 + \tau^2}. \]


  • Pre \(j \neq k\) teda platí, že \[ \frac{1}{2}E(r_{ij} - r_{ik})^2 = \sigma^2 + \tau^2 (1 - g(|t_{ij} - t_{ik}|)) = v(u_{ijk}), \] kde posledný člen sa nazýva variogram (semi-variogram)—t.j., funkcia \(v(\cdot)\), ktorá závisí pouze na časových okamžikoch jednotlivých meraní, t.j. \(t_{i j}\) a \(t_{i k}\). Klesajúca funkcia \(g(\cdot)\) má za následok rastúci priebeh variogramu, pričom \(v(0) = \sigma^2\) a \(\lim_{u \to \infty} v(u) = \sigma^2 + \tau^2\).

  • Príslušný empirický odhad teoretického variogramu (semi-variogramu) môže byť efektívne použitý k exploratívnej analýze variančnej-kovariančnej štruktúry dat.



Samostatne

  • Pre ilustraciu budeme opäť uvažovať datový súbor pacientov so sklerózou multiplex, ale namiesto scatterplot matice využijeme variogram aplikovaný na rezídua. Pomocou SAS helpu pre štatistický program SAS sa podívajte na procedúru PROC VARIOGRAM.

  • Jednoduchá ukážka aplikácie variogramu v programe SAS (na základe odhadnutých rezídui z jednoduchého regresného modelu):

    /* Jednoduchy model pre strednu hodnotu */
proc glm data=sm.data;
model EDSS=time;
output out=out r=residual;
run;

/* Zoradenie podla id a time */
proc sort data=out;
  by id time;
run;

/* Subject-specific variogram*/
proc variogram data=out;
  by id;
  coordinates xc=time yc = id;
  compute robust lagd=1 maxlag=8;
  var residual;
run;

/* Celkovy variogram */
proc variogram data=out;
  coordinates xc=time yc = id;
  compute robust lagd=1 maxlag=8;
  var residual;
run;

    Alternatívny (nie ekvivalentný) spôsob, ako spočítať empirický variogram (z rezídui z daného modelu):

    /* Vytvorenie jednotlivych dvojic pozorovani */
proc sql;
create table pairs as
select 
    a.id,
    a.time as time1,
    b.time as time2,
    a.residual as v1,
    b.residual as v2,
    abs(a.time - b.time) as distance
from out as a
join out as b
on a.id = b.id
where a.time < b.time;
quit;

proc print data = pairs; 
run; 

/* Empiricky variogram */
data variogram;
    set pairs;
    vario = (v1 - v2)**2 / 2;
run;

/* Vytvorenie premennej lag */
data variogram;
    set variogram;
    lag = distance;   /* or group: floor(distance) */
run;

/* Priemer cez jednotlive subjekty */
proc means data=variogram;
    class lag;
    var vario;
    output out=vario_avg mean=gamma;
run;

proc print data = vario_avg; 
run; 

proc sgplot data=vario_avg;
    scatter x=lag y=gamma;
    series x=lag y=gamma;
    xaxis label="Time lag";
    yaxis label="Semivariance";
    title "Empirical Variogram";
run;
  • Namiesto rezídui modelu je možné uvažovať aj namerané hodnoty závislej premennej – teda EDSS. Zmodifikujte vyššie uvedený SAS kód a vytvorte empirický variogram na základe EDSS.



2. Odhadovanie modelu s náhodnými efektami v SAS

Základným nástrojom pre odhadovanie lineárnych regresných modelov s náhodnými efektami v programe SAS je procedúra PROC MIXED. Parametre sú odhadované pomocou metódy maximálnej vierohodnosti, pomocou metódy rezidálnej vierohodnosti (REML – default v programe SAS) a prípadne dalšími iterativnými numerickými algoritmami (užitočné v prípade komplexných modelov s veľkým počtom pozorovaní).



Pre správnu implementáciu modelu sú dôležité následújúce časti procedúry:

  1. CLASS Statement
    Definícia faktorových/diskrétných premenných. Pre program SAS je nutné definovať aj čas opakovaných pozorovaní (v rámci spojitého stochastického procesu \(W_I(t)\)) ako diskrétnú premennú–t.j., pozorovania v rámci diskreétnych časových okamžikoch \(t_{i1} < \dots < t_{i n_i}\)
  2. MODEL Statement
    Špecifikácia podmienenej strednej hodnoty \(\mathbb{X}_i \boldsymbol{\beta}\)–t.j. pevných efektov. Dodatočný parameter solution zabezpečí vypísanie odhadnutých pevných efektov–odhadnutej štruktúry podmienenej strednej hodnoty. Rôzne parametrizácie pre faktorové premenné je možne ovplyvniť voľbou parametra noint.
  3. RANDOM Statement
    Špecifikácia náhodných efektov a to vrátane náhodného interceptu, ktorý musí byť uvedený explicitne. Potrebná je aj špecifikácia subjektov, v rámci ktorých sa predpokladá nezávislosť. Parameter type = upresňuje variančnú-kovariančnú maticu \(\mathbb{G}\) pre vektor náhodných efektov. V návode (help) procedúry pozri aj parametre g, gcorr, v a vcorr.
  4. REPEATED Statement
    Špcifikácia poradia opakovaných pozorovaní v rámci konkrétneho subjektu. Opakované pozorovania sú špecifikované pomocou faktorovej/ordinálnej premennej (viď 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. 



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

proc mixed data = sm.data2 method = reml; 
class id gender timeCls;
model EDSS = gender age time time2 age*time  / solution;
random intercept  / type = vc subject = id;
repeated timeCls / type = simple subject = id;
run;



Samostatne

Všimnite si, že v prípade variogramu spočítanom vyššie, boli použité rezidua z jednoduchého regresného modelu, kde EDSS záviselo pouze na čase (jednotlivé roky od začiatku liečby). Model, ktorý je prezentovaný vyššie (asi nie finálny model, ale určite model, ktorý je zo štatistického hľadiska vhodnejší než model so závislosťou EDSS pouze na čase), používa komplexnejšiu štruktúru pre modelovanie podmienenej strednej hodnoty.

  • Spočítajte príslušný variogram na základe rezídui z komplexného modelu a porovnajte oba variogramy.
  • Aké sú závery plynpce z variogramu na základe rezídui z komplexnejšieho modelu?
  • Pokúste sa model dopracovať a podmienenú strednú hodnotu modelovať čo najlepším modelom (v rámci frameworku lineárnych regresných modelov). Využijte rezídua modelu a opäť spočítajte príslušný variogram (semi-variogram).



3. Inferencia v modeloch s náhodnymi efektami

Otázka inferencie v lineárnom regresnom modeli s náhodnými efektami sa týká jednak štruktúry podmienenej strednej hodnoty (t.j., vektor nezámych parametrov \(\boldsymbol{\beta} \in \mathbb{R}^p\)), ale taktiež aj neznámych parametrov \(\boldsymbol{\alpha} \in \mathbb{R}^q\), ktoré parametrizuju podmienenú variančnú/kovariančnú štruktúru.

Niektoré (najčastejšie používané) základne inferenčné nástroje:



V nasledujúcom SAS kóde si všimnite úlohu jednotlivých parametrov (pomocou návodu zistite ich význam) –inferencia o parametroch \(\boldsymbol{\beta}\in \mathbb{R}^p\) (contrast statement a voľba (napr.) chisc) a tiež inferencia o parametroch \(\boldsymbol{\alpha} \in \mathbb{R}^q\) (voľba covtest v proc mixed statement).

V nasledujúcich častiach je postupne robená inferencia vzhľadomk k \(\boldsymbol{\beta}\in \mathbb{R}^p\) parametrom pomocou približného Waldovho testu, pomocou \(t\)-testu (\(F\)-testu resp.) s Kenward-Roger aproximáciou stupňov voľnosti, robustná inferencia, inferenia ohľadom náhodných efektov a individuaálne profily.

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


/* Model pomocou REML -- inferencia o parametroch pre strednu hodnotu je nepouzitelna */
proc mixed data = sm.data2 method = reml empirical; 
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*time 1, 
                                    age*time2 1 / chisq;
run;

Ohľadom inferencie pre variančnú-kovariančnú štruktúru je možné využiť aj vierohodnosť počítanú z rezidui, pokiaľ je štruktúra podmienenej strednej hodnoty (pevné efekty v modeli) rovnaká. Zameriame sa na jednodlivé variančné komponenty.

/* Model s neobmedzenou (UNstructured) kovariancnou strukturou */
proc mixed data = sm.data2 method = reml covtest cl; 
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 group = gender;
run;


/* Model s VC (variance components) strukturou */
proc mixed data = sm.data2 method = reml covtest cl; 
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 group = gender;
run;


/* Model s homoskedasticitou v ramci pohlavi */
proc mixed data = sm.data2 method = reml covtest cl; 
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;



Rozdiel v logaritmickej (restricted) vierohodnosti teda je \(1579.40 - 1564.79 = 14.61\), čo je hodnota (konkrétna realizácia) náhodnej testovej štatistiky, ktorá má za platnosti nulovej hypotézy \(\chi^2\) rozdelenie s jedným stupňom voľnosti. Významnosť náhodných efektov je možné testovať aj pomocou testu pomerom vierohodnosti (ak je zachovaná štruktúra podmienenej strednej hodnoty – pevné efekty v modeli).

Je nutné si ale uvedomiť, že v podstate testujeme nulovú hypotézu na hranici parametrického priestoru, pretože nulová hypotéze o homogenite medzi mužskými a ženskými pacientami, teda \[ H_0: \sigma_{males}^2 = \sigma_{females}^2 \] môže byť ekvivalentne prepísana aj ako \[ H_0: \delta = 0 \] kde predpokládame, že \(\sigma_{males}^2 = \sigma^2 > 0\), \(\sigma_{females}^2 = \sigma^2 = \delta\) a \(\delta > 0\).

Príslušný štatistoclý test je teda na hranici parametrického priestoru pre \(\delta > 0\). Výsledné rozdelenie testovej štatistiky preto nemá \(\chi^2\) rozdelenie s jedným stupňom voľnosti (ako by tomu bolo v prípade testu mimo hranicu paramtrického priestoru), ale jedná sa o rovnoměrnou změs dvoch rozdelení – \(\chi^2\) rozdelenie s jedným stupňom voľnosti a \(\chi^2\) rozdelenie s nula stupňami voľnosti.

V praxi sa často používa \(\chi^2\) rozdelenie s jedným stupňom voľnosti, avšak výsledná \(p\)-hodnota býva mierne konzervatívna.

/* FPlny model s heteroskedasticitou */
  proc mixed data = sm.data2 method=ml  covtest;
    class id gender timeCls;
    model EDSS = gender age time time2 age*time age*time2;
    random intercept time / type = vc subject = id;
    repeated timeCls / type = simple subject = id group = gender;
    ods output FitStatistics=full_fit;
  run;

  /* Model za platnosti nulovej hypotezy (homoskedasticita) */
  proc mixed data = sm.data2 method=ml  covtest;
    class id gender timeCls;
    model EDSS = gender age time time2 age*time age*time2;
    random intercept time / type = vc subject = id;
    repeated timeCls / type = simple subject = id;
    ods output FitStatistics=reduced_fit;
  run;

  /* Extract -2 Log Likelihoods */
  data lrt_result;
    ll_diff = 1530.05 - 1515.34;
    df = 1;
      
    p0 = (ll_diff = 0);
    p1 = 1 - cdf('CHISQ', ll_diff, 1);

    /* Mixture p-value: 0.5 * p0 + 0.5 * p1 */
    p_value = 0.5 * p0 + 0.5 * p1;

    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;