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

Funkcja varIdent pozwala nam modelować różne wariancje na warstwę. Aby zobaczyć jak to działa, najpierw zasymulujemy dane z tą właściwością. Poniżej używamy set.seed(11) na wypadek, gdyby ktoś chciał zasymulować te same losowe dane. Następnie ustawiamy n równe 400, czyli liczbę obserwacji, które będziemy symulować. x jest generowany tak samo jak poprzednio. W tym przykładzie dołączamy dodatkowy predyktor o nazwie g, który można traktować jako płeć. Próbujemy losowo 400 wartości „m” i „f”. Następnie symulujemy wektor dwóch odchyleń standardowych i zapisujemy jako msd. Jeśli g == "m", to odchylenie standardowe jest \(2 razy 2,5). W przeciwnym razie jest to \(2\) dla g == "f". Używamy tego wektora w następnej linii, aby wygenerować y. Zauważ, gdzie msd jest wstawione do funkcji rnorm. Zauważ też, jak generujemy y. Włączamy interakcję pomiędzy x i y. Gdy g == "f", intercept i nachylenie wynoszą 1,2 i 2,1. Kiedy g == "m", intercept i nachylenie wynoszą odpowiednio (1.2 + 2.8) i (2.1 + 2.8). Na koniec umieszczamy nasze symulowane dane w ramce danych i tworzymy wykres za 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()

Zauważ różne wariancje dla każdego poziomu g. Wariancja dla „m” jest znacznie większa niż wariancja dla „f”. Ma o wiele większy rozrzut niż „f”. W ten sposób symulowaliśmy dane. Ustawiliśmy wariancję dla „m” na 2,5 razy większą niż wariancja dla „f”.

Teraz użyjmy funkcji gls z varIdent w próbie odzyskania tych prawdziwych wartości. Używamy tego samego sposobu co poprzednio: definiujemy naszą funkcję wariancji w argumencie weights. Poniżej określamy w argumencie form, że wzór na modelowanie wariancji jest warunkowy na g. Wyrażenie ~ 1|g jest formułą jednostronną, która mówi, że wariancja różni się pomiędzy poziomami g. Wyrażenie 1 oznacza jedynie, że w naszym modelu wariancji nie mamy żadnych dodatkowych zmiennych kowariancyjnych. Możliwe jest włączenie formuły takiej jak ~ x|g, ale w tym przypadku byłoby to niepoprawne, ponieważ nie użyliśmy x w generowaniu wariancji. Spojrzenie na wykres pokazuje również, że podczas gdy zmienność w y jest różna między grupami, zmienność nie rośnie wraz ze wzrostem x.

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

Pierwsze spostrzeżenie na dole danych wyjściowych, że szacowany resztowy błąd standardowy wynosi 2.0053974, bardzo blisko „prawdziwej” wartości 2. Zauważ również, że w sekcji „Funkcja wariancji” otrzymujemy szacunkową wartość 2.58127 dla grupy „m”, która jest bardzo blisko „prawdziwej” wartości 2.5, której użyliśmy do wygenerowania różnej wariancji dla grupy „m”. Ogólnie rzecz biorąc, kiedy używasz funkcji varIdent do oszacowania różnych wariancji pomiędzy poziomami warstw, jeden z poziomów zostanie ustawiony jako bazowy, a pozostałe zostaną oszacowane jako wielokrotności resztowego błędu standardowego. W tym przypadku, ponieważ „f” występuje przed „m” alfabetycznie, „f” zostało ustawione na wartość podstawową, czyli 1. Oszacowany resztowy błąd standardowy dla grupy „f” wynosi \(2,005397 \ razy 1\). Szacowany resztkowy błąd standardowy dla grupy „m” to \(2.005397 \ razy 2.58127 \).

Ważne jest, aby zauważyć, że są to tylko szacunki. Aby poznać niepewność tych oszacowań, możemy użyć funkcji intervals z pakietu nlme do uzyskania 95% przedziałów ufności. Aby zmniejszyć rozmiar danych wyjściowych, zapisujemy wynik i oglądamy wybrane elementy obiektu.

int <- intervals(vm2)

Element varStruct zawiera 95% przedział ufności dla parametru oszacowanego w funkcji wariancji. Parametrem w tym przypadku jest mnożnik resztowego błędu standardowego dla grupy mężczyzn.

int$varStruct

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

Element sigma zawiera 95% przedział ufności dla resztowego błędu standardowego.

int$sigma

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

Oba przedziały są dość małe i zawierają „prawdziwą” wartość, której użyliśmy do wygenerowania danych.

varPower

Funkcja varPower pozwala nam modelować wariancję jako potęgę wartości bezwzględnej zmiennej. Jeszcze raz, aby to wyjaśnić, zasymulujemy dane z tą właściwością. Poniżej główną rzeczą, na którą należy zwrócić uwagę jest argument sd funkcji rnorm. To właśnie tam bierzemy odchylenie standardowe 2, a następnie mnożymy je przez wartość bezwzględną x podniesioną do potęgi 1.5. Jest to podobne do tego, jak symulowaliśmy dane, aby zademonstrować funkcję varFixed powyżej. W tym przypadku po prostu założyliśmy, że wykładnik wynosi 0,5. (Przypomnijmy, że pierwiastek kwadratowy z liczby jest równoważny podniesieniu jej do potęgi 0.5). Tutaj arbitralnie wybraliśmy potęgę 1.5. Gdy użyjemy gls z varPower będziemy próbowali odzyskać „prawdziwą” wartość 1,5 oprócz 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()

Widzimy, że dane wykazują klasyczną formę wariancji rosnącej wraz ze wzrostem predyktora. Aby poprawnie zamodelować te dane za pomocą gls, zaopatrujemy go w formułę y ~x i używamy argumentu weights z varPower. Zauważ, że określamy formułę jednostronną tak samo, jak w przypadku funkcji varFixed. W tym modelu otrzymamy jednak oszacowanie mocy.

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:"

Moc szacowana jest na 1,52, co jest bardzo bliskie „prawdziwej” wartości 1,5. Oszacowany resztowy błąd standardowy wynoszący 1,8729439 jest również dość bliski wartości 2. Oba przedziały odzwierciedlają wartości, których użyliśmy do symulacji danych. Współczynniki w modelu, z drugiej strony, są nieco słabo oszacowane. Nie jest to zaskakujące, biorąc pod uwagę, jak duża jest zmienność w y, szczególnie dla \(x > 2\).

varExp

Funkcja varExp pozwala nam modelować wariancję jako funkcję wykładniczą zmiennej. Jeszcze raz wyjaśnimy tę funkcję wariancji używając symulowanych danych. Jedyna zmiana dotyczy argumentu sd funkcji rnorm. Mamy stałą wartość 2, którą mnożymy przez x, który sam jest mnożony przez 0.5 i wykładany. Zauważ, jak szybko wariancja wzrasta, gdy zwiększamy x. Wynika to z wykładniczego wzrostu wariancji.

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()

Aby pracować wstecz i odzyskać te wartości, używamy funkcji varExp w argumencie weights funkcji gls. Wzór jednostronny nie ulega w tym przypadku zmianie. Mówi ona o modelowaniu wariancji jako funkcji x. Funkcja varExp mówi, że x została pomnożona przez pewną wartość i wykładnicza, więc gls spróbuje oszacować tę wartość.

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:"

Oszacowanie „expon” wynoszące 0,4845623 w sekcji „funkcja wariancji” jest bardzo bliskie określonej przez nas wartości 0,5. Podobnie szacowany resztkowy błąd standardowy wynoszący 2.1191918 jest bliski „prawdziwej” wartości 2. Oszacowania współczynników modelu są również bliskie wartościom, których użyliśmy do wygenerowania danych, ale zauważ niepewność w (Intercept). Test hipotezy w podsumowaniu nie może wykluczyć ujemnego punktu przecięcia. To znowu nie jest zaskakujące, ponieważ y ma tak dużo niestałej wariancji, jak również fakt, że nie mamy obserwacji x równej 0. Ponieważ punkt przecięcia jest wartością, którą y przyjmuje, gdy x równa się 0, nasz szacowany punkt przecięcia jest ekstrapolacją do zdarzenia, którego nie obserwowaliśmy.

varConstPower

Funkcja varConstPower pozwala nam modelować wariancję jako dodatnią stałą plus potęgę wartości bezwzględnej zmiennej. To trochę skomplikowane, ale jest to w zasadzie to samo, co funkcja varPower, z tą różnicą, że teraz dodajemy do niej dodatnią stałą. Poniższy kod symuluje dane, dla których funkcja varConstPower byłaby odpowiednia do użycia. Zauważ, że jest on identyczny z kodem, którego użyliśmy do symulacji danych dla sekcji varPower, z tą różnicą, że do x w funkcji abs dodajemy 0,7. Dlaczego 0.7? Nie ma specjalnego powodu. To po prostu dodatnia stała, którą wybraliśmy.

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()

Prawidłowym sposobem modelowania tych danych jest użycie gls z funkcją varConstPower. Wzór jednostronny jest taki sam jak poprzednio. Podsumowanie zwraca trzy oszacowania dla modelu wariancji: stałą, moc i resztowy błąd standardowy. Przypomnijmy, że „prawdziwe” wartości wynoszą odpowiednio 0,7, 1,5 i 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:"

Przedziały wychwytują „prawdziwe” wartości, choć przedział dla stałej jest dość szeroki. Gdybyśmy nie znali prawdziwej wartości, wydawałoby się, że stała mogłaby wynosić 2 lub nawet 4.

varComb

Wreszcie, funkcja varComb pozwala nam modelować kombinacje dwóch lub więcej modeli wariancji poprzez mnożenie razem odpowiednich funkcji wariancji. Oczywiście może to obejmować niektóre bardzo złożone modele wariancji. Zasymulujemy kilka podstawowych danych, które byłyby odpowiednie do użycia funkcji varComb.

To, co robimy poniżej, to połączenie dwóch różnych procesów wariancji:

• Jeden, który pozwala na to, aby odchylenie standardowe różniło się w zależności od płci („f” = \(2\), „m” = \(2 \ razy 2.5))
• Jeden, który pozwala na wzrost odchylenia standardowego wraz ze wzrostem x, gdzie x jest mnożone przez 1,5 i wykładane

Aby pomóc w wizualizacji danych ograniczyliśmy x do zakresu od 1 do 3.

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()

Wykres pokazuje rosnącą wariancję wraz ze wzrostem x, ale także różnice w wariancji między płciami. Wariancja y dla grupy „m” jest znacznie większa niż wariancja y w grupie „f”, zwłaszcza gdy x jest większa niż 1.5.

Aby poprawnie zamodelować proces generowania danych, który określiliśmy powyżej i spróbować odzyskać prawdziwe wartości, używamy funkcji varComb jako wrappera wokół dwóch kolejnych funkcji wariancji: varIdent i varExp. Dlaczego akurat te dwie? Ponieważ mamy różne wariancje między płciami i mamy wariancję rosnącą wykładniczo jako funkcja 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

Sekcja podsumowująca zawiera dwie sekcje do modelowania wariancji. Pierwsza z nich szacuje mnożnik dla grupy „m” na 2,437, co jest bardzo bliskie prawdziwej wartości 2,5. Parametr wykładniczy szacowany jest na 1,54, bardzo blisko prawdziwej wartości 1,5, której użyliśmy podczas generowania danych. Wreszcie, resztkowy błąd standardowy jest szacowany na około 1,79, blisko prawdziwej wartości 2. Wywołanie intervals(vm6) pokazuje bardzo wąskie przedziały ufności. Przedrostki A i B w elemencie varStruct to tylko etykiety dla dwóch różnych modeli wariancji.

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:"

Niestety z powodu dużej zmienności wykładniczej oszacowania współczynników modelu są żałośnie złe.

Co to ma na celu?

Po co więc symulowaliśmy wszystkie te fałszywe dane, a następnie próbowaliśmy odzyskać „prawdziwe” wartości? Jaki to ma sens? Słuszne pytania. Odpowiedź brzmi: pomaga nam zrozumieć, jakie założenia przyjmujemy, gdy określamy i dopasowujemy model. Kiedy określamy i dopasowujemy następujący model…

m <- lm(y ~ x1 + x2)

…mówimy, że uważamy, że y jest w przybliżeniu równe ważonej sumie x1 i x2 (plus punkt przecięcia) z błędami będącymi losowymi wyciągnięciami z rozkładu normalnego o średniej 0 i pewnym ustalonym odchyleniu standardowym. Podobnie, jeśli określimy i dopasujemy następujący model…

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

…mówimy, że uważamy, że y jest w przybliżeniu równe ważonej sumie x1 i x2 (plus punkt przecięcia) z błędami będącymi losowymi wyciągnięciami z rozkładu normalnego ze średnią 0 i odchyleniem standardowym, które rośnie jako wielokrotność x1 podniesiona o pewną potęgę.

Jeśli możemy symulować dane odpowiednie dla tych modeli, wtedy mamy lepsze zrozumienie założeń, które robimy, kiedy używamy tych modeli. Mam nadzieję, że teraz lepiej rozumiesz, co możesz zrobić, aby modelować wariancję za pomocą pakietu nlme.

W przypadku pytań lub wyjaśnień dotyczących tego artykułu, skontaktuj się z UVA Library StatLab: [email protected]

Zobacz całą kolekcję artykułów UVA Library StatLab.

Clay Ford
Konsultant ds. badań statystycznych
University of Virginia Library
7 kwietnia 2020