Come rimuovere un parametro di ordine inferiore in un modello quando i parametri di ordine superiore rimangono?

Question

Il problema: Non riesco a rimuovere un parametro di ordine inferiore (ad es. Un parametro degli effetti principali) in un modello fintanto che i parametri di ordine superiore (cioè le interazioni) rimangono nel modello. Anche in questo caso, il modello viene refactored e il nuovo modello non è nidificato nel modello superiore.
Vedere il seguente esempio (come io vengo da ANOVA io uso 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 

Come si può vedere, anche se tolgo un parametro (A), il nuovo modello (m2) viene riscritta ed è non nidificato nel modello più grande (m1). Se trasformo i miei fattori per mano in variabili di contrasto numeriche, posso ottenere i risultati desiderati, ma come faccio a usare le capacità del fattore R?

la domanda: Come posso rimuovere un fattore di ordine inferiore in R e ottenere un modello che manca davvero questo parametro e non è refactoring (vale a dire, il numero di parametri del modello più piccolo deve essere inferiore)?

---

Ma perché? Desidero ottenere "Tipo 3" come valori p per un modello lmer utilizzando la funzione KRmodcomp dal pacchetto pbkrtest. Quindi questo esempio è davvero solo un esempio.

Perché non CrossValidated? Ho la sensazione che questa sia davvero più una domanda di R allora una statistica (ad esempio, so che non dovresti mai adattare un modello con interazioni ma senza uno degli effetti principali, ma voglio comunque farlo).

Answer 1

Ecco una sorta di risposta; non c'è modo che io conosco per formulare questo modello direttamente dalla formula ...

dati Construct come sopra:

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

Confermare risultato originale che semplicemente sottraendo il fattore dalla formula non funziona :

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 

formulare la nuova matrice modello:

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

lm.fit permette l'immissione diretta del modello matrice:

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

Questo metodo funziona solo per alcuni casi speciali che consentono la matrice modello da specificare esplicitamente (ad esempio lm.fit, glm.fit).

Più in generale:

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

Questo approccio ha lo svantaggio che non riconoscerà più variabili di input come derivanti dal medesimo predittore (cioè, più livelli di fattore da un fattore più-che-2-livello).

Answer 2

Penso che la soluzione più semplice sia utilizzare model.matrix.Forse, potresti ottenere quello che vuoi con un gioco di fantasia e contrasti personalizzati. Tuttavia, se si desidera "p 3 valori di tipo 3 esque", probabilmente lo si desidera per ogni termine nel proprio modello, nel qual caso, penso che il mio approccio con model.matrix sia comunque conveniente perché è possibile collegare in modo implicito tutti i modelli lasciando cadere una colonna Al tempo. La previsione di un possibile approccio non è un'approvazione dei meriti statistici di esso, ma penso che tu abbia formulato una domanda chiara e sembri sapere che potrebbe non essere statisticamente statica, quindi non vedo alcun motivo per non rispondere.

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

Questo si traduce in questo output finale:

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

Che è bello finora per modelli con lm. Hai menzionato questo in definitiva per lmer(), quindi ecco un esempio che utilizza modelli misti. Credo che potrebbe diventare più complesso se si ha più di un'intercettazione casuale (vale a dire, gli effetti devono essere eliminati dalle parti fisse e casuali del modello).

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

Il che ci dà ...

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