<table width = 900 border = 0><tr><td>
# <b>NMST539 | Cvičenie 12</b>

## <b>Zmesové (mixture) modely</b>
### <font color = "#880000"><b>LS 2020/2021 |  18/05/21 | (online výuka)</b></font>
###### <a href = "cvicenie2021_12.Rmd">zdrojový Rmd súbor (kódovanie UTF8)</a>

**************


<b>Outline 12. cvičenia:</b> 

<ul>
<li>teoretické základy zmesových distribučných modelov;</li>
<li>niektoré zmesové modely v programe R;</li>
<li>praktické aplikácie na reálných datach;</li>
</ul><br>


<b>Štatisticky program R</b> je k dispozícii (GNU public licence) na adrese <a href = "https://www.r-project.org">https://www.r-project.org</a><br><br>

<b>RStudio</b> -- "user-friendly" interface (jeden z mnohých, ktoré na internete nájdete): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Užitočné návody a manuály pre prácu s programom R:</b><br>

<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i>  (<a href = "NMFM301/Rmanual2.pdf">PDF súbor</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i>  (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF súbor</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (<a href = "NMFM301/uvodr.pdf" >PDF súbor</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) </li>		
</ul><br>


<b>Odborná literatúra:</b><br>
<ul>
<li>Bouveyron, C., Celeux , G., Murphy, T.B., and Raftery A.E.: Model-Based Clustering and Classification for Data Science: With Applications in R. Cambridge University Press, 2019.</li>
<li>Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015</li>
<li>Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.</li>
</ul>

<br><br>

**************

  

### <font color = "#880000"><b>1. Zmesové (mixture) modely</b></font>

Zmesové modely predstavujú kategóriu tzv. "model-based" zhlukovacích metód. Pri zhlukovaní (cluster analysis) sa obecně snažíme využiť informáciu v jednotlivých pozorovaniach (mnohorozmerné pozorovania -- rôzne premenné)
k vytvoreniu konkrétnej miery vzájomnej podobnosti a nepodobnosti, pomocou ktorej sú následne jednotlivé subjekty (pozorovania) rozdelené do niekoľkých skupín. Skutočné začlenenie objektov do jednotlivých skupín je na rozdiel od diskriminačných techník ale neznáme. 
To zároveň umožňuje aj v určitom zmysle nejednoznačné priradenie niektorých subjektov do viac ako jedného zhluku (cluster) --  čo opäť z pohľadu diskriminačných metód nie je možné (každý subjekt totíž obsahuje aj konkrétnu informáciu o príslušnosti do danej skupiny). 

Zhlukové algoritmy, ktoré sme spominali, boli zatiaľ výhradne založené na konkrétnej miere podobnosti/nepodobnosti -- tzv. matici vzdialenosti medzi jednotlivými pozorovaniami. Samotná stochastická povaha dat bola ale v určitom zmysle irelevantná (t.j., zhlukovacie algoritmy 
nevyuživali žiaden konkrétny predpoklad ohľadom distribučnej povahy dat -- konkrétne podkladové rozdelenie generujúce data).

Fundamentálnym zhlukovacím postupom, ktorý berie do úvahy aj distribučnú povahu podkladových dat je tzv. zmesový model -- mixture model (nemýliť si zmesový "mixture" model s modelom s náhodnými efektami -- tzv. "mixed effects model"). 
Zmesový model štandardne predpokladá, že jednotlivé pozorovania sú náhodne vyberané z obecně $K \in \mathbb{N}$ rôznych populácii (počet populácii je opäť apriórny predpoklad), ale namiesto jednoduchémo priradenia pouze príslušnosti do určitej populácie (pridelenie nálepky ku každému pozorovaniu, jak fungujú distribution-free zhlukovacie algoritmy) zmesový model odhaduje konkrétny tvar (distribúciu) príslušných populácii (zhlukov). Výstupom zhlukového "mixture" modelu je preto výrazne komplexnejšia informácia, 
než je tomu v prípade výstupu z ľubovolného distribution-free algoritmu (napr. K-means, alebo hierarhické zhlukovanie -- dendogram). Zároveň ale zmesový model vyžaduje výrazne konkrétnejšie predpoklady na podkladové data -- napr. 
predpoklad určitej parametrickej rodiny pre jednotlivé populácie. <br><br>

Pre ilustráciu využijeme opäť datový súbor `mtcars` o nameranej spotrebe automobilov v Sponejn7ých Štátoch v  80tych rokoch. Data obsahujú 32 pozorovaní (32 automobilov predávaných v danom období na americkom trhu) a 11 rôznych premenných (spojité aj kategorické). 

```{r, warning = F, message = F}
rm(list = ls())
attach(mtcars)
head(mtcars)
```
<br><br>

Vzhľadom k povahe datového súboru budeme apriorne predpokladať existenciu troch rôznych populácii (t.j. $K = 3$) definovaných počtom válcov v danom aute (v datovom súbore sú automobily so štyrma, šiestimi a ôsmimi válcami -- premenná `cyl`).
Z teoretického hľadiska môžeme pre každú populáciu definovať určité kvalitatícne a kvantitatívne charakteristiky -- napr. stredná hodnota, rozptyl, médian, distribučná funkcia, prípadne hustota... a mnohé ďalšie. 
Zároveň je možné s využitím dat odhadnúť proporčné zastúpenie jednotlivých sub-populácii v celkovej populácii (t.j. odhadnúť percentuálne zastúpenie štvorvalcových, šesťvalcových a osemvalcových automobilov). 

S využitím dodatočného predpokladu normality -- normálne rozdelenie napr. premennej `mpg` -- je už ľahké odhadnúť jednotlivé populácie -- t.j. získať výberové stredné hodnoty a výberové rozptyly podmienene pri 
konkrétnom počte válcov. Výsledný zmesový model, ktorý obsahuje tri zmesy -- zhluky -- ktoré sú tvorené troma rôznymi normálnými rozdeleniami, je teda následovný: 

```{r, fig.width=10}
### population mean estimates
mean1 <- mean(mpg[cyl == 4])
mean2 <- mean(mpg[cyl == 6])
mean3 <- mean(mpg[cyl == 8])

### estimated proportions for each population
prop <- table(cyl)/dim(mtcars)[1]

### the overall mixture distribution
mixture <- function(x, prop, means){
  mixValue <- 0
  for (i in 1:length(prop)){
    mixValue <- mixValue + prop[i] * dnorm(x, means[i])
  }
  return(mixValue)
}

xSeq <- seq(10, 34, length = 100)
density <- mixture(xSeq, prop, c(mean1, mean2, mean3))
plot(density ~ xSeq, col = "red", type = "l", xlab = "Miles per Gallon", ylab = "Density")
```

Pre jednoduchosť sme v predchádzajúcom príklade predpokládali jednotkový rozptyl pre každú sub-populáciu -- t.j. pre každé normálne rozdelenie 
spotreby -- premennej `mpg` u štvor, šesť a osemvalcových automobilov. Výsledný zmesový model je postavený na troch zmesiach -- zhlukoch --- 
ktoré sú zmiešané dohromady príslušne podľa odhadnutých proporcií zastúpenia automobilov s konkrétnym počtom válcov. 

```{r, fig.width = 10}
plot(density ~ xSeq, col = "red", type = "l", lwd = 3, xlab = "Miles per Gallon", ylab = "Density")
d1 <- density(rnorm(10000, mean1, 1))
d2 <- density(rnorm(10000, mean2, 1))
d3 <- density(rnorm(10000, mean3, 1))
lines(d1$x, prop[1] * d1$y, col = "black", lwd = 2)
lines(d2$x, prop[2] * d2$y, col = "violet", lwd = 2)
lines(d3$x, prop[3] * d3$y, col = "gray", lwd = 2)
```

<br><br>

Zmesový model -- zmesové rozdelenie --- obsahuje v tomto prípade niekoľko lokálnych maxím -- presnejšie jeden pre každý zhluk -- resp. pre každú subpopuláciu. 
Jednotlivé maxima sú lokalizované v bodov $\widehat{\mu}_k \in \mathbb{R}$, pre $k = 1, 2, 3$, ktoré reprezentujú príslušné odhady podmienených stredných hodnôt 
pri danej informácii o počte válcov. Samozrejme, pre dve subpopulácie, ktoré sú hodne podobné, je možne, že výsledný zmesový model bude unimodálny. 

Z obrázku je tatiež zrejmá nejednoznačnosť priradenia niektorých pozorovaní do príslušnej subpopulácie. Zatiaľ čo automobil so spotrebou 15 míľ na jeden galón 
by sme intuitivne asi jednoznačne zaradili medzi automobily s ôsmimi válcami,  u automobilu so spotrebou 17,5 míľ na jeden galón by sme váhali 
medzí zaradením medzi osemvalcové, alebo šesťvalcové automobily (aposteriórna pravdepodobnosť príslušnosti do jednej, alebo do druhej skupiny sa zdá byť v tomto prípade rovnaká).<br><br>

V praktických situáciach je komplexnosť zmesových modelov posunutá ešte o úroveň vyššie. Napríklad počet populácii $K \in \mathbb{N}$ môže byť taktiež neznámy
a je nutné ho nejakým vhodným spôsobom odhadnúť. Pre odhad jednotlivých parametrov sa štandardne využíva metóda maximálnej věrohodnosti (čo samozrejme predpokladá znalosť 
prametrickej rodiny jednotlivých subpopulácii). Dodatočne podmienky regularity sú ale nutné k získaniu rozumných dohadov.  <br><br>

Rozmyslite si, ako by dopadol zmesový model vyššie (t.j. pre premennú `mpg`, o ktorej predpoládame, že bude mať normálne rozdelenie), ale pre neznámy počet zhlukov (t.j. 
hypoteticky každé pozorovanie by tvorilo jeden samostatný zhluk). Riešenie pomocou metódy maximálnej vierohodnosti je priamočiaré a okamžité. Praktické využitie takéhoto 
zmesového modelu je ale nulové. Regularizačné podmienky by boli v tomto prípade ekvivalentné s obmedzením rozptylu (zdola). <br><br>

Uvedomte si, že v tejto interpretácii je z určitého pohľadu  možné považovať  aj klasický lineárny regresný model za určitý limitný zmesový model.
<br><br>


#### <font color = "#880000">Samostatne</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Podobný problém ale nastane aj v prípade, že rozptyl je obmedzený (pre každú subpopuláciu predpokládame jednotkový rozptyl), ale počet zhlukov je naďalej neobmedzený... 
$$
Maximize \prod_{k = 1}^{K} \prod_{i = 1}^{n_k} f_k(X_{k i}) \equiv Maximize \prod_{k = 1}^{K} \prod_{i = 1}^{n_k} \frac{1}{\sqrt{2 \pi}} \cdot exp \left\{\frac{-1}{2} (X_{k i} - \mu_k)^2\right\}
$$
kde maximalizujeme vzhľadom ku $K \in \{1, \dots, n\}$, a $\mu_1, \dots, \mu_K \in \mathbb{R}$ (kde $n \in \mathbb{N}$ je celkový počet pozorovaní, ktoré sú k dispozícii a $n_k \in \mathbb{N}$ je počet subjektov 
zaradených do $k$-tej subpopulácie -- t.j. pozorovania $X_{k 1}, \dots, X_{k n_k}$). Bez akýchkoľvek dodatočných regularizačných predpokladov bude výsledný zmesový model vyzerať nasledujúco: 
 
```{r, fig.width = 10}
plot(density ~ xSeq, col = "gray", type = "l", lwd = 2, xlab = "Miles per Gallon", ylab = "Density", lty = 2, ylim = c(0,0.5))
for (i in 1:length(mpg)){
  lines(rep(mpg[i], 2), c(dnorm(0), 0), col = "red", lwd = 1, lty = 3)
  points(mpg[i], dnorm(0), pch = 21, cex = 0.5, bg = "red")
}
```
 
 
 </li>
 <li>Ako by ste interpretovali obrázok vyššie? Dokážete obrázok vysvetliť na základe toho, čo viete o normálnom rozdelení (hustote normálneho rozdelenia)?</li>
 z
 <li>Využitie takto vytvoreného zmesového modelu je úplne nepraktické. Ako by ste napr. využili tento zmesový model pre zaradenie nového pozorovania -- napr. automobilu 
 so spotrebou 12 míľ na jeden galón?</li>
 
  </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

Zmesový model je možné samozrejme vylepšiť zohľadnením heteroskedasticity. Modifikácia výsledného zmesového modelu je vpodstate okamžitá. 

```{r, fig.width=10}
### population standard error estimates
sd1 <- sd(mpg[cyl == 4])
sd2 <- sd(mpg[cyl == 6])
sd3 <- sd(mpg[cyl == 8])

### the overall mixture distribution
mixtureHetero <- function(x, prop, means, sd){
  mixValue <- 0
  for (i in 1:length(prop)){
    mixValue <- mixValue + prop[i] * dnorm(x, means[i], sd[i])
  }
  return(mixValue)
}

xSeq <- seq(10, 34, length = 100)
density <- mixtureHetero(xSeq, prop, c(mean1, mean2, mean3), c(sd1, sd2, sd3))
plot(density ~ xSeq, col = "red", type = "l", xlab = "Miles per Gallon", ylab = "Density")
```

Ako by ste v tomto prípade rozhodli o zaradení automobilu so spotrebou 17,5 míľ na jeden galón? Lze nejak 
vhodne zohľadniť variability spotreby osemvalcových automobilov a variabilitu spotreby šesťvalcových automobilov? <br><br>

Zmesový model môže samozrejme obsahovať jednotlivé zmesy z rôznych rozdelení a tieto rozdelenia nemusia byť nutné 
rovnaké. Dôležité je správne naformulovať združenú vierohodnosť, ktorú je následne potrebne maximalizovať za účelom 
získania finálneho zmesového modelu. Alternatívou býva využitie tzv. EM algoritmu. 
Parametrické rozdelenia sú často používané, nakoľko umožňujú explicitné vyjadrenie 
vierohodnosti ako funkcie neznámých parametrov a následna maximalizácia je viacmenej priamočiará. <br><br>

Často používanýcm predpokladom je buď apriórna znalosť výsledného počtu zhlukov, alebo, ak je počet zhlukov neznámy, 
tak sa výsledný zmesový model regularizuje predpokladmi na výsledný tvar -- napr. predpoklad unimodality, alebo log-konkavity a pod. 
Regularizačné predpoklady sú nutné, aby sme zabránili niečomu, čo sa často nazýva "Diracova katastrova". <br><br>

Pre zmesové modely v programe R je k dispozícii celá knižnica <a href = "https://cran.r-project.org/web/packages/mixtools/vignettes/mixtools.pdf">mixtools</a>.<br><br>


### <font color = "#880000"><b>2. Zmesové modely v regresii</b></font>

Zmesové modely sú často aplikované aj v regresnej analýze, v rôznych regresných modeloch. Pre ilustráciu si môžeme predstaviť klasický 
GLM regresný model (s logaritmiským linkom) pre Poissonové počty -- napr. pre modelovanie počtu pozitívnych Covid-19 testov (napr. v priebehu dňa, alebo vrámci nejakej špecifickej lokality a pod.). 
Vhľadom k pomerne nízkej prevalencii je možné predpokladať, že nezanedbateľné množstvo testov bude negatívných (a ich pozorovaný počet nebude korespondovať s predpokladom získaným z klasického 
Poissonového modelu). V takom prípade je možné využiť tzv. "zero-inflated" Poissonov model, ktorý je vlastne zmesovým modelom, ktorý môžeme zapísať následovne: 


$$
Y_i | \boldsymbol{x}_i \sim \left\{ 
\begin{array} 
 00 & \textrm{s pravdepodobnosťou } p \in (0,1);\\
Poiss(\lambda_i) & \textrm{s pravdepodobnosťou  } 1 - p_i.
\end{array}   \right.
$$
Jedná sa o zmesový model z dvoma zmesami -- Diracova miera pre výskyt negativného testu a Poissonov model pre modelovanie celkového počtu pozitívnych testov. 
Nezávisle premenné obsiahnuté v $\boldsymbol{x}_i$ sú tzv. "subject-specific" kovariáty, pomocou ktorých subjekty zaradíme do jednej subpopulácie (zmesy), alebo do druhej. <br><br>

Z teoretického hľadiska je pravdepodobnosť, že budeme sledovať negatívný test (t.j. máme jedinca z prvej subpopulácie negatívnych jedincov) vyjadrená ako 

$$P[Y_i = 0 | \boldsymbol{x}_i] = p + (1 - p) \cdot e^{- \boldsymbol{x}_i^\top \boldsymbol{\beta}}$$ 

a pravdepodobnosť, že budeme sledovať $k > 0$ pozitívnych prípadov, môžeme vyjadriť ako 

$$P[Y_i = k | \boldsymbol{x}_i] = (1 - p) \cdot \frac{(\boldsymbol{x}_i^\top \boldsymbol{\beta})^k}{k!}e^{- \boldsymbol{x}_i^\top \boldsymbol{\beta}}$$

pričom predpokládame, že pre $\lambda_i$ platí GLM regresný model 

$$\lambda_i = E[Y_i | \boldsymbol{x}_i] = \boldsymbol{x}_i^\top \boldsymbol{\beta},$$ 


s vektorom neznámych parametrov $\boldsymbol{\beta} \in \mathbb{R}^p$. Finálny zmesový model získame opǎť pomocou maximalizácie vierohodnostnej funkcie.  
Pre tzv. "zero-inflated" Poissonové procesy je v programe R určená (napr.) funkcia `zeroinfl()` z knižnice 'pscl'. <br><br>

Princip tzv. "zero-inflated" zmesových/regresných modelov je v štatistike hodne častý -- napr. zero-inflated binomial model, zero-inflated negative binomial models, zero-inflated GLM a ďalšie. 
Jedná sa o pomerne bohatú skupinu modelov. S využitím princípu vierohodnosti je vpodstate jednoduché tento "inflation" princíp využiť aj pre iné hodnoty ako nula a v súvislosti 
s ľubovolným iným rozdelením. Napr. špeciálne "zero-inflated" modelom je venovaná aj táto publikácia: <a href = "https://www.nhbs.com/beginners-guide-to-zero-inflated-models-with-r-book">Zuur and Ieno (2016). Beginner's Guide to Zero-Inflated Models with R.</a>  
<br><br>






<br><br>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 

## <font color = "#880000"><b>Domáca (samostatná) úloha </b></font>
### (Deadline: Cvičenie č.13 / 25.05.2021)

<hr size="3" width="100%" noshade style="color:#440000" align="left" />

<ul>
<li>Na vhodnom datovom súbore aplikujte tzv. "model-based clustering" -- tzv. vytvorte vhodný zmesový model. 
Na základe povahy dat uvažujte konkrétny počet zhlukov. Využijte buď analógiu postupu v tomto skripte (ideálne ale pre iné rozdelenie, napr.
častá je  zmes dvoch a viacerých exponenciálnych rozdelení), alebo použijte knižnicu `mixtools`, prípadne aplikujte niektorý z tzv. "zero-inflated" modelov (k dispozícii 
sú rôzne balíčky a knižnice).</li>

<li>Navrhnutý model priamo aplikujte (resp. dostatočnej jasne vysvetlite, ako aplikovať) na zaradenie nových pozorovaní 
do existujúcich zhlukov. Tento princíp zaradenia bol vpodstate veľmi intuitivne zmienený už výššie. Ako by to ale vyzeralo více formálne? </li>


<li> Riešenie umiestnite na svoju webovú stránku, <font color = "red">najneskôr v utorok, 25.05.2021, do 14:00</font>.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

<br><br>

</td></tr></table>
<br><br><br>