伽瑪障礙分佈 | 數據科學進階

哪個結果最重要？

這是一個常見的情境：進行了一個A/B測試，隨機選擇了一些單位（例如顧客）參加一個活動，這些顧客接受了處理A。另一組顧客則接受了處理B。“A”可能是一個促銷或優惠，而“B”則可能沒有任何促銷或優惠。“A”可能是打九折，而“B”可能是打八折。兩組顧客，兩種不同的處理方式，A和B是兩種不同的處理，但這個概念可以擴展到多於兩種的處理方式。

活動結束後，結果出爐。透過我們的後端系統，我們可以追蹤這些顧客中哪些進行了我們關心的行動（例如，購買了商品），哪些則沒有。此外，對於那些有行動的顧客，我們還記錄了他們的行動強度。常見的情況是，我們可以追蹤購買金額，這通常被稱為平均訂單金額或每位顧客的收入指標。或者有一百種不同的名稱，意思都是一樣的——對於那些購買的顧客，他們平均花了多少錢？

對於某些情況，行銷人員會關心前者指標——購買率。例如，我們在獲客活動中使用處理A或B，是否吸引了更多（可能是第一次）購買的顧客？有時，我們更關心提高每位顧客的收入，因此會強調後者。

不過，通常我們更關心以成本效益的方式提高收入，而我們真正關心的是這次活動總共產生了多少收入。處理A或B產生了更多的收入？我們不一定有平衡的樣本數（可能因為成本或風險考量），因此我們將測量的收入除以每組中接受處理的顧客數量（稱為N_A和N_B）。我們想比較這兩組的這個指標，因此標準的對比方式就是：

這只是處理A的平均收入減去處理B的平均收入，這個平均值是針對整個目標單位集計算的，無論他們是否有回應。它的解釋也很簡單——從處理A到處理B，每個推廣單位的平均收入變化是多少？

當然，這最後的指標同時考慮了前兩者：回應率乘以每位回應者的平均收入。

不確定性？

顧客的消費金額變化很大，某一組的幾筆大額購買可能會顯著影響平均值。同樣，樣本變異也可能很大。因此，我們想了解在這個平均值比較中，我們的信心有多大，並量化觀察到的差異的“顯著性”。

所以，你可以將數據放入t檢驗中，然後看p值。但是等等！不幸的是，對於行銷人員來說，大多數情況下，購買率相對較低（有時非常低），因此有很多零收入的值——通常是絕大多數。t檢驗的假設可能會被嚴重違反。非常大的樣本數可能會有所幫助，但還有一種更原則性的方式來分析這些數據，這在多方面都很有用，接下來會解釋。

示例數據集

讓我們開始使用示例數據集，讓事情變得實用。我最喜歡的直接行銷數據集之一來自KDD Cup 98。

url="https://kdd.ics.uci.edu/databases/kddcup98/epsilon_mirror/cup98lrn.zip"

filename="cup98LRN.txt"



r = requests.get(url)

z = zipfile.ZipFile(io.BytesIO(r.content))

z.extractall()



pdf_data = pd.read_csv(filename, sep=',')

pdf_data = pdf_data.query('TARGET_D >=0')

pdf_data['TREATMENT'] =  np.where(pdf_data.RFA_2F >1,'A','B')

pdf_data['TREATED'] =  np.where(pdf_data.RFA_2F >1,1,0)

pdf_data['GT_0'] = np.where(pdf_data.TARGET_D >0,1,0)

pdf_data = pdf_data[['TREATMENT', 'TREATED', 'GT_0', 'TARGET_D']]

在上面的代碼片段中，我們下載了一個zip檔案（特別是學習數據集），提取它並將其讀入Pandas數據框。這個數據集的性質是來自一個非營利組織的活動歷史，該組織通過直接郵寄尋求捐款。這個數據集中沒有處理變量，因此我們假裝並根據過去捐款的頻率對數據集進行分段。我們稱這個指標為TREATMENT（作為類別變量），並創建TREATED作為‘A’的二元指標。可以將這視為隨機對照試驗的結果，其中部分樣本人口接受了優惠，剩下的則沒有。我們追蹤每個個體並累計他們的捐款金額。

所以，如果我們檢查這個數據集，我們會看到大約有95,000名被推廣的個體，通常在兩個處理中均勻分佈：

處理A的回應率較高，但整體數據集的回應率僅約5%。因此，我們有95%的零。

對於那些捐款的人，處理A似乎與較低的平均捐款金額相關。

將所有被針對的人合併在一起，處理A似乎與較高的平均捐款金額相關——較高的回應率超過了回應者的較低捐款金額，但差異不大。

最後，這裡顯示了捐款金額的直方圖，這是針對兩種處理的總和，顯示了零的集中和右偏。

對於這兩個處理組的數值摘要量化了上述現象——雖然處理A似乎驅動了顯著更高的回應，但那些接受A處理的回應者在回應時的平均捐款金額較低。這兩個指標的淨值，即我們最終追求的——每個被針對單位的整體平均捐款——似乎仍然對處理A較高。我們對這一發現的信心是本分析的主題。

伽瑪障礙

一種建模這些數據並回答我們的研究問題的方法是使用伽瑪障礙分佈。這類似於更知名的零膨脹泊松（ZIP）或負二項（ZINB）分佈，這是一種混合分佈，其中一部分與零的集中有關，而另一部分則在隨機變量為正的情況下，使用伽瑪密度函數。

這裡π代表隨機變量y大於0的概率。換句話說，這是伽瑪過程的概率。同樣，(1- π)是隨機變量為零的概率。在我們的問題中，這涉及到捐款是否被做出，如果是，則其值。

讓我們開始使用這個分佈進行回歸的組件部分——邏輯回歸和伽瑪回歸。

邏輯回歸

這裡的logit函數是連結函數，將對數賠率與我們的預測變量的線性組合相關聯，對於像我們的二元處理指標這樣的單一變量，形式如下：

其中π代表結果為“正面”（標記為1）事件的概率，例如購買，而(1-π)代表結果為“負面”（標記為0）事件的概率。此外，上述的π是我們關心的數量，由反logit函數定義：

擬合這個模型非常簡單，我們需要找到兩個beta的值，使數據（結果y）的可能性最大化——假設N是獨立同分佈的觀察值：

我們可以使用多個庫來快速擬合這個模型，但將演示PYMC作為建立簡單貝葉斯邏輯回歸的手段。

在沒有任何貝葉斯工作流程的正常步驟下，我們使用MCMC擬合這個簡單模型。

import pymc as pm

import arviz as az

from scipy.special import expit



with pm.Model() as logistic_model:



    # 非信息性先驗

    intercept = pm.Normal('intercept', 0, sigma=10)

    beta_treat = pm.Normal('beta_treat', 0, sigma=10)



    # 通過反logit將處理變量的線性組合壓縮到0和1之間

    p =  pm.invlogit(intercept + beta_treat * pdf_data.TREATED)



    # 個體級二元變量（回應或不回應）

    pm.Bernoulli(name="logit", p=p, observed=pdf_data.GT_0)



    idata = pm.sample(nuts_sampler="numpyro")

az.summary(idata, var_names=[‘intercept’, ‘beta_treat’])

如果我們構建兩個處理的平均回應率對比，我們發現，正如預期的那樣，處理A的平均回應率比處理B高0.026，94%的可信區間為(0.024 , 0.029)。

# 在後驗中創建一個新列，對比處理A - B

idata.posterior['TREATMENT A - TREATMENT B'] = expit(idata.posterior.intercept + idata.posterior.beta_treat) -  expit(idata.posterior.intercept)



az.plot_posterior(

    idata,

    var_names=['TREATMENT A - TREATMENT B']

)

伽瑪回歸

下一個組件是伽瑪分佈及其概率密度函數的一種參數化，如上所示：

這個分佈是為嚴格正的隨機變量定義的，並且如果用於商業，則用於成本、顧客需求支出和保險索賠金額等值。

由於伽瑪的均值和方差是根據α和β的公式定義的：

對於伽瑪回歸，我們可以通過α和β或μ和σ進行參數化。如果我們將μ定義為預測變量的線性組合，則可以使用μ將伽瑪定義為α和β：

伽瑪回歸模型假設（在這種情況下，反連結是另一個常見選擇）使用對數連結，旨在“線性化”預測變量和結果之間的關係：

幾乎完全按照與回應率相同的方法，我們僅限於回應者的數據集，並使用PYMC擬合伽瑪回歸。

with pm.Model() as gamma_model:



    # 非信息性先驗

    intercept = pm.Normal('intercept', 0, sigma=10)

    beta_treat = pm.Normal('beta_treat', 0, sigma=10)



    shape = pm.HalfNormal('shape', 5)



    # 通過exp確保線性預測器為正

    mu =  pm.Deterministic('mu',pm.math.exp(intercept + beta_treat * pdf_responders.TREATED))



    # 個體級二元變量（回應或不回應）

    pm.Gamma(name="gamma", alpha = shape, beta = shape/mu,  observed=pdf_responders.TARGET_D)



    idata = pm.sample(nuts_sampler="numpyro")

az.summary(idata, var_names=[‘intercept’, ‘beta_treat’])

# 在後驗中創建一個新列，對比處理A - B

idata.posterior['TREATMENT A - TREATMENT B'] = np.exp(idata.posterior.intercept + idata.posterior.beta_treat) -  np.exp(idata.posterior.intercept)



az.plot_posterior(

    idata,

    var_names=['TREATMENT A - TREATMENT B']

)

再次，正如預期的那樣，我們看到處理A的平均值提升預測值等於樣本值的-7.8。94%的可信區間為(-8.3, -7.3)。

這些組件，回應率和每位回應者的平均金額，顯示出來的簡單程度幾乎是我們能做到的最簡單的。但，這是一個直接的擴展，可以添加額外的預測變量，以便1）估計條件平均處理效應（CATE），當我們預期處理效應會因細分而異時，或2）通過對預處理變量進行條件化來減少平均處理效應估計的變異。

障礙模型（伽瑪）回歸

此時，應該很容易看出我們的進展方向。對於障礙模型，我們有一個條件可能性，根據特定觀察值是0還是大於零，如上所示的伽瑪障礙分佈。我們可以同時擬合這兩個組件模型（邏輯回歸和伽瑪回歸）。我們免費獲得它們的乘積，在我們的例子中，這是每個被針對單位的捐款金額的估計。

使用基於結果變量值的開關語句的可能性函數擬合這個模型並不困難，但PYMC已經為我們編碼了這個分佈。

import pymc as pm

import arviz as az



with pm.Model() as hurdle_model:



    ## 非信息性先驗 ##

    # 邏輯回歸

    intercept_lr = pm.Normal('intercept_lr', 0, sigma=5)

    beta_treat_lr = pm.Normal('beta_treat_lr', 0, sigma=1)



    # 伽瑪回歸

    intercept_gr = pm.Normal('intercept_gr', 0, sigma=5)

    beta_treat_gr = pm.Normal('beta_treat_gr', 0, sigma=1)



    # α

    shape = pm.HalfNormal('shape', 1)



    ## 預測變量的均值函數 ##

    p =  pm.Deterministic('p', pm.invlogit(intercept_lr + beta_treat_lr * pdf_data.TREATED))

    mu =  pm.Deterministic('mu',pm.math.exp(intercept_gr + beta_treat_gr * pdf_data.TREATED))

    

    ## 可能性 ##

    # psi是pi

    pm.HurdleGamma(name="hurdlegamma", psi=p, alpha = shape, beta = shape/mu, observed=pdf_data.TARGET_D)



    idata = pm.sample(cores = 10)

如果我們檢查追蹤摘要，我們會看到這兩個組件模型的結果完全相同。

如前所述，伽瑪障礙分佈的均值是π * μ，因此我們可以創建一個對比：

# 在後驗中創建一個新列，對比處理A - B

idata.posterior['TREATMENT A - TREATMENT B'] = ((expit(idata.posterior.intercept_lr + idata.posterior.beta_treat_lr))* np.exp(idata.posterior.intercept_gr + idata.posterior.beta_treat_gr)) - \

                                                    ((expit(idata.posterior.intercept_lr))* np.exp(idata.posterior.intercept_gr))



az.plot_posterior(

    idata,

    var_names=['TREATMENT A - TREATMENT B']

這個模型的均值預期值為0.043，94%的可信區間為(-0.0069, 0.092)。我們可以檢查後驗，看看捐款每位顧客預測為更高的比例，以及任何其他對我們案例有意義的決策函數——包括將更完整的損益表納入估計（即包括利潤和成本）。

註：某些實現對伽瑪障礙模型的參數化不同，其中零的概率是π，因此伽瑪障礙的均值涉及(1-π)。還要注意，在撰寫本文時，PYMC中的nuts取樣器似乎有問題，我們不得不回退到Python的默認實現來運行上述代碼。

總結

通過這種方法，我們獲得了兩個模型的相同推斷，並且額外獲得了第三個指標。使用PYMC擬合這些模型使我們能夠享受貝葉斯分析的所有好處——包括注入先前的領域知識和完整的後驗來回答問題和量化不確定性！

致謝：

所有圖片均為作者所有，除非另有說明。

使用的數據集來自KDD 98 Cup，由Epsilon贊助。https://kdd.ics.uci.edu/databases/kddcup98/kddcup98.html (CC BY 4.0)