Растровые прогнозы не воспроизводят сеанс за сеансом, когда фактор-переменная включена в модель

Question

Этот вопрос связан с вопросом, который я опубликовал полтора года назад: Воспроизводимость результатов функции предиката () - растровый пакет . Но так как у него не было примера, я создал новый вопрос также с обновленной информацией.

У меня несколько неясная проблема с воспроизведением моих прогнозов в растре. Я создаю модель gbm с числовыми переменными и однофакторной переменной. Затем я использую растровый пакет, чтобы прогнозировать растр, используя мою обученную модель. Прогнозы варьируются от сеанса к сеансу, но воспроизводятся в течение одного сеанса 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)

Answer 1

Похоже, что эта проблема не из-за пакета raster, а из-за пакета gbm. После некоторых копаний я обнаружил, что пакет gbm потерян в марте 2017 года, и есть новый пакет gbm, который называется gbm3 на github (еще не доступен в 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)

Answer 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.