2012-07-04 23 views
9

Vấn đề: Tôi không thể xóa thông số thứ tự thấp hơn (ví dụ: thông số hiệu ứng chính) trong mô hình miễn là các tham số thứ tự cao hơn (tức là tương tác) vẫn còn trong mô hình. Ngay cả khi làm như vậy, mô hình được tái cấu trúc và mô hình mới không được lồng trong mô hình cao hơn.
Xem ví dụ sau (như tôi đến từ ANOVAs tôi sử dụng contr.sum):Làm thế nào để loại bỏ một tham số thứ tự thấp hơn trong một mô hình khi các tham số thứ tự cao hơn vẫn còn?

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 

Như có thể thấy, mặc dù tôi loại bỏ một tham số (A), mô hình mới (m2) được refactored và được không lồng nhau trong mô hình lớn hơn (m1). Nếu tôi chuyển đổi các yếu tố của tôi trên tay theo các biến số tương phản, tôi có thể nhận được kết quả mong muốn, nhưng làm thế nào để tôi có được nó bằng cách sử dụng các khả năng yếu tố của R?

Câu hỏi: Làm cách nào để loại bỏ yếu tố thứ tự thấp hơn trong R và lấy mô hình thực sự bỏ tham số này và không được cấu trúc lại (tức là số tham số trong mô hình nhỏ hơn phải thấp hơn)?


Nhưng tại sao? Tôi muốn lấy 'Loại 3' như giá trị p cho một mô hình lmer bằng cách sử dụng chức năng KRmodcomp từ gói pbkrtest. Ví dụ này thực sự là một ví dụ.

Tại sao không được CrossValidated? Tôi có cảm giác rằng đây thực sự là một câu hỏi về số liệu thống kê (nghĩa là, tôi biết rằng bạn không bao giờ nên phù hợp với mô hình có tương tác nhưng không có một trong những hiệu ứng chính, nhưng tôi vẫn muốn làm điều đó).

+1

đọc Bill Venables [http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf](http:// www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf) trên tổng loại III của hình vuông. Đó là một câu hỏi thống kê. – mnel

+3

Một cách để làm điều này là xây dựng ma trận thiết kế đầy đủ (sử dụng 'model.matrix'), xóa các cột bạn không muốn, và sau đó phù hợp với mô hình với các cột còn lại. Tôi sẽ làm ví dụ nếu/khi tôi có cơ hội ... –

+0

Hãy xem gói ['MixMod'] (http://cran.r-project.org/web/packages/MixMod/). Cơ sở 'R' sẽ không hỗ trợ điều này (xem bình luận trước đó của tôi về Bill Venables. – mnel

Trả lời

8

Đây là một loại câu trả lời; không có cách nào mà tôi biết để xây dựng mô hình này trực tiếp theo công thức ...

Construct dữ liệu như trên:

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

Xác nhận việc tìm kiếm ban đầu mà chỉ trừ yếu tố từ công thức không hoạt động :

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 

Xây dựng ma trận mô hình mới:

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

lm.fit cho phép đặc điểm kỹ thuật trực tiếp của ma trận mô hình:

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

Phương pháp này chỉ hoạt động cho một số trường hợp đặc biệt cho phép ma trận mô hình được chỉ định rõ ràng (ví dụ: lm.fit, glm.fit).

Tổng quát hơn:

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

Cách tiếp cận này có những bất lợi mà nó sẽ không nhận ra nhiều biến đầu vào như bắt nguồn từ những yếu tố dự báo tương tự (ví dụ, nồng độ yếu tố nhiều từ một yếu tố hơn-hơn-2 cấp).

+0

Cảm ơn câu trả lời tuyệt vời. Tôi sẽ chấp nhận câu trả lời của bạn (tôi nghĩ cả hai đều tương tự nhau), khi bạn chỉ ra cách xây dựng công thức và đề cập đến vấn đề của những người dự đoán phân loại với hơn hai cấp độ. – Henrik

5

Tôi nghĩ giải pháp đơn giản nhất là sử dụng model.matrix.Có thể, bạn có thể đạt được những gì bạn muốn với một số bước chân ưa thích và tương phản tùy chỉnh. Tuy nhiên, nếu bạn muốn giá trị p "loại 3 esque", bạn có thể muốn nó cho mọi thuật ngữ trong mô hình của bạn, trong trường hợp này, tôi nghĩ rằng cách tiếp cận của tôi với model.matrix là thuận tiện. tại một thời điểm. Việc cung cấp một cách tiếp cận có thể không phải là một sự chứng thực về thành tích thống kê của nó, nhưng tôi nghĩ rằng bạn đã xây dựng một câu hỏi rõ ràng và dường như biết nó có thể không có ý nghĩa thống kê vì vậy tôi thấy không có lý do gì để không trả lời nó.

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

Điều này dẫn đến kết quả cuối cùng này:

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

Đó là tốt đẹp cho đến nay cho các mô hình sử dụng lm. Bạn đã đề cập điều này là cuối cùng cho lmer(), vì vậy, đây là một ví dụ sử dụng các mô hình hỗn hợp. Tôi tin rằng nó có thể trở nên phức tạp hơn nếu bạn có nhiều hơn một đánh chặn ngẫu nhiên (tức là, các hiệu ứng cần phải được giảm từ các phần cố định và ngẫu nhiên của mô hình).

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

nào cho chúng ta ...

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

Cảm ơn rất nhiều vì câu trả lời tuyệt vời. Tôi sẽ trao cho bạn 100 điểm (như bạn đã chỉ ra cách sử dụng 'lmer') nhưng sẽ chấp nhận câu trả lời của Ben Bolker (xem đó là lý do). – Henrik

Các vấn đề liên quan