R

1 什么是时空扫描统计？ 扫描统计学主要用于观测多个位置上时间序列数据或事件数据,检测数据集中的异常聚集， 在scanstatistics包(Benjamin Allévius 2023)主要能实现以下功能： 首先将研究区域划分为多个位置(位置点或网格)，观测每个位置上一段时间内的数据变化或事件发生情况。 然后从现在向过去延伸，在空间上选取若干位置，在时间上选取一定时间段,形成空间-时间聚集。每个聚集对应选定的位置集合和时间间隔。 对每个空间-时间聚集计算一个统计量,反映该聚集内观测到的数据或事件与预期值的差异。 报告统计量最大的几个空间-时间聚集，它们可能对应数据或事件的异常聚集。 重复这个过程,对不同的空间-时间聚集进行检验，找到研究区域内所有潜在的异常聚集。 # 使用扫描空间统计的常见步骤有： 数据资料准备。 如事件数据、时间序列数据等。数据要包含空间位置信息和测量值或事件计数信息。 设置研究区域和空间格网。根据数据空间范围设置研究区域，并在区域内划分网格。网格形状和大小根据具体研究目标确定，通常情况使用行政区划的地图。 选取空间聚集的时间窗口。从当前时间向过去延伸，选择一个或几个时间窗口，如最近1周、2周等。 在研究区域内，选择若干空间窗口和时间窗口组合，形成空间-时间聚集。每个聚集对应一组位置(空间窗口)和一段时间区间(时间窗口)。 计算每个空间-时间聚集的观测值与预期值的差异，得到聚集的统计量。统计量可直接使用观测值与预期值的差异，或基于该差异计算s标准华的差异比值等。 识别统计量最大的k个空间-时间聚集，作为潜在的异常聚集中心。k的选择根据研究目标确定，通常为3到5个。 对k个异常聚集中心，进行进一步分析，如数据动态监测、空间相关探索等，判断聚集的显著性和可能原因。 重复步骤4到步骤7，调整空间-时间聚集的设置，完成对整个研究区域的扫描检验，找到所有潜在异常聚集。 根据扫描检验结果和进一步分析，得出研究结论，识别空间聚集模式以及影响因素。 # 示例：重庆市输入性疟疾时空聚集分析 # scanstatistics包中的主要函数 ## 扫描统计函数 scan_eb_poisson：最常用的函数， 基于Poisson分布模型期望值的扫描统计检测聚集时间(M. Kulldorff et al. 1998)。 scan_pb_poisson: 基于总体值(M. Kulldorff 2010). scan_eb_negbin： 基于负二项分布模型。 scan_eb_zip： 基于零膨胀Poission分布模型期望值。(B. Allévius and Hhle 2017) scan_permutation：基于时间-空间置换模型。(Martin Kulldorff et al....

library(tidyverse) rm(list = ls()) 本系列主要是又来挖坑想对tidyverse(Wickham 2023)这个整合包的各个包的使用说明做一个完整的中文介绍。本章主要介绍purrr包(Wickham and Henry 2023)， purrr包通常被称为“优雅的循环迭代包”，它比for循环的速度快，代码更容读。 # map()函数 0.1 一个列表 map(.x, .f, ...)map系列函数会将一个函数应用到输入的列表或向量的每个元素上面，然后返回一个list。比如： x <- list(1:10, 11:20, 21:30) l1 <- list(x = c("a", "b"), y = c("c", "d")) map(l1, sort, decreasing = T ) # $x # [1] "b" "a" # # $y # [1] "d" "c" map_dbl(.x, .f, ...)返回双精度向量 # map_dbl(x, ~ mean(.x, na.rm = T)) map_dbl(x, mean) # [1] 5.5 15.5 25....

这些教程主要市基于“modrndive”包(Kim and Ismay 2022)的参考手册和自己浅薄的R使用经验编写，主要目的还是引导重庆的公共卫生人使用R语言分析数据，撰写报告，减轻重复劳动的负担。 图0.1: Family of tidy packages. 1 需要使用到的包： library(tidyverse) # tidy系列数据分析套装 library(moderndive) # 统计推断、数据模拟，统计学导学 library(nycflights13) # 航班数据库 library(gapminder) # gap minder 数据库 library(gtsummary) # 超棒的统计表格自动生成 library(skimr) # 描述性统计包 library(infer) # 统计推断 在R里面对于数据可视化当然要用到数据可视化的王者ggplot2(Wickham et al. 2022)。 Gapminder 数据 TED talk “The best stats you’ve ever seen”。其中的数据来源就是，随机抽取其中5行数据看看1.1。 表1.1: Gapminder数据 country continent year lifeExp pop gdpPercap Gambia Africa 1997 55.9 1235767 654 Netherlands Europe 1997 78....

library(tidyverse) library(survival) library(survminer) library(gtsummary) theme_gtsummary_journal("lancet") theme_gtsummary_language("zh-cn") 关于survial包做生存分析的记录 基础知识 生存分析 生存分析（survival analysis）是生物医学研究中常用的分析方法。在队列随访研究中，我们会事先定义一些观察终点，比如肿瘤复发、患者死亡、血压达标等，这些终点称为事件（event）。从研究开始到发生事件的时间间隔称为生存时间（survival time），某些场景下也称为失效时间（failure time）。由于生存时间数据具有以下两个特点，所以提出生存分析这一特殊的分析方法： （1）偏态分布：生存时间通常具有明显的偏态分布，有正态分布假设的统计方法不能适用。 （2）删失（censoring）：研究对象在观察时间内没有发生事件称为删失。一种情况是研究对象在中途失访或退出，导致没有观察到事件；另一种情况是超过了最长的随访时间事件仍未发生。删失数据是一种不完整数据，是生成分析独有的重要组成部分。 包里的常用函数 Surv() 函数主要用于对时间和状态变量进行转换，主要作为模型的自变量，放在模型方程左边，例如： Surv(time, status ) 右截断数据（right censored data）在进行随访观察中，研究对象观察的起始时间已知，但终点事件发生的时间未知，无法获取具体的生存时间，只知道生存时间大于观察时间，这种类型的生存时间称为右删失。 Surv(time, endpoint=='death') 状态变量是因子或字符的时候 Surv(t1, t2, status) 连续的过程变量 Surv(t1, ind, type='left') 左删失 左删失（Left censored）假设研究对象在某一时刻开始进入研究接受观察，但是在该时间点之前，研究所感兴趣的时间点已经发生，但无法明确具体时间，这种类型即为左删失数据。例如，某项关于脑卒中复发危险因素的研究，生存时间规定为从第一次脑卒中发病到下一次脑卒中发病之间的时间间隔。在研究起始时刻对研究对象进行问卷调查，询问是否发生过脑卒中，以及第一次脑卒中发病的时间，如果研究对象回答“发生过脑卒中。 aareg() Aalen’s additive regression model. coxph() Cox比例风险模型 coxph(Surv(time, status) ~ x, data=aml) -标准cox模型 coxph(Surv(t1, t2, stat) ~ (age + surgery) * transplant) -时间相关协变量...

0.1 新增感染者情况 8月17日0-24时，全国（除港澳台外）25个省市报告新增本土感染者3424例（较前日增加504例），新增数排前5位的依次是海南、西藏、新疆、陕西、青海，分别为2018例、889例、233例、65例、60例。详见图1。 Figure 1: 全国新冠肺炎新增感染者情况 0.2 累计报告感染者情况 自2022年7月1日起，全国累计报告感染者人数较多的省份（剔出无症状转确诊）为海南13814例，甘肃4516例，广西3739例，西藏3469例，新疆3193例，安徽2326例，河南1096例等省市,见图2。 Figure 2: 7月1日至今全国累计报告感染者分布图 0.3 近10内各省报告感染者数量变化趋势 对近10日内各省报告数据分析显示，感染者数量呈上升趋势的省份（增速）分别为西藏（0）、海南（0）、陕西（0）、云南（0）。 报告数量呈下降趋势的省份（降速）分别为（） ，见图3。 Figure 3: 近10日全国各省感染者报告情况 Table 1: 近10日各省报告感染者数线性回归分析结果 省份 变化趋势 标准误 p值 显著性 西藏 0.001 0 0.000 <0.05 海南 0.001 0 0.020 <0.05 青海 0.000 0 0.070 / 陕西 0.000 0 0.046 <0.05 兵团 0.000 0 0.084 / 江西 0....

使用到的R包： 1 gghigglight包的使用

初衷 在这次疫情处理过程中，了解到在梳理传播链的过程中，很多时候仍然是手工在powerpoint等软件绘制传播链的。采用这种方式的优点能够在图中根据设计者需要加入较多的信息，比如：人员大致位置分布，接触的途径和强度等信息。不足之处在于，在链条上节点（感染者）较少的时候还能够梳理得很明确，但一旦节点达到一定数量，其中关系复杂度将呈几何倍数增加（比如1人传多个，1人与多个感染者有接触之类）。 在这种情况下。单纯的手工整理，将耗费非常多的脑力。最严重的缺点是，当现场流调信息变更，对链条进行修订的时，其中一个节点或链接的变化，会因连锁作用导致整个链条的变化。节点越多，变化的影响范围越大，越复杂，就像整理线头一样。当感染者人数上升到一定数量时，手动整理已经变成了一件难以完成的事情。由于本人对R的热衷，探索了一下能不能使用软件自动化链就是自己懒嘛的方式绘制传播，使用igraph，ggraph和networkD3最终效果如下面几张图,个人觉得还是networkD3炫酷的互动效果最好。 具体制作过程 参见我使用的3个包的说明。。。。。详细步骤待补充。 数据 节点数据 节点数据里面只需要包含所有感染者的基本信息，比如编号，姓名，类别等等。 边数据 边数据最基础的要求为，节点数据左右感染者的对应关系，简单说就像Excel两列，第一列from， 第二列to，代表每一行两个感染者的关系，从谁传播到谁，当然这些资料需要辛苦在现场的流调专家们提供。 可视化 igraph 首先使用graph_from_data_frame(d =line, vertices = node, directed = T)将节点和边转换成igraph，就可以直接plot第一张图, 参数自己可以调节。 netwokd3 个人最喜欢的效果，使用igraph_to_networkD3命令，将igraph数据转换一下，就可以使用simpleNetwork,forceNetwork，sankeyNetwork(画出交互性网络图了。试验了下，手机浏览器一样可以互动，包括拖动节点，放大，移动等，非常棒的体验。 {"x":{"links":{"source":[2,8,10,12,3,7,56,2,14,13,1,56,0,2,2,2,2,3,4,5,6,7,7,1,56,1,9,2,11,6,24,3,56,7,12,12,1,12,12,6,7,10,15,15,7,28,30,12,7,56,16,14,14,17,37,18,53,2,56,1,5,3,6],"target":[10,11,12,12,13,14,15,16,17,18,19,19,1,20,21,22,23,24,25,26,27,28,29,2,2,30,31,32,33,34,35,35,36,37,38,39,3,40,41,42,43,44,45,46,47,47,47,48,49,4,50,51,52,53,53,54,55,5,6,7,8,9,9],"value":[8,4,8,3,8,8,4,8,8,4,8,3,3,8,8,8,8,8,8,8,8,8,8,8,3,4,8,8,8,8,3,8,4,8,6,8,3,8,8,8,8,8,8,8,8,3,3,8,8,4,8,8,8,8,3,8,4,8,4,4,8,4,2],"colour":["#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666","#666"]},"nodes":{"name":["0号","1号","2号","3号","4号","5号","6号","7号","8号","9号","10号","11号","12号","13号","14号","15号","16号","17号","18号","19号","20号","21号","22号","23号","24号","25号","26号","27号","28号","29号","30号","31号","32号","33号","34号","35号","36号","37号","38号","39号","40号","41号","42号","43号","44号","45号","46号","47号","48号","49号","50号","51号","52号","53号","54号","55号","56号"],"group":[2,2,2,1,7,2,1,6,5,1,2,5,8,4,3,9,2,3,4,2,2,2,2,2,1,7,2,1,6,6,6,1,2,5,1,1,2,3,8,8,8,8,1,6,2,9,9,6,8,6,2,3,3,3,4,3,2],"nodesize":[200,150,30.6,34,35.6,43.1,40.5,10.8,46.9,44.5,42.8,46.8,22.8,27.3,25.4,18.7,11,14.2,23.1,12.4,31.3,11.1,21.7,11.6,11.2,1.3,24,48.2,33.1,2.9,29.2,15.5,45.1,38,23.8,44.4,32.7,15,36,41.2,4.8,22.9,1.6,38.9,41.1,38.6,43.3,9.9,2.1,46.1,12,42.3,1.8,18.9,40.6,27.9,20]},"options":{"NodeID":"name","Group":"group","colourScale":"d3.scaleOrdinal(d3.schemeCategory20);","fontSize":12,"fontFamily":"serif","clickTextSize":30,"linkDistance":50,"linkWidth":"function(d) { return Math.sqrt(d.value); }","charge":-30,"opacity":2,"zoom":true,"legend":false,"arrows":true,"nodesize":true,"radiusCalculation":" Math.sqrt(d.nodesize)+6","bounded":false,"opacityNoHover":1,"clickAction":null}},"evals":[],"jsHooks":[]} {"x":{"links":{"source":[2,8,10,12,3,7,56,2,14,13,1,56,0,2,2,2,2,3,4,5,6,7,7,1,56,1,9,2,11,6,24,3,56,7,12,12,1,12,12,6,7,10,15,15,7,28,30,12,7,56,16,14,14,17,37,18,53,2,56,1,5,3,6],"target":[10,11,12,12,13,14,15,16,17,18,19,19,1,20,21,22,23,24,25,26,27,28,29,2,2,30,31,32,33,34,35,35,36,37,38,39,3,40,41,42,43,44,45,46,47,47,47,48,49,4,50,51,52,53,53,54,55,5,6,7,8,9,9],"value":[8,4,8,3,8,8,4,8,8,4,8,3,3,8,8,8,8,8,8,8,8,8,8,8,3,4,8,8,8,8,3,8,4,8,6,8,3,8,8,8,8,8,8,8,8,3,3,8,8,4,8,8,8,8,3,8,4,8,4,4,8,4,2]},"nodes":{"name":["0号","1号","2号","3号","4号","5号","6号","7号","8号","9号","10号","11号","12号","13号","14号","15号","16号","17号","18号","19号","20号","21号","22号","23号","24号","25号","26号","27号","28号","29号","30号","31号","32号","33号","34号","35号","36号","37号","38号","39号","40号","41号","42号","43号","44号","45号","46号","47号","48号","49号","50号","51号","52号","53号","54号","55号","56号"],"group":["0号","1号","2号","3号","4号","5号","6号","7号","8号","9号","10号","11号","12号","13号","14号","15号","16号","17号","18号","19号","20号","21号","22号","23号","24号","25号","26号","27号","28号","29号","30号","31号","32号","33号","34号","35号","36号","37号","38号","39号","40号","41号","42号","43号","44号","45号","46号","47号","48号","49号","50号","51号","52号","53号","54号","55号","56号"]},"options":{"NodeID":"name","NodeGroup":"name","LinkGroup":null,"colourScale":"d3.scaleOrdinal(d3.schemeCategory20);","fontSize":16,"fontFamily":null,"nodeWidth":15,"nodePadding":10,"units":"Letter(s)","margin":{"top":null,"right":null,"bottom":null,"left":null},"iterations":20,"sinksRight":true}},"evals":[],"jsHooks":[]} ggraph ggraph基本研用了ggplot2绘图的方式，画出来的图也相对更漂亮。首先使用tidygraph包将igraph类型的数据转换为ggraph更合适的元数据。然后可以愉快地使用ggplot2的方式画图了。

今天本地的疫情终于没有增加，抽了点时间关注其他地区的疫情形势。看官方通报的数据，对曾经的模范城市的疫情有兴趣，决定来简单看看。 1 数据获取 1.1 数据来源 要获取准确的数据，当然是上官方网站。打开上海市卫健委的官网(https://wsjkw.sh.gov.cn/xwfb/index.html)，疫情数据公告都在“新闻发布”栏目中疫情信息的标题中就包含了所有新增、确诊数据。真是太方便了。 library(rvest) library(tidyverse) library(lubridate) library(readxl) library(openxlsx) library(ggforce) library(mgcv) library(deSolve) library(FME) 1.2 数据 上海的疫情变化主要从3月开始 1.3 数据清洗 这一步比较麻烦的是对标题中日期字符的整理。使用str_extract_all命令后提取的日期，变成了列表。再合并为向量形势的日期格式数据时出了点麻烦。最后使用了笨办法for循环unlist后再paste0合并。其实直接用标题前的日期-1没什么大的误差，主要是在跟自己较劲搞得这么麻烦。 2 好了，开始分析吧 2.1 先画个简单的图看看大趋势 以报告感染者类型为颜色看，从20年1月到21年12月期间，上海的感染人数几乎在处于一个长期稳定的状态，从2022年3月开始呈直升飞机式的增长。奇怪的是中间咋有个空白区，没有数据。查看原始网页，发现网站从2021年11月6日-2022年1月1日没有更新数据。这是一个奇怪的现象。不过没关系，这不影响我们后面的分析，这对原因分析有很大关系。 2.2 输入感染者的趋势 对上海这波疫情有个合理的猜测是，1月初上海接纳大量某地的航班的，导致上海市输入疫情压力陡增，再加上Omicron变异株超强传播能力，双重压力下导致这个模范城市失守。好吧，我们来看看是否能验证输入压力陡增这个猜测。从图上看新增输入的感染者数量变化并不大，鉴于确诊和无症状都属与感染这，下一步我们把本地和输入的感染者的合计数的变化可视化看看。 2.3 感染者总数 先用二者做个散点图看看,如图，完全看不出啥关系啊。。。。输入感染者较多的时候，反而本地感染处于低水平。这个图像，线性回归暂时也不考虑。我们还是从感染者数量和时间的关系看看。 画出来如下，大致能看出在3月14日前，上海几乎没有本土感染病例。 2.4 2022年3-4月 我们把时间尺度拉大到2月中旬到3月看看。放大后（下图1）到3月14日，输入感染人数开始上升，此时上海本地感染数量仍然无幅度的改变，处于平稳状态。为了更仔细看清楚变化情况，我们将y轴的数量调整至0-1000例，这时候能清楚看到每天的变化，本地感染者（蓝色）数量从3月16日快速上升者。 2.5 重点关注3月的数据 2.6 地区分布 获取了地区分布的疫情数据，结果很明显，浦东新区感染者数量最多，见下表： Local 2022-04-07 2022-04-06 2022-04-05 2022-04-04 2022-03-28 2022-03-31 2022-03-30 2022-03-29 2022-03-27 2022-04-01 2022-03-24 2022-03-26 2022-04-03 2022-03-23 2022-03-20 2022-03-21 2022-03-19 2022-03-22 2022-03-25 浦东新区 9050 8457 8145 7071 2506 2407 2207 2183 1429 NA 193 323 NA 436 220 169 135 NA NA 闵行区 2257 2409 2937 1381 NA 392 780 987 619 1043 980 972 NA 256 NA 122 NA NA NA 徐汇区 2076 1107 920 1229 NA 226 404 1100 277 639 167 331 NA 212 46 130 122 NA NA 黄浦区 1380 1044 658 970 NA 121 361 110 56 260 204 320 824 40 42 98 72 118 NA 松江区 1288 781 796 1106 NA NA 476 426 190 948 184 158 NA 68 18 76 32 NA NA 普陀区 957 1033 483 254 NA 15 146 113 70 245 38 69 NA 68 36 31 50 NA NA 嘉定区 933 NA 481 237 NA 54 158 255 0 0 130 0 NA 69 0 NA 81 NA NA 长宁区 852 350 84 33 115 256 128 30 90 74 118 29 NA 22 13 52 24 38 8 虹口区 594 668 410 608 NA 128 84 73 27 61 58 50 NA 52 24 46 NA NA NA 宝山区 414 660 554 265 311 17 504 363 94 45 87 153 NA 138 26 134 NA 32 8 杨浦区 601 630 623 220 NA 174 100 99 70 134 21 23 NA 22 16 20 28 NA NA 奉贤区 92 556 144 64 230 183 130 96 90 161 48 168 119 12 26 22 12 NA NA 静安区 381 545 302 50 104 164 102 175 70 189 70 108 338 96 64 88 30 NA 4 崇明区 520 126 79 94 NA 44 386 55 472 166 164 54 NA 108 NA NA NA 38 82 青浦区 493 470 384 314 NA 95 77 174 86 174 26 44 NA 20 8 10 10 40 26 金山区 129 79 77 52 NA 42 66 26 0 49 0 0 NA 0 22 NA NA NA NA 画个玫瑰图看看，其实这个图并并能很好反应数据特征。...

对于一个公卫狗来说，很多工作都在分析数据、撰写报告总结的内容。这些文档有周报、月报、季报、年报。在使用R之前，经常会对同样架构的不同时间的的数据进行重复的分析和撰写。Rmarkdwon提供了非常便利的可重复复文档，简单来说就是根据跟新的数据库一键生成同样格式的文档。文档的类型可以是PDF和WORD。非常的便利。 1 Rmarkdown基础 1.1 markdown Markdown是一种轻量级标记语言，排版语法简洁，让人们更多地关注内容本身而非排版。它使用易读易写的纯文本格式编写文档，可与HTML混编，可导出HTML、PDF以及本身的.md格式的文件。因简洁、高效、易读、易写，Markdown被大量使用，如Github、Wikipedia、简书等。 Markdown已经是一种轻量化的语言，他的不足在于当需要转换为PDF或者WORD格式时有点力不从心了。针对这问题，目前供职与Rstudio的谢益辉大神在Markdown的基础上开发了rmarkdwon包。Rmarkdown的优点有： 具备markdown简单的文本标记语言； 相较于markdown有更丰富的文档排版能力； 可以在生成文档时直接运行插入的代码块，并将运行结果输出在最终文档中。 1.2 基本语法 Rmarkdown的基本语法与markdwon相同，几乎所有markdown能用的语法都可以在其中使用。 常用的markdown语法见表1.1，更详细的语法请参考Markdown官方教程 表1.1: Rmarkdown的优点 输入 结果 普通文本 普通文本 *斜体文本* 斜体文本 **加粗文本** 加粗文本 ***斜体加粗文本*** 斜体加粗文本 ~~删除线文本~~ 删除线文本 上标^文本^ 上标文本 `代码文本` 代码文本 # 一级标题 1.3 建立Rmarkdwon文件 在RStudio选择文件-新建-Rmarkdown就可以建立一个最简单的rmarkdown的模板。（图 1.1） 图1.1: 新建rmarkdown文件 新建文件后会弹出新的窗口，如图1.2所示。窗口中可以设置标题、作者和输出文件类型。在左边栏选择”From Template”会有更多的模板更选择。你也可以安装相应的R包扩展更多的模板。 图1....

加载常用的包 library(tidyverse) library(readxl) library(openxlsx) library(lubridate) 1 R基础相关 1.1 基础命令 1.1.1 set_names命令 set_names()命令可以很方便的创建带名称的向量，set_names(x, nm = x, ...)其中x表示需要赋名的变量，nm表示向量的名称，示例： set_names(1:4, letters[1:4]) # a b c d # 1 2 3 4 1.1.2 交集并集 interesect()求交集，union()并集,setdiff()补，结果如图1.1： 图1.1: R集合运算命令 1.2 数学公式测试 行内公式\(x^2+y^2=z^2\) 行间公式:\[f(x)=ax+b\] 1.3 逻辑运算符号 R中基础逻辑运算符号有（=;!=;>;<;>=;<=），返回的值是true和false，也可分别作为为数字1和0进行数学运算： x <- c(1, 2, 3, 4, 5, 6) x > 3 # [1] FALSE FALSE FALSE TRUE TRUE TRUE sum(x > 3) # [1] 3 1....