当因子变量包含在模型中时，栅格预测不会在会话之间重现答案

【问题标题】：Raster predictions do not reproduce session to session when factor variable is included in model当因子变量包含在模型中时，栅格预测不会在会话之间重现
【发布时间】：2018-10-08 20:12:53
【问题描述】：

这个问题与我一年半前发布的这个问题有关：Reproducibility of results from predict() function - raster package。但由于它没有示例，因此我创建了一个新问题，并提供了更新的信息。

在将我的预测复制到栅格时，我遇到了一个有点模糊的问题。我正在创建一个带有数值变量和单因素变量的 gbm 模型。然后，我使用 raster 包使用我训练的模型预测到栅格。预测因会话而异，但在单个 R 会话中重现。如果我删除因子变量，结果会重现会话到会话。此外，在下面的示例中，如果训练数据中的因子级别比栅格变量版本中的多，我可以让它重现会话到会话。是什么原因造成的？如何在包含因子变量的同时重现我的结果会话？

# This code will not reproduce session to session, but does if I leave many many factor levels in newwine with the
# commented out code

library(breakDown)
library(gbm)
library(dplyr)
library(raster)

# leave in many levels and code will reproduce session to session
#newwine <- wine[1:500,c(1:3,6)]

# specify only levels which are in the below raster and code will not reproduce session to session
newwine <- wine[,c(1:3,6)] %>% 
           filter(free.sulfur.dioxide == 3 | free.sulfur.dioxide == 10 | free.sulfur.dioxide == 15 |
                  free.sulfur.dioxide == 37 | free.sulfur.dioxide == 76)

head(newwine)

# make free.sulfur.dioxide as factor variable
newwine$free.sulfur.dioxide <- as.factor(newwine$free.sulfur.dioxide)
levels(newwine$free.sulfur.dioxide)

set.seed(123)
model <- gbm(fixed.acidity ~ ., data = newwine, 
             distribution = "gaussian",
             bag.fraction = 0.50,
             n.trees = 1000, 
             interaction.depth = 16, 
             shrinkage = 0.016, 
             n.minobsinnode = 10, verbose = FALSE)

summary(model)
plot(model, i.var = 3, n.trees = 1000)


# make some rasters for the predictor variables
free.sulfur.dioxide <- c(rep(3,times=10), rep(10, times = 10), 
                         rep(15, times = 10), rep(37, times = 10), 
                         rep(76, times = 10))

free.sulfur.dioxide.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(free.sulfur.dioxide.r) <- free.sulfur.dioxide

set.seed(123)
volatile.acidity <- newwine %>% 
                    dplyr::select(volatile.acidity) %>% 
                    sample_n(50) 
volatile.acidity <- as.vector(volatile.acidity)[,1]
volatile.acidity.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(volatile.acidity.r) <- volatile.acidity

set.seed(123)
citric.acid <- newwine %>% 
               dplyr::select(citric.acid) %>% 
               sample_n(50) 
citric.acid <- as.vector(citric.acid)[,1]
citric.acid.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(citric.acid.r) <- citric.acid

# create a raster stack
r <- stack(free.sulfur.dioxide.r, volatile.acidity.r, citric.acid.r)
names(r) <- c("free.sulfur.dioxide", "volatile.acidity", "citric.acid")

###########################################################################################################################

# predict to a raster with raster predict
pred <- predict(r, model, n.trees = model$n.trees, format="GTiff")
writeRaster(pred, "prediction1.tif", overwrite = TRUE)

###########################################################################################################################

# close the session and reopen, run until line 61, then run below to make a new prediction, called prediction 2
pred <- predict(r, model, n.trees = model$n.trees, format="GTiff")
writeRaster(pred, "prediction2.tif", overwrite = TRUE)

# read in the previous prediction
prediction1 <- raster("prediction1.tif")
prediction2 <- raster("prediction2.tif")

# compare rasters built across sessions
compareRaster(prediction1, prediction2, values = TRUE)
summary(prediction1-prediction2)

# compare rasters built within same session
pred2 <- predict(r, model, n.trees = model$n.trees, format="GTiff")
compareRaster(pred, pred2, values = TRUE)

但是，下面的代码不使用因子变量，并且会在会话之间重现。

### Same exercise but without setting the free sulfur dioxide to factor
## this code will reproduce session to session

library(breakDown)
library(gbm)
library(dplyr)
library(raster)

newwine <- wine[1:500,c(1:3)]

head(newwine)

set.seed(123)
model <- gbm(fixed.acidity ~ ., data = newwine, 
             distribution = "gaussian",
             bag.fraction = 0.50,
             n.trees = 1000, 
             interaction.depth = 16, 
             shrinkage = 0.016, 
             n.minobsinnode = 10, verbose = FALSE)

summary(model)

set.seed(123)
volatile.acidity <- newwine %>% 
  dplyr::select(volatile.acidity) %>% 
  sample_n(50) 
volatile.acidity <- as.vector(volatile.acidity)[,1]
volatile.acidity.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(volatile.acidity.r) <- volatile.acidity

set.seed(123)
citric.acid <- newwine %>% 
  dplyr::select(citric.acid) %>% 
  sample_n(50) 
citric.acid <- as.vector(citric.acid)[,1]
citric.acid.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(citric.acid.r) <- citric.acid

# create a raster stack
r <- stack( volatile.acidity.r, citric.acid.r)
names(r) <- c( "volatile.acidity", "citric.acid")

#######################################################################################################################

# predict to a raster with raster predict
pred <- predict(r, model, n.trees = model$n.trees, format="GTiff")
writeRaster(pred, "prediction1.tif", overwrite = TRUE)

#######################################################################################################################

# close the session and reopen to make a new prediction, called prediction 2
pred <- predict(r, model, n.trees = model$n.trees, format="GTiff")
writeRaster(pred, "prediction2.tif", overwrite = TRUE)

# read in the previous prediction
prediction1 <- raster("prediction1.tif")
prediction2 <- raster("prediction2.tif")

# compare rasters built across sessions
compareRaster(prediction1, prediction2, values = TRUE)
summary(prediction1-prediction2)

# compare rasters built within same session
pred2 <- predict(r, model, n.trees = model$n.trees, format="GTiff")
compareRaster(pred, pred2, values = TRUE)
summary(pred-pred2)

【问题讨论】：

你能用代码生成一些数据吗？如果没有示例数据，很难提供答案。展示结果的不同也是很好的。
您好，数据是内置在breakDown包中的。当您加载库时，“wine”数据集可用，这就是本示例所使用的。
问题可能与此处接受的答案有关：stackoverflow.com/questions/25121725/…。如果我使用 gbm.fit() 拟合模型并明确设置 x 和 y 而不是使用 gbm() 中的公式，我可以在 R 会话中重现结果。
好的，谢谢。我会尝试看看它（尽管我认为永远不应该保存会话）。
我不保存会话，而是进行光栅预测，关闭会话（不保存），然后重新开始进行第二次预测。然后，我通过重新读取两个栅格来比较两个预测。它们应该完全相同，但略有偏差。我将差异与summary(prediction1 - prediction2) 进行比较。

标签： r session prediction r-raster gbm

【解决方案1】：

看来这个问题不是由raster 包引起的，而是由gbm 包引起的。经过一番挖掘，我发现 gbm 包在 2017 年 3 月被孤立，并且有一个新的 gbm 包，在 github 上称为 gbm3（在 CRAN 上尚不可用）https://github.com/gbm-developers/gbm3。当您对栅格进行预测时，您可以使用模型类型要求的任何预测方法（例如，predict.gbm() 用于gbm，predict.GBMFit() 用于gbm3。似乎predict.gbm() 只是不处理来自模型中栅格的因素正确。它可能是一个错误，也可能不是一个错误，但在任何一种情况下，这个包都不再被维护。gbm3 可以解决问题并且是可重现的。

# This code will reproduce session to session for the gbm3 model, but not for old gbm model

library(breakDown)
# install gbm3 from github
library(gbm3)
library(dplyr)
library(raster)

# specify only levels which are in the below raster 
newwine <- wine[,c(1:3,6)] %>% 
           filter(free.sulfur.dioxide == 3 | free.sulfur.dioxide == 10 | free.sulfur.dioxide == 15 |
                  free.sulfur.dioxide == 37 | free.sulfur.dioxide == 76)

head(newwine)

# make free.sulfur.dioxide as factor variable
newwine$free.sulfur.dioxide <- as.factor(newwine$free.sulfur.dioxide)
levels(newwine$free.sulfur.dioxide)

#set.seed(123)
# model <-  gbm(fixed.acidity ~ ., data = newwine, #gbm.fit(x = newwine[,2:4], y = newwine[,1], 
#              distribution = "gaussian",
#              bag.fraction = 0.50,
#              n.trees = 1000, 
#              interaction.depth = 16, 
#              shrinkage = 0.016, 
#              n.minobsinnode = 10, verbose = FALSE)
set.seed(123)
model <- gbmt(fixed.acidity ~ ., data = newwine, distribution = gbm_dist("Gaussian")) 

summary(model)
plot(model, var_index = 3, num_trees = 1000)

# make some rasters for the predictor variables
free.sulfur.dioxide <- c(rep(3,times=10), rep(10, times = 10), 
                         rep(15, times = 10), rep(37, times = 10), 
                         rep(76, times = 10))

free.sulfur.dioxide.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(free.sulfur.dioxide.r) <- free.sulfur.dioxide

set.seed(123)
volatile.acidity <- newwine %>% 
                    dplyr::select(volatile.acidity) %>% 
                    sample_n(50) 
volatile.acidity <- as.vector(volatile.acidity)[,1]
volatile.acidity.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(volatile.acidity.r) <- volatile.acidity

set.seed(123)
citric.acid <- newwine %>% 
               dplyr::select(citric.acid) %>% 
               sample_n(50) 
citric.acid <- as.vector(citric.acid)[,1]
citric.acid.r <- raster(ext = extent(-10, 5, -10, 5), nrows = 5, ncols = 10)
values(citric.acid.r) <- citric.acid

# create a raster stack
r <- stack(free.sulfur.dioxide.r, volatile.acidity.r, citric.acid.r)
names(r) <- c("free.sulfur.dioxide", "volatile.acidity", "citric.acid")

###########################################################################################################################

# predict to a raster with raster predict
pred <- raster::predict(r, model, n.trees = 2000, format="GTiff")
writeRaster(pred, "prediction1.tif", overwrite = TRUE)

# predict to a vector with predict
v <- values(r)
v <- data.frame(v)
v$free.sulfur.dioxide <- as.factor(v$free.sulfur.dioxide)
vpred <- predict(model, v, n.trees = 2000)
write.table(vpred, "vector_predict.txt", row.names = FALSE, col.names = TRUE)

###########################################################################################################################

# close the session and reopen, run until #### line, then run below to make a new prediction, called prediction 2
pred <- raster::predict(r, model, n.trees = 2000, format="GTiff")
writeRaster(pred, "prediction2.tif", overwrite = TRUE)

# predict to a vector with predict
v <- values(r)
v <- data.frame(v)
v$free.sulfur.dioxide <- as.factor(v$free.sulfur.dioxide)
vpred <- predict(model, v, n.trees = 2000)
write.table(vpred, "vector_predict2.txt", row.names = FALSE, col.names = TRUE)

# read in the previous prediction
prediction1 <- raster("prediction1.tif")
prediction2 <- raster("prediction2.tif")

# compare rasters built across sessions
compareRaster(prediction1, prediction2, values = TRUE)
summary(prediction1-prediction2)

# compare rasters built within same session
pred2 <- raster::predict(r, model, n.trees = 2000, format="GTiff", factors = f)
compareRaster(pred, pred2, values = TRUE)

# compare the vector predictions
p1 <- read.delim("vector_predict.txt")
p2 <- read.delim("vector_predict2.txt")

plot(p1$x,p2$x)

summary(p1$x - p2$x)

【讨论】：

【解决方案2】：

这不是解决方案，而是试图解决问题。在我看来，这与raster 无关。

当我这样做时：

v <- values(r)
pred <- predict(model, data.frame(v), n.trees = model$n.trees)
rpred <- predict(r, model, n.trees = model$n.trees)

退出，保存会话，开始新会话并执行以下操作：

library(gbm)
library(raster)
pred2 <- predict(model, data.frame(v), n.trees = model$n.trees )
rpred2 <- predict(r, model, n.trees = model$n.trees)

我看到pred 和pred2 的值不太一样。（请参阅plot(pred, pred2)。但是，pred2 和 rpred2 的值是相同的：plot(values(rpred2), pred2)。

或者，当我保存pred (saveRDS(pred, 'pred.rds')，并将其加载到新会话pred1 <- readRDS(pred.rds) 时，结果并不完全相同。

这表明gbm 的某处发生了一些不受set.seed 控制的随机化。

【讨论】：

你的意思是rpred和rpred2的值是一样的吗？
不，它们不是，但它们与 pred 和 pred2 相同。也就是说，似乎差异源于gbm.predict；而不是在 raster::predict 中。
啊，好吧，我明白了。但是，如果我使用tpred <- predict(model, newwine, n.trees = model$n.trees 预测训练数据，则每次我在新会话中预测时结果都是相同的。 values() 是否对因子变量进行任何转换？
class(data.frame(v)$free.sulfur.dioxide) 返回数字。
如果我这样做：v <- values(r) v <- data.frame(v) v$free.sulfur.dioxide <- as.factor(v$free.sulfur.dioxide) vpred <- predict(model, v, n.trees = model$n.trees) 它将重现会话到会话。