GAM с mrf smooth - ошибки (несоответствие между именами областей nb / polys и именами областей данных - PullRequest
Новая
       45

GAM с mrf smooth - ошибки (несоответствие между именами областей nb / polys и именами областей данных

7 голосов
/ 11 июня 2019

Я пытаюсь соответствовать результатам выборов в местные органы власти Польши за 2015 год, опубликованным в блоге @GavinSimpson.https://www.fromthebottomoftheheap.net/2017/10/19/first-steps-with-mrf-smooths/ Я присоединяюсь к xls с данными shp по 6-значному идентификатору (могут быть первые 0 с).Я сохранил это текстовой переменной.РЕДАКТИРОВАТЬ, я упростил идентификатор и теперь использую последовательность от 1 до nrow, чтобы упростить мой вопрос.

library(tidyverse)
library(sf)
library(mgcv)

# Read data
# From https://www.gis-support.pl/downloads/gminy.zip shp file

boroughs_shp <- st_read("../../_mapy/gminy.shp",options = "ENCODING=WINDOWS-1250",
                     stringsAsFactors = FALSE ) %>% 
  st_transform(crs = 4326)%>% 
  janitor::clean_names() %>% 
# st_simplify(preserveTopology = T, dTolerance = 0.01) %>% 
  mutate(teryt=str_sub(jpt_kod_je, 1, 6)) %>% 
  select(teryt, nazwa=jpt_nazwa, geometry)

# From https://parlament2015.pkw.gov.pl/wyniki_zb/2015-gl-lis-gm.zip data file
elections_xls <-
  readxl::read_excel("data/2015-gl-lis-gm.xls",
             trim_ws = T, col_names = T) %>% 
  janitor::clean_names() %>% 
  select(teryt, liczba_wyborcow, glosy_niewazne)

elections <-
  boroughs_shp %>% fortify() %>% 
  left_join(elections_xls, by = "teryt") %>% 
  arrange(teryt) %>%
  mutate(idx = seq.int(nrow(.)) %>% as.factor(), 
         teryt = as.factor(teryt)) 

# Neighbors

boroughs_nb <-spdep::poly2nb(elections, snap = 0.01, queen = F, row.names = elections$idx )
names(boroughs_nb) <- attr(boroughs_nb, "region.id")

# Model

ctrl <- gam.control(nthreads = 4) 
m1 <- gam(glosy_niewazne ~ s(idx, bs = 'mrf', xt = list(nb = boroughs_nb)), 
          data = elections,
          offset = log(liczba_wyborcow), # number of votes
          method = 'REML', 
          control = ctrl,
          family = betar())

Вот сообщение об ошибке: 

    Error in smooth.construct.mrf.smooth.spec(object, dk$data, dk$knots) : 
  mismatch between nb/polys supplied area names and data area names
In addition: Warning message:
In if (all.equal(sort(a.name), sort(levels(k))) != TRUE) stop("mismatch between nb/polys supplied area names and data area names") :
  the condition has length > 1 and only the first element will be used

выборы $ idx является фактором.Я использую его, чтобы дать имена boroughs_nb, чтобы быть абсолютно уверенным, что у меня одинаковое количество уровней.Что я делаю не так?

РЕДАКТИРОВАТЬ: Условие, указанное в сообщении об ошибке, выполняется: 

> all(sort(names(boroughs_nb)) == sort(levels(elections$idx)))
[1] TRUE

1 Ответ

0 голосов
/ 21 июня 2019

Кажется, что я решил проблему, возможно, не совсем понимая, как она поступила со статусом новичка.

Во-первых, ни один NA не должен присутствовать в смоделированных данных.Был один.После этого mcgv, казалось, запускался, но это заняло много времени (четверть часа) и необъяснимо для меня, только когда я ограничил число узлов k=50, с плохими результатами (меньше или больше, и он не дал никакого результата)и с предупреждением быть осторожным о результатах.Затем я попытался удалить offset=log(liczba_wyborcow), то есть сместить число избирателей, и набрал число недействительных голосов на 1000 мою прогнозируемую переменную. 

elections <-
 boroughs_shp %>%  
 left_join(elections_xls, by = "teryt") %>% na.omit() %>% 
 arrange(teryt) %>% 
 mutate(idx = row_number() %>% as.factor()) %>% 
 mutate(void_ratio=round(glosy_niewazne/liczba_wyborcow,3)*1000)

Теперь, когда это подсчет, почему бы не попробоватьизменить family = betar() в формуле гаммы на poisson() - все еще не очень хороший результат, а затем на отрицательный бином family = nb() Теперь моя формула выглядит как 

m1 <-
gam(
 void_ratio ~ s(
 idx,
 bs = 'mrf',
 k =500,
 xt = list(nb = boroughs_nb),
 fx = TRUE),
 data = elections_df,
 method = 'REML', 
 control = gam.control(nthreads = 4),
 family = nb()
)

Кажется, теперьбыть невероятно быстрым и возвращать достоверные результаты без предупреждений или ошибок.На ноутбуке с 4 ядрами Intel Core I7 6820HQ @ 2,70 ГГц 16 ГБ Win10 сейчас требуется 1-2 минуты, чтобы собрать модель.

Вкратце, я изменил следующее: удалите один NA, удалите смещение из формулы и используйте отрицательное биномиальное распределение .

Вот результат того, что я хотел достичь, слева направо, реальную ставку недействительных голосов, ставку, сглаженную моделью, и остатки, указывающие на расхождения.Код mcgv позволяет мне это сделать.

expected result

...