使用循环，执行回归，预测每个人的新值。

Question

我有一个包含20个变量的数据集。其中10个变量是非常感兴趣的变量，但这些变量需要根据年龄和性别的群体差异进行调整。我用回归的方法来预测年龄和性别的数值。

有很多变量，很多人，所以我想要一个循环或者类似的。

下面是我尝试的一个例子

# Load example data
library(survival)
library(dplyr)
data(lung) # example data

# I want to obtain adjusted values for the following two variables, called "dependents"
dependents <- names(select(lung, 7:8))

new_data <- lung # copies data set
for (i in seq_along(dependents)) {
     eq <- paste(dependents[i],"~ age + sex")
     fit <- lm(as.formula(eq), data= new_data)
     new_data$predicted_value <- predict(fit, newdata=new_data, type='response')
     new_data <- rename(new_data, paste(dependents[i], "_predicted", sep="") = predicted_value)
}
View(new_data)

这未能为我提供经调整(即预测)形式的“受抚养人”。

有什么想法吗？

提前感谢

Answer 1

下面是另一种方法，使用我的扫帚包中的tidyr包和augment函数：

library(tidyr)
library(broom)

new_data <- lung %>%
  gather(dependent, value, ph.karno:pat.karno) %>%
  group_by(dependent) %>%
  do(augment(lm(value ~ age + sex, data = .)))

这将对数据进行重新组织，以便将每个依赖项(ph.karno和pat.karno)叠加在一起，并由dependent列来区分。augment函数将每个模型转换为一个数据框架，其中包含符合值、残差和其他您所关心的值的列(更多信息请参见?lm_tidiers )。然后，.fitted列给出拟合的值：

new_data
#> Source: local data frame [452 x 12]
#> Groups: dependent
#> 
#>    dependent .rownames value age sex  .fitted  .se.fit     .resid
#> 1   ph.karno         1    90  74   1 78.86709 1.406553  11.132915
#> 2   ph.karno         2    90  68   1 80.53347 1.115994   9.466530
#> 3   ph.karno         3    90  56   1 83.86624 1.226463   6.133759
#> 4   ph.karno         4    90  57   1 83.58851 1.181024   6.411490
#> 5   ph.karno         5   100  60   1 82.75532 1.078170  17.244683
#> 6   ph.karno         6    50  74   1 78.86709 1.406553 -28.867085
#> 7   ph.karno         7    70  68   2 80.18860 1.419744 -10.188596
#> 8   ph.karno         8    60  71   2 79.35540 1.555365 -19.355404
#> 9   ph.karno         9    70  53   1 84.69943 1.388600 -14.699433
#> 10  ph.karno        10    70  61   1 82.47759 1.056850 -12.477586
#> ..       ...       ...   ... ... ...      ...      ...        ...
#> Variables not shown: .hat (dbl), .sigma (dbl), .cooksd (dbl), .std.resid
#>   (dbl)

作为使用这些数据的一种方法，您可以绘制对因变量的预测是如何不同的：

ggplot(new_data, aes(age, .fitted, color = dependent, lty = factor(sex))) +
     geom_line()

但是，如果您希望控制年龄和性别，那么您可能希望使用.resid列。

Answer 2

还有第三种方法。

new_data <- na.omit(lung[,c("sex","age",dependents)])
result <- lapply(new_data[,dependents],
                 function(y)predict(lm(y~age+sex,data.frame(y=y,new_data[,c("age","sex")]))))
names(result) <- paste(names(result),"predicted",sep="_")
result <- cbind(new_data,as.data.frame(result))
head(result)
#   sex age ph.karno pat.karno ph.karno_predicted pat.karno_predicted
# 1   1  74       90       100           78.83030            77.34670
# 2   1  68       90        90           80.59974            78.53841
# 3   1  56       90        90           84.13862            80.92183
# 4   1  57       90        60           83.84371            80.72321
# 5   1  60      100        90           82.95899            80.12736
# 6   1  74       50        80           78.83030            77.34670

您的原始代码有一些微妙的问题(除了它没有运行的事实)。响应变量有几个NA，它们由lm(...)自动删除，因此预测比原始数据集有更少的行，并且当您尝试添加新列时，例如

 new_data$predicted_value <- predict(fit, newdata=new_data, type='response') 

你犯了个错误。您必须首先从NA中删除new_data s，如上面的代码所示。

我也想知道，既然你的数据似乎是计数，你是否应该使用泊松glm而不是lm？