¿Cómo eliminar un parámetro de orden inferior en un modelo cuando los parámetros de orden superior permanecen?

Question

Problema: No puedo eliminar un parámetro de orden inferior (por ejemplo, un parámetro de efectos principales) en un modelo, siempre que los parámetros de orden superior (es decir, las interacciones) permanezcan en el modelo. Incluso al hacerlo, el modelo se refactoriza y el nuevo modelo no se anida en el modelo superior.
Véase el siguiente ejemplo (como estoy viniendo ANOVAs utilizo 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 

Como se puede ver, aunque se quita un parámetro (A), el nuevo modelo (m2) es refactorizado y es no anidado en el modelo más grande (m1). Si transformo mis factores por mano en las variables de contraste numérico, puedo obtener los resultados deseados, pero ¿cómo puedo obtenerlo utilizando las capacidades del factor R?

La Pregunta: ¿Cómo puedo eliminar un factor de orden inferior en R y obtener un modelo que realmente pierde este parámetro y no se refactorizado (es decir, el número de parámetros en el modelo más pequeño debe ser menor)?

---

¿Pero por qué? Deseo obtener 'tipo 3' como valores p para un modelo lmer usando la función KRmodcomp del paquete pbkrtest. Entonces este ejemplo es solo un ejemplo.

¿Por qué no CrossValidated? Tengo la sensación de que esto es realmente más una pregunta de R que de estadísticas (es decir, sé que nunca debes encajar en un modelo con interacciones pero sin uno de los efectos principales, pero aún así quiero hacerlo).

Answer 1

Aquí hay un tipo de respuesta; no hay manera que yo sepa para formular este modelo directamente por la fórmula ...

datos Constructo que el anterior:

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

Confirmar constatación inicial de que simplemente restando el factor de la fórmula no funciona :

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 

formular la nueva matriz de modelo:

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

lm.fit permite la especificación directa de la matriz de modelo:

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

Este método sólo funciona para unos pocos casos especiales que permiten que la matriz de modelo que se especifica explícitamente (por ejemplo lm.fit, glm.fit).

más general:

## 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 

Este enfoque tiene la desventaja de que no reconocerá múltiples variables de entrada como derivado de la misma predictor (es decir, múltiples niveles de los factores de un factor de nivel más-que-2)

Answer 2

Creo que la solución más directa es usar model.matrix.Posiblemente, puede lograr lo que quiere con un trabajo de pies elegante y contrastes personalizados. Sin embargo, si desea valores p "tipo 3 esque", probablemente lo desee para cada término de su modelo, en cuyo caso, creo que mi enfoque con model.matrix es conveniente porque puede pasar implícitamente fácilmente todos los modelos que abandonan una columna a la vez La provisión de un posible enfoque no es un respaldo de los méritos estadísticos del mismo, pero creo que formuló una pregunta clara y parece saber que puede no ser estadísticamente correcta, así que no veo ninguna razón para no contestarla.

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

Esto da lugar a este resultado final:

[[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 

Eso es bueno hasta ahora para los modelos usando lm. Usted mencionó que esto es en última instancia para lmer(), así que aquí hay un ejemplo usando modelos mixtos. Creo que puede volverse más complejo si tiene más de una intercepción aleatoria (es decir, los efectos deben eliminarse de las partes fijas y aleatorias del modelo).

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

Lo que nos da ...

[[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