University of Virginia Library Research Data Services + Sciences

Posted by admin Articles

線形モデリングの基本前提の1つは、一定、または均質な分散であることです。 これは具体的にはどういうことでしょうか。

以下では、xという1から10までの数値の並べ替えベクトルを作成し、xの関数であるyという数値のベクトルを作成します。 x と y をプロットすると、切片が 1.2、傾きが 2.1 の直線が得られます。

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

ここで、y が完全に x で決まらないようにデータにいくつかの「ノイズ」を追加してみましょう。 これは平均0、分散を設定した理論的な正規分布からランダムに値を引き、それをyを生成する式に追加することで実現できます。 Rのrnorm関数を使えば、簡単にこれができます。 以下では、平均0、標準偏差2の正規分布からランダムに100個の値を引き、noiseというベクトルとして保存しています。 (標準偏差は単に分散の平方根であることを思い出してください）そして、ノイズを加えたyを生成します。 関数set.seed(222)を使えば、同じように「ランダムな」データを得ることができますので、参考にしてください。 

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

ここで、線形決定性成分（ \(y = 1.2 + 2.1x) ）とランダムノイズ（ \(N(0, 2)\ ）分布から引いたデータを再投影しています。 これは、従来の線形モデルを適用する際のデータに関する基本的な仮定です。 以下では、lm 関数を使用して、データを生成するために使用した「真の」値を「回復」します。1

  • The standard deviation: 2
    • 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

    係数セクションの(Intercept)とx（つまり傾斜）の推定値は、それぞれ1.27と2.09とかなり1.2と2.1に近くなっています。 残差標準誤差1.926も定数2の値にかなり近いです。yの生成方法を知り、lm関数にフィットする「正しい」モデルを与えたので、「良い」モデルが生成されたのです。 実生活でこれを行うことはできませんが、線形モデリングの仮定を理解するのに役立つ演習です。

    さて、分散が一定でなかった場合はどうでしょうか。 標準偏差 2 に x の平方根を掛けたらどうでしょうか。 下のプロットでわかるように、xが大きくなると点の垂直方向の散らばりが大きくなります。

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

    2をsqrt(x)で乗じたのは、標準偏差を指定しているためです。 もし分散を指定できれば、4にxだけをかけたはずです。

    lmを使って同じモデルをあてはめると、Residual Standard errorは4.488になります。

    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

    データをシミュレーションしたので、これが正しくないことは分かっています。 ノイズ」を作るときに、一定の標準偏差はなかったのです。 それぞれのランダムな値は異なる正規分布から引かれ、平均は0ですが、標準偏差はxに従って変化しています。 つまり、分散が一定であるという仮定に反しているのです。

    最も一般的な方法は、残差と適合値をプロットすることです。 これはRで簡単にできます。モデル・オブジェクトに対してplotを呼び出すだけです。 これは従来のモデリング仮定を評価するために4つの異なるプロットを生成します。 詳しくはこのブログ記事をご覧ください。 

    plot(m2, which = 1)

    plot(m2, which = 3)

    最初のプロットは、右側に行くに従って0付近の変動が大きくなり、フィット値が大きくなっていることがわかる。 3番目のプロットでも右に行くにつれて変動が大きくなっているのがわかりますが、今回は残差は標準化され平方根のスケールに変換されています。

    さて、分散が一定でないという仮定が無効であることが確認されたので、私たちは何ができるでしょうか。 1つの方法は、データを対数変換することです。 これは応答変数が正で大きく偏っている場合に有効な場合があります。 しかし、ここではそのようなことはありません。 yはわずかに傾いているだけです。 (確認するためにyでhist()を呼び出してください。) さらに、我々のデータは対数スケールでないことを知っています。

    もうひとつのアプローチは、一定でない分散をモデル化することです。

    これを行うには、Rの基本インストールに含まれるnlmeパッケージの関数を使用します。 この関数はlm関数と同じように使いますが、weightsの引数といくつかの分散関数も使います。 4083>

    以下では、一定でない分散を示すデータに対して「正しい」モデルを当てはめました。 gls関数1が使えるようにnlmeパッケージをロードします。 モデルの構文は前と同じようにy ~ xを指定します。 そして、weights引数で分散関数、この場合はvarFixedを指定しますが、これもnlmeパッケージの一部です。 これは分散関数がパラメータを持たず、共変量が1つのxであり、まさにデータを生成する方法であることを示しています。 ~xの構文は、「xの関数としてのモデル分散」として読むことができる片側式です。

    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

    結果は、残差標準誤差を1.925とし、データを生成するのに我々が使った「真の」値である2に非常に近くなっています。 切片と傾きの値、1.37 と 2.09 も 1.2 と 2.1 に非常に近い値です。

    もう一度、モデルオブジェクトで plot を呼び出すことができます。 

    plot(vm1)

    このプロットは良さそうです。 分散をxの関数としてモデル化する限り、フィットしたモデルは系統的な方法でオーバーフィットもアンダーフィットもしません（lmを使ってモデルをフィットさせ、一定の分散を仮定したときとは異なります）

    関数varFixedは、各オブザベーションに対して重み、つまり \(w_i\) とシンボル化したものを作成します。 これは、ある観測値が高い重みを持っているほど、そのノイズ成分を引いた正規分布の分散が小さくなるという考え方です。 我々の例では、xが増加すると、分散が増加します。 したがって、より高い重みはより小さいx値と関連づけられるべきです。 これは、xの逆数をとることで実現できます。つまり、(w_i = 1/xenta)です。 つまり、”we_i = 1/x “であれば、”we_i = 1/2 “です。 となり、(x = 10) の場合は、 \(w = 1/10) となります。 最後に、重みが大きいと分散が小さくなるように、定数分散を重みで割ります。 数学的に言うと、

    または平方根を取って標準偏差を表現すると、

    つまり分母が大きいと（つまり重みが大きく、xが小さい）、分散が小さく、観測されたyがより正確になるわけですね。 これもgls関数と同じようにweights引数を持ちます。 唯一の違いは、重みを表現する方法についてより明確にしなければならないことです。 この場合、xの逆数を指定しなければなりません。 4083> 

    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

    nlmeパッケージの強力な点は、さまざまな分散関数を使用できることです。 先ほど例示したvarFixed関数は最も単純なもので、lmでできるものです。

    • varIdent
    • varPower
    • varExp
    • varConstPower
    • varComb

    シミュレーションデータで、それぞれの分散を調べてみましょう。 その効果を見るために、まずこの性質を持つデータをシミュレートしてみます。 以下では、誰かが同じランダムなデータをシミュレートしたい場合に備えて、set.seed(11)を使用しています。 次にnを400に設定し、シミュレーションするオブザベーションの数にします。 xは前と同じように生成されます。 この例では、gという予測変数が追加されていますが、これは性別と考えることができます。 我々は、”m “と “f “の400の値を無作為にサンプリングします。 次に、2つの標準偏差のベクトルをシミュレートし、msdとして保存します。 g == "m"の場合、標準偏差は୧⃛(๑⃙⃘◡̈๑⃙⃘)୨⃛になります。 それ以外の場合は、g == "f"の場合は”˶‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾”です。 このベクトルを次の行で使ってyを生成します。 msdがrnorm関数にプラグインされているところに注目してください。 また、yの生成方法にも注目してください。 これは、xとyの間に交互作用があるためです。 g == "f"のとき、切片と傾きは1.2と2.1です。 g == "m"の場合は、切片と傾きはそれぞれ(1.2 + 2.8), (2.1 + 2.8)となる。 最後にシミュレーションしたデータをデータフレームに入れ、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()

    gの各レベルで分散が異なることに注意してください。 m “の分散は “f “の分散よりずっと大きい。 また、”f “よりも多くの散らばりがあります。 そのようにデータをシミュレートしてみました。 mの分散はfの分散の2.5倍としました。

    さて、これらの真の値を復元するためにvarIdentでgls関数を使いましょう。 先ほどと同じ方法で、weights引数で分散関数を定義します。 以下、form引数で分散のモデル化式がgを条件とすることを指定しています。 ~ 1|gという式は、gの水準間で分散が異なるという片寄った式です。 1は、分散のモデルに共変量を追加していないことを意味しているだけです。 4101＞のような式を入れることも可能ですが、この場合は分散を生成する際に＜6204＞を使っていないので不正解になります。 プロットを見ると、yの変動はグループ間で異なっていますが、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

    まず出力の一番下の推定残留標準誤差が2であることに注目してください。また、「分散関数」のセクションでは、グループ「m」の推定値が2.58127であり、グループ「m」の異なる分散を生成するために使用した「真の」値2.5に非常に近いことに注目してください。 一般に、varIdent関数を用いて層別の分散を推定する場合、層の1つを基準値とし、他は残差標準誤差の倍数として推定されます。 この場合、アルファベット的には「f」が「m」の前に来るので、「f」をベースライン、つまり「1」に設定した。 f群の推定残差標準誤差は、Ⓐ（2.005397Ⓑ）である。 グループ “m “の推定残留標準誤差は、 \(2.005397 \times 2.58127).

    It’s important to note these are just estimated. これらの推定値の不確かさを知るには、nlmeパッケージのintervals関数を使って95%信頼区間を求めることができる。 

    int <- intervals(vm2)

    varStruct要素は分散関数で推定されたパラメーターの95%信頼区間を含んでいます。 この場合のパラメータは、男性グループの残留標準誤差乗数です。

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

    sigma要素は、残留標準誤差の95%信頼区間を含んでいます。

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

    両方の区間は非常に小さく、データを生成するために使用した「真の」値を含みます。

    varPower

    関数により、共変数の絶対値のべき乗として、分散をモデル化することが可能になります。 もう一度、これを説明するために、この性質を持つデータをシミュレートしてみましょう。 以下、注目すべきはrnorm関数のsd引数です。 これは標準偏差を2として、それにxの1.5乗の絶対値をかけているところです。 これは、上記のvarFixed関数のデモのためにデータをシミュレートした方法と似ています。 この場合、単純に指数を0.5と仮定した。 (数の平方根を取ることは0.5の累乗に相当することを思い出してください)。 ここでは、任意に1.5の累乗を選びました。 

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

    予測変数が増加すると分散が増加するという典型的な形のデータであることが分かります。 このデータをglsを使って正しくモデル化するために、y ~x式で与え、varPowerでweights引数を使用します。 varFixed関数で行ったのと同じように、片側だけの式を指定していることに注意してください。 4083> 

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

    powerは1.52と見積もられ、「真の」値1.5に非常に近くなっていることがわかるでしょう。 推定残差標準誤差1.8729439も2にかなり近いです。両方の区間は、我々がデータをシミュレートするために使用した値を捕捉しています。 一方、モデルの係数は、やや貧弱に推定されています。 これは、yの変動幅が大きいことを考えると、特に \(x > 2

    ).

    varExp

    関数では、共変量の指数関数として、分散をモデル化することが可能です。 さらにまた、この分散関数をシミュレーションデータを使って説明します。 変更点はrnorm関数のsd引数のみです。 固定値である2をxに掛け、それに0.5を掛けて指数化しています。 xを増加させると分散が急速に増加することに注目してください。 

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

    逆算してこれらの値を回復するには、glsのweights引数にvarExp関数を使用します。 この場合、片側だけの公式は変わりません。 分散をxの関数としてモデル化すると書いてあります。 varExp関数は、xがある値で掛けられ、指数化されているので、glsはその値を推定しようとします。

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

    「分散関数」セクションでのexponの推定値が 0.4845623 は我々の指定値 0.5 に非常に近くなっていることがわかります。 同様に推定残差標準誤差2.1191918は「真の」値である2に近いです。モデル係数推定値もデータを生成するために使用した値に近いですが、(Intercept)の不確実性に注意してください。 サマリーの仮説検定では、負の切片を除外することができません。 切片はxが0になったときのyの値なので、切片の推定値は観測していない事象に対する外挿値となります。

    varConstPower

    この関数により、正の定数に共変量の絶対値のべき乗を加えたものとして分散をモデル化することができます。 というと語弊がありますが、これは基本的にvarPower関数と同じで、今は正の定数を加えています。 次のコードはvarConstPower関数の使用に適したデータをシミュレートしています。 これは、abs関数のxに0.7を加えた以外は、varPowerセクションのデータのシミュレーションに使ったコードと同じであることに注意してください。 なぜ0.7なのか？ 特別な理由はありません。 4083> 

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

    このデータの正しいモデル化方法は、glsとvarConstPower関数を使用することです。 片側式は以前と同じです。 要約は分散モデルに対する3つの推定値、すなわち定数、べき乗、残差標準誤差を返します。 4083> 

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

    定数の間隔はかなり広いですが、間隔は「真の」値をとらえています。

    varComb

    最後に、varComb関数は、対応する分散関数を掛け合わせることにより、2つ以上の分散モデルの組み合わせをモデル化することを可能にします。 明らかにこれは非常に複雑な分散モデルに対応することができます。

    以下では、2つの異なる分散処理を組み合わせています。

    • 一つは、標準偏差が男女で異なるもの（「f」= \(2), 「m」= \(2 \times 2.5))
    • One that allow standard deviation as increasing as x are multiplied by 1.5 and exponentiated

    データの可視化のためにxを1から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()

    プロットは、xが増加するにつれて分散が増加しますが、性別による分散の違いも示しています。 グループ「m」のyの分散はグループ「f」のyの分散よりはるかに大きく、特にxが1.5より大きいときです。

    上で指定したデータ生成過程を正しくモデル化して真の値を復元しようと、varComb関数を使ってさらに二つの分散関数にラップして使用します。 varIdentとvarExpです。 なぜこの2つなのか。 4083> 

    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

    サマリーセクションには、分散をモデル化するための2つのセクションがあります。 1つ目は “m “グループの乗数を2.437と推定し、これは真の値である2.5に非常に近い値である。 指数パラメータは1.54と推定され、データ作成時に使用した真の値1.5に極めて近い値です。 最後に残差の標準誤差は約1.79と推定され、真の値である2に近い値となっています。 varStruct 要素の A と B という接頭辞は、2種類の分散モデルのラベルにすぎません。

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

    残念ながら大きな指数変動により、モデル係数の推定値はひどく悪いものになっています。

    何のために？

    では、なぜこれほどまでに偽のデータをシミュレーションし、「真の」値を回復しようとしたのでしょうか？ それは何の役に立つのでしょうか？ 公平な質問です。 その答えは、モデルを指定し適合させるときに、どのような仮定をしているのかを理解するのに役立つからです。 

    m <- lm(y ~ x1 + x2)

    …と指定し、フィットさせるとき、我々はyがx1とx2の加重和（プラス切片）にほぼ等しく、誤差は平均0、ある決まった標準偏差の正規分布からの無作為抽出であると考えていることを意味しています。 同様に、次のモデルを指定してフィットさせると…

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

    …我々は、yはx1とx2の重み付き和（＋切片）とほぼ等しいと考えており、エラーは、平均0、ある乗で上昇するx1の倍数となる標準偏差の正規分布からランダムに抽出される、と述べているのです。

    これらのモデルに適したデータをシミュレートできれば、それらのモデルを使用する際の前提条件をよりよく理解できるようになります。 nlme パッケージを使用して分散をモデル化するために何ができるのか、理解が深まることを願っています。

    この記事に関する質問または説明については、UVA Library StatLab: [email protected]

    View the entire collection of UVA Library StatLab articles.

    Clay Ford
    Statistical Research Consultant
    University of Virginia Library
    April 7, 2020

    ＜512>1.この記事に関する質問は、以下の連絡先にご連絡ください。 nlme パッケージは、おそらく混合効果モデル（つまり、固定効果とランダム効果の両方を持つモデル）の適合に使用される lme 関数でよりよく知られています。 このブログでは、ランダム効果にフィットしないglsを使った分散関数を紹介します。 しかし、このブログ記事で紹介するものはすべてlmeで使用できます。

    コメントを残す

    メールアドレスが公開されることはありません。