Bodo Winter, Statistics for Linguists: An Introduction Using R (poznámky a příklady)
290, 300 hezké příklady vztahu t-testu, ANOVA a lm atd. Ukazatele t vycházejí zcela stejně. Též chi2 nezávislost.
Zajímavé funkce
36 grid.arrange
82,83 broom tidy glance
140 count() jako alternativa k table, pro více podskupin přehlednější výpis
Transformace proměnných
86 centering, standardizing - jen zmínka, 98 příklad v R
- Centrování mění pouze intercept, nemění slope ani r2 ani p-hodnoty a může vést k lepším interpretacím (viz později hlavní efekty)
- Standardizace mění slope, zkrátka mění měřítka jednotek. Vhodné pro srovnání veličin s velmi rozdílnou dynamikou. P-hodnotu ani r2 nemění, viz 108
91 log a exp + 96 příklad
Demonstrace TB: Poměry reakčních dob, transformace pomocí log(). Ukazuji 4 příklady různých poměrů časů. Vždy počítám několik možných koeficientů, které porovnávají jejich trvání. Logaritmování je určitě nejvhodnější. Z příkladů je snad zřejmý postup různých převodů včetně získání rozdílů v procentech a milisekundách. Důležité je rozumět, jak pro tu samou realitu vycházejí pochopitelně jiné koeficienty, které ovšem znamenají totéž.
Z příkladů 2 a 3 vyplývá, že log(t2) - log(t1) je lineární veličina vhodná pro původně poměrové veličiny.
# Př. 1: t1 = t2
t1 <- 0.100 # s
t2 <- 0.100 # s
t2 - t1 # 0 = stejné trvání, ale nelineární veličina, "rozdíl"
t2/t1 # 1 = stejné trvání, "poměr"
t2/t1 * 100 # "poměr v procentech", 100 = stejné trvání
log(t2) - log(t1) # 0 = stejné trvání, lineární veličina. "Logaritmus poměru = rozdíl logaritmů", doporučuji.
log(t2/t1) # totéž
exp(log(t2) - log(t1)) # 1 = už zase poměr, totéž co t2/t1 (zde ukazuji, jak z doporučené veličiny získat zpět poměr)
t2new <- t1 * exp(log(t2) - log(t1)) # totéž co t2, ale vypočteno na základě inverze rozdílu logaritmů
t2new - t1 # rozdíl v sekundách, 0 = stejné trvání, nelineární veličina
# Př. 2: t2 = 200 % t1
t1 <- 0.100 # s
t2 <- 0.200 # s
t2 - t1 # 0.1 => t2 trvá o 100 ms více
t2/t1 # 2 = dvojnásobné trvání, poměr
t2/t1 * 100 # poměr v procentech, 200 = dvojnásobné trvání, o 100 % víc než t1
log(t2) - log(t1) # 0.69 = kladné číslo znamená delší trvání, lineární veličina
log(t2/t1) # totéž
exp(log(t2) - log(t1)) # 2 = už zase poměr, totéž co t2/t1
t2new <- t1 * exp(log(t2) - log(t1)) # totéž co t2, ale vypočteno na základě inverze rozdílu logaritmů
t2new - t1 # rozdíl v sekundách, 0.1 => t2 trvá o 100 ms více, nelineární veličina
# Př. 3: t2 = 50 % t1
t1 <- 0.100 # s
t2 <- 0.050 # s
t2 - t1 # -0.05 => t2 je o 50 ms kratší
t2/t1 # 0.5 = poloviční trvání, poměr
t2/t1 * 100 # poměr v procentech, 50 = poloviční trvání, o 100 - 50 = 50 % kratší než t1
log(t2) - log(t1) # -0.69 = záporné číslo znamená kratší trvání, lineární veličina (srovnej viz předchozí příklad +0.69)
log(t2/t1) # totéž
exp(log(t2) - log(t1)) # 0.5 = už zase poměr, totéž co t2/t1
t2new <- t1 * exp(log(t2) - log(t1)) # totéž co t2, ale vypočteno na základě inverze rozdílu logaritmů
t2new - t1 # rozdíl v sekundách, -0.05 => t2 je o 50 ms kratší, nelineární veličina
# Př. 4: t2 je 80 % t1, tedy o 20 % kratší než t1
t1 <- 0.200 # s
t2 <- 0.160 # s
t2 - t1 # -0.04 => t2 je o 40 ms kratší
t2/t1 # 0.8 = poměr trvání
t2/t1 * 100 # poměr v procentech, t2 je 80 % t1, je tedy o 100 - 80 = 20 % kratší než t1
log(t2) - log(t1) # -0.22 = záporné číslo znamená kratší trvání, lineární veličina
log(t2/t1) # totéž
exp(log(t2) - log(t1)) # 0.8 = už zase poměr, totéž co t2/t1
t2new <- t1 * exp(log(t2) - log(t1)) # totéž co t2, ale vypočteno na základě inverze rozdílu logaritmů
t2new - t1 # rozdíl v sekundách, -0.04 => t2 je o 40 ms kratší, nelineární veličinalm
90 r2 = Pearsonova korelace na druhou
114 collinearity a vif() pro její detekci
122-123 elegantní automatická predikce středů pro jednotlivé hodnoty factor
124 relevel() pro změnu reference ve factor (intercept)
126 contrasts() a contr.sum() pro sum-coding (vystředění) - pro 2 hodnoty analogie k centering - smysl použít jen pro main effects (viz dále 142), ale i tam opatrně, protože mezi kategoriemi -1 a +1 je krok 2. Rozhodně člověk zvládne main effects i bez sum-coding. Tak i tak, musí se přemýšlet a dávat pozor.
134 interakce kategorická * spojitá 137 Už nelze interpretovat samostatné intercepts. Centrování se vyplatí - koeficient ukazuje rozdíl pro průměrné hodnoty, nikoliv pro 0.
149 interakce kategorická * kategorická
142 simple effects vs main effects
146 interakce spojitá * spojitá 147 připomenutí: každý koeficient je pro nulové ostatní koeficienty. Centrováním získáme efekty pro průměrné hodnoty ostatních veličin. A standardizací srovnáme metriky proměnných. Užitečný obrázek 8.4a. Domyslet rovnici.
153 nelineární - parabolická regrese. Centrovaný prediktor i ten samý na druhou jako vstup.
159 Cohenovo d - effect size mezi dvěma středy
162 effsize, cohen.d vč. conf. int.
182-183 graf intervalových odhadů regresních koeficientů
184-187 post-hoc párové testy pomoci emmeans
188-194 graf nejistoty kategoriálních prediktorů
194-197 nejistota kontinuálních prediktorů, dvě tabulky v jednom grafu
lmer
235 fig 14.1 varying intercepts vs. int. + slopes a 237 fig 14.2
Pozn. v této knize B. Winter preferuje termín “varying” intercept či slope oproti běžnějšímu “random” intercept či slope. Jedná se ale o totéž.
236 fixed effects - vše, co bylo doteď. Opakovatelné při nových verzích experimentu, jejich vliv se předpokládá konstantní napříč těmito experimenty. Mohou být kategoriální i spojité. Random effects vždy jen kategoriální, často participant a item. Většinou nekonečně možností, je to vzorek z populace, zajímá nás spíše jen variabilita.
238-239 vysvětlení výstupů lme4 pro fixed i random effects effects (intercept a slope variation mezi účastníky). A též random effect correlations: trial -0.38 znamená, že vyšší intercepty mívají nižší trial slopes (trial zde znamená index opakovaných pokusů, kde ve fixed obecně vyšlo, že s dalšími pokusy se čas zkracuje), tedy kdo začne na začátku pomalu (vyšší čas interceptu) má větší možnost později zrychlit (snížit čas). Tato korelace interceptu se slopes je dalším parametrem, který je v modelu odhadován a obsahuje užitečnou informaci, proto by neměla být normálně vynechávána (může ukázat tzv. ceiling effect - saturaci - když je pro někoho úkol příliš lehký, je pak málo prostoru pro objevení vlivu nějakého efektu, viz 240 fig 14.3).
241 tab 14.1 syntaxe lme4, jak se případně vyhnout slope/intercept korelaci
Příklad 1: jak přemýšlet o rovnici modelu (242)
Pro které proměnné modelovat varying slopes (jiné označení termínu “random slopes”)?
Obecně vynechání důležitých věcí z modelu, tedy např. varying slopes, vede na příliš optimistické p-hodnoty a je to tedy chyba.
Ne vše lze ale z našich dat odhadnout a krédo “keep it maximal” není vždy reálné.
Příklad - má trénování vliv na reakční dobu? Podmínka pro srovnání je tedy “před” a “po”.
Náhled začátku dat
Výzkumník se rozhodl modelovat vztahem
RT ~ condition + age + (1 + condition|participant) + (1|item)Zásadní otázka: mění se v mých datech konkrétní proměnná napříč jednotlivci?
Zde: condition - účastníci byli vystavěni oběma podmínkám -> můžeme modelovat varying slopes.
A dává to teoreticky smysl? Lidé jsou různí a je velmi pravděpodobné, že na každého to bude působit jinak -> ano, dává to smysl.
Age - mění se v rámci účastníka? Ne. Zde varying slopes nepůjdou. Kdybychom ale měli longitudinální (dlouhodobou) studii, kde by účastníci stárli a my to měli zaznamenané, tak je to jiná situace.
Item - položky opakované napříč účastníky -> závislost, něco unikátního v rámci položky může ovlivnit reakce napříč účastníky. Varying intercepts tedy určitě ano. A co varying slopes? Pro podmínku “před” a “po” byly využité jiné položky (tzv. between-items design na rozdíl od within-items design), tudíž se condition napříč items nemění a varying slopes nelze modelovat.
Age napříč items se mění, takže varying slope (1+age|item) lze modelovat. Ale dává to teoreticky smysl? To bude závislet na konkrétním studii, kterou provádíme. Takováto rozhodnutí by rozhodně měla být prováděna dopředu, ideálně ještě před sběrem dat (aby jich případně bylo dost), rozhodně alespoň před tím, než je začneme analyzovat (pak třeba dostaneme odpověď, že z daných dat něco zjistit nelze a máme pořídit nová data). Rozhodně je nesmysl vymýšlet otázky na základě toho, co data dovolí. Otázky by měly být dány předem a musíme se pak smířit s tím, že model občas řekne, že na něco odpovědět nelze.
Příklad 2: model a umělá data pro lme4 (252-260)
Dopředu víme správné výsledky. Opět můžeme vyhodnocovat funkčnost výpočtu p-hodnoty v souvislosti s hladinou významnosti a schopnost modelu odhalit hypotézy, podobně jako u konfidenčních intervalů a testování hypotéz na začátku semestru.
Navíc pochopení vztahu mezi procesem generování umělých dat a následnou konstrukcí vzorce pro model je nezbytné pro uvědomění si toho, co je schopen model s daným předpisem v datech modelovat.
Vowel durations v závislosti na word frequency. Očekáváme, že běžnější slova budou mít kratší samohlásky. 6 participants x 20 items (stejné položky napříč účastníky)
logfreqs - logaritmus frekvence = transformace linearizující poměrovou veličinu. Násobení *5 je jen z důvodu, aby čísla připomínala běžně se vyskytující hodnoty, není to jinak nic zásadního.
# 15.2. Simulating vowel durations for a mixed model analysis:
library(tidyverse)
set.seed(666) # Set random seed value
ppt_ids <- gl(6, 20) # Generate participant identifiers for 6 participants & 20 items
it_ids <- gl(20, 1) # 20 item identifiers
it_ids <- rep(it_ids, 6) # Repeat each 6 times for participant numbers
logfreqs <- round(rexp(20) * 5 , 2) # Create predictor values (simulated frequencies)
logfreqs <- rep(logfreqs, 6) # Repeat for 6 participants
xdata <- tibble(ppt = ppt_ids, item = it_ids, freq = logfreqs) # Put predictors into tibbleIntercept for vowel duration: 300 ms, napříč účastníky sd = 40 ms normálně rozdělené. Položky mají také každá svoji odchylku, sd = 20 ms. A náhodný šum (zbytková chyba) sd = 20 ms.
Frequency effect: s růstem log frekvence o +1 se změní vowel duration o -5 ms.
Dur je součtem všech těchto dílčích složek, které budeme chtít později v modelu odhalit.
xdata$int <- 300 # Add intercept
ppt_ints <- rnorm(6, sd = 40) # Create participant deviation scores
xdata$ppt_ints <- rep(ppt_ints, each = 20) # Add them, each repeat 20 times to main tibble
item_ints <- rnorm(20, sd = 20) # Same for items
item_ints <- rep(item_ints, times = 6)
xdata$item_ints <- item_ints # Add to tibble
xdata$error <- rnorm(120, sd = 20) # Add general error term
xdata$effect <- -5 * xdata$freq # Add effect
# Put all the different effects into the main response:
xdata <- mutate(xdata,
dur = int + ppt_ints + item_ints +
error + effect)
# For pedagogical purposes, get rid of what's been used to create stuff:
xreal <- dplyr::select(xdata, -(int:effect))
# This is what you would have as the researcher:
xrealxreal[1:6, ]## # A tibble: 6 × 4
## ppt item freq dur
## <fct> <fct> <dbl> <dbl>
## 1 1 1 2.74 314.
## 2 1 2 9.82 238.
## 3 1 3 0.7 366.
## 4 1 4 0.05 245.
## 5 1 5 2.76 304.
## 6 1 6 17.6 198.xreal[21:26, ]## # A tibble: 6 × 4
## ppt item freq dur
## <fct> <fct> <dbl> <dbl>
## 1 2 1 2.74 306.
## 2 2 2 9.82 209.
## 3 2 3 0.7 317.
## 4 2 4 0.05 276.
## 5 2 5 2.76 332.
## 6 2 6 17.6 215.1. model: odpovídající přesně tomu, co jsme generovali
Intercept int a effect -5*freq (v generujícím modelu) jsou systematické složky (fixed) - ty, co jsme doposud zkoumali pomocí lm(). Ostatní jsou náhodné (random). Mimochodem: freq je spojitá proměnná a náhodné efekty jsou vždy faktorové.
# --------------------------------------------------------
# 15.3. Analyzing the simulated vowel durations with mixed models:
library(lme4)
# Create varying intercept model:
xmdl <- lmer(dur ~ freq + (1|ppt) + (1|item),
data = xreal, REML = FALSE)
# REML = FALSE přepíná model do nastavení pro maximum likelihood estimation.
# Inspect:
summary(xmdl)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: dur ~ freq + (1 | ppt) + (1 | item)
## Data: xreal
##
## AIC BIC logLik deviance df.resid
## 1105.1 1119.0 -547.5 1095.1 115
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09699 -0.60947 0.06482 0.60761 2.39757
##
## Random effects:
## Groups Name Variance Std.Dev.
## item (Intercept) 589.3 24.28
## ppt (Intercept) 1296.5 36.01
## Residual 284.0 16.85
## Number of obs: 120, groups: item, 20; ppt, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 337.973 16.735 20.196
## freq -5.460 1.004 -5.438
##
## Correlation of Fixed Effects:
## (Intr)
## freq -0.339Model skutečně odpovídá přesně tomu, co jsme generovali. Fixní itercept = 338, freq = -5.46. Random ppt = 36, item = 24, residual = 17 (Std.Dev.). Docela podobné našim číslům (nemůžeme očekávat zázraky ze 6 účastníků). Samozřejmě coef(model) neobsahuje varying slopes, ty jsme negenerovali ani nemodelovali.
# --------------------------------------------------------
# 15.4. Extracting information out of lme4 objects:
# Fixed effects coefficients:
fixef(xmdl)
# Works as usual, e.g., extract slope:
fixef(xmdl)[2]
# Retrieve coefficient table:
summary(xmdl)$coefficients
# Random effects coefficients:
coef(xmdl)
# Participant random effects:
coef(xmdl)$ppt
# Participant deviation scores:
ranef(xmdl)$ppt2. model: špatný, bez (1|item)
# --------------------------------------------------------
# 15.5. Messing up the model:
# Dropping item random intercepts:
xmdl_bad <- lmer(dur ~ freq + (1|ppt),
data = xdata, REML = FALSE)
summary(xmdl_bad)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: dur ~ freq + (1 | ppt)
## Data: xdata
##
## AIC BIC logLik deviance df.resid
## 1182.1 1193.2 -587.0 1174.1 116
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6270 -0.6397 0.1089 0.7462 1.9778
##
## Random effects:
## Groups Name Variance Std.Dev.
## ppt (Intercept) 1240.8 35.23
## Residual 877.6 29.62
## Number of obs: 120, groups: ppt, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 337.9730 14.8829 22.71
## freq -5.4601 0.4813 -11.34
##
## Correlation of Fixed Effects:
## (Intr)
## freq -0.183Residual se skoro zdvojnásobil ze 16.85 na 29.62. Std. Error u freq je poloviční (0.48 vs. 1), což vypadá sice dobře, ale v tom je právě ten problém - ona to totiž není pravda. Model neví o důležitém zdroji rozptylu a je více anti-conservative, tedy přehnaně optimistický.
3. model: včetně proměnných slopes napříč participants (které jsme ovšem negenerovali)
# Adding random slopes even though we didn't generate them:
xmdl_slope <- lmer(dur ~ freq + (1 + freq|ppt) + (1|item),
data = xdata, REML = FALSE)## boundary (singular) fit: see help('isSingular')summary(xmdl_slope)$varcor## Groups Name Std.Dev. Corr
## item (Intercept) 24.27780
## ppt (Intercept) 35.32440
## freq 0.12079 1.000
## Residual 16.83644# Extract participant slopes:
coef(xmdl_slope)$ppt## (Intercept) freq
## 1 315.7902 -5.535971
## 2 302.3449 -5.581948
## 3 362.9170 -5.374818
## 4 318.3843 -5.527100
## 5 324.7279 -5.505408
## 6 403.6741 -5.235448Varování boundary (singular) fit: problém - modelu nelze věřit. Std.Dev napříč varying slopes pro freq (0.12079) je velmi malá. Vysvětlení: v umělých datech jsme nic takového nevytvářeli, mělo by to být 0. Další nápověda průšvihu: Corr varying slope / intercept je přesně rovna 1. Perfektní korelace v lingvistice nejsou. Interpretovali bychom to takto: větší intercept jde přesně beze zbytku s vyšší frekvencí; jenže díky tomu, že je to singular fit, tak to vůbec interpretovat nemáme.
4. model: ponechme freq|ppt a odstraňme “zlobivou” korelaci
# De-correlate random effects structure:
xmdl_slope_nocorr <- lmer(dur ~ freq +
(1|ppt) + (1|item) +
(0 + freq|ppt),
data = xdata, REML = FALSE)## boundary (singular) fit: see help('isSingular')summary(xmdl_slope_nocorr)$varcor## Groups Name Std.Dev.
## item (Intercept) 24.276
## ppt freq 0.000
## ppt.1 (Intercept) 36.006
## Residual 16.853# Notice that the random slopes are all the same:
coef(xmdl_slope_nocorr)$ppt## (Intercept) freq
## 1 315.3049 -5.460115
## 2 301.6253 -5.460115
## 3 363.4343 -5.460115
## 4 318.0924 -5.460115
## 5 324.4672 -5.460115
## 6 404.9142 -5.460115Boundary (singular) fit, lme4() neodhalilo žádné varying slopes napříč participants. Pokud bychom se detailně nepodívali na random strukturu, neodhalili bychom to. Klidně by to někdo publikoval jako random slope model (což vypadá hezky), ovšem random slopes vůbec nebyly odhadnuty.
Výpočet p-hodnoty
260 Probability = pravděpodobnost výsledku (outcome) při nějaké hodnotě parametrů. Likelihood = věrohodnost, jak věrohodné jsou určité parametry pro pozorovaný outcome. Maximum likelihood estimation = hledání nejvěrohodnějších parametrů modelu pro pozorovaná data.
261 likelihood ratio test pro model s 1 fixním modelem: porovnání anova s modelem jen s interceptem
262 jak referovat o p-hodnotě v textu
262 p-hodnota pro významnost random effect
263 p-hodnoty pro mnoho fixed effects library(afex), mixed(…, method = “LRT”)
264 R2 pro smíšené a generalizované modely
264 konfidenční intervaly pro predikce - jen odkaz
Problémy s konvergencí
265-267 problémy s konvergencí. Naprosto nutné kontrolovat strukturu random effects, random effects correlation nesmí být přesně 0, 1 nebo -1. Dobrý zdroj informací je ?convergence
argument lmerControl může něco vyřešit, ?lmerControl vysvětluje ledasco
library(afex) all_fit(proměnná_s_lmer_modelem) zkusí různé optimalizační algoritmy a řekne, které konvergují. A vždy domyslet, zda má vůbec model v datech k dispozici dané kombinace hodnot efektů, aby se měl čeho chytit, viz Příklad 1 (242) výše.
Příklad 2: pokračování (p-hodnoty)
# --------------------------------------------------------
# 15.6. Likelihood ratio tests:
# Construct null model without fixed effect:
xmdl_nofreq <- lmer(dur ~ 1 + (1|ppt) + (1|item),
data = xreal, REML = FALSE)
# Model comparison:
anova(xmdl_nofreq, xmdl, test = 'Chisq')## Data: xreal
## Models:
## xmdl_nofreq: dur ~ 1 + (1 | ppt) + (1 | item)
## xmdl: dur ~ freq + (1 | ppt) + (1 | item)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## xmdl_nofreq 4 1121.0 1132.2 -556.51 1113.0
## xmdl 5 1105.1 1119.0 -547.55 1095.1 17.933 1 2.288e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Full model:
x_REML <- lmer(dur ~ freq + (1|ppt) + (1|item),
data = xdata, REML = TRUE)
# Reduced model (no item random intercepts):
x_red <- lmer(dur ~ freq + (1|ppt),
data = xdata, REML = TRUE)
anova(x_red, x_REML, test = 'Chisq', refit = FALSE)## Data: xdata
## Models:
## x_red: dur ~ freq + (1 | ppt)
## x_REML: dur ~ freq + (1 | ppt) + (1 | item)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## x_red 4 1174.4 1185.5 -583.2 1166.4
## x_REML 5 1095.8 1109.7 -542.9 1085.8 80.608 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Use afex for testing fixed effects:
library(afex)
xmdl_afex <- mixed(dur ~ freq + (1|ppt) + (1|item),
data = xdata, method = 'LRT')
xmdl_afex## Mixed Model Anova Table (Type 3 tests, LRT-method)
##
## Model: dur ~ freq + (1 | ppt) + (1 | item)
## Data: xdata
## Df full model: 5
## Effect df Chisq p.value
## 1 freq 1 17.93 *** <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1fixef(xmdl_afex$full_model)## (Intercept) freq
## 337.973044 -5.460115Závěrem
271 fig 15.1 shrinkage effect lmer oproti individuálním modelům - lmer vidí věci dohromady a může si dovolit snížit rozptyl odhadů, i když dovolíme individuální intercept a slope random effects
272 poslední odstavec 15.9 random effect structure mohou občas obsahovat skoro zajímavější závěry než fixed effects, vyplatí se neignorovat, jak mnoho lidí činí
Exploratory vs confirmatory statistics
275 exploratorní analýza generuje hypotézy, confirmatory statistics je testuje. Není správné dělat obojí na stejných datech, obdobně jako u klasifikátorů rozdělujeme data na trénovací a testovací, jinak hrozí přetrénování a ztráta generalizace - příliš se soustředíme na VS a je vyšší riziko chyb 1. druhu. Stejnou chybou je “ladění” modelu, tj. přidávání či ubírání prediktorů, dokud nevyjde hezké p. Pokud se neúspěšné modely nereportují, tak je to klasický postup hledání chyby 1. druhu (dokud ji nenajdu, hledám dál - tak dlouho budu generovat rnorm se středem 0, až mi jednou vyjde, že to může mít střed 12… A to jediné pak opublikuji) a hladina významnosti přestává dávat smysl. Hypotézy i rovnice modelu mají být dány před pořízením dat. Rovnice a hladina významnosti určuje, jak mají data vypadat (kolik od čeho). V seriozním výzkumu (medicína) by se měly tyto informace předregistrovat, aby bylo transparentní, kolikrát se co povedlo a kolikrát ne. A pozor, skoro každý matematický výsledek dokáže vědec zdůvodnit. Jeden i druhý se zcela opačným závěrem. Takže ani na “teoretické předpoklady” nelze spoléhat.
276 odstavec When following this cookbook approach…
278 poslední odstavec - ve studii by mělo být jednoznačně řečeno, jestli je to exploratory nebo confirmatory. Obojí je velmi důležité, ale nelze to vágně míchat - dělat jedno a tvářit se, že je to v podstatě i to druhé.
Logistická regrese pro výstupní binární faktorovou proměnnou
198 logistická regrese, 199 fig 12.1
201 plogis() převod přímky (-oo -> -oo, 0 -> 0, +oo -> +oo) na logistickou s-křivku od 0 do 1 (-oo -> 0, 0 -> 0.5, +oo -> +1), viz 202 fig 12.4
202 log odds (logit) = log(p / (1-p) ), kde log je přirozený. Odds = šance, je to poměr pravděpodobnosti, že věc nastane, ku pravděpodobnosti, že nenastane. Padesát na padesát (p=0.5) => 1, log odds = 0. p=0.1 => 0.11 (1/9), log = -2.2. p=0.9 => 9, log = +2.2. Log odds je kontinuální škála od -oo do +oo, pro naše účely (odds je poměr) krásně symetrická. Přehledně: 0 odpovídá p=0.5, kladné odpovídá p>0.5, záporné p<0.5. Převod p na kontinuální lineární škálu umožní napasování úlohy na lm… Logistická funkce plogis() fig 12.4 je inverze k log odds - vstupem log odds, dostáváme p. A logit() je inverze k plogis() - převod p na log odds.
204 jasný příklad od začátku do konce (výstup binární kategorie, vstup kontinuální) včetně tvorby fig 12.1. Ke koeficientům a jejich SE, z a p se dostane přes tidy(model), viz 205
207 příklad (vstup je binární kategoriální proměnná)
210 příklad (vstup je diskrétní intervalová 9bodová škála). A predikce pravděpodobností včetně confint pro tuto škálu. Ovšem dělá to přes +-1.96*SE, neumí predict kromě SE vrátit rovnou i confint?
267 příklad na smíšený model glmer logistický včetně automaticky nalezeného optimizer
Poissonova a NB regrese pro modelování četností
220 glm Poisson pro modelování počtů u procesů vznikajících dle Poissonova rozdělení (rozptyl je rovný střední hodnotě)
225 přidání exposure variable, která má vliv a je asi přímo úměrná modelované četnostní veličině
227 Negative binomial regression: pro počty neřídící se Poissonovým rozdělením (většinou: rozptyl je větší než očekávaný - overdispersion). Zde není nutná rovnost rozptylu se střední hodnotou.
228 test, zda je nutný glm.nb oproti Poisson