Smíšené modely
Bi7491 Regresní modelování
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Co již znáte z minulých hodin?
Užití lineárního regresního modelu – spojité výsledky
(definice, předpoklady – analýza reziduí, praktické užití)
Binární výsledky (např. onemocnění) – práce s logistickou regresí,
specifika, interpretace výsledků – poměr šancí
Analýza deviance, Poissonova regrese, nadměrný rozptyl
Kauzální vztahy – zkreslující faktory, interakce, kauzální diagramy
Příprava dat: kategoriální proměnné, spojité proměnné – transformace,
centrování, škálování
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Co byste po dnešní hodině měli
vědět a umět?
popsat problém s klasickými statistickými metodami v případě úloh
s opakovanými pozorováními u stejných subjektů (skupin)
znát rozdíl mezi pevnými a náhodnými efekty
znát definici smíšeného a longitudinálního modelu
provést základní hodnocení dat se shluky prostřednictvím
popsaných metod
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• Předchozí příklady: průřezový (cross-sectional) design
Měření jsou provedena v jediném časovém okamžiku
Např. srovnání podílu tělesného tuku u 10letých a 15letých
dívek – dvě kohorty, nepárový t-test, možná zkreslení
• Longitudinální design
Stejné dívky, měřeny v 10 a 15 letech – párový t-test
Nyní již všechna měření nejsou nezávislá (předpoklad
standardních statistických technik)
Dívky tvoří „shluky“ v datech – pozorování uvnitř shluku
budou zřejmě podobnější, než v různých shlucích
Motivace
Smíšené modely
Opakování – analýza rozptylu
(viz Biostatistika pro MB)
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad – CHOPN
Maximální inspirační tlak dle stadií
kPa5,3
kPa9,8
9
1
1
1



s
y
n
kPa9,2
kPa6,6
12
2
2
2



s
y
n
kPa5,2
kPa4,5
27
3
3
3



s
y
n
Stadium
III IVII
Celkový průměr („grand mean“)
kPa0,3
kPa4,6
48



s
y
n
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Značení
Součty:
Průměry:
Celková variabilita v souboru:
Variabilita v rámci skupin (reziduální součet čtverců):
Variabilita mezi skupinami (příslušná sledovanému vlivu = proměnné):
   
k
i
n
j iije
i
yYS 1 1
2
)(
  
k
i iiA yynS 1
2
)(
   in
j iji YY 1    
k
i
n
j ij
i
YY 1 1
iii nYy /  nYy / 
Skupinový průměr
(„population mean“)
Celkový průměr
(„grand mean“)
   
k
i
n
j ijT
i
yYS 1 1
2
)( Stupně volnosti: 1 ndfT
Stupně volnosti: kndfe 
Stupně volnosti: 1 kdfA
k skupin, v i-té skupině ni pozorování
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Vztahy mezi odhady variability
Platí:
Dále se dá ukázat, že platí:
Tedy platí, že celková variabilita se dá
rozložit na variabilitu v rámci skupin a
variabilitu mezi skupinami:
)()(   yyyYyY iiijij
Stadium
III IVII
AeT SSS 
 
 
   
  


k
i ii
k
i
n
j iij
k
i
n
j ij
yynyY
yY
i
i
1
2
1 1
2
1 1
2
)()(
)(
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Umělý příklad
Léčba
Pozorovaná
hodnota
Skupinový
průměr
Skupinový průměr
– celkový průměr
Pozorovaná hodnota
– skupinový průměr
Pozorovaná hodnota
– celkový průměr
1 10 12 -4 -2 -6
1 12 12 -4 0 -4
1 14 12 -4 2 -2
2 19 20 4 -1 3
2 20 20 4 0 4
2 21 20 4 1 5
3 14 16 0 -2 -2
3 16 16 0 0 0
3 18 16 0 2 2
Celkový průměr = 16 Součet čtverců = 96 Součet čtverců = 18 Součet čtverců = 114
Stupně volnosti = 2 Stupně volnosti = 6 Stupně volnosti = 8
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Princip analýzy rozptylu
Testová statistika analýzy rozptylu:
ee
AA
k
i
n
j iij
k
i ii
dfS
dfS
kn
yY
k
yyn
F i
/
/
)(
1
)(
1 1
2
1
2






 

  
 
Za platnosti H0 platí:
),1(~ knkFF 
vysvětlená variabilita
reziduální variabilita
(pokud faktor A nic nevysvětluje – výběrové průměry ve skupinách jsou blízké –
testové kritérium blízké 0 – nezamítáme nulovou hypotézu)
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Výsledek analýzy rozptylu
Výsledné počty se standardně zaznamenávají do tzv. tabulky analýzy
rozptylu:
Nulovou hypotézu zamítneme/nezamítneme buď na základě srovnání
výsledné p-hodnoty se zvolenou hladinou významnosti testu α, nebo
srovnáním výsledné F statistiky s kritickou hodnotou (kvantilem) rozdělení
F(k – 1, n – k) příslušnou zvolené hladině významnosti testu α.
Variabilita
Součet
čtverců
Počet stupňů
volnosti
Průměrný
čtverec
F statistika p-hodnota
Mezi skupinami SA = 96 dfA = k – 1 = 2 MSA = 48 F = 16 0,004
Uvnitř skupin Se = 18 dfe = n – k = 6 MSe = 3
Celkem ST = 114 dfT = n – 1 = 8
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Výsledek umělého příkladu
16F14,5)6,2(
95,0
),1(
1 
 FF knk

),1( knkFf 
Na hladině významnosti
α =0,05 zamítáme H0 o
rovnosti středních
hodnot.
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Kontrolní dotaz:
Jak vypadá matice plánu pro uvedený model?
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Kontrolní dotaz:
Jak vypadá matice plánu pro uvedený model?
• jedná se o lineární model s jedním kategoriálním
prediktorem s k hodnotami

































101
101
011
011
001
001










X
iEY
1 iEY
H0:
H1:
1 kiEY 





















 0
0
1
1

k






















 0
0
1
1

k

Smíšené modely
Pevné a náhodné efekty
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• pevný efekt
neznámá konstanta, kterou se snažíme odhadnout
konkrétní parametr – těmi jsme se zabývali doteď
• náhodný efekt
– představuje náhodnou veličinu
– jeho přidání do lineárního prediktoru umožňuje zavést
korelaci mezi pozorováními – např. v podobě „náhodného
interceptu“, který reprezentuje nepozorovatelnou
individuální charakteristiku
– neodhadujeme náhodný efekt (což ani nelze), ale
parametry popisující jeho rozdělení – nezajímají nás efekty
konkrétních jedinců, ale informace o cílové populaci jako
celku
Pevné a náhodné efekty
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• obsahuje pevné i náhodné efekty
• jednoduchý příklad: analýza rozptylu (two-way)
Smíšený model
ijkjiijkY  
pevné efekty náhodné efekty
reziduum
),0(~ 2
 Nijk
),0(~ 2
 Nj
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Smíšený model
ijkjiijkY  
pevné efekty
náhodné efekty
iH i  ,0:0 
),0(~ 2
 Nj
0: 2
0 H
jediný parametr
několik
parametrů
Smíšené modely
Odhad parametrů
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• analýza rozptylu (one-way), vyvážený design
Nejjednodušší model
ijiijY  
náhodné efekty
ai ,...,1
nj ,...,1
skupiny
pozorování
ve skupině
),0(~ 2
 Nij ),0(~ 2
 Ni
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Náhodné efekty vytvářejí korelaci
22
2







• koeficient vnitrotřídní korelace
(intraclass correlation coefficient, ICC)
22
)(   ijYD)( ijYE
),cov( '' jiij YY
',',22
jjii   
',',2
jjii 
',',0 jjii 
rozptyl
v jedné skupině
v různých skupinách
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Součet čtverců – odhad ANOVA
    

a
i
n
j
i
a
i
n
j
iij
a
i
n
j
ij YYYYYY
1 1
2
...
1 1
2
.
1 1
2
.. )()()(
celkový součet čtverců
(SST)
reziduální součet
čtverců (SSE)
součet čtverců
efekt alfa (SSA)
2
)1()(  naSSEE ))(1()( 22
   naSSAE
MSEnaSSE  ))1(/(ˆ2

n
MSEMSA
n
aSSA 



2
2 ˆ)1/(
ˆ 



vyvážený design – n je počet ve skupině
stupně volnosti: a-1an-1 an-a
průměrné čtverce MST MSE MSA
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Odhad metodou maximální
věrohodnosti
• pevné efekty s normálními chybami
• doplnění náhodných efektů – smíšený model
  ZXY
  XY X je matice n x p
matice plánu pro pevné efekty
Z je matice n x q
matice plánu pro náhodné efekty
• v praxi není ANOVA estimátor příliš vhodný
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Odhad metodou maximální
věrohodnosti
• v praxi se využívá tzv. REML (restricted maximum likelihood)
obecně dostáváme méně zkreslené odhady parametrů
),0(~ 2
DN 
))(,(~ 2 T
ZDZIXNY 
  ZXY
),0(~ 2
IN 
Máme věrohodnostní funkci odhad parametrů β, σ2, D
Smíšené modely
Užití smíšených modelů
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Bloky jako náhodné efekty
• blok – experimentální jednotka
– definovány podmínkami experimentu nebo úsudkem
– např. konkrétní jedinec, laboratoř vyhodnocující vzorky, zdravotnické
zařízení, vesnice, rodina, vrh, ...
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• snažíme se porovnat čtyři výrobní procesy A, B, C, D
– definovány a stanoveny, zajímají nás – PEVNÉ EFEKTY
• médium (kukuřičný výluh) lze vždy vytvořit v množství pro
čtyři experimenty
– náhodně utvořeny, nejsou předmětem výzkumu – NÁHODNÉ EFEKTY
• výsledkem je množství získaného penicilinu
Příklad:
Produkce penicilinu
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Produkce penicilinu
• ANOVA s pevnými efekty:
> lmod <- aov(yield ~ blend + treat, penicillin)
> summary(lmod)
Df Sum Sq Mean Sq F value Pr(>F)
blend 4 264 66.00 3.504 0.0407 *
treat 3 70 23.33 1.239 0.3387
Residuals 12 226 18.83
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> coef(lmod)
(Intercept) blend1 blend2 blend3 blend4 treat1
86 6 -3 -1 2 -2
treat2 treat3
-1 3
• POZOR: Jsou využity odlišné typy kontrastů, které srovnávají
výsledek s průměrem (interceptem) – součet všech je 0
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Kategoriální prediktory
Součtové kontrasty
Původní Nové proměnné
treat treat1 treat2 treat3
A 1 0 0
B 0 1 0
C 0 0 1
D -1 -1 -1
0 i
)( 321  iEY
2 iEY
3 iEY
• Stanovena dodatečná podmínka:
1 iEY
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Produkce penicilinu
• Smíšený model:
> mmod <- lmer(yield ~ treat + (1|blend), penicillin)
pevný efekt náhodný efekt
data seskupena podle blend
1 ... jen intercept
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Produkce penicilinu
> summary(mmod)
Linear mixed model fit by REML
Formula: yield ~ treat + (1 | blend)
Data: penicillin
AIC BIC logLik deviance REMLdev
118.6 124.6 -53.3 117.3 106.6
Random effects:
Groups Name Variance Std.Dev.
blend (Intercept) 11.792 3.4339
Residual 18.833 4.3397
Number of obs: 20, groups: blend, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 86.000 1.817 47.34
treat1 -2.000 1.681 -1.19
treat2 -1.000 1.681 -0.59
treat3 3.000 1.681 1.78
Correlation of Fixed Effects:
(Intr) treat1 treat2
treat1 0.000
treat2 0.000 -0.333
treat3 0.000 -0.333 -0.333
pevný efekt
náhodný efekt
pro orientační test použijeme normální
aproximaci t-statistiky
obdobně i test vnořených modelů pomocí
ANOVA (nutno použít ML místo REML)
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Produkce penicilinu
• Odhad náhodných efektů??
Už to nejsou parametry, ale náhodné veličiny
(s nulovou střední hodnotou)
• Lze však spočítat tzv. posteriorní střední hodnotu
• Odhad hodnoty pro některý výluh
Kombinace pevných efektů a posteriorní střední hodnoty:
BEST LINEAR UNBIASED PREDICTOR (BLUP)
> ranef(mmod)$blend
(Intercept)
Blend1 4.2878788
Blend2 -2.1439394
Blend3 -0.7146465
Blend4 1.4292929
Blend5 -2.8585859
Blend1 6
Blend2 -3
Blend3 -1
Blend4 2
Blend5 -4
Z modelu s pevnými efekty:
Odhady se ve srovnání s
původním modelem „scvrkly“:
SHRINKAGE ESTIMATOR
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Produkce penicilinu
• Analýza reziduí
Srovnání s predikovanými hodnotami
-2 -1 0 1 2
-6-4-20246
Theoretical Quantiles
SampleQuantiles
82 84 86 88 90 92
-6-4-20246
Fitted
Residuals
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Vnořené efekty
• Příklad: Laboratorní testy
Měření obsahu tuku ve vaječném prášku
• 6 laboratoří
– v každé dva laboranti
• každý dva vzorky
– u každého dvě měření
ijklijkijiijkl STLY  
cmod <- lmer(Fat ~ 1
+ (1|Lab)
+ (1|Lab:Technician)
+ (1|Lab:Technician:Sample),
data=eggs)
Smíšené modely
Longitudinální data
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• proměnná je u jedince měřena opakovaně
• longitudinální studie se zabývají změnou příslušného
výsledku v čase
• cílem je charakterizovat změnu a faktory které ji ovlivňují
• měření u jedince jsou korelovaná!!!
Longitudinální data
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Pro každého jedince...
náhodné efekty
každý jedinec si „vylosuje“ vektor
rozdělení společné pro celou populaci
),0(~ 2
DNi 
),(~| 2
iiiiii ZXNY  
• u každého jedince (i) je provedeno ni měření ve vektoru Yi
),(~ iii XNY  )(2 T
iiii DZZ 
možná explicitní
autokorelace
pevné efekty
vektor společný pro celou populaci
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• americká studie, vývoj příjmů
• 85 osob, alespoň 11 záznamů mezi 1968-1990
Příklad:
Panelová studie příjmové dynamiky
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Panelová studie příjmové dynamiky
year
income
0
20000
40000
60000
80000
70 75 80 85 90 70 75 80 85 90 70 75 80 85 90
0
20000
40000
60000
80000
0
20000
40000
60000
80000
70 75 80 85 90 70 75 80 85 90
0
20000
40000
60000
80000
• vývoj u 20 osob:
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Panelová studie příjmové dynamiky
• srovnání ženy x muži
• zdá se, že muži mají vyšší plat, u žen však rychleji roste
year
log(income+100)
6
8
10
12
70 75 80 85 90
F
70 75 80 85 90
M
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Panelová studie příjmové dynamiky
• pro každého jedince lze sestavit vlastní model
– vlastní absolutní člen (interpretovatelnost – v roce 1978) i sklon přímky
7 8 9 10 11
0.00.10.20.3
Intercept
Slope
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Panelová studie příjmové dynamiky
• Smíšený model:
> mmod <- lmer(log(income) ~ cyear*sex+age+educ+(cyear|person),psid)
pevný efekt náhodný efekt
data seskupena podle osob
Jaké parametry model bude obsahovat?
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Panelová studie příjmové dynamiky
• Smíšený model:
> mmod <- lmer(log(income) ~ cyear*sex+age+educ+(cyear|person),psid)
pevný efekt náhodný efekt
data seskupena podle osob
Jaké parametry model bude obsahovat?
ijijj
jajeijysjsiyijincome




year
ageeducyearsexsexyear)log(
10
... i-tý rok u j-té osoby
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Příklad:
Panelová studie příjmové dynamiky
Linear mixed model fit by REML
Formula: log(income) ~ cyear * sex + age + educ +
(cyear | person)
Data: psid
AIC BIC logLik deviance REMLdev
3840 3894 -1910 3786 3820
Random effects:
Groups Name Variance Std.Dev. Corr
person (Intercept) 0.28166 0.53071
cyear 0.00240 0.04899 0.187
Residual 0.46727 0.68357
Number of obs: 1661, groups: person, 85
Fixed effects:
Estimate Std. Error t value
(Intercept) 6.674178 0.543334 12.284
cyear 0.085312 0.008999 9.480
sexM 1.150315 0.121293 9.484
age 0.010932 0.013524 0.808
educ 0.104212 0.021437 4.861
cyear:sexM -0.026307 0.012238 -2.150
pevný efekt
náhodný efekt
pro velký vzorek lze použít
normální aproximaci t-statistiky
(cyear – centrovaný)
Smíšené modely
Závěr
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Co byste po dnešní hodině měli
vědět a umět?
popsat problém s klasickými statistickými metodami v případě úloh
s opakovanými pozorováními u stejných subjektů (skupin)
znát rozdíl mezi pevnými a náhodnými efekty
znát definici smíšeného a longitudinálního modelu
provést základní hodnocení dat se shluky prostřednictvím
popsaných metod
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
• MIT Growth and Development Study
• Průvodní text a data k dispozici na adrese:
http://www.biostat.harvard.edu/~fitzmaur/ala/fat.txt
1. Prostudujte si text ke studii
2. Načtěte příslušná data do software R
3. Proveďte základní popis longitudinálních dat
4. Sestavte smíšený model pro vývoj podílu tělesného tuku
5. Odpovězte na otázku, zda růst podílu tělesného tuku je
stejný před a po menarche
Praktický úkol
Ondřej Májek, 2018
Bi7491 Smíšené modely
Lékařská fakulta Masarykovy univerzity
Institut biostatistiky a analýz
Použitá literatura
Julian J. Faraway: Extending the Linear Model with R
Garrett M. Fitzmaurice a kol.: Applied Longitudinal Analysis
viz také
http://support.sas.com/documentation/cdl/en/statug/63033/HTML/de
fault/viewer.htm#statug_mixed_sect022.htm
https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html