9

問題：高次のパラメーター（つまり、交互作用）がモデルに残っている限り、モデルの低次のパラメーター（たとえば、主効果パラメーター）を削除できません。その場合でも、モデルはリファクタリングされ、新しいモデルは上位モデルにネストされません。
次の例を参照してください（私が使用しているANOVAから来ているためcontr.sum）：

d <- data.frame(A = rep(c("a1", "a2"), each = 50), B = c("b1", "b2"), value = rnorm(100))
options(contrasts=c('contr.sum','contr.poly'))
m1 <- lm(value ~ A * B, data = d)
m1

## Call:
## lm(formula = value ~ A * B, data = d)
## 
## Coefficients:
## (Intercept)           A1           B1        A1:B1  
##   -0.005645    -0.160379    -0.163848     0.035523  

m2 <- update(m1, .~. - A)
m2

## Call:
## lm(formula = value ~ B + A:B, data = d)

## Coefficients:
## (Intercept)           B1       Bb1:A1       Bb2:A1  
##   -0.005645    -0.163848    -0.124855    -0.195902

ご覧のとおり、1つのパラメーター（A）を削除しましたが、新しいモデル（m2）はリファクタリングされ、より大きなモデル（ ）にネストされていませんm1。手あたりの因子を数値コントラスト変数に変換すると、目的の結果を得ることができますが、Rの因子機能を使用してどのように得るのですか？

質問： Rの低次因子を削除して、このパラメーターを実際に見逃し、リファクタリングされていないモデルを取得するにはどうすればよいですか（つまり、小さいモデルのパラメーターの数を少なくする必要があります）。

しかし、なぜ？パッケージの関数lmerを使用して、モデルのp値のような「タイプ3」を取得したいと思います。したがって、この例は実際には単なる例です。KRmodcomppbkrtest

なぜCrossValidatedではないのですか？これは実際には統計の質問よりもRの方が多いと感じています（つまり、モデルを相互作用に適合させるべきではないが、主要な効果の1つを持たないようにする必要があることはわかっていますが、それでもやりたいと思います）。

4

2 に答える 2

8

これが一種の答えです。このモデルを式で直接定式化する方法はありません...

上記のようにデータを作成します。

d <- data.frame(A = rep(c("a1", "a2"), each = 50),
                B = c("b1", "b2"), value = rnorm(100))
options(contrasts=c('contr.sum','contr.poly'))

数式から係数を引くだけでは機能しないという元の結果を確認します。

m1 <- lm(value ~ A * B, data = d)
coef(m1)
## (Intercept)          A1          B1       A1:B1 
## -0.23766309  0.04651298 -0.13019317 -0.06421580 

m2 <- update(m1, .~. - A)
coef(m2)
## (Intercept)          B1      Bb1:A1      Bb2:A1 
## -0.23766309 -0.13019317 -0.01770282  0.11072877

新しいモデル行列を作成します。

X0 <- model.matrix(m1)
## drop Intercept column *and* A from model matrix
X1 <- X0[,!colnames(X0) %in% "A1"]

lm.fitモデル行列を直接指定できます。

m3 <- lm.fit(x=X1,y=d$value)
coef(m3)
## (Intercept)          B1       A1:B1 
## -0.2376631  -0.1301932  -0.0642158

この方法は、モデル行列を明示的に指定できるいくつかの特殊なケースでのみ機能します（例lm.fit：）glm.fit。

より一般的には：

## need to drop intercept column (or use -1 in the formula)
X1 <- X1[,!colnames(X1) %in% "(Intercept)"]
## : will confuse things -- substitute something inert
colnames(X1) <- gsub(":","_int_",colnames(X1))
newf <- reformulate(colnames(X1),response="value")
m4 <- lm(newf,data=data.frame(value=d$value,X1))
coef(m4)
## (Intercept)          B1   A1_int_B1 
##  -0.2376631  -0.1301932  -0.0642158

このアプローチには、複数の入力変数が同じ予測変数に由来するものとして認識されないという欠点があります（つまり、2レベルを超える因子からの複数の因子レベル）。

于 2012-07-08T20:18:03.597 に答える
6

最も簡単な解決策はを使用することだと思いますmodel.matrix。おそらく、いくつかの派手なフットワークとカスタムコントラストであなたが望むものを達成することができます。ただし、「タイプ3のような」p値が必要な場合は、モデル内のすべての項に必要になる可能性があります。その場合、model.matrix1つの列をドロップするすべてのモデルを暗黙的にループできるため、とにかく便利だと思います。時間。可能なアプローチの提供は、その統計的メリットを支持するものではありませんが、明確な質問を作成し、統計的に不健全である可能性があることを知っているように思われるので、答えない理由はありません。

## initial data
set.seed(10)
d <- data.frame(
  A = rep(c("a1", "a2"), each = 50),
  B = c("b1", "b2"),
  value = rnorm(100))

options(contrasts=c('contr.sum','contr.poly'))

## create design matrix
X <- model.matrix(~ A * B, data = d)

## fit models dropping one effect at a time
## change from 1:ncol(X) to 2:ncol(X)
## to avoid a no intercept model
m <- lapply(1:ncol(X), function(i) {
  lm(value ~ 0 + X[, -i], data = d)
})
## fit (and store) the full model
m$full <- lm(value ~ 0 + X, data = d)
## fit the full model in usual way to compare
## full and regular should be equivalent
m$regular <- lm(value ~ A * B, data = d)
## extract and view coefficients
lapply(m, coef)

これにより、次の最終出力が得られます。

[[1]]
   X[, -i]A1    X[, -i]B1 X[, -i]A1:B1 
  -0.2047465   -0.1330705    0.1133502 

[[2]]
X[, -i](Intercept)          X[, -i]B1       X[, -i]A1:B1 
        -0.1365489         -0.1330705          0.1133502 

[[3]]
X[, -i](Intercept)          X[, -i]A1       X[, -i]A1:B1 
        -0.1365489         -0.2047465          0.1133502 

[[4]]
X[, -i](Intercept)          X[, -i]A1          X[, -i]B1 
        -0.1365489         -0.2047465         -0.1330705 

$full
X(Intercept)          XA1          XB1       XA1:B1 
  -0.1365489   -0.2047465   -0.1330705    0.1133502 

$regular
(Intercept)          A1          B1       A1:B1 
 -0.1365489  -0.2047465  -0.1330705   0.1133502

これは、を使用するモデルにとってこれまでのところ素晴らしいことですlm。これは最終的にはのためであるとおっしゃってlmer()いたので、混合モデルを使用した例を次に示します。ランダムな切片以上のものがある場合（つまり、モデルの固定部分とランダム部分から効果を削除する必要がある場合）、より複雑になる可能性があると思います。

## mixed example
require(lme4)

## data is a bit trickier
set.seed(10)
mixed <- data.frame(
  ID = factor(ID <- rep(seq_along(n <- sample(3:8, 60, TRUE)), n)),
  A = sample(c("a1", "a2"), length(ID), TRUE),
  B = sample(c("b1", "b2"), length(ID), TRUE),
  value = rnorm(length(ID), 3) + rep(rnorm(length(n)), n))

## model matrix as before
X <- model.matrix(~ A * B, data = mixed)

## as before but allowing a random intercept by ID
## becomes trickier if you need to add/drop random effects too
## and I do not show an example of this
mm <- lapply(1:ncol(X), function(i) {
  lmer(value ~ 0 + X[, -i] + (1 | ID), data = mixed)
})

## full model
mm$full <- lmer(value ~ 0 + X + (1 | ID), data = mixed)
## full model regular way
mm$regular <- lmer(value ~ A * B + (1 | ID), data = mixed)

## view all the fixed effects
lapply(mm, fixef)

それは私たちに...

[[1]]
   X[, -i]A1    X[, -i]B1 X[, -i]A1:B1 
 0.009202554  0.028834041  0.054651770 

[[2]]
X[, -i](Intercept)          X[, -i]B1       X[, -i]A1:B1 
        2.83379928         0.03007969         0.05992235 

[[3]]
X[, -i](Intercept)          X[, -i]A1       X[, -i]A1:B1 
        2.83317191         0.02058800         0.05862495 

[[4]]
X[, -i](Intercept)          X[, -i]A1          X[, -i]B1 
        2.83680235         0.01738798         0.02482256 

$full
X(Intercept)          XA1          XB1       XA1:B1 
  2.83440919   0.01947658   0.02928676   0.06057778 

$regular
(Intercept)          A1          B1       A1:B1 
 2.83440919  0.01947658  0.02928676  0.06057778
于 2012-07-08T20:35:19.627 に答える