有没有办法改变 R 在模型输出中标记交互参数的方式?
Is there a way to change the way R labels the interaction parameter in model output?
我遇到了一个看似简单但非常令人沮丧的问题。当您 运行 模型在 R 中具有交互项时,R 将生成的参数命名为 "var1:var2" 等。不幸的是,这种命名约定使我无法计算预测值和 CI 的新数据所在的位置必需,因为“:”不是可以包含在列 header 中的字符,并且原始数据框中的名称必须与 newdata 中的名称完全匹配。还有其他人遇到过这个问题吗?
这是我的代码示例:
wemedist2.exp = glm(survive/trials ~ sitedist + type + sitedist*type + roaddist, family = binomial(logexp(wemedata$expos)), data=wemedata)
summary(wemedist2.exp)
wemepredict3 = with(wemedata, data.frame(sitedist=mean(sitedist),roaddist=mean(roaddist), type=factor(1:2)))
wemepredict3 = cbind(wemepredict3, predict(wemedist2.exp, newdata = wemepredict3, type = "link", se = TRUE))
这会生成一个 table,其中包含每个变量在指定水平的预测值,但没有交互作用。
对于您的 newdata
数据框,您不应包含交互列。调用 predict
.
时,将为您计算交互变量的乘积(并乘以估计系数)
例如:
创建一些虚拟数据:
set.seed(1)
n <- 10000
X <- data.frame(x1=runif(n), x2=runif(n))
X$x1x2 <- X$x1 * X$x2
head(X)
# x1 x2 x1x2
# 1 0.2655087 0.06471249 0.017181728
# 2 0.3721239 0.67661240 0.251783646
# 3 0.5728534 0.73537169 0.421260147
# 4 0.9082078 0.11129967 0.101083225
# 5 0.2016819 0.04665462 0.009409393
# 6 0.8983897 0.13091031 0.117608474
b <- runif(4)
y <- b[1] + c(as.matrix(X) %*% b[-1]) + rnorm(n, sd=0.1)
拟合模型并比较估计系数与真实系数:
M <- lm(y ~ x1 * x2, X)
summary(M)
# Call:
# lm(formula = y ~ x1 * x2, data = X)
#
# Residuals:
# Min 1Q Median 3Q Max
# -0.43208 -0.06743 -0.00170 0.06601 0.37197
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.202040 0.003906 51.72 <2e-16 ***
# x1 0.128237 0.006809 18.83 <2e-16 ***
# x2 0.156942 0.006763 23.21 <2e-16 ***
# x1:x2 0.292582 0.011773 24.85 <2e-16 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.09906 on 9996 degrees of freedom
# Multiple R-squared: 0.5997, Adjusted R-squared: 0.5996
# F-statistic: 4992 on 3 and 9996 DF, p-value: < 2.2e-16
b
# [1] 0.2106027 0.1147864 0.1453641 0.3099322
创建要预测的示例数据,并进行预测。请注意,我们只创建 x1
和 x2
,并且 而不是 创建 x1:x2
:
X.predict <- data.frame(x1=runif(10), x2=runif(10))
head(X.predict)
# x1 x2
# 1 0.26037592 0.7652155
# 2 0.73988333 0.3352932
# 3 0.02650689 0.9788743
# 4 0.84083874 0.1446228
# 5 0.85052685 0.7674547
# 6 0.13568509 0.9612156
predict(M, newdata=X.predict)
# 1 2 3 4 5 6 7
# 0.4138194 0.4221251 0.3666572 0.3681432 0.6225354 0.4084543 0.4711018
# 8 9 10
# 0.7092744 0.3401867 0.2320834
或者...
另一种方法是通过计算交互项的乘积将交互包含在模型拟合数据中,然后将其也包含在新数据中。我们已经完成了上面第 1 点中的第一步,我们在其中创建了一个名为 x1x2
的列。
然后我们将拟合模型:lm(y ~ x1 + x2 + x1x2, X)
并对以下数据进行预测:
X.predict <- data.frame(x1=runif(10), x2=runif(10), x1x2=runif(10)
如果交互中涉及分类变量...
当您有涉及分类变量的交互时,模型会估计系数来描述属于每个级别相对于属于参考级别的影响。因此,例如,如果我们有一个连续预测变量 (x1
) 和一个分类预测变量 (x2
,水平为 a
、b
和 c
),则模型 y ~ x1 * x2
将估计六个系数,描述:
截距(即x1
为零时的预测y
且观测属于x2
的参考水平);
当观测值属于x2
参考水平时改变x1
的效果(即斜率,对于x2
参考水平);
属于二级的影响(即属于二级相对于属于参考级截距的变化);
属于第三级的影响(即由于属于第三级相对于属于参考级截距的变化);
x1
的效果因属于第二层级相对于属于参考层级的变化(即斜率的变化);和
由于属于第三级相对于属于参考级x1
的影响变化(即斜率变化)
如果您想拟合和预测模型 with/to 描述交互的预计算数据,您可以创建一个包含列的数据框:x1
; x2b
(二进制,表示观察是否属于级别b
); x2c
(二进制,表示观察是否属于级别c
); x1x2b
(x1
和 x2b
的乘积);和 x1x2c
(x1
和 x2c
的乘积)。
一个快速的方法是 model.matrix
:
set.seed(1)
n <- 1000
d <- data.frame(x1=runif(n), x2=sample(letters[1:3], n, replace=TRUE))
head(d)
# x1 x2
# 1 0.2655087 b
# 2 0.3721239 c
# 3 0.5728534 b
# 4 0.9082078 c
# 5 0.2016819 a
# 6 0.8983897 a
X <- model.matrix(~x1*x2, d)
head(X)
# (Intercept) x1 x2b x2c x1:x2b x1:x2c
# 1 1 0.2655087 1 0 0.2655087 0.0000000
# 2 1 0.3721239 0 1 0.0000000 0.3721239
# 3 1 0.5728534 1 0 0.5728534 0.0000000
# 4 1 0.9082078 0 1 0.0000000 0.9082078
# 5 1 0.2016819 0 0 0.0000000 0.0000000
# 6 1 0.8983897 0 0 0.0000000 0.0000000
b <- rnorm(6) # coefficients
y <- X %*% b + rnorm(n, sd=0.1)
您可以将 X
的列重命名为您想要的任何名称,只要您在稍后 predict
将模型转换为新数据时使用一致的命名即可。
现在拟合模型。在这里,我告诉 lm
不要计算截距(使用 -1
),因为变量 (Intercept)
已经存在于 X
中,并且将为它计算一个系数。我们也可以通过拟合数据 as.data.frame(X[, -1])
:
来做到这一点
(M <- lm(y ~ . - 1, as.data.frame(X)))
# Call:
# lm(formula = y ~ . - 1, data = as.data.frame(X))
#
# Coefficients:
# `(Intercept)` x1 x2b x2c `x1:x2b` `x1:x2c`
# 1.14389 1.09168 -0.88879 0.20405 0.09085 -1.63769
创建一些新数据进行预测,并执行预测:
d.predict <- expand.grid(x1=seq(0, 1, 0.1), x2=letters[1:3])
X.predict <- model.matrix(~x1*x2, d.predict)
y.predict <- predict(M, as.data.frame(X.predict))
我遇到了一个看似简单但非常令人沮丧的问题。当您 运行 模型在 R 中具有交互项时,R 将生成的参数命名为 "var1:var2" 等。不幸的是,这种命名约定使我无法计算预测值和 CI 的新数据所在的位置必需,因为“:”不是可以包含在列 header 中的字符,并且原始数据框中的名称必须与 newdata 中的名称完全匹配。还有其他人遇到过这个问题吗?
这是我的代码示例:
wemedist2.exp = glm(survive/trials ~ sitedist + type + sitedist*type + roaddist, family = binomial(logexp(wemedata$expos)), data=wemedata)
summary(wemedist2.exp)
wemepredict3 = with(wemedata, data.frame(sitedist=mean(sitedist),roaddist=mean(roaddist), type=factor(1:2)))
wemepredict3 = cbind(wemepredict3, predict(wemedist2.exp, newdata = wemepredict3, type = "link", se = TRUE))
这会生成一个 table,其中包含每个变量在指定水平的预测值,但没有交互作用。
对于您的 newdata
数据框,您不应包含交互列。调用 predict
.
例如:
创建一些虚拟数据:
set.seed(1) n <- 10000 X <- data.frame(x1=runif(n), x2=runif(n)) X$x1x2 <- X$x1 * X$x2 head(X) # x1 x2 x1x2 # 1 0.2655087 0.06471249 0.017181728 # 2 0.3721239 0.67661240 0.251783646 # 3 0.5728534 0.73537169 0.421260147 # 4 0.9082078 0.11129967 0.101083225 # 5 0.2016819 0.04665462 0.009409393 # 6 0.8983897 0.13091031 0.117608474 b <- runif(4) y <- b[1] + c(as.matrix(X) %*% b[-1]) + rnorm(n, sd=0.1)
拟合模型并比较估计系数与真实系数:
M <- lm(y ~ x1 * x2, X) summary(M) # Call: # lm(formula = y ~ x1 * x2, data = X) # # Residuals: # Min 1Q Median 3Q Max # -0.43208 -0.06743 -0.00170 0.06601 0.37197 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 0.202040 0.003906 51.72 <2e-16 *** # x1 0.128237 0.006809 18.83 <2e-16 *** # x2 0.156942 0.006763 23.21 <2e-16 *** # x1:x2 0.292582 0.011773 24.85 <2e-16 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 0.09906 on 9996 degrees of freedom # Multiple R-squared: 0.5997, Adjusted R-squared: 0.5996 # F-statistic: 4992 on 3 and 9996 DF, p-value: < 2.2e-16 b # [1] 0.2106027 0.1147864 0.1453641 0.3099322
创建要预测的示例数据,并进行预测。请注意,我们只创建
x1
和x2
,并且 而不是 创建x1:x2
:X.predict <- data.frame(x1=runif(10), x2=runif(10)) head(X.predict) # x1 x2 # 1 0.26037592 0.7652155 # 2 0.73988333 0.3352932 # 3 0.02650689 0.9788743 # 4 0.84083874 0.1446228 # 5 0.85052685 0.7674547 # 6 0.13568509 0.9612156 predict(M, newdata=X.predict) # 1 2 3 4 5 6 7 # 0.4138194 0.4221251 0.3666572 0.3681432 0.6225354 0.4084543 0.4711018 # 8 9 10 # 0.7092744 0.3401867 0.2320834
或者...
另一种方法是通过计算交互项的乘积将交互包含在模型拟合数据中,然后将其也包含在新数据中。我们已经完成了上面第 1 点中的第一步,我们在其中创建了一个名为 x1x2
的列。
然后我们将拟合模型:lm(y ~ x1 + x2 + x1x2, X)
并对以下数据进行预测:
X.predict <- data.frame(x1=runif(10), x2=runif(10), x1x2=runif(10)
如果交互中涉及分类变量...
当您有涉及分类变量的交互时,模型会估计系数来描述属于每个级别相对于属于参考级别的影响。因此,例如,如果我们有一个连续预测变量 (x1
) 和一个分类预测变量 (x2
,水平为 a
、b
和 c
),则模型 y ~ x1 * x2
将估计六个系数,描述:
截距(即
x1
为零时的预测y
且观测属于x2
的参考水平);当观测值属于
x2
参考水平时改变x1
的效果(即斜率,对于x2
参考水平);属于二级的影响(即属于二级相对于属于参考级截距的变化);
属于第三级的影响(即由于属于第三级相对于属于参考级截距的变化);
x1
的效果因属于第二层级相对于属于参考层级的变化(即斜率的变化);和由于属于第三级相对于属于参考级
x1
的影响变化(即斜率变化)
如果您想拟合和预测模型 with/to 描述交互的预计算数据,您可以创建一个包含列的数据框:x1
; x2b
(二进制,表示观察是否属于级别b
); x2c
(二进制,表示观察是否属于级别c
); x1x2b
(x1
和 x2b
的乘积);和 x1x2c
(x1
和 x2c
的乘积)。
一个快速的方法是 model.matrix
:
set.seed(1)
n <- 1000
d <- data.frame(x1=runif(n), x2=sample(letters[1:3], n, replace=TRUE))
head(d)
# x1 x2
# 1 0.2655087 b
# 2 0.3721239 c
# 3 0.5728534 b
# 4 0.9082078 c
# 5 0.2016819 a
# 6 0.8983897 a
X <- model.matrix(~x1*x2, d)
head(X)
# (Intercept) x1 x2b x2c x1:x2b x1:x2c
# 1 1 0.2655087 1 0 0.2655087 0.0000000
# 2 1 0.3721239 0 1 0.0000000 0.3721239
# 3 1 0.5728534 1 0 0.5728534 0.0000000
# 4 1 0.9082078 0 1 0.0000000 0.9082078
# 5 1 0.2016819 0 0 0.0000000 0.0000000
# 6 1 0.8983897 0 0 0.0000000 0.0000000
b <- rnorm(6) # coefficients
y <- X %*% b + rnorm(n, sd=0.1)
您可以将 X
的列重命名为您想要的任何名称,只要您在稍后 predict
将模型转换为新数据时使用一致的命名即可。
现在拟合模型。在这里,我告诉 lm
不要计算截距(使用 -1
),因为变量 (Intercept)
已经存在于 X
中,并且将为它计算一个系数。我们也可以通过拟合数据 as.data.frame(X[, -1])
:
(M <- lm(y ~ . - 1, as.data.frame(X)))
# Call:
# lm(formula = y ~ . - 1, data = as.data.frame(X))
#
# Coefficients:
# `(Intercept)` x1 x2b x2c `x1:x2b` `x1:x2c`
# 1.14389 1.09168 -0.88879 0.20405 0.09085 -1.63769
创建一些新数据进行预测,并执行预测:
d.predict <- expand.grid(x1=seq(0, 1, 0.1), x2=letters[1:3])
X.predict <- model.matrix(~x1*x2, d.predict)
y.predict <- predict(M, as.data.frame(X.predict))