问题描述
这个问题是上一个问题 (Linear Regression prediction in R using Leave One out Approach) 的第二部分。
我正在尝试为每个国家/地区构建模型并使用留一法生成线性回归预测。换句话说，在下面的代码中，在构建模型 1 和模型 2 时，使用的“数据”不应是整个数据集。相反，它应该是数据集（国家/地区）的子集。每个国家/地区的数据都应使用基于该国家/地区特定数据构建的模型进行评估。
下面的代码返回一个错误。我如何修改/修复下面的代码来做到这一点？或者有更好的方法吗？
library(modelr)
install.packages("gapminder")
library(gapminder)
data(gapminder)
#CASE 1
model1 <- lm(lifeExp ~ pop,data = gapminder,subset = country)
model2 <- lm(lifeExp ~ pop + gdpPercap,subset = country)
models <- list(fit_model1 = model1,fit_model2 = model2)
gapminder %>% nest_by(continent,country) %>%
bind_cols(
map(1:nrow(gapminder),function(i) {
map_dfc(models,function(model) {
training <- data[-i,]
fit <- lm(model,data = training)
validation <- data[i,]
predict(fit,newdata = validation)
})
}) %>%
bind_rows()
)
解决方法
最简洁直接的解决方案是嵌套 for 循环方法，其中外循环是两个模型公式，而内循环是我们想要省略的统一体。这也可以使用 outer 来完成，我也会在后面展示。
为了清楚起见，我首先展示了如何在每次迭代中省略一个观察（即一行）（第一部分）。我稍后将展示如何省略一个集群（例如国家/地区）（第二部分）。我还使用了内置的 iris 数据集，它更小，因此更容易处理。它包含一个 "Species" 列，用于对应您数据中的 "countries"。
第一部分
首先，我们将两个公式放入一个列表中，并按照我们希望它们稍后出现在结果列中的方式命名它们。
FOAE <- list(fit1=Petal.Length ~ Sepal.Length,fit2=Petal.Length ~ Sepal.Length + Petal.Width)
对于循环，我们要初始化一个矩阵 im，其行对应于我们要省略的行数，列对应于模型公式的数量。
im <- matrix(NA,nrow=nrow(iris),ncol=length(FOAE),dimnames=list(NULL,names(FOAE)))
这看起来像这样：
head(im,n=3)
# fit1 fit2
# [1,] NA NA
# [2,] NA NA
# [3,] NA NA
现在我们如上所述循环公式和行。
for (i in seq(FOAE)) {
for(j in seq(nrow(iris))) {
train <- iris[-j,]
test <- iris[j,]
fit <- lm(FOAE[[i]],data=train)
im[j,i] <- predict(fit,newdata=test)
}
}
im 现在已经被填充了，我们可以将它cbind到原来的iris数据集来得到我们的结果res1。
res1 <- cbind(iris,im)
head(res1)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit1 fit2
# 1 5.1 3.5 1.4 0.2 setosa 2.388501 1.611976
# 2 4.9 3.0 1.4 0.2 setosa 2.014324 1.501389
# 3 4.7 3.2 1.3 0.2 setosa 1.639805 1.392955
# 4 4.6 3.1 1.5 0.2 setosa 1.446175 1.333199
# 5 5.0 3.6 1.4 0.2 setosa 2.201646 1.556620
# 6 5.4 3.9 1.7 0.4 setosa 2.944788 2.127184
为了替代 outer 方法，我们将 for 循环内的代码放入 formula 的 Vectorize 中，以便它可以处理矩阵列（即向量）。
FUN1 <- Vectorize(function(x,y) {
train <- iris[-x,]
test <- iris[x,]
fit <- lm(y,data=train)
predict(fit,newdata=test)
})
现在我们将 FOAE 和行 1:nrow(iris) 随后省略，连同 FUN1 放入 outer()。这已经为我们提供了结果，我们可以以与上述相同的方式 cbind 到 iris 得到我们的结果 res2。
o1 <- outer(FOAE,1:nrow(iris),FUN1)
res2 <- cbind(iris,o1)
head(res2)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit1 fit2
# 1 5.1 3.5 1.4 0.2 setosa 2.388501 1.611976
# 2 4.9 3.0 1.4 0.2 setosa 2.014324 1.501389
# 3 4.7 3.2 1.3 0.2 setosa 1.639805 1.392955
# 4 4.6 3.1 1.5 0.2 setosa 1.446175 1.333199
# 5 5.0 3.6 1.4 0.2 setosa 2.201646 1.556620
# 6 5.4 3.9 1.7 0.4 setosa 2.944788 2.127184
## test if results are different is negative
stopifnot(all.equal(res1,res2))
第二部分
我们可能会在遗漏集群（即物种或国家/地区）时采用类似的方法。我在这里展示了 outer 方法。我们想要改变的是，我们现在想要忽略属于特定集群的观察，这里 "Species"（在您的情况下为 "countries"），我们将这些 unique 值放入向量中Species.u 。由于值采用 "character" 或 "factor" 格式，我们使用 data[!data$cluster %in% x,] 而不是 data[-x,] 对数据进行子集化。因为 predict 会在集群中产生多个值，但我们希望在各自的集群中具有相同的值，所以我们可能想要使用统计数据，例如每个集群的 mean 预测。我们根据集群使用 rownames。
FUN2 <- Vectorize(function(x,y) {
train <- iris[!iris$Species %in% x,]
test <- iris[iris$Species %in% x,data=train)
mean(predict(fit,newdata=test))
})
Species.u <- unique(iris$Species)
o2 <- `rownames<-`(outer(Species.u,FOAE,FUN2),Species.u)
这现在给了我们一个小于我们的数据集的矩阵。多亏了 rownames，我们可以match 预测它们所属的集群。
o2
# fit1 fit2
# setosa 3.609943 2.662609
# versicolor 3.785760 3.909919
# virginica 4.911009 5.976922
res3 <- cbind(iris,o2[match(iris$Species,rownames(o2)),])
head(res3)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit1 fit2
# setosa 5.1 3.5 1.4 0.2 setosa 3.609943 2.662609
# setosa.1 4.9 3.0 1.4 0.2 setosa 3.609943 2.662609
# setosa.2 4.7 3.2 1.3 0.2 setosa 3.609943 2.662609
# setosa.3 4.6 3.1 1.5 0.2 setosa 3.609943 2.662609
# setosa.4 5.0 3.6 1.4 0.2 setosa 3.609943 2.662609
# setosa.5 5.4 3.9 1.7 0.4 setosa 3.609943 2.662609
tail(res3)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit1 fit2
# virginica.44 6.7 3.3 5.7 2.5 virginica 4.911009 5.976922
# virginica.45 6.7 3.0 5.2 2.3 virginica 4.911009 5.976922
# virginica.46 6.3 2.5 5.0 1.9 virginica 4.911009 5.976922
# virginica.47 6.5 3.0 5.2 2.0 virginica 4.911009 5.976922
# virginica.48 6.2 3.4 5.4 2.3 virginica 4.911009 5.976922
# virginica.49 5.9 3.0 5.1 1.8 virginica 4.911009 5.976922
编辑
在这个版本的 FUN2,FUN3 中，每个集群的模型的输出都经过rbinded（当然是两列，因为有两个模型）。
FUN3 <- Vectorize(function(x,data=train)
(predict(fit,newdata=test))
},SIMPLIFY=F)
Species.u <- unique(iris$Species)
o3 <- `rownames<-`(outer(Species.u,FUN3),Species.u)
res32 <- cbind(iris,apply(o3,2,unlist))
head(res32)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit1 fit2
# setosa.1 5.1 3.5 1.4 0.2 setosa 3.706940 2.678255
# setosa.2 4.9 3.0 1.4 0.2 setosa 3.500562 2.547587
# setosa.3 4.7 3.2 1.3 0.2 setosa 3.294183 2.416919
# setosa.4 4.6 3.1 1.5 0.2 setosa 3.190994 2.351586
# setosa.5 5.0 3.6 1.4 0.2 setosa 3.603751 2.612921
# setosa.6 5.4 3.9 1.7 0.4 setosa 4.016508 3.073249
编辑 2
正如我在您的评论中了解到的，您需要 1. 沿着集群的数据子集。这将是下面 ss 中的 FUN4。然后 ss 也通过在子集 z 的行上省略一行 ss 进行子集化。
FUN4 <- Vectorize(function(x,y) {
## subsets first by cluster then by row
ss <- iris[iris$Species %in% x,] ## cluster subset
sapply(1:nrow(ss),function(z) { ## subset rows using `sapply`
train <- ss[-z,] ## train data w/o row z
test <- ss[z,] ## test data for `predict`,just row z
fit <- lm(y,data=train)
predict(fit,newdata=test)
})
},SIMPLIFY=F)
## the two models
FOAE <- list(fit1=Petal.Length ~ Sepal.Length,fit2=Petal.Length ~ Sepal.Length + Petal.Width)
## unique cluster names
Species.u <- unique(iris$Species)
## with the `outer` we iterate over all the permutations of clusters and models `FOAE`.
o4 <- `rownames<-`(outer(Species.u,FUN4),Species.u)
## `unlist`ed result is directly `cbind`able to original data
res4 <- cbind(iris,apply(o4,unlist))
## result
head(res4)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit1 fit2
# setosa.1 5.1 3.5 1.4 0.2 setosa 1.476004 1.451029
# setosa.2 4.9 3.0 1.4 0.2 setosa 1.449120 1.431737
# setosa.3 4.7 3.2 1.3 0.2 setosa 1.426185 1.416492
# setosa.4 4.6 3.1 1.5 0.2 setosa 1.404040 1.398103
# setosa.5 5.0 3.6 1.4 0.2 setosa 1.462460 1.441295
# setosa.6 5.4 3.9 1.7 0.4 setosa 1.504990 1.559045