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

Funktio varIdent mahdollistaa erilaisten varianssien mallintamisen ositteittain. Nähdäksemme, miten se toimii, simuloimme ensin dataa, jolla on tämä ominaisuus. Alla käytämme funktiota set.seed(11) siltä varalta, että joku haluaa simuloida samaa satunnaisdataa. Sitten asetamme n arvoksi 400, joka on simuloitavien havaintojen määrä. x generoidaan samalla tavalla kuin aiemmin. Tässä esimerkissä otamme mukaan ylimääräisen ennustajan nimeltä g, jota voidaan ajatella sukupuoleksi. Otamme satunnaisesti 400 arvoa ”m” ja ”f”. Seuraavaksi simuloimme kahden keskihajonnan vektorin ja tallennamme sen nimellä msd. Jos g == "m", niin keskihajonta on \(2 \ kertaa 2,5\). Muussa tapauksessa se on \(2\) g == "f":lle. Käytämme tätä vektoria seuraavalla rivillä luodaksemme y. Huomaa, missä msd on liitetty rnorm-funktioon. Huomaa myös, miten generoimme y. Sisällytämme vuorovaikutuksen x ja y välille. Kun g == "f", leikkauspiste ja kaltevuus ovat 1,2 ja 2,1. Kun g == "m", leikkauspiste ja kaltevuus ovat (1,2 + 2,8) ja (2,1 + 2,8). Lopuksi sijoitamme simuloidut datamme datakehykseen ja piirretään ggplot2:lla.

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

Huomaa erilaiset varianssit kullakin tasolla g. Varianssi ”m”:lle on paljon suurempi kuin varianssi ”f:lle”. Siinä on paljon enemmän hajontaa kuin ”f”:ssä. Simuloimme datan näin. Asetimme ”m:n” varianssin 2,5 kertaa suuremmaksi kuin ”f:n”.

Käytetään nyt gls-funktiota varIdent:n kanssa yrittäen palauttaa nämä todelliset arvot. Käytämme samaa tapaa kuin aiemmin: määrittelemme varianssifunktiomme weights-argumentissa. Alla määrittelemme form-argumentissa, että varianssin mallintamisen kaava on ehdollinen g:lle. Lauseke ~ 1|g on yksipuolinen kaava, joka sanoo, että varianssi eroaa g:n tasojen välillä. 1 tarkoittaa vain sitä, että meillä ei ole ylimääräisiä kovariaatteja varianssin mallissa. On mahdollista sisällyttää kaava kuten ~ x|g, mutta se olisi tässä tapauksessa virheellinen, koska emme käyttäneet x:ää varianssin muodostamisessa. Kuvaajan tarkastelu osoittaa myös, että vaikka y:n vaihteluväli on erilainen ryhmien välillä, vaihteluväli ei kasva x:n kasvaessa.

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

Huomaa ensin tulosteen alareunasta, että estimoitu jäännösstandardivirhe on 2.0053974, joka on hyvin lähellä ”todellista” arvoa 2. Huomaa myös kohdassa ”Varianssifunktio”, että saamme ryhmälle ”m” estimoiduksi arvoksi 2,58127, joka on hyvin lähellä ”todellista” arvoa 2,5, jota käytimme luodaksemme ryhmän ”m” erilaisen varianssin. Yleensä, kun käytät varIdent-funktiota estimoidaksesi erilaisia variansseja ositteiden tasojen välillä, yksi tasoista asetetaan perustasoksi ja muut estimoidaan jäännösstandardivirheen kertoimina. Tässä tapauksessa, koska ”f” on aakkosjärjestyksessä ennen ”m”:ää, ”f” asetettiin perustasoksi eli 1:ksi. Ryhmän ”f” estimoitu jäännösstandardivirhe on \(2,005397 \ kertaa 1\). Ryhmän ”m” estimoitu jäännösstandardivirhe on \(2.005397 \times 2.58127\).

On tärkeää huomata, että nämä ovat vain arvioita. Saadaksemme käsityksen näiden estimaattien epävarmuudesta, voimme käyttää nlme-paketin intervals-funktiota saadaksemme 95 prosentin luottamusvälit. Tuloksen pienentämiseksi tallennamme tuloksen ja tarkastelemme kohteen valittuja elementtejä.

int <- intervals(vm2)

Elementti varStruct sisältää varianssifunktiossa estimoidun parametrin 95 %:n luottamusvälin. Parametri on tässä tapauksessa miesryhmän jäännösstandardivirheen kerroin.

int$varStruct

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

Elementti sigma sisältää jäännösstandardivirheen 95 prosentin luottamusvälin.

int$sigma

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

Kumpikin intervalli on melko pieni ja sisältää sen ”todellisen” arvon, jota käytimme aineiston tuottamiseen.

varPower

Funktion varPower avulla voimme mallintaa varianssin kovariaatin absoluuttisen arvon potenssina. Jälleen kerran tämän selittämisen helpottamiseksi simuloimme dataa, jolla on tämä ominaisuus. Alla tärkein huomioitava asia on rnorm-funktion sd-argumentti. Siinä otamme keskihajonnan 2 ja sitten kerromme sen absoluuttisella arvolla x korotettuna potenssiin 1,5. Tämä on samanlaista kuin se, miten simuloimme dataa edellä esitetyn varFixed-funktion demonstroimiseksi. Siinä tapauksessa oletimme yksinkertaisesti, että eksponentti oli 0,5. (Muistathan, että luvun neliöjuuren ottaminen vastaa sen korottamista potenssiin 0,5). Tässä valitsimme mielivaltaisesti potenssin 1,5. Kun käytämme gls yhdessä varPower:n kanssa, yritämme saada takaisin ”todellisen” arvon 1,5 lisäksi 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()

Näemme, että datassa näkyy klassinen muoto, jossa varianssi kasvaa ennustajan kasvaessa. Jotta voimme mallintaa tämän datan oikein käyttämällä gls, annamme sille kaavan y ~x ja käytämme weights-argumenttia varPower kanssa. Huomaa, että määrittelemme yksipuolisen kaavan aivan kuten teimme varFixed-funktion kanssa. Tässä mallissa saamme kuitenkin estimaatin potenssille.

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

Potenssiksi estimoidaan 1,52, joka on hyvin lähellä ”todellista” arvoa 1,5. Arvioitu jäännösstandardivirhe 1,8729439 on myös melko lähellä arvoa 2. Molemmat väliarvot kuvaavat niitä arvoja, joita käytimme datan simuloinnissa. Mallin kertoimet ovat sen sijaan hieman heikosti estimoituja. Tämä ei ole yllättävää, kun otetaan huomioon, kuinka paljon vaihtelua y:ssä on, erityisesti \(x > 2\).

varExp

Funktion varExp avulla voimme mallintaa varianssia kovariaatin eksponentiaalisena funktiona. Selitämme jälleen kerran tätä varianssifunktiota simuloitujen tietojen avulla. Ainoa muutos on rnorm-funktion sd-argumentissa. Meillä on kiinteä arvo 2, jonka kerromme x:llä, joka puolestaan kerrotaan 0,5:llä ja eksponentoidaan. Huomaa, kuinka nopeasti varianssi kasvaa, kun kasvatamme x. Tämä johtuu varianssin eksponentiaalisesta kasvusta.

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

Työstääksemme taaksepäin ja saadaksemme nämä arvot takaisin käytämme varExp-funktiota gls:n weights-argumentissa. Yksipuolinen kaava ei muutu tässä tapauksessa. Sen mukaan varianssi mallinnetaan x:n funktiona. Funktio varExp kertoo, että x on kerrottu jollakin arvolla ja eksponenttiarvoistettu, joten gls yrittää arvioida tätä arvoa.

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

Varanssifunktio-kohdan ”expon”-estimaatti 0,4845623 on hyvin lähellä määrittelemäämme arvoa 0,5. Samoin estimoitu jäännösstandardivirhe 2,1191918 on lähellä ”todellista” arvoa 2. Mallin kerroinestimaatit ovat myös lähellä niitä arvoja, joita käytimme tietojen tuottamiseen, mutta huomaa epävarmuus (Intercept). Yhteenvedossa oleva hypoteesitesti ei voi sulkea pois negatiivista interceptiä. Tämäkään ei ole yllättävää, koska y:llä on niin paljon ei-vakioivaa varianssia ja koska meillä ei ole yhtään havaintoa siitä, että x olisi 0. Koska leikkauspiste on arvo, jonka y saa, kun x on 0, estimoimamme leikkauspiste on ekstrapolointi tapahtumaan, jota emme ole havainneet.

varConstPower

Funktion varConstPower avulla voimme mallintaa varianssin positiivisena vakiona plus kovariaatin absoluuttisen arvon potenssina. Se on suusanallinen, mutta tämä on periaatteessa sama kuin varPower-funktio, paitsi että nyt lisäämme siihen positiivisen vakion. Seuraavassa koodissa simuloidaan dataa, johon varConstPower-funktio sopisi käytettäväksi. Huomaa, että se on identtinen koodin kanssa, jota käytimme simuloidaksemme dataa varPower-jaksoa varten, paitsi että lisäämme 0,7 x:ään abs-funktiossa. Miksi 0,7? Ei erityistä syytä. Se on vain positiivinen vakio, jonka valitsimme.

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

Oikein tapa mallintaa tätä dataa on käyttää gls yhdessä varConstPower-funktion kanssa. Yksipuolinen kaava on sama kuin aiemmin. Yhteenveto palauttaa varianssimallille kolme estimaattia: vakio, potenssi ja residuaalin keskivirhe. Muistutetaan, että ”todelliset” arvot ovat 0.7, 1.5 ja 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:"

Välit kaappaavat ”todelliset” arvot, vaikkakin vakion intervalli on melko laaja. Jos emme tietäisi todellista arvoa, näyttäisi siltä, että vakio voisi uskottavasti olla 2 tai jopa 4.

varComb

Funktion varComb avulla voimme lopuksi mallintaa kahden tai useamman varianssimallin kombinaatioita kertomalla vastaavat varianssifunktiot yhteen. On selvää, että näin voidaan toteuttaa joitakin hyvin monimutkaisia varianssimalleja. Simuloimme joitakin perusdatoja, joita olisi tarkoituksenmukaista käyttää varComb-funktion kanssa.

Alhaalla yhdistämme kaksi erilaista varianssiprosessia:

• Yksi, joka sallii keskihajonnan vaihdella sukupuolten välillä (”f” = \(2\), ”m” = \(2 \times 2.5\))
• Yksi, joka sallii keskihajonnan kasvaa x:n kasvaessa, jolloin x kerrotaan 1.5:llä ja eksponentoidaan

Datan visualisoinnin helpottamiseksi olemme rajoittaneet x:n vaihteluväliksi 1:stä 3:een.

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

Kuvio osoittaa, että varianssi kasvaa x:n kasvaessa, mutta myös sukupuolten väliset erot varianssissa. Ryhmän ”m” varianssi y on paljon suurempi kuin ryhmän ”f” varianssi y, erityisesti kun x on suurempi kuin 1,5.

Mutta jotta voimme mallintaa oikein edellä määrittelemämme tiedonmuodostusprosessin ja pyrkiä saamaan takaisin todelliset arvot, käytämme varComb-funktiota kääreenä kahden muun varianssifunktion ympärillä: varIdent ja varExp. Miksi nämä kaksi? Koska meillä on erilaisia variansseja sukupuolten välillä, ja meillä varianssi kasvaa eksponentiaalisesti funktiona 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

Yhteenveto-osio sisältää kaksi osiota varianssin mallintamista varten. Ensimmäisessä arvioidaan ”m”-ryhmän kertoimeksi 2,437, joka on hyvin lähellä todellista arvoa 2,5. Eksponentiaaliparametrin arvoksi arvioidaan 1,54, mikä on erittäin lähellä todellista arvoa 1,5, jota käytimme aineistoa luodessamme. Lopuksi jäännösstandardivirhe arvioidaan noin 1,79:ksi, joka on lähellä todellista arvoa 2. Kutsu intervals(vm6) osoittaa erittäin tiukat luottamusvälit. Elementin varStruct etuliitteet A ja B ovat vain merkintöjä kahdelle eri varianssimallille.

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

Suuren eksponentiaalisen vaihtelun vuoksi mallin kertoimien estimaatit ovat valitettavan huonoja.

Mitä järkeä siinä on?

Miksi siis simuloimme kaiken tämän väärennetyn datan ja yritimme sitten palauttaa ”todelliset” arvot? Mitä hyötyä siitä on? Reiluja kysymyksiä. Vastaus on, että se auttaa meitä ymmärtämään, mitä oletuksia teemme, kun määrittelemme ja sovitamme mallia. Kun määrittelemme ja sovitamme seuraavan mallin…

m <- lm(y ~ x1 + x2)

… sanomme, että uskomme y olevan suunnilleen yhtä suuri kuin x1:n ja x2:n painotettu summa (sekä leikkauspiste) virheiden ollessa satunnaisotoksia normaalijakaumasta, jonka keskiarvo on 0 ja kiinteä keskihajonta. Vastaavasti, jos määrittelemme ja sovitamme seuraavan mallin…

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

…sanomme, että mielestämme y on suunnilleen yhtä suuri kuin x1:n ja x2:n (sekä leikkauspisteen) painotettu summa, jossa virheet ovat satunnaisotoksia normaalijakaumasta, jonka keskiarvo on 0 ja keskihajonta, joka kasvaa moninkertaisena jollakin potenssilla korotetulle arvoasteikolla.

Jos pystymme simuloimaan dataa, joka soveltuu näihin malleihin, ymmärrämme paremmin, mitä oletuksia teemme, kun käytämme näitä malleja. Toivottavasti sinulla on nyt parempi käsitys siitä, mitä voit tehdä varianssin mallintamiseksi nlme-paketilla.

Jos sinulla on tätä artikkelia koskevia kysymyksiä tai selvennyksiä, ota yhteyttä UVA:n kirjaston StatLabiin: [email protected]

Katsele koko UVA:n kirjaston StatLab-artikkelikokoelmaa.

Clay Ford
Statistinen tutkimuskonsultti
Virginiankirjaston yliopiston kirjasto
Huhtikuun 7. päivä, 2020