R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化...

全文链接：https://tecdat.cn/?p=35518

在统计建模过程中，经常会遇到空间自相关性的问题。空间自相关性是指相近位置的观测值往往比远离位置的观测值更相似。在尝试估计参数或进行预测时，空间自相关性可能会导致结果产生偏差（点击文末“阅读原文”获取完整代码数据）。

相关视频

INLA（Integrated Nested Laplace Approximation，集成嵌套拉普拉斯近似）是一种在潜在高斯模型中，包括具有空间成分的模型，进行贝叶斯推断的流行工具。本文中我们将帮助客户探讨如何使用INLA处理统计建模中的空间自相关性。

本文将带您逐步进行一项分析。这将涉及：

• 在INLA中进行模型选择。

• INLA模型的组成部分（基础）。

• 设置空间分析。

• 空间模型的修改和规格。

• 时空分析（简要涉及）。

数据

本文将使用一组捕获的野生动物的数据集。结合了单独的抗蠕虫治疗和营养补充，以研究它们如何影响寄生虫强度。

研究问题

不同的治疗方法如何影响寄生虫的活动，这种活动是否受到空间模式的影响？

研究人员在四个网格中捕获宿主，其中两个网格补充了高质量的食物。一些个体接受了抗寄生虫化合物的治疗，而另一些则没有。在每次捕获时，都会记录诸如身体状况等表型数据，并计算寄生虫数量。

1.导入数据

让我们导入数据。

Hosts <- read.csv(paste0(Root, "/Hostures.csv"), header = T)

检查数据，查看各列内容。

phen <- c("Grid", "ID", "Easting", "Northing") # 包含我们需要的空间信息的基础列  
 
	  
 
	resp <- "Parasite.count" # 响应变量  
 
	  
 
	covar <- c("Month", # 采样的月份  
 
	           "Sex", # 性别  
 
	           "Smi", # 身体状况  
 
	TestHosts <- na.omit(Hosts[, c(phen, resp, covar)]) # 去除NA值，选择成年人  
 
	# 我们使用[]进行子集划分，并仅提取特定的列  
 
	  
 
	# 将变量转换为因子  
 
	TestHosts$Month <- as.factor(TestHosts$Month)  
 
	TestHosts$Grid <- as.factor(TestHosts$Grid)

INLA的工作原理与其他许多统计分析包类似，如lme4或MCMCglmm。如果您在这些包中运行相同的简单模型，应该会得到类似的结果。

在空间中绘制采样位置。

# 绘制采样位置图  
 
	(samp_locations <- ggplot(TestHosts, aes(Easting, Northing)) +   
 
	
	  labs(colour = "Grid"))

不同个体在不同网格中被捕获的频率如何？

table(with(TestHosts, tapply(Grid, ID, function(x) length(unique(x)))))

似乎个体倾向于停留在同一个网格中。

---

点击标题查阅往期内容

R语言用贝叶斯层次模型进行空间数据分析

左右滑动查看更多

01

02

03

04

2. 在INLA中进行模型选择

首先，我们将使用所有我们认为会影响数据的协变量进行完整的分析。正如我之前所说，您可以将INLA像其他建模包一样使用，但在这里我将使用公式规范来指定模型。

# 运行模型  
 
	IM0.1  <- inla(Parasite.count ~ Month + Sex + Smi + Supp.corrected + Treated,   
 
	               family = "nbinomial", # 指定分布族。
 
 
 
	  
 
	# 然后添加ID随机效应 ####  
 
	  
 
	
	                          paste(covar, collapse = " + "),   
 
	                          " +  f(ID, model = 'iid')")) # 这是如何包含典型的随机效应的。

接下来，我们将可视化模型的结果。我们将绘制效应大小以及围绕它们的可信区间。

 这个模型可能包含了过多的解释变量。让我们进行模型选择，以移除那些不重要的协变量。

这涉及到逐一移除协变量，并观察这如何根据模型的偏差信息准则（DIC，一种类似于赤池信息准则（AIC）的贝叶斯度量）来改变模型拟合度。

我们可以将这个函数应用于我们的数据，看看应该在模型中包含哪些变量。

最终我们移除了身体条件和食物补充，而治疗、性别和月份保留在最终模型中。提醒一下，INLA中没有P值。变量的重要性或显著性可以通过检查它们的2.5%和97.5%后验估计与零的重叠程度来推断。通过绘图可以更容易地进行这种检查。比起仅查看模型估计，我更喜欢使用DIC来比较变量对模型拟合的贡献。

关于模型选择的详细阐述

为了检查空间自相关的重要性，我们接下来查看具有不同随机效应结构的一系列竞争模型的DIC。鉴于我的采样地点布局，我决定有几种可能的方式可以在这个数据集中编码空间自相关。

1. 在整个研究期间和研究区域内，空间自相关保持恒定（空间，1个网格）。

2. 在整个研究区域内，空间自相关保持恒定，但在研究期间有所变化（时空，X个网格）。

3. 空间自相关在每个网格内变化，以忽略网格之间的空间模式（空间，4个网格）。

3. INLA模型的组成部分

到目前为止的设置涉及使用相当简单的模型公式。下一步通常是人们感到困惑的地方，因为它涉及更独特于INLA且难以分解的模型设置。

INLA之所以计算效率高，是因为它使用SPDE（随机偏微分方程）来估计数据的空间自相关。这涉及使用离散采样位置的“网格”进行插值，以估计空间中的连续过程。

4. 设置空间分析

设置网格

MeshA <- inla.mesh.2d(jitter(Locations), max.edge = c(20, 40))

网格有几个重要的方面。三角形的大小（由max.edge和cutoff的组合决定）决定了方程将如何精确地适应数据。使用较小的三角形会增加精度，但也会成倍增加计算量。通常，网格函数会自动创建一个类似于网格A的网格，其中更靠近的采样位置会产生较小的三角形。在这个数据集中，采样位置分布得如此均匀，以至于我不得不通过抖动它们来在网格A中显示这一点。在探索/设置初步分析时，使用类似网格B的网格。在报告中用于论文的分析时，使用类似网格C的网格。请注意边缘，并尝试在采样区域周围留出一些空间，以便INLA进行估计。可以将边缘的三角形做得更大一些，以减少计算量。

# 创建A矩阵  
 
	  
 
	HostsA <- inla.sc = Locations) # 创建A矩阵

A矩阵与模型矩阵和随机效应组合成一种称为堆栈（stack）的格式。

# 创建模型矩阵 ####  
 
	  
 
	X0 <- model.m -1 + ", paste(Finalcovar, collapse = " + "))), data = TestHosts) # 使用最终模型选择公式（无响应变量）创建模型矩阵。
	  
	X <- as.data.frame(X0[,-which(colnames(X0)%in%c("Month7"))]) # 转换为数据框。如果适用，则消除第一个分类变量的基础水平（您将在下面手动指定截距）  
	  
	head(X)  
	  data = list(y = TestHosts[,resp]), # 指定响应变量  
	    
	  A = list(1, 1, 1, HostsA), # 随机效应和固定效应的乘法因子向量  
	
	    Intercept = rep(1, N), # 指定手动截距！
	      
	    X = X, # 附加模型矩阵  
	      
	    ID = TestHosts$ID, # 插入任何随机效应的向量

堆栈（stack）包括以下内容：

1. 响应变量（编码为“y”）

2. 乘法因子向量。这通常是一系列1（用于截距、随机效应和固定效应），后跟您之前指定的空间A矩阵。

3. 效应。您需要分别指定截距、随机效应、模型矩阵和spde。要记住的是，堆栈的第二部分（乘法因子）的组件与第三部分（效应）的组件相关。添加效应需要在乘法因子中添加另一个1（在正确的位置） 。

添加随机效应？将其添加到效应中，并在A向量中添加一个1。

假设我试图添加网格的随机效应：

data = list(y = TestHosts[,resp]), # 指定响应变量  
 
	    
 
	  A = list(1, 1, 1, HostsA), # 随机效应和固定效应的乘法因子向量  
 
	    
 
 
	    Intercept = rep(1, N), # 指定手动截距！
 
	      
 
	    X = X, # 附加模型矩阵  
 
	      
 
	    ID = TestHosts$ID, # 插入任何随机效应的向量  
 
	
	    w = w.Host)) # 省略其他内容
 
 
 
	  data = list(y = TestHosts[,resp]), # 指定响应变量  
 
	    
 
	  A = list(1, 1, 1, 1, HostsA), # 随机效应和固定效应的乘法因子向量  
 
	    
 
	      
 
	    Intercept = rep(1, N), # 指定手动截距！
 
	      
 
	    X = X, # 附加模型矩阵  
 
	      
 
	    ID = TestHosts$ID, # 插入任何随机效应的向量  
 
 
	    w = w.Host)) # 省略其他内容

运行模型

在运行模型之前，请确保所有需要的组件都已正确添加到堆栈中，并且乘法因子向量与效应列表中的组件相匹配。这确保了模型能够正确地解释和处理您的数据。为了完整性，我们尝试三个不同的模型：

• 仅包含固定效应，

• 固定效应 + ID 随机效应，

• 固定效应 + ID + SPDE 随机效应。

", paste0(colnames(X), collapse = " + "), " +  f(ID, model = 'iid') + f(w, model = Hosts.spde)"))  
	  
	IM1 <- inla(f1, # 基础模型（无随机效应）  
	       
	            control.predictor = list(A = inla.stack.A(StackHost))  
	)  
	  
	IM2 <- inla(f2, # f1 + ID随机效应  
	        
	            control.compute = list(dic = TRUE),  
	            control.predictor = list(A = inla.stack.A(StackHost))  
	)  
	  
	IM3 <- inla(f3, # f2 + SPDE随机效应  
	            family = "nbinomial",  
	            data = inla.stack.data(StackHost),

绘制空间场

r复制代码
	ggFielette = "Blues")

空间中的自相关在什么范围内衰减？具有较大kappa（逆范围）参数的INLA模型在空间上变化非常快。那些模型能够更精细地捕捉空间结构的细节，但这也意味着它们对数据的拟合更加敏感，可能会更容易过拟合。较小的kappa值意味着自相关在空间上衰减得更慢，模型会更加平滑，但可能无法充分捕捉数据的局部变化。选择适当的kappa值需要在模型的拟合度和复杂度之间找到平衡。

查看范围

# 该函数接收（一系列）模型，并在用户定义的范围内绘制空间自相关的衰减情况  
 
	  
 
	# 让我们在我们的模型上试试这个函数 ###  
 
	  
 
	# 定义合理的最大范围：研究区域在东西方向上是80个单位宽，所以让我们设定：
 
	  
 
	Maxrange <- 40  
 
	  
 
	INLARange(axrange)

然而，能够可视化空间模式并不一定意味着空间自相关对模型有实质性的影响，而且范围并不对应于自相关的重要性！为了调查这一点，我们需要查看模型拟合度。这些模型的DIC是如何比较的？

sapply(Spat(f) f$dic$dic)

这很难直观地表示

# 让我们在我们的数据上试试这个函数 ####  
 
	  
 
	INLADICFi c("Base", "IID", "SPDE"))

看起来空间自相关并没有按照我们编码的方式影响这些数据！进行这项研究的任何人都可以继续他们的工作，而无需再担心空间自相关。但是，我们有一些预期，认为这里可能存在其他类型的空间自相关。

5. 修改和指定空间INLA模型

季节模型

现在，如果空间场随季节变化怎么办？我们需要以不同的方式指定A矩阵、SPDE和模型，以产生几个不同的组。

# 指定一组新的SPDE组件 ####  
 
	  
 
	Groups <- "Month"  
 
	  
 
 
	HostA2 <- inla.spde.make.A(Mesh, # 保持不变  
 
	                           loc = Locations, # 保持不变  
 
	                        [,Groups])), # 这必须是一个从1开始计数的数值。如果组变量是因子，这将在默认情况下发生。
 
	              
 
	  data = list(y = TestHosts[,resp]), # 保持不变  
 
	    
 
	  A = list(1, 1, 1, HostA2), # 将A矩阵更改为新的矩阵  
 
	  
	    Intercept = rep(1, N), # 保持不变  
 
	
	      
 
	    w = w.Host2)) # 修改处

现在既然已经指定了这些，让我们运行模型。

w, model = Hosts.spde,   
 
	group = w.group,                           # 这部分是新的
 
 
 
	            family = "nbinomial",  
 
	            data = inla.stack.data(StackHost2), # 别忘了更改堆栈  
 
	     
 
	SpatialHostList[[4]] <- IM4

这段代码展示了如何根据季节变化来修改空间模型，并运行新的模型。既然模型已经运行完毕，让我们来绘制结果图吧

ggField(IM4, Mesh, Groups = NGroups) + # 注意groups参数，使用唯一月份的数量。
 
 
	  facls), ncol = 3) # 手动设置这个以更改分面标签

6. 时空分析

有一种更快的方法可以将空间场分割成组，使用repl而不是将它们分成组并通过iid模型连接。在上面的模型中，我们假设每月的空间场彼此之间是完全不相关的。然而，我们可以使用“可交换”模型来强制它们之间建立相关性，并推导出字段之间的rho相关性。

control.predictor = list(A = inla.stack.A(StackHost2))
)
 
SpatialHostList[[5]] <- IM5
 
	# 与上面的函数相同
 
	  
 
	INLADICFig(SpatiID", "SPDE", "SPDE2", "SPDE3"))

ggField(IM5, Mesh, Groups = NGroups) + # 注意groups参数，使用唯一月份的数量。
 
	  scal")

这段内容主要介绍了时空分析中的不同模型，并强调了使用AR1模型时，时间上更近的字段将有更高的相关性。同时，提供了绘制不同模型比较图（通过DIC，即偏差信息准则）和绘制空间场地图的R代码示例。通过这些分析，可以更好地理解和比较不同模型在时空数据上的表现，并选择最适合的模型进行后续研究。如果你对AR1模型感兴趣，并认为它适用于你的数据，你可以尝试在自己的数据上运行相关代码，以获得更准确的结果。

7. 格内模型

为了完整性起见，让我们尝试使用repl而不是group。简而言之，这稍微快一些，但只能在你不指定字段之间的链接时使用。

我们将研究将研究区域限制为四个形状相同的网格是否会比在乡村的大量空白空间（从未捕获到宿主）中提高拟合度。

Group2 = "Grid"  
 
	  
 
	# 重新编码数据，将每个网格内的坐标转换为相对于网格最小坐标的相对坐标  
 
	TestHosts$Easting2 <- TestHosrthing - with(TestHosts, tapply(Northing, Grid, min))[TestHosts$Grid]  
 
	  
 
	# 创建新的位置矩阵  
 
	Locations2 = cbind(TestH$Northing2)  
 
	  
 
	# 使用新的位置矩阵创建网格  
 
	Mesh2 <- inla.m = c(20, 40))  
 
	  
 
	# 计算网格数量  
 
	NGroup2 <- les[,Group2]))  
 
	  
 
	# 定义SPDE模型  
 
	Hosts.spde2 , prior.range = c(10, 0.5), prior.sigma = c(.5, .5))  
 
	  
 
	# 创建A矩阵  
 
	HostA3 <- inla.sp
	                           n.repl = NGroup2)  
 
	  
 
	# 创建权重索引  
 
	w.Host3 <- inlapde,  
 
	  n.repl = NGroup2)  
 
	  
 
	# 创建堆叠数据  
 
	StackHo
	  A = list(1, 1, 1, HostA3), # 更新A矩阵  
 
	  effects = list(  
 
	  
 
	# 定义模型公式  
 
	f6 = as.formula(pastde2, replicate = w.repl)"))  
	  
	# 拟合模型  
	IM6 <- inla(f6,tack.A(StackHost3))  
	)  
	  
	# 将模型添加到列表中  
	SpatialHostList[[6]] <- IM6

这个模型是否更好地拟合了数据？

要确定模型是否更好地拟合了数据，我们可以查看模型的偏差信息准则（DIC）和其他拟合统计量，比如对数似然值。此外，我们还可以通过查看残差图、预测值与实际观测值的对比等来进行可视化检查。但是，仅凭代码本身无法直接判断模型是否拟合得更好，需要进一步的计算和可视化分析。

INLADICFiase", "IID", "SPDE", "SPDE2", "SPDE3", "GridSPDE"))

没有。

TestHost
  ggsave("Fields6.png", un120, height = 100, dpi = 300)

一样的。

总结

拟合度最好的模型是SPDE 3（模型5）。该模型具有每个月不同的空间字段，字段之间具有相关性。然而，与非空间模型相比，这种模型仅稍微提高了拟合度。

EfxpE", "SPDE2", "SPDE3", "GridSPDE"))

---

点击文末“阅读原文”

获取全文完整代码数据资料。

本文选自《R语言贝叶斯INLA空间自相关、混合效应、季节空间模型、SPDE、时空分析野生动物数据可视化》。

点击标题查阅往期内容

R语言贝叶斯广义线性混合（多层次/水平/嵌套）模型GLMM、逻辑回归分析教育留级影响因素数据

R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例

用SPSS估计HLM多层（层次）线性模型模型

R语言用线性混合效应（多水平/层次/嵌套）模型分析声调高低与礼貌态度的关系

R语言LME4混合效应模型研究教师的受欢迎程度

R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model分析藻类数据实例

R语言混合线性模型、多层次模型、回归模型分析学生平均成绩GPA和可视化

R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例

R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据

R语言 线性混合效应模型实战案例

R语言混合效应逻辑回归（mixed effects logistic）模型分析肺癌数据

R语言如何用潜类别混合效应模型（LCMM）分析抑郁症状

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言建立和可视化混合效应模型mixed effect model

R语言LME4混合效应模型研究教师的受欢迎程度

R语言 线性混合效应模型实战案例

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言如何解决线性混合模型中畸形拟合(Singular fit)的问题

基于R语言的lmer混合线性回归模型

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

R语言分层线性模型案例

R语言用WinBUGS 软件对学术能力测验（SAT）建立分层模型

使用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

SPSS中的多层（等级）线性模型Multilevel linear models研究整容手术数据

用SPSS估计HLM多层（层次）线性模型模型