由于相似性和 pvalues 而未定义 Coef
Coef not defined because of similarity and pvalues
我有一个用于实验的 CRM 数据集,其中虚拟 W 对应于 treatment/control 组(参见下面的代码)。当我测试 W 与其他特征的独立性时,我意识到两件事:
- 使用 model.matrix 时,一些系数(此虚拟数据集中的 1)由于相似性而未定义。直接将 DT 馈送到 lm()
时不会发生这种情况
- 两种情况下获得的模型产生不同的结果,即各个特征的 p 值发生变化
我(认为我)理解多重共线性的概念,但在这种特殊情况下我不太理解 a) 为什么会出现 b) 为什么它对 model.matrix 和lm
我错过了什么?
非常感谢!
set.seed(1)
n = 302
DT = data.table(
zipcode = factor(sample(seq(1,52), n, replace=TRUE)),
gender = factor(sample(c("M","F"), n, replace=TRUE)),
age = sample(seq(1,95), n, replace=TRUE),
days_since_last_purchase = sample(seq(1,259), n, replace=TRUE),
W = sample(c(0,1), n, replace=TRUE)
)
summary(DT)
m = model.matrix(W ~ . +0, DT)
f1 = lm(DT$W ~ m)
f2= lm(W~ ., DT)
p_value_ratio <- function(lm)
{
summary_randomization = summary(lm)
p_values_randomization = summary_randomization$coefficients[, 4]
L = length(p_values_randomization)
return(sum(p_values_randomization <= 0.05)/(L-1))
}
all.equal(p_value_ratio(f1), p_value_ratio(f2))
alias(f1)
alias(f2)
您的问题是 model.matrix
中的 + 0
。第二次拟合包括模型矩阵中的截距。如果排除它,则会排除较少的因子水平(通常由截距表示):
colnames(model.matrix(W ~ ., DT))
#excludes zipcode1 and genderf since these define the intercept
colnames(model.matrix(W ~ . + 0, DT))
#excludes only genderf
请注意 f1
包含一个拦截,它由 lm
添加(我相信是通过对 model.matrix
的内部调用,但尚未检查):
m = model.matrix(W ~ . + 0, DT);
f1 = lm(DT$W ~ m );
model.matrix(f1)
你可能想要这个:
m = model.matrix(W ~ ., DT);
f1 = lm(DT$W ~ m[,-1]);
(如果想直接使用lm.fit
,通常只能手动构建模型矩阵。)
f2= lm(W~ ., DT);
all.equal(unname(coef(f1)), unname(coef(f2)))
#[1] TRUE
最后,这归结为您对治疗对比的理解。通常,您不应从模型矩阵中排除截距。
我有一个用于实验的 CRM 数据集,其中虚拟 W 对应于 treatment/control 组(参见下面的代码)。当我测试 W 与其他特征的独立性时,我意识到两件事:
- 使用 model.matrix 时,一些系数(此虚拟数据集中的 1)由于相似性而未定义。直接将 DT 馈送到 lm() 时不会发生这种情况
- 两种情况下获得的模型产生不同的结果,即各个特征的 p 值发生变化
我(认为我)理解多重共线性的概念,但在这种特殊情况下我不太理解 a) 为什么会出现 b) 为什么它对 model.matrix 和lm
我错过了什么?
非常感谢!
set.seed(1)
n = 302
DT = data.table(
zipcode = factor(sample(seq(1,52), n, replace=TRUE)),
gender = factor(sample(c("M","F"), n, replace=TRUE)),
age = sample(seq(1,95), n, replace=TRUE),
days_since_last_purchase = sample(seq(1,259), n, replace=TRUE),
W = sample(c(0,1), n, replace=TRUE)
)
summary(DT)
m = model.matrix(W ~ . +0, DT)
f1 = lm(DT$W ~ m)
f2= lm(W~ ., DT)
p_value_ratio <- function(lm)
{
summary_randomization = summary(lm)
p_values_randomization = summary_randomization$coefficients[, 4]
L = length(p_values_randomization)
return(sum(p_values_randomization <= 0.05)/(L-1))
}
all.equal(p_value_ratio(f1), p_value_ratio(f2))
alias(f1)
alias(f2)
您的问题是 model.matrix
中的 + 0
。第二次拟合包括模型矩阵中的截距。如果排除它,则会排除较少的因子水平(通常由截距表示):
colnames(model.matrix(W ~ ., DT))
#excludes zipcode1 and genderf since these define the intercept
colnames(model.matrix(W ~ . + 0, DT))
#excludes only genderf
请注意 f1
包含一个拦截,它由 lm
添加(我相信是通过对 model.matrix
的内部调用,但尚未检查):
m = model.matrix(W ~ . + 0, DT);
f1 = lm(DT$W ~ m );
model.matrix(f1)
你可能想要这个:
m = model.matrix(W ~ ., DT);
f1 = lm(DT$W ~ m[,-1]);
(如果想直接使用lm.fit
,通常只能手动构建模型矩阵。)
f2= lm(W~ ., DT);
all.equal(unname(coef(f1)), unname(coef(f2)))
#[1] TRUE
最后,这归结为您对治疗对比的理解。通常,您不应从模型矩阵中排除截距。