7

在@GavinSimpson 的 supperb 博客之后，我正在尝试适应 2015 年波兰地方政府选举结果。 https://www.fromthebottomoftheheap.net/2017/10/19/first-steps-with-mrf-smooths/ 我用 6 位标识符的 shp 数据加入 xls（可能有前导 0 s）。我将其保留为文本变量。编辑，我简化了标识符，现在使用从 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
4

1 回答 1

1

似乎我解决了这个问题，也许没有完全意识到它是如何成为 stat 初学者的。

首先，建模数据中不应存在单个 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()gam 公式更改为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.70GHZ 16GB Win10 的笔记本电脑上，现在需要 1-2 分钟来构建模型。

简而言之，我改变的是：删除单个 NA，从公式中删除偏移量并使用负二项分布。

这是我想要实现的结果，从左到右，实际无效投票率、模型平滑率和表示差异的残差。mcgv 代码让我这样做。

预期结果

于 2019-06-21T13:14:19.353 回答