2016-07-01 20 views
2

假设我们有一个玩具数据帧:model.matrix():为什么我失去的对比度控制在这种情况下

x <- data.frame(x1 = gl(3, 2, labels = letters[1:3]), 
       x2 = gl(3, 2, labels = LETTERS[1:3])) 

我想构建一个模型矩阵

# x1b x1c x2B x2C 
# 1 0 0 0 0 
# 2 0 0 0 0 
# 3 1 0 1 0 
# 4 1 0 1 0 
# 5 0 1 0 1 
# 6 0 1 0 1 

model.matrix(~ x1 + x2 - 1, data = x, 
      contrasts.arg = list(x1 = contr.treatment(letters[1:3]), 
            x2 = contr.treatment(LETTERS[1:3]))) 

但实际上我得到:

# x1a x1b x1c x2B x2C 
# 1 1 0 0 0 0 
# 2 1 0 0 0 0 
# 3 0 1 0 1 0 
# 4 0 1 0 1 0 
# 5 0 0 1 0 1 
# 6 0 0 1 0 1 
# attr(,"assign") 
# [1] 1 1 1 2 2 
# attr(,"contrasts") 
# attr(,"contrasts")$x1 
# b c 
# a 0 0 
# b 1 0 
# c 0 1 

# attr(,"contrasts")$x2 
# B C 
# A 0 0 
# B 1 0 
# C 0 1 

我有点困惑在这里:

  • 我已经明确的对比矩阵传递给第一滴因子水平;
  • 我要求丢弃拦截。

那么为什么我会得到5列的模型矩阵?我怎样才能得到我想要的模型矩阵?

+0

我认为这样做是正确的。你不能只有只有零的行。 –

+0

正如你在答案中所说的那样,总是可以自己创建虚拟变量。用这种方法,你可以得到你想要的精确模型矩阵。然而,如果你有只包含零的行,矩阵将是单数的(用统计学的话来说,你有“完美的共线性”),这意味着无法获得参数估计值。请参阅:https://stats.stackexchange.com/questions/70899/ –

+0

如果您确实需要破碎的模型矩阵,让它有一个截距,然后放下列:'model.matrix(〜x1 + x2,data = x)[,-1]'。虽然我无法真正想象这会对...有用... – Gregor

回答

0

每当我们失去对R级别的控制时,在C级别上必定存在一些默认的,不可更改的行为。为model.matrix.default() C代码可以在R源之中包被发现在:

R-<release_number>/src/library/stats/src/model.c 

我们可以在这里找到的解释:

/* If there is no intercept we look through the factor pattern */ 
/* matrix and adjust the code for the first factor found so that */ 
/* it will be coded by dummy variables rather than contrasts. */ 

让我们对这个小测试,用数据帧

x <- data.frame(x1 = gl(2, 2, labels = letters[1:2]), x2 = sin(1:4)) 
  1. 如果我们在RHS上只有x2,我们可以放下拦截器成功:

    model.matrix(~ x2 - 1, data = x) 
    #   x2 
    #1 0.8414710 
    #2 0.9092974 
    #3 0.1411200 
    #4 -0.7568025 
    
  2. 如果我们只在RHS x1,对比度不适用:

    model.matrix(~ x1 - 1, data = x) 
    # x1a x1b 
    #1 1 0 
    #2 1 0 
    #3 0 1 
    #4 0 1 
    
  3. 当我们都x1x2,对比度不适用:

    model.matrix(~ x1 + x2 - 1, data = x) 
    # x1a x1b   x2 
    #1 1 0 0.8414710 
    #2 1 0 0.9092974 
    #3 0 1 0.1411200 
    #4 0 1 -0.7568025 
    

这意味着虽然有差异NCE之间:

lm(y ~ x2, data = x) 
lm(y ~ x2 - 1, data = x) 

存在不能确保数值稳定性,但要保证灵敏度

lm(y ~ x1, data = x) 
lm(y ~ x1 - 1, data = x) 

lm(y ~ x1 + x2, data = x) 
lm(y ~ x1 + x2 - 1, data = x) 

的原因这样的行为之间没有差异的估计/预测。如果我们真的降了拦截,同时将对比x1,我们结束了一个模型矩阵:

# x1b 
    #1 0 
    #2 0 
    #3 1 
    #4 1 

的效果是,我们限制估计a级别设置为0

在这篇文章:How can I force dropping intercept or equivalent in this linear model?,我们有一个数据集:

#   Y X1 X2 
#1 1.8376852 TRUE TRUE 
#2 -2.1173739 TRUE FALSE 
#3 1.3054450 FALSE TRUE 
#4 -0.3476706 TRUE FALSE 
#5 1.3219099 FALSE TRUE 
#6 0.6781750 FALSE TRUE 

这个数据集中没有联合存在(X1 = FALSE, X2 = FALSE)。但从广义上讲,model.matrix()必须做一些安全和明智的事情。假设在训练数据集中不存在两个因子水平的联合存在意味着它们不需要被预测,这是有偏差的。如果我们在应用对比的时候真的放弃了拦截,那么这种联合存在被限制为0.然而,这个帖子的OP故意要这样的非标准行为(出于某种原因),在这种情况下,在我的答案中给出了可能的解决方法。