R语言入门与实践笔记(第四章)

R语言入门与实践笔记(第四章),第1张

在R中有 6中索引编写方式 ,包括 正整数、负整数、零、空格、逻辑值、名称

与正整数索引相反,它的含义是 不包含 负整数索引所对应的元素。

说实话,零索引并没有多大用处。这里就不介绍了

代表选取该索引位置所代表维度的所有元素。

当索引提供一个包含TRUE和FALSE逻辑值的向量时,R会匹配索引值为TRUE的元素。 此索引方式非常重要

编写一个可以返回第一行所有元素的函数

问题:这样每次发牌都是黑桃K,所以我们要在每次发完牌后进行洗牌,然后再发,现在写一个洗牌的函数

下面写一个输入进去deck输出一个洗牌后的数据框的函数

$ 可以提取数据框或列表对象中的值。

列表提取元素

掌握R语言的索引,最基本 *** 作为 写出对象名字,并在随后中括号里写出对应的索引即可 。若对象是一维的,如向量,只需要提供一个位置索引;若对象是二维的,如数据框,则提供两个位置索引,中间用逗号隔开。n维则用n个索引。另外数据框和列表还可用 $ 来索引。

地理加权回归(GWR)在R里面怎么实现?
121 人关注0 条评论
写回答
查看全部 5 个回答
写回答
叶山Shan Ye
GIS/地质/人文地理/可持续发展
A2A 谢邀,
我和我认识的一些人,刚开始用R做空间分析的时候,也遇到过这个问题。R这种开源的东西,优点是各种包很丰富,缺点是有些包的说明写得很乱,地理加权回归(GWR)的R包其实功能很强大,但大部分说明都不大靠谱。
GWR在R里面可以用好几个不同的包来实现,其中步骤最简单的是spgwr。思路就两步:建立窗口、用窗口扫全局。这其实就是GWR本质上的两步。比如我要在全美国范围内统计某两个(或多个)变量之间的回归关系,我可以做一个全局回归(global regression),但因为这些变量在空间分布上或许会有异质性(heterogeneity),表现在统计结果上就是空间不稳定性(nonstationarity),因此只看全局的统计,可能看不出什么结果来。举个不完全恰当但是很容易领会精神的例子,你比如说,我要分析亚洲范围内,经济发展程度与牛肉销量之间的关系,经济越发达的地方,人们就越吃得起牛肉。可是等我统计到印度的时候,坏了,印度大部分人不吃牛肉,这不是经济状况导致的,这一下就影响了全局统计的参考价值,那怎么办呢?我们可以建立一个窗口(正规说法是带宽窗口,bandwidth window),每次只统计窗口范围内的经济与牛肉销量的关系,然后用这个窗口去扫过全局的范围。等统计到印度的时候,印度内部的各地和印度自己比,吃牛肉的人的比例就不会突然减少,这样就能减少这种空间不稳定性对全局统计的影响。
所以,第一步就是要建立这样一个『窗口』。当然了,首先要安装包,我们要用到的R包有:
library(spgwr)
library(rgdal)
library(sf)
library(spData)
library(sp)
library(lattice)
library(ggplot2)
library(ggthemes)
其中,spgwr是做GWR的包,rgdal是用来读取矢量要素的,sf,sp和spData都是用来处理矢量数据的,别的基本都是画图用。
以下默认你会R和GWR的基本 *** 作。并且,以下只展现方法,不要纠结我的数据和结果,我随便找的数据,这个数据本身没有什么意义,所以做出的统计看起来很『壮观』。
我们先导入数据。这里我用的是美国本土48州各个县(county,也有翻译成郡的)的人口普查数据和农业数据,来源是ESRI Online数据库。为啥用这个数据呢?因为我电脑里面就存了这么个可以用来做GWR的数据
我们用rgdal读取数据,然后把它画出来看看
require(rgdal)
usa_agri <- readOGR(dsn = "~/Documents/Spatial", layer = "usa_counties")
plot(usa_agri)
会得到这个东西:
readOGR里面,dsn后面加储存shp的路径(加到文件夹为止),layer后面写shp的文件名(不加shp)。不喜欢rgdal的同学可以不用,用maptools或者spData等别的处理shp的R包代替。不过如果用maptools,要注意处理一下参考系。
我们看一下这个shp里面的列联表都有什么:
可见,shp里面有3108个县的数据,数据有61种。然后再看data下面有什么:
总之就是各种人口普查的数据,后面截不完图,还有经济、房地产和农业之类的数据。那我们就随便选两个来当变量。我就随便挑了,因变量选AVESIZE12,即2012年各个县农场的平均占地面积。自变量选POP_SQMI,也就是人口密度(每平方英里的人口)。
现在正式建立窗口,调用的是spgwr里面的gwrsel函数:
bw <- gwrsel( AVE_SIZE12 ~ POP_SQMI, data = usa_agri, gweight = gwrGauss,
verbose = FALSE, method = "cv")
其中~前后分别是因变量和自变量。GWR里因变量只能有1个,但自变量可以选多个,如果需要多个自变量的话,就在代码POP_SQMI之后用+号连接就行。gweight是你的空间加权的函数(随空间距离增大而不断衰减的函数,衰减率由下面要提到的带宽控制),这里用的是比较常用的高斯函数,其余的还有gwrbisquare等函数可以调用。verbose决定是否汇报制定窗口的过程。method是决定构建带宽窗口模型的方法,这里用的cv指的是cross validation,即交叉验证法,也是最常用的方法,简单说就是把数据分成不同的组,分别用不同的方法来做回归计算,计算完了之后记录下结果,然后打乱重新分组,再回归计算,再看结果,周而复始,最后看哪种计算方法的结果最靠谱,这种方法就是最优解。还有一种很常见的选择最佳拟合模型的方法是AIC optimisation法,把method后面的cv改成aic就可以用。具体AIC optimisation是什么:AIC(赤池信息准则)_百度百科。总之,空间加权函数和带宽窗口构建方法的选择是GWR里面十分重要的步骤。
以上便是固定带宽窗口的示意图。比如我在对佐治亚做GWR,这一轮的regression target是红色的这个县,根据做出来的窗口,圆圈以内的县都要被算为红色县的邻县,其权重根据高斯函数等空间权重函数来赋值,而圆圈以外的县,空间权重都赋为0。
不喜欢固定带宽窗口的同学也可以不用它,而是用符合Tobler地理学第一定律的非固定带宽邻域统计, *** 作方法是在gwrsel里面加一个命令adapt = TRUE,这样的情况下,根据你设置的k邻居数,每一轮统计的时候,和本轮对象在k以内相邻的多边形的权重参数会被赋值为0到1之间的一个数,比如下图:
我在对佐治亚做GWR,这一轮的regression target是红色的这个县,那么图上标为1的县就是红色县的1阶邻县,标为2的是2阶(邻县的邻县),标为3的是3阶(邻县的邻县的邻县)。如果用非固定带宽邻域统计,k为3,那么1、2、3都被定义为红色县的邻县,它们的权重从3到1依次增加,会按比例被赋上0和1之间的值,而其它没有标注的县,权重为0。
下一步就是用前一步做出的窗口去扫过全局区域:
gwr_result <- gwr(AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwrGauss, hatmatrix = TRUE)
这一步如果数据量大,可能会要跑一阵,跑完之后我们看看结果里面有什么:
Call:
gwr(formula = AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwrGauss, hatmatrix = TRUE)
Kernel function: gwrGauss
Fixed bandwidth: 2058803
Summary of GWR coefficient estimates at data points:
Min 1st Qu Median 3rd Qu Max Global
XIntercept 73883e+01 21081e+02 32802e+02 66691e+02 85705e+03 6255656
POP_SQMI -80085e+01 -45983e-01 -14704e-01 -73703e-02 -21859e-03 -00426
Number of data points: 3108
Effective number of parameters (residual: 2traceS - traceS'S): 1196193
Effective degrees of freedom (residual: 2traceS - traceS'S): 2988381
Sigma (residual: 2traceS - traceS'S): 104878
Effective number of parameters (model: traceS): 8490185
Effective degrees of freedom (model: traceS): 3023098
Sigma (model: traceS): 1042741
Sigma (ML): 10284
AICc (GWR p 61, eq 233; p 96, eq 421): 5210955
AIC (GWR p 96, eq 422): 520177
Residual sum of squares: 3287040139
Quasi-global R2: 04829366
基本上你做GWR该需要的结果这里都有了。比如窗口大小(Fixed bandwidth)是2058803,意思是前一步构建的带宽窗口是半径20588千米的圆。Effective number of parameters显示的是你带宽窗口的大小合不合适。Sigma是残差的标准差,这个值要尽量小。Residual sum of squares(RSS)也是对拟合程度的一个评估值。最重要的是最后那个R2,越靠近1说明统计的拟合度越好。我这里面Sigma很大,R2也不是很大,因为我这里只是呈现方法,用的数据本来就是互不相干、没什么太大意义的,所以不用太纠结。如果你是真正的统计数据要来做GWR,就需要注意这些值了。
然后,我们就可以把每个县的R2画在地图上。首先,前面报告里的这些数据,比如R2,要先自己去生成的GWR结果里面去找,然后自己再算一下每个县的local R2,并把它们赋值到shp里面去:

你好! r如何简便将pdf文件转换成word文件,这个问题提得很好,在这里我做一下回答!
pdf转换成word,其实很简单,只要使用软件就可以搞定了!我推荐的是迅捷pdf转换成word软件,功能韩强大,效果很好的!
迅捷pdf转换成word转换器软件特点:

1支持设置输出的文档中是否保留图像。

2支持输出MS Word文档(doc)和富文本格式(rtf)两种格式。

3支持自定义转换页面范围。

4支持转换加密后的PDF文件(需要手动输入PDF文档密码)。

5支持批量添加PDF文件。

6转换速度快,效果良好,可较好地保留PDF文档中的、超链接、布局。

7生成的Word文档可直接应用于编辑。

8不依赖于Adobe Acrobat,Acrobat Reader 软件。
迅捷pdf转换成word转换器软件的使用步骤:
第一步:添加PDF文件。点击软件界面的“添加PDF文件”,将需要转换的PDF文件加入到文件列表当中。
第二步:输出选项。这里指的是转换后的Word文件内容样式。“保留原始版面”可以确保转换前后的文件内容保持一致;“仅文本(无图像)”则意味着转换后的Word文件内容不含文件,可以减少转换后的文件体积。
第三步:输出格式。用户可以选择转换成为Word文件格式或者富文本文件格式。通常来说,选择Word文件格式较为常用。

迅捷pdf转换成word转换器下载地址

>

数据准备

1、通过表达式赋值创建

变量名←表达式
以上语句中的"表达式"部分可以包含多种运算符和函数。

2、通过transform函数创建

为原数据框添加新的列,可以改变原变量列的值,也可以赋值NULL删除列变量
transform( _data , )
data:要修改的数据;
:进行修改的内容。

1、variable[condition] <- expression
语句variable[condition] <- expression将仅在condition的值为TRUE时执行赋值。

2、使用within函数进行转化
within(data, expr, )
data:要处理的数据;
expr:计算表达式。

1、fix()函数

使用fix()函数调用交互式编辑器修改变量名。例如,要修改df数据集p8列的变量名称为v5,运行fix(df)结果如下:

edit和fix的区别
edit()函数也可以调出交互式编辑器,修改数据后关闭窗口发现数据还是原来的值,所以需要进行赋值 *** 作才能保存修改结果,比如我的数据修改可以写为df <- edit(df)。 fix()函数调出的交互式编辑器,修改数据后关闭窗口发现已经保存了修改后的值,不需要赋值 *** 作。

2、reshape包 rename()函数
rename(dataframe, c(oldname="newname", oldname="newname", ))
dataframe:要修改的数据框;
oldname="newname":指定修改前变量名和修改后变量名。

3、names()函数重命名变量名

参考资料:

目录

一般的VCF文件都很大,用手动提取里面的信息肯定不大现实。用 vcfR 就可以轻松实现。
vcfR 自带测试文件 vcfR_test 。就用这个文件来 *** 作一下吧。

在分区 Genotype 里,通过观察 FORMAT 列可以看到一共有四种类型的数据 GT:GQ:DP:HQ ,至于这四种类型的数据个各自代表什么意思大家可以查阅知乎百度谷歌。我们可以提取出我们想要的数据类型。比方说最重要的 GT (genotype)。

同样,我们也可以提取例如 DP (测序深度Read Depth)的数字矩阵。

值的注意的是这里用到了参数 asnumeric = TRUE 使得数据自动转换成了数字。但是并不是对所有类型的数据都有效,比方说我们重复一下提取 gt 。

在没有任何报错的情况下 gt 变成了一堆毫无意义的数字,很明显不合理,不要用这些经过错误转换的数据进行下一步分析,比方说喜闻乐见的主成分分析。

在一些类型的数据里可能会出现一个以上的结果,比方说上面的 HQ 数据。

一般情况下我们只需要每一列的第一个数字

不需要samtools之类的软件我们也可以实现vcf数据读取自由,关键是可以直接写入内存进行下一步的统计分析和数据可视化,个人感觉是很有效的提高了生产力。值得花时间学习一下这个工具。


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/13347773.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-07-19
下一篇 2023-07-19

发表评论

登录后才能评论

评论列表(0条)

保存