उच्च क्रम पैरामीटर बने रहने पर मॉडल में निचले क्रम पैरामीटर को कैसे निकालें?

Question

समस्या: जब तक उच्च क्रम पैरामीटर (यानी, इंटरैक्शन) मॉडल में रहते हैं, तब तक मैं मॉडल में निम्न ऑर्डर पैरामीटर (उदा।, मुख्य प्रभाव पैरामीटर) को हटा नहीं सकता। ऐसा करने पर भी, मॉडल को दोहराया जाता है और नए मॉडल को उच्च मॉडल में घोंसला नहीं दिया जाता है।
निम्न उदाहरण देखें (के रूप में मैं 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) को हटाने रिफैक्टर्ड और है बड़े मॉडल में नेस्टेड नहीं किया (m1)। यदि मैं संख्यात्मक विपरीत चर में प्रति हाथ अपने कारकों को बदलता हूं तो मुझे वांछित परिणाम मिल सकते हैं, लेकिन मैं आर की कारक क्षमताओं का उपयोग करके इसे कैसे प्राप्त करूं?

प्रश्न: मैं कैसे आर में एक निचले क्रम के कारक को हटाने और एक मॉडल है कि वास्तव में इस पैरामीटर याद करते हैं और पुनर्संशोधित नहीं है (अर्थात, छोटे मॉडल में पैरामीटर की संख्या कम होना चाहिए) प्राप्त कर सकते हैं?

---

लेकिन क्यों? मैं pbkrtest पैकेज से KRmodcomp फ़ंक्शन का उपयोग करके lmer मॉडल के लिए पी-मानों जैसे 'टाइप 3' प्राप्त करना चाहता हूं। तो यह उदाहरण वास्तव में सिर्फ एक उदाहरण है।

क्रॉस वैलिडेटेड क्यों नहीं? मुझे यह एहसास है कि यह वास्तव में एक आर के आंकड़े हैं और फिर आंकड़े प्रश्न (यानी, मुझे पता है कि आपको कभी भी इंटरैक्शन के साथ मॉडल को फिट नहीं करना चाहिए, लेकिन मुख्य प्रभावों में से एक के बिना, लेकिन मैं अभी भी इसे करना चाहता हूं)।

Answer 1

यहां एक प्रकार का उत्तर है;

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 स्तर कारक से एक ही भविष्यवक्ता से उत्पन्न के रूप में कई इनपुट चर पहचान नहीं होगा)।

Answer 2

मुझे लगता है कि सबसे सरल समाधान model.matrix का उपयोग करना है।संभवतः, आप कुछ फैंसी फुटवर्क और कस्टम विरोधाभासों के साथ जो चाहते हैं उसे प्राप्त कर सकते हैं। हालांकि, अगर आप "टाइप 3 एस्क्यू" पी-वैल्यू चाहते हैं, तो शायद आप इसे अपने मॉडल में हर शब्द के लिए चाहते हैं, इस मामले में, मुझे लगता है कि 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