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

在空间中绘制采样位置。

	  

	# 绘制采样位置图  

	(samp_locations <- ggplot(TestHosts, aes(Easting, Northing)) +   

	
	  labs(colour = "Grid"))

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)) # 省略其他内容

 ", 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值意味着自相关在空间上衰减得更慢，模型会更加平滑，但可能无法充分捕捉数据的局部变化。

查看范围

	# 该函数接收（一系列）模型，并在用户定义的范围内绘制空间自相关的衰减情况  

	  

	# 让我们在我们的模型上试试这个函数 ###  

	  

	# 定义合理的最大范围：研究区域在东西方向上是80个单位宽，所以让我们设定：  

	  

	Maxrange <- 40  

	  

	INLARange(axrange)

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

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"))