加权后的数据怎么用r转换出来

Posted 2023-05-13

地理加权回归（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.33; p. 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里面去： 参考技术A R语言学习笔记总结

R语言初步-用dplyr进行数据转换

install.packages("tidyverse")
install.packages("nycflights13")#仍然记得要先安装
library(nycflights13)#航班信息文件
library(tidyverse)

?flights#查看数据信息的说明书
flights#查看航班信息
5、其他常用的摘要函数

之前使用了均值、求和和计数

5.1、位置度量：median()函数

median()用法和mean()类似，只不过是中位数而已
注意：本例与median()无关
Not_cancelled <- flights %>%
filter(!is.na(dep_delay),!is.na(arr_delay))
Not_cancelled %>%
group_by(year,month,day)%>%
summarise( #平均延误时间
avg_delay1=mean(arr_delay), #平均正延误时间
avg_delay2=mean(arr_delay[arr_delay>0]), )
运行：
A tibble: 365 x 5
Groups: year, month [12]
year month day avg_delay1 avg_delay2
1 2013 1 1 12.7 32.5
2 2013 1 2 12.7 32.0
3 2013 1 3 5.73 27.7
4 2013 1 4 -1.93 28.3
5 2013 1 5 -1.53 22.6
6 2013 1 6 4.24 24.4
7 2013 1 7 -4.95 27.8
8 2013 1 8 -3.23 20.8
9 2013 1 9 -0.264 25.6
10 2013 1 10 -5.90 27.3
... with 355 more rows
5.2、分散程度度量：sd()、IQR()、mad()函数

* sd():标准误差函数：standard deviation，分散程度的标准度量方式
* IQR():四分位距
* mad()：绝对中位差

注：mad()与IQR()基本等价，但是IQR()更适合有离群点的情况。
Not_cancelled %>%
group_by(dest)%>%
summarise(
distance_sd=sd(distance))%>% #计算distance列的标准误差
arrange(desc(distance_sd) #降序排序此行
)
运行：
A tibble: 104 x 2 dest distance_sd
1 EGE 10.5
2 SAN 10.4
3 SFO 10.2
4 HNL 10.0
5 SEA 9.98
6 LAS 9.91
7 PDX 9.87
8 PHX 9.86
9 LAX 9.66
10 IND 9.46
... with 94 more rows
5.3、秩的度量：min()、quantile()、max()函数

quantile():分位数函数，是中位数函数的拓展
使用说明：quantile(x,0.25)是指将x按从小到大顺序排列，找到大于前25%，小于后75%的值。

每天最早和最晚的航班是是什么时候：
Not_cancelled %>%
group_by(year,month,day)%>%
#先按时间分组
summarise(
first=min(dep_time), #最小值
last=max(dep_time) #最大值
)
运行：
A tibble: 365 x 5
Groups: year, month [12]
year month day first last
1 2013 1 1 517 2356
2 2013 1 2 42 2354
3 2013 1 3 32 2349
4 2013 1 4 25 2358
5 2013 1 5 14 2357
6 2013 1 6 16 2355
7 2013 1 7 49 2359
8 2013 1 8 454 2351
9 2013 1 9 2 2252
10 2013 1 10 3 2320
... with 355 more rows
5.4、定位度量：first()、nth()、last()函数

这三个函数的作用相当于x[1]、x[2]、x[length(x)]
通过此函数也可以找出最早和最晚出发的航班

Not_cancelled %>%
group_by(year,month,day)%>%
summarise(
first_dep=first(dep_time),
last_dep=last(dep_time)
)

运行：
A tibble: 365 x 5
Groups: year, month [12]
year month day first_dep last_dep
1 2013 1 1 517 2356
2 2013 1 2 42 2354
3 2013 1 3 32 2349
4 2013 1 4 25 2358
5 2013 1 5 14 2357
6 2013 1 6 16 2355
7 2013 1 7 49 2359
8 2013 1 8 454 2351
9 2013 1 9 2 2252
10 2013 1 10 3 2320
... with 355 more rows
5.5、计数n(),count()

n():不需要任何参数，返回当前分组的大小
sum(!is.na(x)):计算非缺失值的数量
n_distinct(x):计算唯一值的数量
count()函数：用于只需要计数的情况

例如：
计算哪个目的地有最多的航空公司？
Not_cancelled %>%
group_by(dest)%>%
summarise(
carriers=n_distinct(carrier))%>%
arrange(desc(carriers))
运行：
A tibble: 104 x 2
dest carriers
1 ATL 7
2 BOS 7
3 CLT 7
4 ORD 7
5 TPA 7
6 AUS 6
7 DCA 6
8 DTW 6
9 IAD 6
10 MSP 6
... with 94 more rows
count()函数用法举例：计算目的地不同的飞机数量
Not_cancelled %>%
count(dest)
运行：
A tibble: 104 x 2
dest n
1 ABQ 254
2 ACK 264
3 ALB 418
4 ANC 8
5 ATL 16837
6 AUS 2411
7 AVL 261
8 BDL 412
9 BGR 358
10 BHM 269
... with 94 more rows
count()函数中可以添加加权变量，例如distance，用于计算飞机飞行里程（相当于求和）
Not_cancelled %>%
count(tailnum,wt=distance)

运行：
A tibble: 4,037 x 2
tailnum n
1 D942DN 3418
2 N0EGMQ 239143
3 N10156 109664
4 N102UW 25722
5 N103US 24619
6 N104UW 24616
7 N10575 139903
8 N105UW 23618
9 N107US 21677
10 N108UW 32070
... with 4,027 more rows
5.6、逻辑值的计数和比例

当需要用数值表示结果，TRUE=1，FALSE=0。
sum():可以找出TRUE的数量
mean():可以找出比例

以下一例：找出出发时间小于5:00的航班总数
Not_cancelled %>%
group_by(year,month,day)%>%
summarise(
n_nearly=sum(dep_time<500) #出发时间小于5:00的航班总数
)
运行：
A tibble: 365 x 4
Groups: year, month [12]
year month day n_nearly
1 2013 1 1 0
2 2013 1 2 3
3 2013 1 3 4
4 2013 1 4 3
5 2013 1 5 3
6 2013 1 6 2
7 2013 1 7 2
8 2013 1 8 1
9 2013 1 9 3
10 2013 1 10 3
... with 355 more rows
sum(dep_time<500)换成count(dep_time<500)是没有用的，sum相当于计算了返回值1，而dep_time<50这样的逻辑表达，count()函数是不支持的，其中牵扯到数据的逻辑。

以下一例：找出延误超过一小时的航班比例
Not_cancelled %>%
group_by(year,month,day)%>%
summarise(
hour_perc=mean(arr_delay>60) #延误超过一小时的航班
)
运行：
A tibble: 365 x 4
Groups: year, month [12]
year month day hour_perc

1 2013 1 1 0.0722
2 2013 1 2 0.0851
3 2013 1 3 0.0567
4 2013 1 4 0.0396
5 2013 1 5 0.0349
6 2013 1 6 0.0470
7 2013 1 7 0.0333
8 2013 1 8 0.0213
9 2013 1 9 0.0202
10 2013 1 10 0.0183
... with 355 more rows 参考技术B 第一步就是要建立这样一个『窗口』。当然了，首先要安装包，我们要用到的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里面十分重要的步骤。 参考技术C 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)
会得到这个东西： 参考技术D 这种开源的东西，优点是各种包很丰富，缺点是有些包的说明写得很乱，地理加权回归（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读取数据

R中的AIC：使用加权数据时手动值与内部值的差异

【中文标题】R中的AIC：使用加权数据时手动值与内部值的差异

【英文标题】：AIC in R: differences in manual vs. internal value when using weighted data

【发布时间】：2018-06-22 09:12:31

【问题描述】：

我正在尝试使用 R 进行基于 AIC 统计的模型选择。在比较带或不带加权的线性模型时，我在 R 中的代码告诉我，与不加权相比，加权更可取，并且这些结果在其他软件 (GraphPad Prism) 中得到了证实。我有使用来自标准曲线的真实数据的示例代码：

#Linear Curve Fitting
a <- c(0.137, 0.412, 1.23, 3.7, 11.1 ,33.3)
b <- c(0.00198, 0.00359, 0.00816, 0.0220, 0.0582, 0.184)
m1 <- lm(b ~ poly(a,1))
m2 <- lm(b ~ poly(a,1), weight=1/a)
n1 <- 6 #Number of observations
k1 <- 2 #Number of parameters

当我使用 R 中的内部函数或通过手动计算计算 AIC 时：

AIC = n + n log 2π + n log(RSS/n) + 2(k + 1) 有 n 个观测值和 k 参数

我得到了非加权模型的等效 AIC 值。当我分析加权的影响时，手动 AIC 值较低，但最终结果是内部和手动 AIC 都表明优先考虑加权。

> AIC(m1); n1+(n1*log(2*pi))+n1*(log(deviance(m1)/n1))+(2*(k1+1))
[1] -54.83171
[1] -54.83171
> AIC(m2); n1+(n1*log(2*pi))+n1*(log(deviance(m2)/n1))+(2*(k1+1))
[1] -64.57691
[1] -69.13025

当我尝试使用非线性模型进行相同的分析时，内部函数和手动计算之间的 AIC 差异更加深刻。以下是示例 Michaelis-Menten 动力学数据的代码：

c <- c(0.5, 1, 5, 10, 30, 100, 300)
d <- c(3, 5, 20, 50, 75, 200, 250)
m3 <- nls(d ~ (V * c)/(K + c), start=list(V=10, K=1))
m4 <- nls(d ~ (V * c)/(K + c), start=list(V=10, K=1), weight=1/d^2)
n2 <- 7
k2 <- 2

按照前两个模型的说明计算 AIC：

> AIC(m3); n2+(n2*log(2*pi))+n2*(log(deviance(m3)/n2))+(2*(k2+1))
[1] 58.48839
[1] 58.48839
> AIC(m4); n2+(n2*log(2*pi))+n2*(log(deviance(m4)/n2))+(2*(k2+1))
[1] 320.7105
[1] 0.1538546

与线性示例类似，当数据未加权 (m3) 时，内部 AIC 和手动 AIC 值相同。权重 (m4) 会出现问题，因为手动 AIC 估计要低得多。这种情况类似于在相关问题AIC with weighted nonlinear regression (nls) 中提出的问题。

我之前提到过 GraphPad Prism，对于上面给出的模型和数据集，当使用加权时，它的 AIC 较低。那么我的问题是，为什么在对数据进行加权时，R 中的内部与手动 AIC 估计存在如此差异（非线性模型与线性模型的结果不同）？最终，我应该认为内部 AIC 值还是手动值更正确，还是我使用了错误的公式？

【问题讨论】：

查看另一个问题中的 cmets，您的权重总和是否为 1？就“正确”而言，您链接的问题中有一个交叉验证的链接表明AIC is not a good metric for nonlinear models，因此统计上的正确答案可能是“都不正确”。

【参考方案1】：

您看到的差异是由于在加权模型的手动计算中使用了未加权对数似然公式。例如，您可以通过以下调整复制 m2 和 m4 的 AIC 结果：

对于m2，您只需从计算中减去sum(log(m2$weights))：

AIC(m2); n1+(n1*log(2*pi))+n1*(log(deviance(m2)/n1))+(2*(k1+1)) - sum(log(m2$weights))
[1] -64.57691
[1] -64.57691

对于m4，您必须将deviance 调用与加权残差计算交换，并从结果中减去n2 * sum(log(m4$weights))：

AIC(m4); n2+(n2*log(2*pi))+n2*(log(sum(m4$weights * m4$m$resid()^2)/n2))+(2*(k2+1)) - n2 * sum(log(m4$weights))
[1] 320.7105
[1] 320.7105

我相信logLikin m2 使用的公式的推导是非常直接和正确的，但我不确定m4。从阅读有关logLik.nls()（example 1，example 2）的其他一些线程来看，似乎对 nls 估计的正确方法有些混淆。总而言之，我相信AIC 对m2 是正确的；我无法验证加权 nls 模型的数学运算，并且在这种情况下倾向于再次使用 m2 公式（但用加权残差替换 deviance 计算），或者（可能更好）不使用 @987654342 @ 为nls 模型

【讨论】：

感谢丹尼尔森的回复！您对加权 AIC 方程的渲染与我正在研究的另一个模型的内部 AIC 匹配。