R 中可以用 = 或者 <- 来进行赋值 ， <-的快捷键是 alt + - 。
R的下标是从1开始的，和python等不同（python四从0开始的）
当然我们也可以用逻辑进行筛选，例如
负数下标表示不选这个这些下标，例如：
c() 可以合并向量，例如
向量有个比较有趣的性质，当两个向量进行操作时，如果长度不等， 长度比较短的一个会复制自己直到自己和长的一样长。
a 自动变成了 c(3,4,3,4) 然后与b相加 ， 得到了下面的结果。
遇到不懂得函数，可以用help("函数")查看函数用法。
矩阵，从本质上来说就是多维的向量，我们来看一看 我们如何新建一个矩阵。
可以看到向量元素变为矩阵元素的方式是按列的，从第一列 到第二列，如果我们想按行输入元素，那么需要加入 byrow = TRUE 的参数：
与向量相似，我们可以用下标来筛选矩阵， 例如：
a[行,列]
当我们对两个矩阵相乘，我们得到的结果是 对应元素两两相乘的结果，例如：
而这不是我们想要的矩阵乘法，在 R 中我们在乘法旁边加两个 百分号来做矩阵乘法：
此外，我们可以用 t() 来求矩阵的转置 ， 用 solve() 来求矩阵的逆。
数据框类似矩阵，与矩阵不同的是，数据框可以有不同的数据类型。 一般做数据分析，我们把一个类似 excel 的表格读入 R ，默认的格式 就是数据框 ， 可见数据框是一个非常重要的数据结构。
一般来说我们需要分析的数据，每一行代表一个样本，每一列代表一个 变量。
下面我们用 R 内置的数据集 iris 来看一看数据框的使用。
我们用 data 函数调入了 iris 这个数据集 ， 然后用 head 函数来看一看这个数据 的前几行 ， 可以看到有 sepal 的长度，宽度，petal 的长度和宽度，还有一个变量 Species 来描述样本的类别。
我们可以用 summary 函数来对数据集做大致的了解。
可以直观地看到每个变量的信息，对于几个数值变量，我们可以看到最小值，中位数等等统计信息。而对于 Species 这个分类变量，我们看到的是计数信息。
筛选数据框与矩阵相似,都可以通过数字下标来获取子集，不同地是因为数据框有不同的列名,我们也可以通过列名来获取某一特定列,例如：
我们可以用 names() 函数来获取数据框的列名
并可以通过为其赋值改变列的名字。
列表是一种递归式的向量,我们可以用列表来存储不同类型的数据，比如：
列表有多种索引方式，可以用如下方式获取。
今天我们实验的对象就是一组从原始 R 进化出来的工具链 Tidyverse , 它是由 Hadley Wickham 主导开发的一系列 R 包的集合。 Tidyverse 继承了R语言进行快速统计分析的优势 ， 并实现了一些新的理念 ， 例如 magrittr 包中的管道操作 ， 让线性嵌套的函数组合变得更加清晰易懂；可视化方面中的 ggplot ，使绘图变成搭积木式的图层叠加。
这样的小发明有的改变了分析的运作方式 ， 有的改变了使用者的认知方式 ， 聚在一起形成了一种新的数据分析的生态链 。具体来看 ， Tidyverse 有如下核心组件:
mpg 数据集是刻画不同汽车的排放状况的一个数据集， 总过有 234 个样本 ， 11 个变量 。 这 11 个变量分别是:
manufacture: 制造商
model: 车型
dispel： 汽车排放量
year: 制造年度
cyl: 排气管数量
trans: 排放类型
drv: 驱动方式
cty: 每公里耗油量（城市道路）
hwy: 每公里耗油量（高速路）
fl: 油的种类
class: 车的类型
更多数据相关信息可以通过 help(mpg) 指令获取。
在属性映射中加入 color=class 参数后 ， 我们可以看到每个点的汽车对应的类型被用 不同颜色表现了出来 ， 对于散点图 ， 还有 size(大小) , shape（形状） 等等参数 可以用于确定点的属性。
对于条形图的y轴就是数据框中原本的数值时，必须将geom_bar()函数中stat(统计转换)参数设置为’identity’，即对原始数据集不作任何统计变换，而该参数的默认值为’count’，即观测数量。
R语言-v1-基础知识
Iretara 12-17 21:18
以例题的形式简述R语言基础知识
# 读取文件
setwd(" 文件链接的时候，用 / ")
install.packages(" readxl ")
library(readxl)
library (tidyverse)
hw1_a<- read_excel ("hw1_a.xlsx", col_types=c("numeric", "numeric", "numeric", "numeric", "numeric") )
hw1_b<- read_excel ("hw1_b.xlsx")
#读取csv
library(readr)
hw1_a<- read_csv ("/")
View(hw1_a)
# 描述型函数
hw1_a + hw1_b 表
#描述最小值，最大值，中值，均值，标准差
Str (hw1_a) #查看数据并指出各个 变量的形式
summary (hw1_a) #指出各个变量的形式， 最小值，最大值，中值，均值
library(psych)
describe (hw1_a) #比summary更简便的方法， 可以直接读取标准差等；但是，使用describe不可读取 NA值, 可以尝试使用 Hmisc包中 describe
描述型函数-R
# 连接
hw1_a %>% inner_join (hw1_b, by ="ID")
hw1_a %>% left_join (hw1_b, by ="ID")
hw1_a %>% right_join (hw1_b, by ="ID")
hw1_a %>% full_join (hw1_b, by ="ID")
inner_join<- inner_join (hw1_a,hw1_b, by =“ID”) #报告合并后的 总行数 ，178行
full_join<- full_join (hw1_a,hw1_b, by ="ID")
( nrow (full_join)) #报告合并后的 总行数 ，200行
> length (full_join$ID)
#找出各个列的 缺失值
i<-NA
a<-NA
for(i in 1:length(full_join[1,])){ a[i]<- sum(is.na( full_join[,i] ) ) }
paste("缺失值是",a)
#缺失值总数
sum(is.na(full_join))
#删除缺失值 na.omit()
full_join1=filter(full_join,!is.na(full_join[2]))
full_join1=filter(full_join1,!is.na(full_join1[3]))
full_join1=filter(full_join1,!is.na(full_join1[4]))
full_join1=filter(full_join1,!is.na(full_join1[5]))
full_join1=filter(full_join1,!is.na(full_join1[6]))
full_join1=filter(full_join1,!is.na(full_join1[7]))
full_join1=filter(full_join1,!is.na(full_join1[8]))
sum(is.na(full_join1))
找出Income中的 极端值 并滤掉对应行的数据
quantile (hw1_a$Income,c(0.025,0.975))
hw1_a2= filter (hw1_a,Income>14168.81 &Income<173030.92)
#使用dplyr进行数据转换
arrange()
>arrange (hw1_a,Income) #默认升序
>arrange(hw1_a, desc (Income)) #desc降序，NA排序一般最后
select()
>select (hw1_a, - (Years_at_Address:Income)) #不要变量
>rename (hw1_a, In_come=Income) #改名
>select（hw1_a,Income, exerything ()） #把Income放在前面
拓例题1：
library(nycflights13)
view(flights)
#counts
（1）
not_cancelled <- flights %>%
filter(! is.na(dep_delay), !is.na(arr_delay))
（2）
not_cancelled %>%
group_by (year,month,day) %>%
summarize (mean=mean(dep_delay))
（3）
delays <- not_cancelled %>%
group_by (tailnum) %>%
summarize (delay=mean(arr_delay))
ggplot (data=delays,mapping=aes(x= delay))+
geom_freqpoly (binwidth=10) #freqpoly
（4）
delays <- not_cancelled %>%
group_by(tailnum) %>%
summarize(delay=mean(arr_delay,na.rm=TRUE), n=n() ) #tailnum的次数
ggplot(data=delays,mapping=aes(x= n, y=delay))+
geom_point(alpha=1/10)
拓例题2：
#请按照价格的均值，产生新的变量price_new, 低于均值为“低价格”，高于均值为“高价格”。 同样对市场份额也是，产生变量marketshare_new, 数值为“低市场份额”和“高市场份额”
price=data1$price
pricebar=mean(price)
price_new= ifelse (price>pricebar,“高价格”,”低价格”)
marketshare=data1$marketshare
marketsharebar=mean(marketshare)
marketshare_new=ifelse(marketshare>marketsharebar ,“高市场份额”,”低市场份额”)
data1= mutate (data1,price_new,marketshare_new)
#可视化
#将Income 对数化
lninc<- log (hw1_a$Income)
#画出直方图和 density curve密度曲线
hist (lninc,prob=T)
lines ( density (lninc),col="blue")
# 添加额外变量 的办法，在 aes()中添加 样式 (color、size、alpha、shape)
ggplot(data=inner_join)+
geom_point(mapping = aes(x=Years_at_Employer,y= Income, alpha= Is_Default))
# 按照Is_Default 增加一个维度，使用明暗程度作为区分方式
ggplot(data=inner_join)+
geom_point(mapping = aes(x=Years_at_Employer,y= Income,
alpha=factor( Is_Default ) ))
#使用形状作为另外一种区分方式
ggplot(data=inner_join)+
geom_point(mapping = aes(x=Years_at_Employer,y= Income,
shape=factor( Is_Default)))
可视化-R
拓展：
#将 flight1 表和 weather1 表根据共同变量进行内连接，随机抽取 100000 行数据， 将生产的结果保存为 flight_weather。 (提示：sample_n()函数，不用重复抽取)
flight_weather <- inner_join(flight1, weather1) %>% sample_n(100000)
# 从 flight_weather表中对三个出发机场按照平均出发延误时间排降序，并将结果保留在 longest_delay表中。把结果展示出来
longest_delay<- flight_weather %>%
group_by(origin) %>%
summarize(delay=mean(dep_delay, na.rm=TRUE )) %>%
arrange(desc(delay))
#根据不同出发地（origin）在平行的 3 个图中画出风速 wind_speed（x 轴）和出发 延误时间 dep_delay（y 轴）的散点图。
ggplot(data= flight_weather) +
geom_point(mapping=aes(x=wind_speed,y=dep_delay))+
facet_grid(.~origin, nrow = 3 ) # 按照class分类，分成3行
#根据 flight_weather 表，画出每个月航班数的直方分布图，x 轴为月份，y 轴是每个 月份航班数所占的比例。
ggplot(data=flight_weather)+
geom_bar(mapping=aes(x=month, y=..prop .., group=1))
#根据 flight_weather 表，画出每个月航班距离的 boxplot 图，x 轴为月份，y 轴为 航行距离, 根据的航行距离的中位数从低到高对 x 轴的月份进行重新排序
ggplot(data=flight_weather)+
geom_boxplot(mapping=aes(x= reorder (month,distance,FUN=median),y=distance))
线性回归
# 以Income作为因变量，Years at Employer作为自变量，进行 OLS回归
m1<- lm (Income ~ Years_at_Employer,data=hw1_a)
#通过***判断显著性
summary (m1)
#画出拟合直线
ggplot(data= hw1_a)+
geom_point(aes(x=Income,y=Years_at_Employer))+
geom_abline(data= m1,col= "blue")
#证明拟合直线是最优的
b0=runif(20000,-5,5)
b1=runif(20000,-5,5)
d<-NA
sum<-NA
n<-1
while(n<=20000){
for(i in 1:24){
d[i]<-(hw1_a $ Income[i]-b0[n]-b1[n]*hw2$ Years_at_Employer[i])^2}
sum[n]<-sum(d)
n<-n+1
}
resi=m1$residuals
resi2=sum(resi^2)
check=sum(as.numeric(sum<resi2))
check
地理加权回归（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里面的gwr.sel函数：
bw <- gwr.sel( AVE_SIZE12 ~ POP_SQMI, data = usa_agri, gweight = gwr.Gauss,
verbose = FALSE, method = "cv")
其中~前后分别是因变量和自变量。GWR里因变量只能有1个，但自变量可以选多个，如果需要多个自变量的话，就在代码POP_SQMI之后用+号连接就行。gweight是你的空间加权的函数（随空间距离增大而不断衰减的函数，衰减率由下面要提到的带宽控制），这里用的是比较常用的高斯函数，其余的还有gwr.bisquare等函数可以调用。verbose决定是否汇报制定窗口的过程。method是决定构建带宽窗口模型的方法，这里用的cv指的是cross validation，即交叉验证法，也是最常用的方法，简单说就是把数据分成不同的组，分别用不同的方法来做回归计算，计算完了之后记录下结果，然后打乱重新分组，再回归计算，再看结果，周而复始，最后看哪种计算方法的结果最靠谱，这种方法就是最优解。还有一种很常见的选择最佳拟合模型的方法是AIC optimisation法，把method后面的cv改成aic就可以用。具体AIC optimisation是什么：AIC（赤池信息准则）_百度百科。总之，空间加权函数和带宽窗口构建方法的选择是GWR里面十分重要的步骤。
以上便是固定带宽窗口的示意图。比如我在对佐治亚做GWR，这一轮的regression target是红色的这个县，根据做出来的窗口，圆圈以内的县都要被算为红色县的邻县，其权重根据高斯函数等空间权重函数来赋值，而圆圈以外的县，空间权重都赋为0。
不喜欢固定带宽窗口的同学也可以不用它，而是用符合Tobler地理学第一定律的非固定带宽邻域统计，操作方法是在gwr.sel里面加一个命令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 = gwr.Gauss, hatmatrix = TRUE)
这一步如果数据量大，可能会要跑一阵，跑完之后我们看看结果里面有什么：
Call:
gwr(formula = AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 205880.3
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu.Max. Global
X.Intercept. 7.3883e+01 2.1081e+02 3.2802e+02 6.6691e+02 8.5705e+03 625.5656
POP_SQMI -8.0085e+01 -4.5983e-01 -1.4704e-01 -7.3703e-02 -2.1859e-03 -0.0426
Number of data points: 3108
Effective number of parameters (residual: 2traceS - traceS'S): 119.6193
Effective degrees of freedom (residual: 2traceS - traceS'S): 2988.381
Sigma (residual: 2traceS - traceS'S): 1048.78
Effective number of parameters (model: traceS): 84.90185
Effective degrees of freedom (model: traceS): 3023.098
Sigma (model: traceS): 1042.741
Sigma (ML): 1028.4
AICc (GWR p. 61, eq 2.33p. 96, eq. 4.21): 52109.55
AIC (GWR p. 96, eq. 4.22): 52017.7
Residual sum of squares: 3287040139
Quasi-global R2: 0.4829366
基本上你做GWR该需要的结果这里都有了。比如窗口大小（Fixed bandwidth）是205880.3，意思是前一步构建的带宽窗口是半径205.88千米的圆。Effective number of parameters显示的是你带宽窗口的大小合不合适。Sigma是残差的标准差，这个值要尽量小。Residual sum of squares（RSS）也是对拟合程度的一个评估值。最重要的是最后那个R2，越靠近1说明统计的拟合度越好。我这里面Sigma很大，R2也不是很大，因为我这里只是呈现方法，用的数据本来就是互不相干、没什么太大意义的，所以不用太纠结。如果你是真正的统计数据要来做GWR，就需要注意这些值了。
然后，我们就可以把每个县的R2画在地图上。首先，前面报告里的这些数据，比如R2，要先自己去生成的GWR结果里面去找，然后自己再算一下每个县的local R2，并把它们赋值到shp里面去：
