R中rhohat函数分析异常：金属斧钺分布与环境因子关联咨询

针对你用R语言rhohat分析丹麦、英格兰/威尔士及荷兰金属斧钺遗存与环境协变量关联时遇到的图形异常、结果和QGIS提取信息不符的问题，我整理了以下排查方向、概念解释和替代工具：

一、图形异常与结果不符的可能原因 #

1. 栅格化分辨率不合理
   你创建栅格时用了raster(ncol=180, nrow=180)，这个分辨率可能过于粗糙，把原本精细的空间信息“抹平”了，导致和QGIS里的精细提取结果不一致。建议根据你的数据投影单位（比如米）设置更合理的分辨率，比如：

   r <- raster(extent(soil))
   res(r) <- 500  # 假设投影单位是米，设置500米分辨率
   

2. 协变量数据预处理不到位
   • 土壤质地里的0（无信息）和9（无矿物质地）属于无效值，rhohat会把它们当作有效类别计算，拉偏结果。建议先过滤这些值：

     soil <- soil[soil$质地字段 != 0 & soil$质地字段 != 9, ]
     

   • 土地覆盖是三位编码（如111），类别过多可能导致每个类别的样本量不足，统计结果不稳定。可以考虑按用地类型合并类别（比如把111、112归为“城市用地”），再做分析。
3. 点模式与协变量的范围不匹配
   检查axeheads.pp这个点模式对象的窗口（window）是否和soiltype_im等协变量栅格的范围完全一致。如果窗口范围更小，会截断部分数据；如果更大，会引入无数据区域，都可能导致结果偏差。可以用window(axeheads.pp) <- extent(soiltype_raster)来统一范围。
4. rasterize字段匹配错误
   确认'FAO85LV1'是soil shapefile里的正确字段名，有没有大小写或者拼写错误？如果字段名不对，栅格化后的值会全部为空或者错误，自然和QGIS结果不符。

二、rhohat图Y轴数值的含义 #

rhohat计算的是相对风险（relative risk），Y轴数值代表：在该协变量类别下，金属斧钺遗存的密度是全局平均密度的多少倍。比如Y=1.8意味着这个类别里的遗存密度是全局平均的1.8倍，Y=0.6则是平均密度的60%。

三、其他探究遗存与环境因子关联的函数 #

除了rhohat，这些工具能帮你更深入分析：

• ppm（点过程模型）：直接把环境协变量作为解释变量，量化每个协变量对遗存分布的影响，还能做显著性检验，适合建立统计模型：

  # 假设landcover_im、ph_im都是处理好的im对象
  fit <- ppm(axeheads.pp ~ soiltype_im + landcover_im + ph_im)
  summary(fit)
  

• density.ppp：计算遗存的核密度，结合协变量分层展示，直观对比不同环境类别下的密度差异。
• quadrat.test：通过划分空间方格，检验不同协变量类别下的遗存分布是否符合随机分布，适合初步的显著性验证。
• lppm（对数高斯Cox过程）：适合处理存在空间自相关的点模式数据，能同时考虑环境协变量和空间自相关的影响，结果更可靠。

内容的提问来源于stack exchange，提问作者Linda Marie