x <- seq(1,10, length.out = 100)y <- 1.2 + 2.1 * xplot(x, y)

set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2)y <- 1.2 + 2.1 * x + noiseplot(x, y)

m <- lm(y ~ x)summary(m)

## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -5.5831 -1.2165 0.3288 1.3022 4.3714 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.27426 0.44720 2.849 0.00534 ** ## x 2.09449 0.07338 28.541 < 2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.926 on 98 degrees of freedom## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8915 ## F-statistic: 814.6 on 1 and 98 DF, p-value: < 2.2e-16

set.seed(222)noise <- rnorm(n = 100, mean = 0, sd = 2 * sqrt(x))y <- 1.2 + 2.1 * x + noiseplot(x, y)

m2 <- lm(y ~ x)summary(m2)

## ## Call:## lm(formula = y ~ x)## ## Residuals:## Min 1Q Median 3Q Max ## -15.0460 -2.4013 0.5638 2.8734 10.2996 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.313 1.042 1.26 0.211 ## x 2.096 0.171 12.26 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 4.488 on 98 degrees of freedom## Multiple R-squared: 0.6051, Adjusted R-squared: 0.6011 ## F-statistic: 150.2 on 1 and 98 DF, p-value: < 2.2e-16

plot(m2, which = 1)

plot(m2, which = 3)

library(nlme)vm1 <- gls(y ~ x, weights = varFixed(~x))summary(vm1)

## Generalized least squares fit by REML## Model: y ~ x ## Data: NULL ## AIC BIC logLik## 576.2928 584.0477 -285.1464## ## Variance function:## Structure: fixed weights## Formula: ~x ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.369583 0.6936599 1.974431 0.0511## x 2.085425 0.1504863 13.857900 0.0000## ## Correlation: ## (Intr)## x -0.838## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.8942967 -0.6293867 0.1551594 0.6758773 2.2722755 ## ## Residual standard error: 1.925342 ## Degrees of freedom: 100 total; 98 residual

plot(vm1)

m3 <- lm(y ~ x, weights = 1/x)summary(m3)

## ## Call:## lm(formula = y ~ x, weights = 1/x)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -5.5725 -1.2118 0.2987 1.3013 4.3749 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.3696 0.6937 1.974 0.0511 . ## x 2.0854 0.1505 13.858 <2e-16 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 1.925 on 98 degrees of freedom## Multiple R-squared: 0.6621, Adjusted R-squared: 0.6587 ## F-statistic: 192 on 1 and 98 DF, p-value: < 2.2e-16

varIdent

Funkce varIdent nám umožňuje modelovat různé rozptyly pro jednotlivé vrstvy. Abychom viděli, jak funguje, budeme nejprve simulovat data s touto vlastností. Níže použijeme set.seed(11) pro případ, že by někdo chtěl simulovat stejná náhodná data. Poté nastavíme n na hodnotu 400, což je počet pozorování, která budeme simulovat. Hodnotu x vygenerujeme stejně jako předtím. V tomto příkladu zahrneme další prediktor s názvem g, který si můžeme představit jako pohlaví. Náhodně vybereme 400 hodnot „m“ a „f“. Dále nasimulujeme vektor dvou směrodatných odchylek a uložíme jako msd. Je-li g == "m", pak směrodatná odchylka je \(2 \krát 2,5\). V opačném případě je to \(2\) pro g == "f". Tento vektor použijeme v dalším řádku pro vygenerování y. Všimněte si, kde je msd zapojeno do funkce rnorm. Všimněte si také, jak generujeme y. Zahrnujeme interakci mezi x a y. Když je g == "f", je intercept a sklon 1,2 a 2,1. Při g == "m" je průsečík (1,2 + 2,8) a sklon (2,1 + 2,8). Nakonec umístíme naše simulovaná data do datového rámce a vykreslíme pomocí ggplot2.

set.seed(11)n <- 400x <- seq(1,10, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)library(ggplot2)ggplot(d, aes(x, y, color = g)) + geom_point()

Všimněte si různých rozptylů pro každou úroveň g. Rozptyl pro „m“ je mnohem větší než rozptyl pro „f“. Má mnohem větší rozptyl než „f“. Takto jsme data simulovali. Nastavili jsme rozptyl pro „m“ na 2,5násobek rozptylu pro „f“.

Nyní použijeme funkci gls s varIdent ve snaze získat tyto skutečné hodnoty. Použijeme stejný způsob jako dříve: definujeme naši funkci rozptylu v argumentu weights. Dále v argumentu form uvedeme, že vzorec pro modelování rozptylu je podmíněn g. Výraz ~ 1|g je jednostranný vzorec, který říká, že rozptyl se liší mezi úrovněmi g. Výraz 1 pouze znamená, že v našem modelu pro rozptyl nemáme žádné další kovariáty. Je možné zahrnout vzorec jako ~ x|g, ale to by v tomto případě nebylo správné, protože jsme při generování rozptylu nepoužili x. Pohled na graf také ukazuje, že zatímco variabilita y se mezi skupinami liší, variabilita se s nárůstem x nezvyšuje.

vm2 <- gls(y ~ x*g, data = d, weights = varIdent(form = ~ 1|g))summary(vm2)

## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 2081.752 2105.64 -1034.876## ## Variance function:## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.00000 2.58127 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 0.9686349 0.3237724 2.99172 0.0029## x 2.1222707 0.0525024 40.42239 0.0000## gm -1.9765090 0.9352152 -2.11343 0.0352## x:gm 2.9957974 0.1553551 19.28355 0.0000## ## Correlation: ## (Intr) x gm ## x -0.901 ## gm -0.346 0.312 ## x:gm 0.304 -0.338 -0.907## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74115039 -0.67013954 0.01619031 0.69793776 2.72985748 ## ## Residual standard error: 2.005397 ## Degrees of freedom: 400 total; 396 residual

Nejprve si všimněte v dolní části výstupu, že odhadnutá směrodatná chyba reziduí je 2. To znamená, že v tomto případě se jedná o 2 %.0053974, což je velmi blízko „skutečné“ hodnotě 2. Všimněte si také, že v části „Funkce rozptylu“ dostáváme odhadovanou hodnotu 2,58127 pro skupinu „m“, což je velmi blízko „skutečné“ hodnotě 2,5, kterou jsme použili k vytvoření rozdílného rozptylu pro skupinu „m“. Obecně platí, že když použijete funkci varIdent k odhadu různých rozptylů mezi úrovněmi vrstev, jedna z úrovní bude nastavena jako základní a ostatní budou odhadnuty jako násobky reziduální směrodatné chyby. V tomto případě, protože písmeno „f“ je abecedně před písmenem „m“, bylo písmeno „f“ nastaveno na základní hodnotu neboli 1. Odhadnutá reziduální standardní chyba pro skupinu „f“ je \(2,005397 \krát 1\). Odhadovaná zbytková standardní chyba pro skupinu „m“ je \(2,005397 \krát 2,58127\).

Je důležité si uvědomit, že se jedná pouze o odhady. Abychom získali představu o nejistotě těchto odhadů, můžeme použít funkci intervals z balíčku nlme a získat 95% intervaly spolehlivosti. Pro redukci výstupu uložíme výsledek a zobrazíme vybrané prvky objektu.

int <- intervals(vm2)

Prvek varStruct obsahuje 95% interval spolehlivosti pro parametr odhadnutý ve funkci rozptylu. Parametrem je v tomto případě násobek reziduální standardní chyby pro skupinu mužů.

int$varStruct

## lower est. upper## m 2.245615 2.58127 2.967095## attr(,"label")## "Variance function:"

Prvek sigma obsahuje 95% interval spolehlivosti pro reziduální standardní chybu.

int$sigma

## lower est. upper ## 1.818639 2.005397 2.211335 ## attr(,"label")## "Residual standard error:"

Oba intervaly jsou poměrně malé a obsahují „skutečnou“ hodnotu, kterou jsme použili při generování dat.

varPower

Funkce varPower nám umožňuje modelovat rozptyl jako mocninu absolutní hodnoty kovariáty. Abychom si to pomohli vysvětlit, budeme opět simulovat data s touto vlastností. Níže si všimněte především argumentu sd funkce rnorm. To je místo, kde vezmeme směrodatnou odchylku 2 a pak ji vynásobíme absolutní hodnotou x zvýšenou na mocninu 1,5. To je místo, kde vezmeme směrodatnou odchylku 1. Je to podobné, jako když jsme výše simulovali data pro demonstraci funkce varFixed. V tomto případě jsme jednoduše předpokládali, že exponent je 0,5. (Připomeňme, že odmocnina z čísla je ekvivalentní jeho zvýšení na mocninu 0,5.) Zde jsme libovolně zvolili mocninu 1,5. Když použijeme gls s varPower, pokusíme se získat „skutečnou“ hodnotu 1,5 navíc k hodnotě 2.

set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()

Vidíme, že data vykazují klasický tvar rozptylu, který roste s rostoucím prediktorem. Abychom tato data správně modelovali pomocí gls, dodáme mu vzorec y ~x a použijeme argument weights s varPower. Všimněte si, že zadáváme jednostranný vzorec stejně jako u funkce varFixed. V tomto modelu však získáme odhad mocniny.

vm3 <- gls(y ~ x, data = d, weights = varPower(form = ~x))summary(vm3)

## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 8840.188 8859.811 -4416.094## ## Variance function:## Structure: Power of variance covariate## Formula: ~x ## Parameter estimates:## power ## 1.520915 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.321502 0.4750149 4.887218 0## x 1.582272 0.2241367 7.059404 0## ## Correlation: ## (Intr)## x -0.845## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.892319514 -0.655237190 0.001653178 0.691241690 3.346273816 ## ## Residual standard error: 1.872944 ## Degrees of freedom: 1000 total; 998 residual

int <- intervals(vm3)int$varStruct

## lower est. upper## power 1.445064 1.520915 1.596765## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.650987 1.872944 2.124740 ## attr(,"label")## "Residual standard error:"

Mocnina je odhadnuta na 1,52, což je velmi blízko „skutečné“ hodnotě 1,5. V případě, že je mocnina odhadnuta na 1,5, dostaneme odhad na 1,5. Odhadovaná reziduální směrodatná chyba 1,8729439 je také poměrně blízko hodnotě 2. Oba intervaly vystihují hodnoty, které jsme použili při simulaci dat. Na druhou stranu jsou koeficienty v modelu odhadnuty poněkud špatně. To není překvapivé vzhledem k tomu, jak velká variabilita je v y, zejména pro \(x > 2\).

varExp

Funkce varExp nám umožňuje modelovat rozptyl jako exponenciální funkci kovariáty. Přesto si tuto funkci rozptylu opět vysvětlíme pomocí simulovaných dat. Jediná změna je v argumentu sd funkce rnorm. Máme pevnou hodnotu 2, kterou vynásobíme x, který sám vynásobíme 0,5 a exponenciálně vynásobíme. Všimněte si, jak rychle roste rozptyl, když zvyšujeme x. To je způsobeno exponenciálním růstem rozptylu.

set.seed(55)n <- 400x <- seq(1, 10, length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*exp(0.5*x))d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()

Pro práci zpět a získání těchto hodnot použijeme funkci varExp v argumentu weights gls. Jednostranný vzorec se v tomto případě nemění. Říká, že modelujeme rozptyl jako funkci x. Funkce varExp říká, že x byla vynásobena nějakou hodnotou a exponována, takže gls se pokusí tuto hodnotu odhadnout.

vm4 <- gls(y ~ x, data = d, weights = varExp(form = ~x))summary(vm4)

## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 3873.153 3889.098 -1932.576## ## Variance function:## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 0.4845623 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 1.198525 1.1172841 1.072713 0.284## x 2.073297 0.4933502 4.202486 0.000## ## Correlation: ## (Intr)## x -0.892## ## Standardized residuals:## Min Q1 Med Q3 Max ## -3.16976627 -0.69525696 0.00614222 0.63718022 2.90423217 ## ## Residual standard error: 2.119192 ## Degrees of freedom: 400 total; 398 residual

int <- intervals(vm4)int$varStruct

## lower est. upper## expon 0.4562871 0.4845623 0.5128374## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.786737 2.119192 2.513505 ## attr(,"label")## "Residual standard error:"

Odhad „expon“ 0,4845623 v části „funkce rozptylu“ se velmi blíží námi zadané hodnotě 0,5. To znamená, že funkce varExp se bude snažit odhadnout tuto hodnotu. Stejně tak odhadovaná směrodatná chyba reziduí 2,1191918 je blízká „skutečné“ hodnotě 2. Odhady koeficientů modelu jsou také blízké hodnotám, které jsme použili při generování dat, ale všimněte si nejistoty u (Intercept). Test hypotézy v souhrnu nemůže vyloučit záporný intercept. To opět není překvapivé, protože y má tak velký nekonstantní rozptyl, stejně jako skutečnost, že nemáme žádná pozorování x rovného 0. Protože intercept je hodnota, kterou y nabývá, když se x rovná 0, náš odhadovaný intercept je extrapolací na událost, kterou jsme nepozorovali.

varConstPower

Funkce varConstPower nám umožňuje modelovat rozptyl jako kladnou konstantu plus mocninu absolutní hodnoty kovariáty. To je sice úsměvné, ale je to v podstatě totéž jako funkce varPower, jen k ní nyní přidáváme kladnou konstantu. Následující kód simuluje data, pro která by bylo vhodné použít funkci varConstPower. Všimněte si, že je totožný s kódem, který jsme použili k simulaci dat pro část varPower, až na to, že k x ve funkci abs přidáváme hodnotu 0,7. Proč 0,7? Žádný zvláštní důvod. Je to prostě kladná konstanta, kterou jsme zvolili.

set.seed(4)n <- 1000x <- seq(1,10,length.out = n)y <- 1.2 + 2.1*x + rnorm(n, sd = 2*abs(0.7 + x)^1.5)d <- data.frame(y, x)ggplot(d, aes(x, y)) + geom_point()

Správný způsob modelování těchto dat je použití gls s funkcí varConstPower. Jednostranný vzorec je stejný jako dříve. Souhrn vrátí tři odhady pro model rozptylu: konstantu, mocninu a reziduální směrodatnou chybu. Připomeňme, že „pravé“ hodnoty jsou 0,7, 1,5 a 2.

vm5 <- gls(y ~ x, data = d, weights = varConstPower(form = ~x))summary(vm5)

## Generalized least squares fit by REML## Model: y ~ x ## Data: d ## AIC BIC logLik## 9319.794 9344.323 -4654.897## ## Variance function:## Structure: Constant plus power of variance covariate## Formula: ~x ## Parameter estimates:## const power ## 0.3554921 1.3593016 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) 2.763788 0.7607449 3.633003 3e-04## x 1.456677 0.2915621 4.996112 0e+00## ## Correlation: ## (Intr)## x -0.824## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.846999539 -0.664551194 0.006709051 0.679895874 3.287366683 ## ## Residual standard error: 2.887761 ## Degrees of freedom: 1000 total; 998 residual

int <- intervals(vm5)int$varStruct

## lower est. upper## const 0.02658804 0.3554921 4.753063## power 1.14258677 1.3593016 1.576016## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.814757 2.887761 4.595195 ## attr(,"label")## "Residual standard error:"

Intervaly zachycují „pravé“ hodnoty, i když interval pro konstantu je poměrně široký. Kdybychom neznali skutečnou hodnotu, zdálo by se, že konstanta by mohla být pravděpodobně 2 nebo dokonce 4.

varComb

Nakonec funkce varComb nám umožňuje modelovat kombinace dvou nebo více rozptylových modelů tak, že vynásobíme dohromady odpovídající rozptylové funkce. Je zřejmé, že tímto způsobem lze pojmout některé velmi složité rozptylové modely. Nasimulujeme si několik základních dat, která by bylo vhodné použít s funkcí varComb.

Následující postup spočívá v kombinaci dvou různých rozptylových procesů:

• Jeden, který umožňuje, aby se směrodatná odchylka lišila mezi pohlavími („f“ = \(2\), „m“ = \(2 \krát 2.5\))
• Jeden, který umožňuje, aby se směrodatná odchylka zvyšovala s růstem x, přičemž x se vynásobí 1,5 a exponenciálně vynásobí

Pro usnadnění vizualizace dat jsme omezili rozsah x na 1 až 3. V případě, že je x větší než x, můžeme použít

.

set.seed(77)n <- 400x <- seq(1, 3, length.out = n)g <- sample(c("m","f"), size = n, replace = TRUE)msd <- ifelse(g=="m", 2*2.5, 2) * exp(1.5*x)y <- 1.2 + 2.1*x - 1.5*(g == "m") + 2.8*x*(g == "m") + rnorm(n, sd = msd)d <- data.frame(y, x, g)ggplot(d, aes(x, y, color = g)) + geom_point()

Kresba ukazuje rostoucí rozptyl s rostoucím x, ale také rozdíly v rozptylu mezi pohlavími. Rozptyl y u skupiny „m“ je mnohem větší než rozptyl y u skupiny „f“, zejména když je x větší než 1,5.

Abychom správně modelovali proces generování dat, který jsme uvedli výše, a pokusili se získat skutečné hodnoty, použijeme funkci varComb jako obal kolem dalších dvou funkcí rozptylu: varIdent a varExp. Proč právě tyto dvě? Protože máme různé rozptyly mezi pohlavími a máme rozptyl rostoucí exponenciálně jako funkci x.

vm6 <- gls(y ~ x*g, data = d, weights = varComb(varIdent(form = ~ 1|g), varExp(form = ~x)))summary(vm6)

## Generalized least squares fit by REML## Model: y ~ x * g ## Data: d ## AIC BIC logLik## 4431.45 4459.32 -2208.725## ## Combination of variance functions: ## Structure: Different standard deviations per stratum## Formula: ~1 | g ## Parameter estimates:## f m ## 1.000000 2.437046 ## Structure: Exponential of variance covariate## Formula: ~x ## Parameter estimates:## expon ## 1.540608 ## ## Coefficients:## Value Std.Error t-value p-value## (Intercept) -5.996337 6.539552 -0.9169339 0.3597## x 8.829971 4.849272 1.8208860 0.0694## gm -17.572238 16.548540 -1.0618603 0.2889## x:gm 16.932938 12.192793 1.3887661 0.1657## ## Correlation: ## (Intr) x gm ## x -0.972 ## gm -0.395 0.384 ## x:gm 0.387 -0.398 -0.974## ## Standardized residuals:## Min Q1 Med Q3 Max ## -2.74441479 -0.67478610 -0.00262221 0.66079254 3.14975288 ## ## Residual standard error: 1.794185 ## Degrees of freedom: 400 total; 396 residual

Souhrnná část obsahuje dvě části pro modelování rozptylu. První z nich odhaduje multiplikátor pro skupinu „m“ na hodnotu 2,437, která je velmi blízká skutečné hodnotě 2,5. V druhé části je odhadnuta hodnota multiplikátoru pro skupinu „m“. Exponenciální parametr je odhadnut na 1,54, což je velmi blízko skutečné hodnotě 1,5, kterou jsme použili při generování dat. A konečně standardní chyba reziduí se odhaduje na přibližně 1,79, což se blíží skutečné hodnotě 2. Volání intervals(vm6) ukazuje velmi úzké intervaly spolehlivosti. Předpony A a B v prvku varStruct jsou pouze označení pro dva různé modely rozptylu.

int <- intervals(vm6)int$varStruct

## lower est. upper## A.m 2.119874 2.437046 2.801673## B.expon 1.419366 1.540608 1.661850## attr(,"label")## "Variance function:"

int$sigma

## lower est. upper ## 1.380008 1.794185 2.332669 ## attr(,"label")## "Residual standard error:"

Naneštěstí kvůli velké exponenciální variabilitě jsou odhady koeficientů modelu žalostně špatné.

Jaký to má smysl?“

Proč jsme tedy simulovali všechna tato falešná data a pak se pokoušeli získat „pravé“ hodnoty? K čemu je to dobré? Spravedlivé otázky. Odpověď zní: Pomáhá nám to pochopit, jaké předpoklady přijímáme, když specifikujeme a fitujeme model. Když specifikujeme a fitujeme následující model…

m <- lm(y ~ x1 + x2)

… říkáme, že si myslíme, že y se přibližně rovná váženému součtu x1 a x2 (plus intercepce), přičemž chyby jsou náhodné výběry z normálního rozdělení se střední hodnotou 0 a nějakou pevnou směrodatnou odchylkou. Podobně, pokud zadáme a přizpůsobíme následující model…

m <- gls(y ~ x1 + x2, data = d, weights = varPower(form = ~x1))

… říkáme, že si myslíme, že y se přibližně rovná váženému součtu x1 a x2 (plus intercepce), přičemž chyby jsou náhodné výběry z normálního rozdělení se střední hodnotou 0 a směrodatnou odchylkou, která roste jako násobek x1 zvýšený o určitou mocninu.

Pokud můžeme simulovat data vhodná pro tyto modely, pak lépe rozumíme předpokladům, které přijímáme, když tyto modely používáme. Doufejme, že nyní lépe rozumíte tomu, co můžete udělat pro modelování rozptylu pomocí balíčku nlme.

V případě dotazů nebo vysvětlení týkajících se tohoto článku se obraťte na StatLab knihovny UVA: [email protected]

Podívejte se na celou sbírku článků StatLab knihovny UVA.

Clay Ford
Statistical Research Consultant
University of Virginia Library
Dne 7. dubna 2020