简介这份资源面向古生物学、地球科学与气候学领域的科研人员、学生及数据分析从业者围绕有孔虫化石形态测量数据的气候重建与群落结构评估展开。内容基于377个地质时间点、每点约2000个个体的测量记录涵盖大小、面积、形状因子、伸长率、球形度、周长与灰度值等参数用Python实现数据加载与预处理、基本统计、时间序列分解、多变量分析及环境关联分析并扩展年龄模型验证、自动化参数提取、环境数据整合与分类单元变化分析帮助读者掌握从形态特征到环境变化的完整分析链路。资源包共1个PDF文件约516KB内含可运行代码与逐段解释便于对照复现。目前已有59人学习。读者可借此理解群落结构变化、物种数量波动、不同参数的选择压力差异、数据周期性及形态与环境关联适合具备Pandas、NumPy、Matplotlib基础者进阶实践。1. 有孔虫化石形态测量做气候重建从壳体尺寸到古温度的一整条链路有孔虫是海洋微体古生物里最常被拿来做古气候重建的载体之一。它的壳体保存量大、演化快、对环境敏感壳径、房室数、螺旋比例、孔隙密度这些形态指标往往和海水温度、盐度、溶解氧存在可量化的关联。把「有孔虫化石形态测量」和「气候重建」放在一起本质是干一件事用 pandas 整理测量表、用 numpy 做矩阵运算、用 scipy 做统计检验、用 matplotlib 和 seaborn 把群落结构和环境因子的关系画出来最终得到一条能解释古温度的曲线。这套流程适合做微体古生物、第四纪地质、古海洋学的从业者也适合手里已经有一批壳体测量数据、却卡在「怎么从 Excel 走到可发表图件」的人。下面按数据清洗、形态指标构建、群落结构评估、环境关联建模、避坑、进阶验证六段推进代码全部可复现。2. 数据清洗与形态指标构建把测量表变成可分析的 DataFrame2.1 测量数据的典型结构与 pandas 读取有孔虫形态测量的原始数据通常来自体视显微镜下的逐个体测量字段包括样品编号、深度、物种名、壳长、壳宽、壳厚、房室数、螺旋方向、孔隙数等。常见格式是 Excel 或制表符分隔的 txt。用 pandas 读取时第一件事是确认分隔符和缺失值标记很多测量表用-999或空字符串表示未测。import pandas as pd import numpy as np # 读取 Excel 测量表指定缺失值标记 df pd.read_excel( foram_measurements.xlsx, sheet_nameraw, na_values[, -999, NA, na], dtype{sample_id: str, species: str} ) # 统一列名去掉空格和大小写差异 df.columns [c.strip().lower().replace( , _) for c in df.columns] # 查看基本结构和缺失情况 print(df.shape) print(df.dtypes) print(df.isna().sum().sort_values(ascendingFalse).head(10))这段代码的逻辑是先解决「读进来」的问题。na_values把多种缺失标记统一成 NaN避免后续统计把 -999 当成真实壳长。列名规范化是为了后面写公式时不用反复处理空格。isna().sum()按缺失数量排序能快速定位哪些形态指标测得不全。参数上dtype强制 sample_id 为字符串防止样品编号被 pandas 推断成数字后丢失前导零。2.2 形态指标的派生与单位统一原始测量往往是像素或目镜格值需要先换算成毫米。更关键的是派生指标壳长宽比、壳面积近似、房室密度、孔隙密度。这些派生指标比原始尺寸更能反映环境压力因为不同个体的绝对大小受生长阶段影响而比例指标相对稳定。# 假设原始单位为目镜格换算系数 0.01 mm/grid GRID_TO_MM 0.01 df[length_mm] df[length_grid] * GRID_TO_MM df[width_mm] df[width_grid] * GRID_TO_MM df[thickness_mm] df[thickness_grid] * GRID_TO_MM # 派生形态指标 df[lw_ratio] df[length_mm] / df[width_mm] df[aspect_ratio] df[length_mm] / df[thickness_mm] df[chamber_density] df[chamber_count] / df[length_mm] df[pore_density] df[pore_count] / (df[length_mm] * df[width_mm]) # 用 numpy 做对数变换压缩量纲差异 df[log_length] np.log1p(df[length_mm]) df[log_lw_ratio] np.log1p(df[lw_ratio]) # 检查派生指标是否有无穷值或异常 print(df[[lw_ratio, aspect_ratio, chamber_density, pore_density]].describe())逻辑说明np.log1p对 1x 取对数避免壳长接近 0 时出现负无穷。派生指标里pore_density用面积做分母比用长度更合理因为孔隙分布在二维壳面上。参数上换算系数必须和显微镜标定一致不同实验室的目镜格值不同这一步错了后面全错。describe()用来快速看极值如果 lw_ratio 出现 50 以上多半是宽度测量错误或单位没统一。2.3 物种名标准化与测量偏差剔除物种名拼写不一致是微体古生物数据里的经典问题。同一个种可能有异名、缩写、大小写差异。用 pandas 的映射表统一再用 numpy 做基于分位数的异常剔除比直接删 3σ 更稳健因为形态数据往往右偏。# 物种名映射表 species_map { g. ruber: Globigerinoides ruber, g.ruber: Globigerinoides ruber, globigerinoides ruber: Globigerinoides ruber, n. dutertrei: Neogloboquadrina dutertrei, p. obliquiloculata: Pulleniatina obliquiloculata } df[species_std] ( df[species].str.strip().str.lower().map(species_map) ) # 对每个物种分别做 1% 和 99% 分位剔除 def trim_by_species(group): lo group[length_mm].quantile(0.01) hi group[length_mm].quantile(0.99) return group[(group[length_mm] lo) (group[length_mm] hi)] df_clean df.groupby(species_std, group_keysFalse).apply(trim_by_species) print(清洗前:, df.shape, 清洗后:, df_clean.shape)这里用groupby加apply按物种分别剔除极端值是因为不同物种的壳长量级差异很大混在一起做分位会把大个体的正常值误删。参数上1% 和 99% 是常用阈值如果样品量少于 30 个建议放宽到 2% 和 98%否则会删掉太多有效数据。清洗后要检查每个物种剩余数量少于 10 个的物种在后续群落分析里要谨慎处理。3. 群落结构评估从丰度矩阵到多样性指数3.1 构建样品-物种丰度矩阵群落结构评估的第一步是把长表转成宽表行是样品列是物种值是丰度或相对丰度。pandas 的pivot_table能直接做但要注意重复测量需要先聚合。# 按样品和物种统计个体数 abundance ( df_clean.groupby([sample_id, species_std]) .size() .reset_index(namecount) ) # 转成宽表缺失填 0 abund_matrix abundance.pivot_table( indexsample_id, columnsspecies_std, valuescount, fill_value0 ) # 计算相对丰度 rel_abund abund_matrix.div(abund_matrix.sum(axis1), axis0) print(abund_matrix.shape) print(rel_abund.head())逻辑上groupby().size()统计每个样品里每个物种的个体数pivot_table把物种变成列。fill_value0很重要因为没出现的物种在群落分析里就是零丰度不能留 NaN。div(..., axis0)按行归一化得到相对丰度。参数上如果原始数据已经是相对丰度跳过 count 步骤直接 pivot 即可。3.2 多样性指数与群落参数计算常用指数包括 Shannon、Simpson、Pielou 均匀度以及物种丰富度。scipy 本身没有直接算这些的函数但用 numpy 可以几行实现。这些指数反映群落对环境的响应比如温度升高时浮游有孔虫多样性可能下降。def shannon_index(row): p row[row 0] / row.sum() return -np.sum(p * np.log(p)) def simpson_index(row): p row[row 0] / row.sum() return 1 - np.sum(p ** 2) def pielou_evenness(row): s (row 0).sum() if s 1: return np.nan return shannon_index(row) / np.log(s) diversity pd.DataFrame(indexabund_matrix.index) diversity[richness] (abund_matrix 0).sum(axis1) diversity[shannon] abund_matrix.apply(shannon_index, axis1) diversity[simpson] abund_matrix.apply(simpson_index, axis1) diversity[pielou] abund_matrix.apply(pielou_evenness, axis1) print(diversity.describe())shannon_index里先过滤零丰度再归一化避免 log(0)。pielou_evenness在物种数小于等于 1 时返回 NaN因为均匀度在单物种群落里没有定义。参数上Shannon 用自然对数如果要用 log2 需在公式里改。这些指数对样品个体数敏感建议在比较前做个体数标准化或者用 rarefaction但 rarefaction 超出本文范围常见做法是只保留个体数大于 50 的样品。3.3 用 seaborn 画群落结构热图与箱线图群落结构可视化首选热图看样品-物种相对丰度箱线图看多样性指数沿深度或温度的变化。seaborn 的heatmap和boxplot能直接接 DataFrame。import matplotlib.pyplot as plt import seaborn as sns fig, axes plt.subplots(1, 2, figsize(14, 6)) # 热图相对丰度 sns.heatmap( rel_abund.T, cmapYlGnBu, cbar_kws{label: Relative abundance}, axaxes[0] ) axes[0].set_xlabel(Sample) axes[0].set_ylabel(Species) # 箱线图Shannon 指数按深度分组 diversity[depth_bin] pd.cut( df_clean.groupby(sample_id)[depth_m].mean(), bins5 ) sns.boxplot( datadiversity, xdepth_bin, yshannon, paletteSet2, axaxes[1] ) axes[1].set_xlabel(Depth bin (m)) axes[1].set_ylabel(Shannon index) axes[1].tick_params(axisx, rotation30) plt.tight_layout() plt.savefig(community_structure.png, dpi300)逻辑说明热图转置是为了让样品在 x 轴、物种在 y 轴符合阅读习惯。pd.cut把连续深度分成 5 段方便箱线图分组。参数上cmapYlGnBu对丰度数据视觉友好dpi300满足投稿要求。注意 seaborn 的boxplot在新版本里palette不加hue会警告可以改成color或加hue并设legendFalse。4. 环境关联研究相关分析、回归与古温度转换函数4.1 形态指标与环境因子的相关矩阵把形态指标、多样性指数和实测环境因子温度、盐度、溶解氧放在一起算 Spearman 相关比 Pearson 更稳健因为形态数据不一定是正态。scipy 的spearmanr可以逐对算但更高效的是用 pandas 的corr(methodspearman)。# 假设 env_df 是样品级环境因子表index 为 sample_id env_df pd.read_csv(environment.csv, index_colsample_id) # 合并形态均值、多样性、环境 morph_mean df_clean.groupby(sample_id)[ [length_mm, lw_ratio, chamber_density, pore_density] ].mean() combined morph_mean.join(diversity[[shannon, richness]]).join(env_df) corr_matrix combined.corr(methodspearman) # 只保留和温度相关的列便于阅读 temp_corr corr_matrix[[SST_annual, SST_summer, salinity]].dropna(howall) print(temp_corr.sort_values(SST_annual, ascendingFalse))逻辑上先按样品聚合形态均值再和多样性、环境表做 join。corr(methodspearman)一次性算所有数值列的秩相关。参数上环境因子列名要和实际数据一致SST 通常指海表温度。如果相关矩阵里出现 NaN说明某些样品环境数据缺失需要在 join 时用howinner或后续 dropna。4.2 用 scipy 做线性回归与显著性检验相关不等于因果要建立形态指标到温度的转换函数需要用scipy.stats.linregress或statsmodels。这里用 scipy 做单变量回归输出斜率、截距、R² 和 p 值。from scipy import stats # 以 pore_density 预测 SST_annual 为例 x combined[pore_density].values y combined[SST_annual].values mask ~np.isnan(x) ~np.isnan(y) slope, intercept, r_value, p_value, std_err stats.linregress(x[mask], y[mask]) print(f斜率: {slope:.3f}, 截距: {intercept:.3f}) print(fR²: {r_value**2:.3f}, p值: {p_value:.4f}, 标准误: {std_err:.3f}) # 画散点加回归线 fig, ax plt.subplots(figsize(7, 5)) ax.scatter(x[mask], y[mask], alpha0.6, edgecolork, labelSamples) x_line np.linspace(x[mask].min(), x[mask].max(), 100) ax.plot(x_line, slope * x_line intercept, r-, labelRegression) ax.set_xlabel(Pore density (pores/mm²)) ax.set_ylabel(Annual SST (°C)) ax.legend() plt.savefig(pore_density_sst.png, dpi300)逻辑说明linregress返回五个值R² 是 r_value 的平方。mask 用来排除缺失值。参数上alpha0.6让散点有透明度重叠点可见。如果 p 值大于 0.05说明这个形态指标和温度的关系不显著需要换指标或考虑多变量模型。注意module numpy has no attribute trapz这类报错通常出现在 numpy 2.0 之后trapz改名为trapezoid如果代码里用了旧函数需要改。4.3 多变量转换函数与验证单变量往往不够古温度转换函数常用多变量线性模型或主成分回归。这里用 numpy 做最小二乘手动实现多元回归便于理解参数含义。# 构造设计矩阵加入截距项 feature_cols [pore_density, chamber_density, lw_ratio, shannon] X combined[feature_cols].values y combined[SST_annual].values mask ~np.isnan(X).any(axis1) ~np.isnan(y) X, y X[mask], y[mask] X_design np.column_stack([np.ones(X.shape[0]), X]) # 最小二乘解 beta, residuals, rank, sv np.linalg.lstsq(X_design, y, rcondNone) y_pred X_design beta ss_res np.sum((y - y_pred) ** 2) ss_tot np.sum((y - y.mean()) ** 2) r2 1 - ss_res / ss_tot print(系数:, beta) print(R²:, round(r2, 3))逻辑说明np.column_stack加一列 1 作为截距。np.linalg.lstsq解最小二乘返回的 beta 第一个是截距后面是各特征系数。R² 手动算便于理解。参数上特征列之间如果高度相关系数会不稳定建议先看相关矩阵或者用 PCA 降维。验证时最好留出独立样品做交叉验证避免过拟合。5. 避坑与排查有孔虫形态测量分析里最容易翻车的五件事5.1 现象pandas 读取 Excel 后数值列变成 object原因Excel 里混有文本备注或空格pandas 推断为 object。解决读取时加dtype指定或读完后用pd.to_numeric(errorscoerce)强制转换再检查 NaN 数量。5.2 现象numpy 报module numpy has no attribute trapz原因numpy 2.0 移除了trapz改名为trapezoid。解决升级代码里的函数名或者降级 numpy 到 1.26。检查环境用pip show numpy。5.3 现象seaborn 热图物种顺序混乱看不出梯度原因默认按列名排序不是按丰度或生态意义排序。解决先算物种平均丰度用reindex按丰度排序或者用linkage做聚类排序。5.4 现象回归 R² 很高但预测新样品偏差大原因过拟合特征太多或样品太少。解决减少特征做留一交叉验证看预测残差是否随温度系统偏移。5.5 现象多样性指数在个体数少的样品里异常低原因Shannon 指数对个体数敏感少于 30 个个体的样品不可靠。解决过滤个体数小于 50 的样品或做 rarefaction 标准化。6. 进阶验证用交叉验证和残差诊断判断转换函数能不能用走到这一步手里已经有一条从测量表到古温度转换函数的链路。但能不能用取决于验证。我一般会做两件事留一交叉验证和残差诊断。留一交叉验证用 numpy 手动实现不依赖 sklearn便于理解每一步。def loo_cv(X, y): n len(y) preds np.zeros(n) for i in range(n): mask np.ones(n, dtypebool) mask[i] False X_train np.column_stack([np.ones(mask.sum()), X[mask]]) beta np.linalg.lstsq(X_train, y[mask], rcondNone)[0] preds[i] np.dot(np.append(1, X[i]), beta) return preds preds loo_cv(X, y) rmse np.sqrt(np.mean((y - preds) ** 2)) print(LOO RMSE:, round(rmse, 3)) # 残差诊断图 residuals y - preds fig, axes plt.subplots(1, 2, figsize(12, 5)) axes[0].scatter(preds, residuals, alpha0.6) axes[0].axhline(0, colorr, linestyle--) axes[0].set_xlabel(Predicted SST) axes[0].set_ylabel(Residual) axes[1].hist(residuals, bins15, edgecolork) axes[1].set_xlabel(Residual) plt.savefig(loo_diagnostics.png, dpi300)逻辑说明每次留一个样品用其余样品拟合预测留出的那个。RMSE 衡量整体误差残差图看是否有系统偏差。参数上如果残差随预测值呈喇叭形说明方差不齐可能需要对数变换。如果残差有趋势说明模型漏了非线性项。验证指标可接受范围说明LOO RMSE小于 1.5°C古温度重建常用阈值R²大于 0.6低于此值解释力不足残差均值接近 0系统偏差检查残差与预测相关不显著p 0.05 最好最后说个血泪经验有孔虫形态测量里物种鉴定的一致性比任何统计方法都重要。我见过同一批样品两个人鉴定物种名差 20%后面所有多样性指数和转换函数全废。所以做分析前先花时间统一物种名和测量标准再跑 pandas 和 numpy。希望帮到你。本文还有配套的精品资源点击获取 有孔虫化石形态测量与气候重建：Python数据分析全流程
简介这份资源面向古生物学、地球科学与气候学领域的科研人员、学生及数据分析从业者围绕有孔虫化石形态测量数据的气候重建与群落结构评估展开。内容基于377个地质时间点、每点约2000个个体的测量记录涵盖大小、面积、形状因子、伸长率、球形度、周长与灰度值等参数用Python实现数据加载与预处理、基本统计、时间序列分解、多变量分析及环境关联分析并扩展年龄模型验证、自动化参数提取、环境数据整合与分类单元变化分析帮助读者掌握从形态特征到环境变化的完整分析链路。资源包共1个PDF文件约516KB内含可运行代码与逐段解释便于对照复现。目前已有59人学习。读者可借此理解群落结构变化、物种数量波动、不同参数的选择压力差异、数据周期性及形态与环境关联适合具备Pandas、NumPy、Matplotlib基础者进阶实践。1. 有孔虫化石形态测量做气候重建从壳体尺寸到古温度的一整条链路有孔虫是海洋微体古生物里最常被拿来做古气候重建的载体之一。它的壳体保存量大、演化快、对环境敏感壳径、房室数、螺旋比例、孔隙密度这些形态指标往往和海水温度、盐度、溶解氧存在可量化的关联。把「有孔虫化石形态测量」和「气候重建」放在一起本质是干一件事用 pandas 整理测量表、用 numpy 做矩阵运算、用 scipy 做统计检验、用 matplotlib 和 seaborn 把群落结构和环境因子的关系画出来最终得到一条能解释古温度的曲线。这套流程适合做微体古生物、第四纪地质、古海洋学的从业者也适合手里已经有一批壳体测量数据、却卡在「怎么从 Excel 走到可发表图件」的人。下面按数据清洗、形态指标构建、群落结构评估、环境关联建模、避坑、进阶验证六段推进代码全部可复现。2. 数据清洗与形态指标构建把测量表变成可分析的 DataFrame2.1 测量数据的典型结构与 pandas 读取有孔虫形态测量的原始数据通常来自体视显微镜下的逐个体测量字段包括样品编号、深度、物种名、壳长、壳宽、壳厚、房室数、螺旋方向、孔隙数等。常见格式是 Excel 或制表符分隔的 txt。用 pandas 读取时第一件事是确认分隔符和缺失值标记很多测量表用-999或空字符串表示未测。import pandas as pd import numpy as np # 读取 Excel 测量表指定缺失值标记 df pd.read_excel( foram_measurements.xlsx, sheet_nameraw, na_values[, -999, NA, na], dtype{sample_id: str, species: str} ) # 统一列名去掉空格和大小写差异 df.columns [c.strip().lower().replace( , _) for c in df.columns] # 查看基本结构和缺失情况 print(df.shape) print(df.dtypes) print(df.isna().sum().sort_values(ascendingFalse).head(10))这段代码的逻辑是先解决「读进来」的问题。na_values把多种缺失标记统一成 NaN避免后续统计把 -999 当成真实壳长。列名规范化是为了后面写公式时不用反复处理空格。isna().sum()按缺失数量排序能快速定位哪些形态指标测得不全。参数上dtype强制 sample_id 为字符串防止样品编号被 pandas 推断成数字后丢失前导零。2.2 形态指标的派生与单位统一原始测量往往是像素或目镜格值需要先换算成毫米。更关键的是派生指标壳长宽比、壳面积近似、房室密度、孔隙密度。这些派生指标比原始尺寸更能反映环境压力因为不同个体的绝对大小受生长阶段影响而比例指标相对稳定。# 假设原始单位为目镜格换算系数 0.01 mm/grid GRID_TO_MM 0.01 df[length_mm] df[length_grid] * GRID_TO_MM df[width_mm] df[width_grid] * GRID_TO_MM df[thickness_mm] df[thickness_grid] * GRID_TO_MM # 派生形态指标 df[lw_ratio] df[length_mm] / df[width_mm] df[aspect_ratio] df[length_mm] / df[thickness_mm] df[chamber_density] df[chamber_count] / df[length_mm] df[pore_density] df[pore_count] / (df[length_mm] * df[width_mm]) # 用 numpy 做对数变换压缩量纲差异 df[log_length] np.log1p(df[length_mm]) df[log_lw_ratio] np.log1p(df[lw_ratio]) # 检查派生指标是否有无穷值或异常 print(df[[lw_ratio, aspect_ratio, chamber_density, pore_density]].describe())逻辑说明np.log1p对 1x 取对数避免壳长接近 0 时出现负无穷。派生指标里pore_density用面积做分母比用长度更合理因为孔隙分布在二维壳面上。参数上换算系数必须和显微镜标定一致不同实验室的目镜格值不同这一步错了后面全错。describe()用来快速看极值如果 lw_ratio 出现 50 以上多半是宽度测量错误或单位没统一。2.3 物种名标准化与测量偏差剔除物种名拼写不一致是微体古生物数据里的经典问题。同一个种可能有异名、缩写、大小写差异。用 pandas 的映射表统一再用 numpy 做基于分位数的异常剔除比直接删 3σ 更稳健因为形态数据往往右偏。# 物种名映射表 species_map { g. ruber: Globigerinoides ruber, g.ruber: Globigerinoides ruber, globigerinoides ruber: Globigerinoides ruber, n. dutertrei: Neogloboquadrina dutertrei, p. obliquiloculata: Pulleniatina obliquiloculata } df[species_std] ( df[species].str.strip().str.lower().map(species_map) ) # 对每个物种分别做 1% 和 99% 分位剔除 def trim_by_species(group): lo group[length_mm].quantile(0.01) hi group[length_mm].quantile(0.99) return group[(group[length_mm] lo) (group[length_mm] hi)] df_clean df.groupby(species_std, group_keysFalse).apply(trim_by_species) print(清洗前:, df.shape, 清洗后:, df_clean.shape)这里用groupby加apply按物种分别剔除极端值是因为不同物种的壳长量级差异很大混在一起做分位会把大个体的正常值误删。参数上1% 和 99% 是常用阈值如果样品量少于 30 个建议放宽到 2% 和 98%否则会删掉太多有效数据。清洗后要检查每个物种剩余数量少于 10 个的物种在后续群落分析里要谨慎处理。3. 群落结构评估从丰度矩阵到多样性指数3.1 构建样品-物种丰度矩阵群落结构评估的第一步是把长表转成宽表行是样品列是物种值是丰度或相对丰度。pandas 的pivot_table能直接做但要注意重复测量需要先聚合。# 按样品和物种统计个体数 abundance ( df_clean.groupby([sample_id, species_std]) .size() .reset_index(namecount) ) # 转成宽表缺失填 0 abund_matrix abundance.pivot_table( indexsample_id, columnsspecies_std, valuescount, fill_value0 ) # 计算相对丰度 rel_abund abund_matrix.div(abund_matrix.sum(axis1), axis0) print(abund_matrix.shape) print(rel_abund.head())逻辑上groupby().size()统计每个样品里每个物种的个体数pivot_table把物种变成列。fill_value0很重要因为没出现的物种在群落分析里就是零丰度不能留 NaN。div(..., axis0)按行归一化得到相对丰度。参数上如果原始数据已经是相对丰度跳过 count 步骤直接 pivot 即可。3.2 多样性指数与群落参数计算常用指数包括 Shannon、Simpson、Pielou 均匀度以及物种丰富度。scipy 本身没有直接算这些的函数但用 numpy 可以几行实现。这些指数反映群落对环境的响应比如温度升高时浮游有孔虫多样性可能下降。def shannon_index(row): p row[row 0] / row.sum() return -np.sum(p * np.log(p)) def simpson_index(row): p row[row 0] / row.sum() return 1 - np.sum(p ** 2) def pielou_evenness(row): s (row 0).sum() if s 1: return np.nan return shannon_index(row) / np.log(s) diversity pd.DataFrame(indexabund_matrix.index) diversity[richness] (abund_matrix 0).sum(axis1) diversity[shannon] abund_matrix.apply(shannon_index, axis1) diversity[simpson] abund_matrix.apply(simpson_index, axis1) diversity[pielou] abund_matrix.apply(pielou_evenness, axis1) print(diversity.describe())shannon_index里先过滤零丰度再归一化避免 log(0)。pielou_evenness在物种数小于等于 1 时返回 NaN因为均匀度在单物种群落里没有定义。参数上Shannon 用自然对数如果要用 log2 需在公式里改。这些指数对样品个体数敏感建议在比较前做个体数标准化或者用 rarefaction但 rarefaction 超出本文范围常见做法是只保留个体数大于 50 的样品。3.3 用 seaborn 画群落结构热图与箱线图群落结构可视化首选热图看样品-物种相对丰度箱线图看多样性指数沿深度或温度的变化。seaborn 的heatmap和boxplot能直接接 DataFrame。import matplotlib.pyplot as plt import seaborn as sns fig, axes plt.subplots(1, 2, figsize(14, 6)) # 热图相对丰度 sns.heatmap( rel_abund.T, cmapYlGnBu, cbar_kws{label: Relative abundance}, axaxes[0] ) axes[0].set_xlabel(Sample) axes[0].set_ylabel(Species) # 箱线图Shannon 指数按深度分组 diversity[depth_bin] pd.cut( df_clean.groupby(sample_id)[depth_m].mean(), bins5 ) sns.boxplot( datadiversity, xdepth_bin, yshannon, paletteSet2, axaxes[1] ) axes[1].set_xlabel(Depth bin (m)) axes[1].set_ylabel(Shannon index) axes[1].tick_params(axisx, rotation30) plt.tight_layout() plt.savefig(community_structure.png, dpi300)逻辑说明热图转置是为了让样品在 x 轴、物种在 y 轴符合阅读习惯。pd.cut把连续深度分成 5 段方便箱线图分组。参数上cmapYlGnBu对丰度数据视觉友好dpi300满足投稿要求。注意 seaborn 的boxplot在新版本里palette不加hue会警告可以改成color或加hue并设legendFalse。4. 环境关联研究相关分析、回归与古温度转换函数4.1 形态指标与环境因子的相关矩阵把形态指标、多样性指数和实测环境因子温度、盐度、溶解氧放在一起算 Spearman 相关比 Pearson 更稳健因为形态数据不一定是正态。scipy 的spearmanr可以逐对算但更高效的是用 pandas 的corr(methodspearman)。# 假设 env_df 是样品级环境因子表index 为 sample_id env_df pd.read_csv(environment.csv, index_colsample_id) # 合并形态均值、多样性、环境 morph_mean df_clean.groupby(sample_id)[ [length_mm, lw_ratio, chamber_density, pore_density] ].mean() combined morph_mean.join(diversity[[shannon, richness]]).join(env_df) corr_matrix combined.corr(methodspearman) # 只保留和温度相关的列便于阅读 temp_corr corr_matrix[[SST_annual, SST_summer, salinity]].dropna(howall) print(temp_corr.sort_values(SST_annual, ascendingFalse))逻辑上先按样品聚合形态均值再和多样性、环境表做 join。corr(methodspearman)一次性算所有数值列的秩相关。参数上环境因子列名要和实际数据一致SST 通常指海表温度。如果相关矩阵里出现 NaN说明某些样品环境数据缺失需要在 join 时用howinner或后续 dropna。4.2 用 scipy 做线性回归与显著性检验相关不等于因果要建立形态指标到温度的转换函数需要用scipy.stats.linregress或statsmodels。这里用 scipy 做单变量回归输出斜率、截距、R² 和 p 值。from scipy import stats # 以 pore_density 预测 SST_annual 为例 x combined[pore_density].values y combined[SST_annual].values mask ~np.isnan(x) ~np.isnan(y) slope, intercept, r_value, p_value, std_err stats.linregress(x[mask], y[mask]) print(f斜率: {slope:.3f}, 截距: {intercept:.3f}) print(fR²: {r_value**2:.3f}, p值: {p_value:.4f}, 标准误: {std_err:.3f}) # 画散点加回归线 fig, ax plt.subplots(figsize(7, 5)) ax.scatter(x[mask], y[mask], alpha0.6, edgecolork, labelSamples) x_line np.linspace(x[mask].min(), x[mask].max(), 100) ax.plot(x_line, slope * x_line intercept, r-, labelRegression) ax.set_xlabel(Pore density (pores/mm²)) ax.set_ylabel(Annual SST (°C)) ax.legend() plt.savefig(pore_density_sst.png, dpi300)逻辑说明linregress返回五个值R² 是 r_value 的平方。mask 用来排除缺失值。参数上alpha0.6让散点有透明度重叠点可见。如果 p 值大于 0.05说明这个形态指标和温度的关系不显著需要换指标或考虑多变量模型。注意module numpy has no attribute trapz这类报错通常出现在 numpy 2.0 之后trapz改名为trapezoid如果代码里用了旧函数需要改。4.3 多变量转换函数与验证单变量往往不够古温度转换函数常用多变量线性模型或主成分回归。这里用 numpy 做最小二乘手动实现多元回归便于理解参数含义。# 构造设计矩阵加入截距项 feature_cols [pore_density, chamber_density, lw_ratio, shannon] X combined[feature_cols].values y combined[SST_annual].values mask ~np.isnan(X).any(axis1) ~np.isnan(y) X, y X[mask], y[mask] X_design np.column_stack([np.ones(X.shape[0]), X]) # 最小二乘解 beta, residuals, rank, sv np.linalg.lstsq(X_design, y, rcondNone) y_pred X_design beta ss_res np.sum((y - y_pred) ** 2) ss_tot np.sum((y - y.mean()) ** 2) r2 1 - ss_res / ss_tot print(系数:, beta) print(R²:, round(r2, 3))逻辑说明np.column_stack加一列 1 作为截距。np.linalg.lstsq解最小二乘返回的 beta 第一个是截距后面是各特征系数。R² 手动算便于理解。参数上特征列之间如果高度相关系数会不稳定建议先看相关矩阵或者用 PCA 降维。验证时最好留出独立样品做交叉验证避免过拟合。5. 避坑与排查有孔虫形态测量分析里最容易翻车的五件事5.1 现象pandas 读取 Excel 后数值列变成 object原因Excel 里混有文本备注或空格pandas 推断为 object。解决读取时加dtype指定或读完后用pd.to_numeric(errorscoerce)强制转换再检查 NaN 数量。5.2 现象numpy 报module numpy has no attribute trapz原因numpy 2.0 移除了trapz改名为trapezoid。解决升级代码里的函数名或者降级 numpy 到 1.26。检查环境用pip show numpy。5.3 现象seaborn 热图物种顺序混乱看不出梯度原因默认按列名排序不是按丰度或生态意义排序。解决先算物种平均丰度用reindex按丰度排序或者用linkage做聚类排序。5.4 现象回归 R² 很高但预测新样品偏差大原因过拟合特征太多或样品太少。解决减少特征做留一交叉验证看预测残差是否随温度系统偏移。5.5 现象多样性指数在个体数少的样品里异常低原因Shannon 指数对个体数敏感少于 30 个个体的样品不可靠。解决过滤个体数小于 50 的样品或做 rarefaction 标准化。6. 进阶验证用交叉验证和残差诊断判断转换函数能不能用走到这一步手里已经有一条从测量表到古温度转换函数的链路。但能不能用取决于验证。我一般会做两件事留一交叉验证和残差诊断。留一交叉验证用 numpy 手动实现不依赖 sklearn便于理解每一步。def loo_cv(X, y): n len(y) preds np.zeros(n) for i in range(n): mask np.ones(n, dtypebool) mask[i] False X_train np.column_stack([np.ones(mask.sum()), X[mask]]) beta np.linalg.lstsq(X_train, y[mask], rcondNone)[0] preds[i] np.dot(np.append(1, X[i]), beta) return preds preds loo_cv(X, y) rmse np.sqrt(np.mean((y - preds) ** 2)) print(LOO RMSE:, round(rmse, 3)) # 残差诊断图 residuals y - preds fig, axes plt.subplots(1, 2, figsize(12, 5)) axes[0].scatter(preds, residuals, alpha0.6) axes[0].axhline(0, colorr, linestyle--) axes[0].set_xlabel(Predicted SST) axes[0].set_ylabel(Residual) axes[1].hist(residuals, bins15, edgecolork) axes[1].set_xlabel(Residual) plt.savefig(loo_diagnostics.png, dpi300)逻辑说明每次留一个样品用其余样品拟合预测留出的那个。RMSE 衡量整体误差残差图看是否有系统偏差。参数上如果残差随预测值呈喇叭形说明方差不齐可能需要对数变换。如果残差有趋势说明模型漏了非线性项。验证指标可接受范围说明LOO RMSE小于 1.5°C古温度重建常用阈值R²大于 0.6低于此值解释力不足残差均值接近 0系统偏差检查残差与预测相关不显著p 0.05 最好最后说个血泪经验有孔虫形态测量里物种鉴定的一致性比任何统计方法都重要。我见过同一批样品两个人鉴定物种名差 20%后面所有多样性指数和转换函数全废。所以做分析前先花时间统一物种名和测量标准再跑 pandas 和 numpy。希望帮到你。本文还有配套的精品资源点击获取