Как удалить параметр более низкого порядка в модели, когда параметры более высокого порядка остаются?

The problem: Я не могу удалить параметр более низкого порядка (например, параметр основных эффектов) в модели, пока параметры более высокого порядка (то есть взаимодействия) остаются в модели. Даже при этом модель подвергается рефакторингу, и новая модель не вкладывается в более высокую модель.
Смотрите следующий пример (поскольку я исхожу из ANOVAs, я использую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  

Как видно, хотя я убираю один параметр (A), новая модель (m2) подвергается рефакторингу и являетсяnot nested в большей модели (m1). Если я преобразую свои коэффициенты на руку в числовые переменные контрастности, я могу получить желаемые результаты, но как мне получить их, используя возможности фактора R?

The Question: Как я могу удалить фактор более низкого порядка в R и получить модель, которая действительно пропускает этот параметр и не подвергается рефакторингу (т. Е. Число параметров в меньшей модели должно быть меньше)?

But why? Я хочу получить «тип 3»; как р-значения дляlmer модель с использованиемKRmodcomp функция отpbkrtest пакет. Так что этот пример действительно просто пример.

Why not CrossValidated? У меня такое ощущение, что это действительно скорее вопрос R, чем вопрос статистики (то есть я знаю, что вам никогда не следует подбирать модель с взаимодействиями, но без одного из основных эффектов, но я все еще хочу это сделать).

 Gavin Simpson05 июл. 2012 г., 10:46
@ Генрик, я понимаю твою точку зрения (тип III ожидается в твоей области), но для меня это почти единственная наихудшая причина что-то делать; только потому, что другие делают это. Вы читали Exegeses Venable, на которые ссылается mnel? Прошло много времени с тех пор, как я его прочитал, но IIRC вышел далеко за рамки простой неприязни к типу III SS.
 mnel05 июл. 2012 г., 07:28
Посмотрите наMixMod package, БазаR не поддержит это (см. мой предыдущий комментарий о Билле Венейбле.
 Ben Bolker05 июл. 2012 г., 04:03
Один из способов сделать это - построить полную матрицу дизайна (используяmodel.matrix), удалите ненужные столбцы, а затем подгоните модель к оставшимся столбцам. Я сделаю пример, если / когда у меня будет шанс ...
 mnel05 июл. 2012 г., 02:22
Читайте Билла Венейблсаhttp://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf по типу III суммы квадратов. Это вопрос статистики.
 Henrik05 июл. 2012 г., 08:59
Я знаю, что тема сумм квадратов типа III является чувствительной для сообщества статистики. Однако я специально не просил никаких статистических советов. У меня технический вопрос о том, как сделать определенную вещь в R. Вы должны знать, что мое поле (психология) ожидает тип 3, поэтому я хотел бы получить результаты типа 3. Кроме того, почему отрицательный голос? Мой вопрос сформулирован разумно, содержит воспроизводимый пример и является конкретным. Поскольку объяснения нет, похоже, что оно основано на неприязни к результатам типа 3. Это довольно неразумно.

Ответы на вопрос(2)

Решение Вопроса

Вот своего рода ответ; я не знаю, как сформулировать эту модель непосредственно по формуле ...

Построить данные, как указано выше:

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 

Недостаток этого подхода состоит в том, что он не распознает несколько входных переменных как вытекающие из одного и того же предиктора (то есть многофакторные уровни из фактора более двух уровней).

 Henrik13 июл. 2012 г., 16:13
Спасибо за отличный ответ. Я приму ваш ответ (я думаю, что оба похожи), так как вы покажете, как построить формулу, и упомянете проблему категориальных предикторов с более чем двумя уровнями.

Я думаю, что самым простым решением является использованиеmodel.matrix, Возможно, вы могли бы достичь того, что вы хотите, с какой-то модной работой ног и нестандартными контрастами. Однако, если вы хотите "введите 3 esque" р-значения, вы, вероятно, хотите это для каждого термина в вашей модели, и в этом случае, я думаю, мой подход сmodel.matrix в любом случае это удобно, потому что вы можете легко неявно проходить по всем моделям, опуская один столбец за раз. Предоставление возможного подхода не является подтверждением его статистических достоинств, но я думаю, что вы сформулировали четкий вопрос и, кажется, знаете, что он может быть статистически несостоятельным, поэтому я не вижу причин не отвечать на него.

## 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 
 Henrik13 июл. 2012 г., 16:09
Большое спасибо за отличный ответ. Я награжу вас 100 баллами (как вы конкретно показали, как использоватьlmer) но примет ответ Бена Болкера (см. обоснование).

Ваш ответ на вопрос