Как предсказать марковскую модель случайного поля гаммы в сетку?

Question

Я работаю над аналогичной проблемой, как в примере в пакете mgcv.

Я создал аналогичную модель и хотел бы прогнозировать до уровня сетки вместо районов, указанных в данных обучения. Это возможно? Если да, то как?

В качестве примера я создал набор данных прогнозирования.

library(mgcv)
data(columb)
data(columb.polys)

b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")
plot(b,scheme=1)


#Create prediction dataset

df <- data.frame(district=numeric(0),x=numeric(0),y= numeric(0)) #Create empty df to store x, y and IDs for each polygon

# Extract x and y coordinates from each polygon and assign district ID
for (i in 1:length(columb.polys)) {
  district <- i-1
  x <- columb.polys[[i]][,1]
  y <- columb.polys[[i]][,2]
  df <- rbind(df,cbind(district,x,y)) 
}

sp <- df %>%
  group_by(district) %>%
  do(poly=dplyr::select(., x, y) %>%Polygon()) %>%
  rowwise() %>%
  do(polys=Polygons(list(.$poly),.$district)) %>%
  {SpatialPolygons(.$polys)}


#Convert sp to sf object
sp_sf <- st_as_sf(sp, proj4string='+proj=utm +zone=48N +ellps=WGS84 +units=m')
sp_sf <- st_set_crs(sp_sf, '+proj=utm +zone=48N +ellps=WGS84 +units=m')

#Make grid
grid <- st_make_grid(sp_sf, 
                           cellsize = 0.05, 
                              what = "centers", crs = '+proj=utm +zone=48N +ellps=WGS84 +units=m')

#intersection of grid and Columbus, Ohio
pred_grid <- st_intersection(grid, st_buffer(sp_sf, dist=0))

pred_b <- predict(b, pred_data)

Error in data[[txt]] : subscript out of bounds

Спасибо!

Answer 1

Это невозможно с MRF, по крайней мере, насколько мне известно. В модели, которую вы показываете, у вас есть только сглаживание MRF, поэтому условное распределение отклика моделируется как оценочные значения для каждогоdistrictпри условии, что близлежащие районы похожи друг на друга. Чтобы использовать модель, вам нужно сделать прогноз для районов, а не для каких-то других данных.

Я предполагаю в этой строке

pred_b <- predict(b, pred_data)

Вы имели в виду

pred_b <- predict(b, pred_grid)

Если ваша реальная проблема представляет собой смесь данных в сетке и на уровне региона, вы можете предсказать по сетке, но вам нужно будет назначить каждое местоположение сетки региону, чтобы вы могли передать в столбце, например District что входит в термин MRF.