线性回归预测中的个别项

Individual terms in prediction of linear regression

我在 R 中对某些数据集进行了回归分析,并尝试预测每个独立变量对数据集中每一行的因变量的贡献。

所以像这样:

set.seed(123)                                              
y <- rnorm(10)                                           
m <- data.frame(v1=rnorm(10), v2=rnorm(10), v3=rnorm(10))
regr <- lm(formula=y~v1+v2+v3, data=m)  
summary(regr)
terms <- predict.lm(regr,m, type="terms")

简而言之:运行 回归并使用预测函数计算数据集 m 中 v1、v2 和 v3 的项。但是我很难理解预测函数在计算什么。我希望它将回归结果的系数与可变数据相乘。所以像这样的 v1:

coefficients(regr)[2]*m$v1

但是与预测函数相比,这给出了不同的结果。

自己计算:

0.55293884  0.16253411  0.18103537  0.04999729 -0.25108302  0.80717945  0.22488764 -0.88835486  0.31681455 -0.21356803

和预测函数计算:

0.45870070  0.06829597  0.08679724 -0.04424084 -0.34532115  0.71294132  0.13064950 -0.98259299  0.22257641 -0.30780616

prediciton 函数是 0.1 左右此外,如果您将预测函数中的所有项与常量相加,它不会加起来成为总预测(使用 type=”response”)。预测函数在这里计算什么,我如何告诉它计算我用系数(regr)[2]*m$v1 做了什么?

以下所有行都会产生相同的预测:

# our computed predictions
coefficients(regr)[1] + coefficients(regr)[2]*m$v1 +
  coefficients(regr)[3]*m$v2 + coefficients(regr)[4]*m$v3

# prediction using predict function
predict.lm(regr,m)

# prediction using terms matrix, note that we have to add the constant.
terms_predict = predict.lm(regr,m, type="terms")
terms_predict[,1]+terms_predict[,2]+terms_predict[,3]+attr(terms_predict,'constant')

您可以阅读有关使用 type="terms" here 的更多信息。

你自己的计算(coefficients(regr)[2]*m$v1)和predict函数计算(terms_predict[,1])不同的原因是因为terms矩阵中的列以均值为中心,所以它们的均值变成零:

# this is equal to terms_predict[,1]
coefficients(regr)[2]*m$v1-mean(coefficients(regr)[2]*m$v1)

# indeed, all columns are centered; i.e. have a mean of 0.
round(sapply(as.data.frame(terms_predict),mean),10)

希望对您有所帮助。

函数 predict(...,type="terms") 以每个变量的均值为中心。因此,输出有点难以解释。这是一个替代方案,其中每个变量(constantx1x2)都乘以它的系数。

TLDR:pred_terms <- model.matrix(formula(mod$terms), testData) %*% diag(coef(mod))


library(tidyverse)

### simulate data

set.seed(123)

nobs <- 50

x1 <- cumsum(rnorm(nobs) + 3)
x2 <- cumsum(rnorm(nobs) * 3)

y <- 2 + 2*x1 -0.5*x2 + rnorm(nobs,0,50)

df <- data.frame(t=1:nobs, y=y, x1=x1, x2=x2)

train <- 1:round(0.7*nobs,0)

rm(x1, x2, y)

trainData <- df[train,]
testData <- df[-train,]

### linear model

mod <- lm(y ~ x1 + x2 , data=trainData)

summary(mod)


### predict test set

test_preds <- predict(mod, newdata=testData)

head(test_preds)

### contribution by predictor

test_contribution <- model.matrix(formula(mod$terms), testData) %*% diag(coef(mod))

colnames(test_contribution) <- names(coef(mod))

head(test_contribution)

all(round(apply(test_contribution, 1, sum),5) == round(test_preds,5))  ## should be true

### Visualize each contribution

test_contribution_df <- as.data.frame(test_contribution)
test_contribution_df$pred <- test_preds
test_contribution_df$t <- row.names(test_contribution_df)
test_contribution_df$actual <- df[-train,"y"]

test_contribution_df_long <- pivot_longer(test_contribution_df, -t, names_to="variable")

names(test_contribution_df_long)

ggplot(test_contribution_df_long, aes(x=t, y=value, group=variable, color=variable)) +
  geom_line() +
  theme_bw()