NMST539 | Cvičenie 12 Zmesové (mixture) modely LS 2020/2021 | 18/05/21 | (online výuka) zdrojový Rmd súbor (kódovanie UTF8) Outline 12. cvičenia: teoretické základy zmesových distribučných modelov; niektoré zmesové modely v programe R; praktické aplikácie na reálných datach; Štatisticky program R je k dispozícii (GNU public licence) na adrese https://www.r-project.org RStudio – “user-friendly” interface (jeden z mnohých, ktoré na internete nájdete): RStudio. Užitočné návody a manuály pre prácu s programom R: Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor) Komárek, A.: Základy práce s R. (PDF súbor) Kulich, M.: Velmi stručný úvod do R. (PDF súbor) De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808) Odborná literatúra: 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. Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015 Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979. 1. Zmesové (mixture) modely 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. 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é). rm(list = ls()) attach(mtcars) head(mtcars) ## mpg cyl disp hp drat wt qsec vs am gear carb ## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 ## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 ## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 ## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 ## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 ## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1 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ý: ### 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. 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) 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á). 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. 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). 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. Samostatne 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: 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") } 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)? z 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? Zmesový model je možné samozrejme vylepšiť zohľadnením heteroskedasticity. Modifikácia výsledného zmesového modelu je vpodstate okamžitá. ### 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? 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á. Č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”. Pre zmesové modely v programe R je k dispozícii celá knižnica mixtools. 2. Zmesové modely v regresii 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. 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’. 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: Zuur and Ieno (2016). Beginner’s Guide to Zero-Inflated Models with R. Domáca (samostatná) úloha (Deadline: Cvičenie č.13 / 25.05.2021) 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). 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? Riešenie umiestnite na svoju webovú stránku, najneskôr v utorok, 25.05.2021, do 14:00.