用GeoBUGS做简单的空间数据分析
第一节 实例介绍基本的空间模型
GeoBUGS是WinBUGS的一个模块,专门用来分析空间数据(spatial data)。由于和WinBUGS的基本模型结合得比较好,所以被广泛地使用。目前的GeoBUGS除了自身的地图格式外,还支持Splus, ArcInfo 以及 EpiMap的地图格式。当然了,在使用的时候需要做适当的转化才行。
下面是一个简单的例子,大家也可以在GeoBUGS的Manual中找到它。模型假设为条件自回归模型 Conditional Autoregressive(CAR)。数据为苏格兰唇癌疾病数据,反映的是苏格兰56个郡的唇癌发病率。这个数据比较经典,Clayton and Kaldor (1987) 和 Breslow and Clayton (1993)都曾在他们的论著中分析过该数据。
County Observed Expected Percentage SMR Adjacent
cases (Oi) cases (Ei) in agric.(xi) counties
1 9 1.4 16 652.2 5,9,11,19
2 39 8.7 16 450.3 7,10
… … … … … …
56 0 1.8 10 0.0 18,24,30,33,45,55
County 为郡的编号。
Observed cases(记作:Oi)为实际患病人数。
Expected cases(记作:Ei)为预计患病人数,这个人数基于当地的人口,对象的年龄、性别分布。
Percentage in agric.(记作:xi)为当地农业、渔业、林业人口所占当地总人口的比例。
Adjacent counties 指的是与当前郡相毗邻的郡的编号。
SMR(Standardised Mortality Ratios)为标准死亡率。
通过观察数据,我们可以发现SMR在某些时候(比如Oi和Ei较小时)出现奇异的值(如 0.0),所以我们需要通过smooth方法来调整SMR的值。这里我们采用的方法是在条件自相关(CAR)的先验假定下,拟合具有空间相关的随机混效Poisson模型。模型如下:
$$O_i \sim Poisson (\mu_i) $$
$$log \mu_i = log E_i + \alpha_0 + \alpha_1 x_i / 10 + b_i$$
其中$\alpha_0$
为intercept项反映的是各个区域间患病的相对基准风险。
$b_i$
反映的是与地域相关的潜在的患病风险因子。其他项不言自明。
需要重点提出的是这里的$b_i$
,在GeoBUGS中可以通过car.normal先验分布来描述。在贝叶斯统计中任河变量都可以通过一个分布来描述。
b[1:N] ~ car.normal(adj[], weights[], num[], tau)
adj[] 为邻接郡的编号
weights[] 为描述各个郡之间重要性差异的权因子
num[] 每个郡的相邻郡的个数
tau 反映的是精度,因为不知道,所以在模型设定时要将其放到先验参数中去。
通过前两次介绍的方法,我们很容易就可以得到模型的结果。下面我们来看看如何将结果反映到地图上去。
第二节 GeoBUGS的界面操作
GeoBUGS的地图工具配置界面
第一步,打开Map-> Map Tool菜单,选择Scotland这张地图
第二步,在variable中填O或者E或者b等等模型参数
第三步,设置分割点和地图模板
第四步,点击plot画图
当然还可以在quantity中设置不同的需要反映的量的类型。
很简单吧。
GeoBUGS生成的地图
GeoBUGS还提供了一些小工具,比如Adjacency Map来查看邻接图。
用GeoBUGS显示邻接地图
附录
以下是WinBUGS用到的模型代码:
#Model
model
{
for (i in 1:N) {
O[i] ~ dpois(mu[i])
log(mu[i]) <- log(E[i]) + alpha0 + alpha1 * X[i]/10 +
b[i]
RR[i] <- exp(alpha0 + alpha1 * X[i]/10 + b[i])
# Area-specific relative risk (for maps)
}
# CAR prior distribution for random effects:
b[1:N] ~ car.normal(adj[], weights[], num[], tau)
for (k in 1:sumNumNeigh) {
weights[k] <- 1
}
# Other priors:
alpha0 ~ dflat()
alpha1 ~ dnorm(0, 1e-05)
tau ~ dgamma(0.5, 5e-04)
# prior on precision
sigma <- sqrt(1/tau)
# standard deviation
}
#Data
list(N = 56, O = c(9, 39, 11, 9, 15, 8, 26, 7, 6,
20, 13, 5, 3, 8, 17, 9, 2, 7, 9, 7, 16, 31, 11, 7, 19, 15,
7, 10, 16, 11, 5, 3, 7, 8, 11, 9, 11, 8, 6, 4, 10, 8, 2,
6, 19, 3, 2, 3, 28, 6, 1, 1, 1, 1, 0, 0), E = c(1.4, 8.7,
3, 2.5, 4.3, 2.4, 8.1, 2.3, 2, 6.6, 4.4, 1.8, 1.1, 3.3, 7.8,
4.6, 1.1, 4.2, 5.5, 4.4, 10.5, 22.7, 8.8, 5.6, 15.5, 12.5,
6, 9, 14.4, 10.2, 4.8, 2.9, 7, 8.5, 12.3, 10.1, 12.7, 9.4,
7.2, 5.3, 18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7,
19.6, 3.4, 3.6, 5.7, 7, 4.2, 1.8), X = c(16, 16, 10, 24,
10, 24, 10, 7, 7, 16, 7, 16, 10, 24, 7, 16, 10, 7, 7, 10,
7, 16, 10, 7, 1, 1, 7, 7, 10, 10, 7, 24, 10, 7, 7, 0, 10,
1, 16, 0, 1, 16, 16, 0, 1, 7, 1, 1, 0, 1, 1, 0, 1, 1, 16,
10), num = c(3, 2, 1, 3, 3, 0, 5, 0, 5, 4, 0, 2, 3, 3, 2,
6, 6, 6, 5, 3, 3, 2, 4, 8, 3, 3, 4, 4, 11, 6, 7, 3, 4, 9,
4, 2, 4, 6, 3, 4, 5, 5, 4, 5, 4, 6, 6, 4, 9, 2, 4, 4, 4,
5, 6, 5), adj = c(19, 9, 5, 10, 7, 12, 28, 20, 18, 19, 12,
1, 17, 16, 13, 10, 2, 29, 23, 19, 17, 1, 22, 16, 7, 2, 5,
3, 19, 17, 7, 35, 32, 31, 29, 25, 29, 22, 21, 17, 10, 7,
29, 19, 16, 13, 9, 7, 56, 55, 33, 28, 20, 4, 17, 13, 9, 5,
1, 56, 18, 4, 50, 29, 16, 16, 10, 39, 34, 29, 9, 56, 55,
48, 47, 44, 31, 30, 27, 29, 26, 15, 43, 29, 25, 56, 32, 31,
24, 45, 33, 18, 4, 50, 43, 34, 26, 25, 23, 21, 17, 16, 15,
9, 55, 45, 44, 42, 38, 24, 47, 46, 35, 32, 27, 24, 14, 31,
27, 14, 55, 45, 28, 18, 54, 52, 51, 43, 42, 40, 39, 29, 23,
46, 37, 31, 14, 41, 37, 46, 41, 36, 35, 54, 51, 49, 44, 42,
30, 40, 34, 23, 52, 49, 39, 34, 53, 49, 46, 37, 36, 51, 43,
38, 34, 30, 42, 34, 29, 26, 49, 48, 38, 30, 24, 55, 33, 30,
28, 53, 47, 41, 37, 35, 31, 53, 49, 48, 46, 31, 24, 49, 47,
44, 24, 54, 53, 52, 48, 47, 44, 41, 40, 38, 29, 21, 54, 42,
38, 34, 54, 49, 40, 34, 49, 47, 46, 41, 52, 51, 49, 38, 34,
56, 45, 33, 30, 24, 18, 55, 27, 24, 20, 18), sumNumNeigh = 234)
#Inits
list(tau = 1, alpha0 = 0, alpha1 = 0, b = c(0, 0,
0, 0, 0, NA, 0, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
WinBUGS在统计分析中的应用 第三部分完
发表/查看评论