目录
估计农业保护措施（CP）对减少营养物质流失的影响
加载数据
数据预处理
因果推断
评估与调整
总结
估计农业保护措施（CP）对减少营养物质流失的影响
我们将使用美国农业部（USDA）发布的Measured Annual Nutrient loads from AGricultural Environments (MANAGE)数据，该数据由Qian和Harmel发布，他们在研究中使用了因果推断来量化水土保护措施在减少农田磷流失方面的作用。
数据是以Rdata格式提供的，所以我们首先需要加载并转换它。我们使用rpy2包，该包不是causallib的要求之一，因此如果你要复现分析的话，需要安装它。
加载数据
import rpy2.robjects as objects
import pandas as pd
url = objects.r['url']('https://onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1111%2F1752-1688.12377&file=jawr12377-sup-0001-manage.RData')
try:objects.r['load'](url)
except Exception: # fallback - if the above doesn't work please download the file directlyobjects.r['load']('jawr12377-sup-0001-manage.rdata')# normally we would use pandas to transform the R-datafram into a pandas dataframe,
# from rpy2.robjects import pandas2ri
# pandas2ri.activate()
# but in many versions of pandas this is broken with categorical data as we have here
# so instead we will use R to save the DataFrame to csv and load it in pandas
objects.r['write.csv'](x=objects.r["LoadApr2007"], file="MANAGE.csv", fileEncoding="UTF-8", **{"row.names": False})
manage = pd.read_csv("MANAGE.csv", index_col='RecordID')
print(manage.shape)
manage.head()
(274, 71)
数据预处理
我们将遵循Qian和Harmel使用的程序来预处理数据，其中包括创建一些对数变换和清理数据。
我们感兴趣的干预措施是等高耕作、梯田、缓冲带或建立水道中的任意一项。这些统称为保护措施（CP）。感兴趣的结果是测量到的径流水中磷的含量（AvgTPLoad）。关于为何这是一个问题，请参阅美国地质调查局（USGS）和环境保护署（EPA）的解释。
import numpy as np
#manage.fillna(manage.median(axis=0), inplace=True)
manage["AvgTPLoad"].clip(0.002, inplace=True)
manage["AvgTNLoad"].clip(0.02, inplace=True)
manage["F.AppMethod1"].replace("Other", "Surface Applied", inplace=True)
manage["F.AppMethod1"].replace("Other2", "Unknown", inplace=True)manage["CFarm"] = (manage[["CP1", "CP2", "CP3"]] == "Contour Farming").any(axis=1)
manage["Buffer"] = ((manage[["CP1", "CP2", "CP3"]] == "Filter Strip").any(axis=1)) | \((manage[["CP1", "CP2", "CP3"]] == "Riparian Buffer").any(axis=1))
manage["Terrace"] = (manage[["CP1", "CP2", "CP3"]] == "Terrace").any(axis=1)
manage["Waterway"] = (manage[["CP1", "CP2", "CP3"]] == "Waterway").any(axis=1)
manage["NumCP"] = manage[["CFarm", "Buffer", "Terrace", "Waterway"]].sum(axis=1)
manage["AnyCP"] = manage[["CFarm", "Buffer", "Terrace", "Waterway"]].any(axis=1)manage.dropna(axis=1, thresh=50, inplace=True)
%matplotlib inline
import matplotlib.pyplot as pltdef boxplot(df:pd.DataFrame, name:str, by="AnyCP", ax=None):"""A helper function to plot boxplots in log scale and annotate the median"""ax = ax or plt.gca()_ = df.boxplot(name, by=by, ax=ax)ax.set_yscale('log')ax.autoscale()ax.grid(False)medians = df[name].groupby(df[by]).median().valuesticks = ax.get_xticklabels()for median, tick in zip(medians, ticks):ax.text(x=tick.get_position()[0], y=median, s=median, horizontalalignment="center")
fig, axes = plt.subplots(2, 2)
fig.set_size_inches(14,14)
boxplot(manage, "AvgTPLoad", ax=axes[0][0])
boxplot(manage, "AvgTNLoad", ax=axes[0][1])
boxplot(manage, "Avg.N.Applied", ax=axes[1][0])
boxplot(manage, "Avg.P.Applied", ax=axes[1][1])
plt.tight_layout()
因果推断
上面的箱形图显示结果和几个特征与治疗有关联，这表明数据中存在一定的偏倚。我们现在将使用causallib来通过纠正这些偏倚来估计治疗的效果。
作为特征，我们使用施用的磷和氮的数量、总的径流水量、土壤体积、磷的应用方法以及土地利用类型（LU）
from causallib.estimation import IPW
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(solver='lbfgs', max_iter=500)
ipw = IPW(lr)
y = np.log(manage["AvgTPLoad"])a = manage["AnyCP"]X = pd.DataFrame({"Runoff": np.log(manage["Avg.Runoff"]),"Soil": np.log(manage["Avg.SoilL"] + 1),"PApplied": np.log(manage["Avg.P.Applied"] + 1),"NApplied": np.log(manage["Avg.N.Applied"] + 1)})
X = pd.concat([X, pd.get_dummies(manage["F.AppMethod1"], prefix="Method"), pd.get_dummies(manage["ShortLU"], prefix="LU"),pd.get_dummies(manage["Tillage"], prefix="Tillage")], axis=1)
X.fillna(X.median(axis=0), inplace=True)y = y.dropna() # remove the data where there is no outcome
a, X = a[y.index], X.loc[y.index]print(X.shape)
X.head()
(222, 23)
ipw.fit(X, a) # train the propensity model to apply causal inference
IPW(clip_min=None, clip_max=None, use_stabilized=False,learner=LogisticRegression(max_iter=500))
pd.DataFrame({"causal": np.exp(ipw.estimate_population_outcome(X, a, y)),"marginal": np.exp(y.groupby(a).mean())})
我们可以看到，因果推断识别出使用任何CP都会减少磷负荷，而边缘观察（仅仅检查治疗区域与未治疗区域的情况）却显示出相反的趋势。为什么会这样呢？
评估与调整
理解这种情况的一种方式是评估倾向性模型的表现并检查其一致性。causallib提供了进行这种评估的工具。
from causallib.evaluation import evaluate
evaluation_results = evaluate(ipw, X, a, y)
evaluation_results.plot_all()
plt.tight_layout()
乍一看，模型在ROC曲线（左上角）中的预测性能似乎很好。然而，加权ROC曲线远离对角线的事实提示模型实际上并不擅长平衡数据。这一点也可以从协变量平衡图（右下角）中看出，其中我们可以看到Pasture变量在数据中极端偏斜，即使加权后仍然相当偏斜。因此，按其进行分层分析或完全从分析中去除是有道理的。
为了检查是更好地进行分层还是限制我们的分析到特定的分层，让我们检查Pasture与治疗之间的重叠。我们不妨查看所有的类别特征并检查：
pd.concat({col: pd.crosstab(X[col], a) for col in X.select_dtypes('uint8').columns.values}, axis=0)
我们看到实际上只有一个是牧草地被进行了处理的情况。我们应该将分析限制在非牧草地。
但情况比这还要复杂一点，因为Alfalfa（羽扁豆）、Cotton（棉花）、Peanuts（花生）和其他类别从未观察到过处理。对于Other类别，差异相当大，而其他类别可能仍然可以归因于小样本波动。虽然这可以通过某种统计检验来进行测试，但这不是本例的重点，我们将简单地移除所有这些样本（然后删除列以避免出现全零列）。
X = X.query('LU_Pasture == 0 & LU_Other == 0 & LU_Cotton ==0 & LU_Alfalfa == 0 & LU_Peanuts ==0 & LU_Soybeans == 0')
X = X.drop(columns=['LU_Pasture', 'LU_Other', 'LU_Cotton', 'LU_Alfalfa', 'LU_Peanuts', 'LU_Soybeans'])
a, y = a[X.index], y[X.index]
print(f'We are left with {len(y)} samples')
ipw.fit(X, a)
pd.DataFrame({"causal": np.exp(ipw.estimate_population_outcome(X, a, y)),"marginal": np.exp(y.groupby(a).mean())})
We are left with 107 samples
evaluation_results = evaluate(ipw, X, a, y)
evaluation_results.plot_all()
plt.tight_layout()
现在的性能评估要好得多：
- 所有的特征都得到了充分的平衡；
- 考虑到我们现在的小样本量，模型的校准程度达到了预期的程度；
- 倾向性分布没有表现出任何一个治疗组主导的区域；
- ROC曲线表现良好。
最重要的是，在限制数据之后，因果分析与边缘分析一致，并且仅对估计的因果效应大小进行了小幅调整。
总结
总之，因果分析使我们能够关注数据允许我们检查治疗组间差异的地方。确实，与最初的分析不同，我们现在了解到，在非牧草地中，梯田、缓冲带、水道或等高耕作等保护措施实际上提高了土壤中的磷保留，并减少了种植燕麦、小麦、玉米、大豆和高粱等作物时磷向流域的泄漏。虽然我们可以推测这对其他作物也可能成立，但我们仅仅从数据中无法回答这个问题。
