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

A varIdent függvény lehetővé teszi, hogy rétegenként különböző varianciákat modellezzünk. Hogy lássuk, hogyan működik, először szimuláljuk az adatokat ezzel a tulajdonsággal. Az alábbiakban a set.seed(11)-t használjuk arra az esetre, ha valaki ugyanazt a véletlenszerű adatot szeretné szimulálni. Ezután a n értékét 400-ra állítjuk, a szimulálni kívánt megfigyelések számára. A x ugyanúgy generálódik, mint korábban. Ebben a példában egy további prediktort veszünk fel g néven, amelyet a nemnek tekinthetünk. Véletlenszerűen 400 “m” és “f” értéket veszünk mintát. Ezután szimulálunk egy két standard eltérésből álló vektort, és msd néven elmentjük. Ha g == "m", akkor a szórás \(2 \szer 2,5\). Ellenkező esetben g == "f" esetén \(2\). Ezt a vektort használjuk a következő sorban a y létrehozásához. Figyeljük meg, hogy a msd hol van beillesztve a rnorm függvénybe. Figyeljük meg azt is, hogyan generáljuk y. A x és y között interakciót veszünk fel. Ha g == "f", akkor a metszéspont és a meredekség 1,2 és 2,1. Ha g == "m", akkor a metszéspont és a meredekség (1,2 + 2,8), illetve (2,1 + 2,8). Végül a szimulált adatainkat egy adatkeretbe helyezzük, és ggplot2-vel ábrázoljuk.

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

Figyeljük meg a különböző szórásokat az egyes g szinteknél. Az “m” szórása sokkal nagyobb, mint az “f” szórása. Sokkal nagyobb szórással rendelkezik, mint az “f”. Így szimuláltuk az adatokat. Az “m” varianciáját úgy állítottuk be, hogy az “f” varianciájának 2,5-szerese legyen.

Most használjuk a gls függvényt a varIdent segítségével, hogy megpróbáljuk visszanyerni ezeket a valódi értékeket. Ugyanúgy használjuk, mint korábban: definiáljuk a varianciafüggvényünket a weights argumentumban. Az alábbiakban a form argumentumban megadjuk, hogy a variancia modellezésére szolgáló képlet a g függvénytől függjön. A ~ 1|g kifejezés egy egyoldalú formula, amely azt mondja ki, hogy a variancia a g szintjei között különbözik. A 1 csak azt jelenti, hogy a variancia modellünkben nincsenek további kovariánsok. Lehetséges egy olyan képletet is felvenni, mint a ~ x|g, de ez ebben az esetben helytelen lenne, mivel nem használtuk a x-et a variancia előállításához. A grafikonra pillantva az is látható, hogy míg a y variabilitása különbözik a csoportok között, a variabilitás nem nő a x növekedésével.

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

Először is figyeljük meg a kimenet alján, hogy a becsült reziduális standard hiba 2.0053974, ami nagyon közel van a “valódi” 2-es értékhez. Vegyük észre azt is, hogy a “Varianciafüggvény” szakaszban az “m” csoportra 2,58127-es becsült értéket kapunk, ami nagyon közel van a “valódi” 2,5-ös értékhez, amit az “m” csoport eltérő varianciájának létrehozásához használtunk. Általában, amikor a varIdent függvényt használjuk a rétegek szintjei közötti különböző varianciák becslésére, az egyik szint alapértéknek lesz beállítva, a többit pedig a reziduális standard hiba többszöröseként becsüljük. Ebben az esetben, mivel az “f” betűrendben az “m” előtt áll, az “f” alapértékre, azaz 1-re lett beállítva. Az “f” csoport becsült maradék standard hibája \(2,005397 \szor 1\). Az “m” csoport becsült maradék standard hibája \(2,005397 \times 2,58127\).

Fontos megjegyezni, hogy ezek csak becslések. Ahhoz, hogy megérezzük e becslések bizonytalanságát, a nlme csomag intervals függvényével 95%-os konfidenciaintervallumokat kaphatunk. A kimenet csökkentéséhez mentsük el az eredményt, és tekintsük meg az objektum kiválasztott elemeit.

int <- intervals(vm2)

A varStruct elem tartalmazza a varianciafüggvényben becsült paraméter 95%-os konfidenciaintervallumát. A paraméter ebben az esetben a férfi csoport reziduális standard hiba szorzója.

int$varStruct

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

A sigma elem tartalmazza a reziduális standard hiba 95%-os konfidenciaintervallumát.

int$sigma

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

Mindkét intervallum meglehetősen kicsi, és azt a “valódi” értéket tartalmazza, amelyet az adatok előállításához használtunk.

varPower

A varPower függvény lehetővé teszi, hogy a varianciát egy kovariáns abszolút értékének hatványaként modellezzük. Ennek magyarázatához ismét szimuláljuk az adatokat ezzel a tulajdonsággal. Az alábbiakban a rnorm függvény sd argumentuma a legfontosabb, amit észre kell vennünk. Itt vesszük a 2 standard eltérést, majd megszorozzuk a x abszolút értékének 1,5-ös hatványra emelt értékével. Ez hasonló ahhoz, ahogyan a fenti varFixed függvény bemutatásához szimuláltuk az adatokat. Ebben az esetben egyszerűen feltételeztük, hogy az exponens 0,5. (Emlékezzünk vissza, hogy egy szám négyzetgyökének kivétele egyenértékű a 0,5 hatványára való emeléssel). Itt önkényesen választottunk egy 1,5-es hatványt. Ha a gls-et a varPower-val használjuk, akkor a 2 mellett megpróbáljuk visszanyerni az 1,5 “valódi” értékét is.

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

Az adatokon a klasszikus formáját látjuk annak, hogy a variancia a prediktor növekedésével nő. Ahhoz, hogy ezt az adatot helyesen modellezzük a gls segítségével, a y ~x képletet adjuk meg, és a weights argumentumot a varPower-val használjuk. Vegyük észre, hogy az egyoldalú képletet ugyanúgy adjuk meg, mint a varFixed függvény esetében. Ebben a modellben azonban egy becslést kapunk a teljesítményre.

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

A teljesítmény becsült értéke 1,52, ami nagyon közel van a “valódi” 1,5 értékhez. Az 1,8729439 becsült maradék standard hiba szintén elég közel van a 2-hez. Mindkét intervallum azokat az értékeket ragadja meg, amelyeket az adatok szimulálásához használtunk. A modellben szereplő együtthatók viszont kissé rosszul becsültek. Ez nem meglepő, tekintve, hogy mennyi variabilitás van a y-ben, különösen \(x > 2\) esetében.

varExp

A varExp függvény lehetővé teszi, hogy a varianciát egy kovariáns exponenciális függvényeként modellezzük. Még egyszer elmagyarázzuk ezt a varianciafüggvényt szimulált adatok segítségével. Az egyetlen változás a rnorm függvény sd argumentumában van. Van egy fix 2-es értékünk, amelyet megszorozunk a x értékkel, amelyet maga is megszoroz 0,5-tel és exponenciál. Figyeljük meg, milyen gyorsan nő a variancia, ahogy növeljük a x értéket. Ez a variancia exponenciális növekedésének köszönhető.

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

Azért, hogy visszafelé haladva visszanyerjük ezeket az értékeket, a gls weights argumentumában a varExp függvényt használjuk. Az egyoldalú képlet ebben az esetben nem változik. Azt mondja ki, hogy a varianciát a x függvényeként modellezzük. A varExp függvény azt mondja, hogy x-et megszoroztuk valamilyen értékkel és exponenciáltuk, így a gls megpróbálja megbecsülni ezt az értéket.

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

A “expon” becslés 0,4845623 a “varianciafüggvény” részben nagyon közel van a megadott 0,5 értékünkhöz. Hasonlóképpen a 2,1191918 becsült maradék standard hiba közel van a “valódi” 2 értékhez. A modell együttható becslései szintén közel vannak az adatok generálásához használt értékeinkhez, de vegyük észre a (Intercept) bizonytalanságát. Az összegzésben szereplő hipotézisteszt nem tudja kizárni a negatív interceptet. Ez ismét nem meglepő, mivel y-nek olyan nagy a nem konstans varianciája, valamint az a tény, hogy nincs olyan megfigyelésünk, amikor x 0-val egyenlő. Mivel a metszéspont az az érték, amelyet y akkor vesz fel, amikor x 0-val egyenlő, a becsült metszéspontunk egy olyan eseményre való extrapoláció, amelyet nem figyeltünk meg.

varConstPower

A varConstPower függvény lehetővé teszi számunkra, hogy a varianciát pozitív konstansként és egy kovariáns abszolút értékének hatványaként modellezzük. Ez elég nagy szó, de ez lényegében ugyanaz, mint a varPower függvény, csak most egy pozitív állandót adunk hozzá. Az alábbi kód olyan adatokat szimulál, amelyekhez a varConstPower függvényt alkalmas lenne használni. Vegyük észre, hogy megegyezik azzal a kóddal, amelyet a varPower szakasz adatainak szimulálásához használtunk, kivéve, hogy a x-hez 0,7-et adunk hozzá a abs függvényben. Miért 0,7? Nincs különösebb ok. Ez csak egy pozitív konstans, amit kiválasztottunk.

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

Az adatok modellezésének helyes módja a gls használata a varConstPower függvénnyel. Az egyoldalú képlet ugyanaz, mint korábban. Az összegzés három becslést ad vissza a varianciamodellre: a konstans, a teljesítmény és a reziduális standard hiba. Emlékezzünk vissza, hogy az “igaz” értékek 0,7, 1,5 és 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:"

Az intervallumok az “igaz” értékeket rögzítik, bár a konstans intervalluma meglehetősen széles. Ha nem ismernénk az igazi értéket, úgy tűnik, hogy a konstans hihetően lehet 2 vagy akár 4.

varComb

Végül a varComb függvény lehetővé teszi, hogy két vagy több varianciamodell kombinációit modellezzük a megfelelő varianciafüggvények összeszorzásával. Nyilvánvalóan ez nagyon összetett varianciamodelleket is képes befogadni. Néhány olyan alapadatot fogunk szimulálni, amelyek alkalmasak lennének a varComb függvény használatára.

Az alábbiakban két különböző varianciafolyamatot kombinálunk:

• Az egyik lehetővé teszi, hogy a szórás eltérjen a nemek között (“f” = \(2\), “m” = \(2 \times 2.5\))
• Egyik, amelyik lehetővé teszi, hogy a szórás növekedjen a x növekedésével, ahol x megszorozzuk 1,5-tel és exponenciáljuk

Az adatok vizualizálásának megkönnyítése érdekében a x-t 1 és 3 közötti tartományra korlátoztuk.

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

A grafikon a x növekedésével növekvő varianciát, de a nemek közötti eltéréseket is mutatja. Az y varianciája az “m” csoportban sokkal nagyobb, mint az y varianciája az “f” csoportban, különösen, ha x nagyobb, mint 1,5.

Hogy helyesen modellezzük a fent megadott adatgenerálási folyamatot, és megpróbáljuk visszanyerni a valódi értékeket, a varComb függvényt két további varianciafüggvény köré burkolva használjuk: varIdent és varExp. Miért pont ez a kettő? Mert a nemek között különböző varianciák vannak, és a variancia exponenciálisan növekszik a x függvényében.

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

Az összefoglaló szakasz két szakaszt tartalmaz a variancia modellezésére. Az első az “m” csoportra vonatkozó szorzót 2,437-re becsüli, ami nagyon közel van a 2,5 valódi értékhez. Az exponenciális paramétert 1,54-re becsüljük, ami rendkívül közel van az adatok generálásakor használt 1,5-ös valódi értékhez. Végül a maradék standard hiba becsült értéke körülbelül 1,79, ami közel van a 2 valós értékhez. A intervals(vm6) hívás nagyon szűk konfidenciaintervallumokat mutat. A A és B előtagok a varStruct elemben csak a két különböző varianciamodell címkéi.

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

A nagy exponenciális változékonyság miatt sajnos a modell együtthatóinak becslése siralmasan rossz.

Mire jó ez?

Miért szimuláltuk tehát ezt a sok hamis adatot, majd próbáltuk visszanyerni a “valódi” értékeket? Mire jó ez? Tisztességes kérdések. A válasz az, hogy segít megérteni, milyen feltételezéseket teszünk, amikor egy modellt specifikálunk és illesztünk. Amikor megadjuk és illesztjük a következő modellt…

m <- lm(y ~ x1 + x2)

…azt mondjuk, hogy szerintünk y megközelítőleg egyenlő x1 és x2 súlyozott összegével (plusz egy metszéspont), a hibák pedig véletlenszerű húzások egy 0 átlagú és fix szórású normál eloszlásból. Hasonlóképpen, ha megadjuk és illesztjük a következő modellt…

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

…azt mondjuk, hogy szerintünk y megközelítőleg egyenlő x1 és x2 súlyozott összegével (plusz egy intercept), ahol a hibák véletlenszerű húzások egy 0 átlagú normáleloszlásból és egy olyan szórással, amely a x1 bizonyos hatványértékkel növelt többszöröseként növekszik.

Ha ezeknek a modelleknek megfelelő adatokat tudunk szimulálni, akkor jobban megértjük, hogy milyen feltételezésekkel élünk, amikor ezeket a modelleket használjuk. Remélhetőleg most már jobban érti, hogy mit tehet a variancia modellezésére a nlme csomag segítségével.

A cikkel kapcsolatos kérdések vagy pontosítások esetén forduljon az UVA Library StatLab-hoz: [email protected]

Nézze meg az UVA Library StatLab teljes cikkgyűjteményét.

Clay Ford
Statisztikai kutatási tanácsadó
University of Virginia Library
Aprilis 7, 2020