原文鏈接：http://tecdat.cn/?p=23038

原文出處：拓端數(shù)據(jù)部落公眾號

最近我們被客戶要求撰寫關(guān)于k-means聚類的研究報告，包括一些圖形和統(tǒng)計輸出。

簡介

假設(shè)我們需要設(shè)計一個抽樣調(diào)查，有一個完整的框架，包含目標(biāo)人群的信息（識別信息和輔助信息）。如果我們的樣本設(shè)計是分層的，我們需要選擇如何在總體中形成分層，以便從現(xiàn)有的輔助信息中獲得最大的優(yōu)勢。

換句話說，我們必須決定以何種方式來組合輔助變量（從現(xiàn)在開始是 "X "變量）的值，來確定一個新的變量，稱為 "分層"。

為此，我們必須考慮到抽樣調(diào)查的目標(biāo)變量"Y "變量：如果為了形成分層，我們選擇與Y變量最相關(guān)的X變量，那么由此產(chǎn)生的分層框架所抽取樣本的效率就會大大增加。

每個活動變量的數(shù)值組合都決定了目標(biāo)人群的特定分層，也就是 "最佳 "分層問題的可能解決方案。在這里，我們所說的最佳分層，是指能夠確保最小樣本成本的分層，足以滿足對調(diào)查目標(biāo)變量Y's的估計精度的約束（約束表示為不同興趣領(lǐng)域的最大允許變異系數(shù)）。

當(dāng)數(shù)據(jù)收集的成本在各分層中是統(tǒng)一的，那么總成本就與總體樣本量成正比。一般來說，對于一個給定的總體來說，可能的替代分層的數(shù)量可能非常多，這取決于變量的數(shù)量和它們的值的數(shù)量，在這些情況下，不可能為了評估最佳分層而枚舉它們。一個非常方便的解決方案是采用進(jìn)化方法，包括應(yīng)用遺傳算法，在有限的迭代次數(shù)后可能收斂到一個接近最佳的解決方案。

步驟

抽樣設(shè)計的優(yōu)化首先是提供抽樣框架，確定調(diào)查的目標(biāo)估計值，并確定對其的精度限制。然后，在分析兩組變量（分層和目標(biāo)）之間的相關(guān)性的基礎(chǔ)上，必須在框架中選擇哪些分層變量。當(dāng)所選的分層變量既是分類變量又是連續(xù)變量時，為了使它們具有同質(zhì)性，應(yīng)該對連續(xù)變量進(jìn)行分類（例如使用聚類的K-means算法）。反之，如果分層變量都是連續(xù)類型的，則可以利用 "連續(xù) "方法直接執(zhí)行優(yōu)化步驟。也可以執(zhí)行兩種優(yōu)化，比較結(jié)果并選擇更方便的方法。

在使用遺傳算法進(jìn)行優(yōu)化之前，最好在使用k-means算法的基礎(chǔ)上運(yùn)行一個不同的快速優(yōu)化任務(wù)，其目的有兩個。

• 為最終分層的合適數(shù)量提供提示。

• 獲得一個初始的 "好 "解決方案，作為遺傳算法的 "建議"，以加速其向最終解決方案的收斂。

在優(yōu)化步驟中，可以指出必須選擇的抽樣單位集合（"全取 "層）。優(yōu)化之后，可以通過模擬從框架中選擇大量的樣本來評估解決方案的質(zhì)量，并計算所有目標(biāo)變量的抽樣差異和偏差。還可以根據(jù)可用預(yù)算 "調(diào)整 "優(yōu)化方案的樣本量：如果允許更大的樣本量，則按比例增加各層的抽樣率，直到達(dá)到新的總樣本量；如果我們不得不減少樣本量，則采取相反的做法。?

最后，我們開始選擇樣本。

在下文中，我們將從一個真實的抽樣框架開始說明每個步驟。

優(yōu)化步驟所需的輸入準(zhǔn)備

框架

為簡單起見，讓我們考慮數(shù)據(jù)集的一個子集。

head(mun)

為了限制處理時間，我們只選擇了前三個地區(qū)和我們例子中感興趣的變量。該數(shù)據(jù)集的每一行都包含一個城市的信息，由市政編號和市政名稱標(biāo)識，并屬于三個選定的地區(qū)之一。

假設(shè)我們要計劃一個抽樣調(diào)查，目標(biāo)估計值Ys是3個地區(qū)（感興趣的區(qū)域）中每個地區(qū)的樹林面積和建筑物面積的總數(shù)。假設(shè)每個市鎮(zhèn)的總面積和總?cè)丝诘闹悼偸潜桓隆？聪嚓P(guān)矩陣。?

cor(mun[,c(4:8)])

我們看到，樹林面積和建筑物面積之間的相關(guān)性，以及"有建筑物的區(qū)域"和"總?cè)丝?#34;之間的相關(guān)性都很高（分別為0.77和0.87），因此我們決定選擇"有建筑物的區(qū)域","總?cè)丝?#34;作為我們的框架中的分層變量X。

首先，我們決定將分層變量視為分類變量，所以我們必須對它們進(jìn)行聚類。一個合適的方法是應(yīng)用k-means聚類方法。?

我們現(xiàn)在可以按照要求的格式定義框架數(shù)據(jù)幀。以合適的模型組織數(shù)據(jù)，以便進(jìn)行下一步處理。?

Frame(df = mun,value = "REG")head(frame1)

Strata分層數(shù)據(jù)框

這個數(shù)據(jù)框架不是必需的，因為它是由從數(shù)據(jù)框架中自動生成的。不過，我們需要使用它來分析框架的初始分層，和在沒有優(yōu)化的情況下可能出現(xiàn)相關(guān)樣本量。

Strata(frameF)

該數(shù)據(jù)框架中的每一行都輸出了與給定分層有關(guān)的信息（通過對每個單元與X變量的值進(jìn)行交叉分類獲得），包括:

• 分層的標(biāo)識符（名為 "strato"）。

• 與框架中的變量相對應(yīng)的m個輔助變量（從X1到Xm命名）的值。

• 人口中的單位總數(shù)（名為 "N"）。

• 標(biāo)志（名為'cens'），表示該層是要進(jìn)行普查（=1）還是抽樣調(diào)查（=0）。

• 成本變量，表示該分層中每個單位的訪談成本。

• 每個目標(biāo)變量y的平均數(shù)和標(biāo)準(zhǔn)差，分別命名為 "Mi "和 "Si"）。

• 分層所屬的關(guān)注域的值（'DOM1'）。

精度約束

誤差數(shù)據(jù)框包含對目標(biāo)估計值設(shè)置的精度約束。這意味著要為每個目標(biāo)變量和每個域值定義一個最大的變異系數(shù)。這個框架的每一行都與感興趣的特定子域中的精度約束有關(guān)，由domainvalue值確定。在我們的案例中，我們選擇定義以下約束:

• 分層的標(biāo)識符。

• 與框架中的變量相對應(yīng)的m個輔助變量（從X1到Xm命名）的值。

• 人口中的單位總數(shù)（名為 "N"）。

• 標(biāo)志（名為'cens'），表示該層是要進(jìn)行普查（=1）還是抽樣調(diào)查（=0）。

• 成本變量，表示該分層中每個單位的訪談成本。

• 每個目標(biāo)變量y的平均數(shù)和標(biāo)準(zhǔn)差，分別命名為 "Mi "和 "Si"）。

• 分層所屬的關(guān)注域的值（'DOM1'）。

ndom <- length(unique(REG))cv

這個例子報告了變量Y1和Y2的精度約束（允許的最大CV等于10%），這些約束對于域級DOM1的所有3個不同的子域（都是一樣的。當(dāng)然，我們可以按地區(qū)區(qū)分精度約束。需要強(qiáng)調(diào)的是，'domainvalue'的值與數(shù)據(jù)框中的值相同，并且與分層數(shù)據(jù)框中的變量'DOM1'的值對應(yīng)。

現(xiàn)在我們想用函數(shù)bethel(Bethel(1989))來確定總的樣本量，以及在給定分層下的相關(guān)分配。?

hel(strata1,cv)

這是優(yōu)化前，在當(dāng)前分層下滿足精度約束所需的總樣本量（570）。?

優(yōu)化

執(zhí)行優(yōu)化步驟的函數(shù)。實際上，調(diào)用三個不同的優(yōu)化函數(shù):

1. 優(yōu)化分層（方法=原子，分層變量是分類的）;

2. 優(yōu)化分層2（方法=連續(xù)，分層變量是連續(xù)的）;

3. 優(yōu)化空間分層（方法=空間，分層變量是連續(xù)的，或為分類的）。

這里我們報告與三種方法有關(guān)的最重要的參數(shù):

參數(shù)說明frame數(shù)據(jù)框架的名稱，包含與抽樣框架有關(guān)的信息cens包含在任何情況下要選擇的單位的數(shù)據(jù)框名稱n要獲得的最終優(yōu)化分層的數(shù)量err包含精度等級的數(shù)據(jù)幀min每個層必須分配的最低單位數(shù)

?‘a(chǎn)tomic’方法

第一次運(yùn)行，我們執(zhí)行優(yōu)化步驟（由于分層變量是分類類型的，所以需要使用atomic方法）。

Strata(method = "atomic", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?pops = 10)

執(zhí)行產(chǎn)生了3個不同的優(yōu)化問題的解決方案。圖中說明了從初始解開始向最終解收斂的情況。在X軸上報告了已執(zhí)行的迭代，從1到最大，而在Y軸上報告了為滿足精度約束所需的樣本大小。上部（紅色）線代表每次迭代的平均樣本大小，而下部（黑色）線代表直到第i次迭代所發(fā)現(xiàn)的最佳解決方案。

我們可以通過執(zhí)行函數(shù)來計算（分析）預(yù)期的CVs：

所得到的滿足精度約束的樣本總規(guī)模比簡單應(yīng)用Bethel算法對初始分層所得到的要低得多，但還不能令人滿意。

為了探索其他解決方案，我們可能希望將抽樣框架中的每個單元都視為一個原子分層，并讓優(yōu)化步驟根據(jù)Y變量的值對其進(jìn)行匯總。在任何情況下，由于我們必須指出至少一個X變量，我們可以為此使用一個簡單的遞增數(shù)字。

frame2 <- FrameDF( mun)head(frame2)

我們可以使用這種方法，因為框架中的單位數(shù)量很少。

即便如此，對1823個層的處理也可能很慢。

為了加快向最優(yōu)解收斂的速度，可以給一個初始解作為 "建議"。通過考慮所有目標(biāo)變量Y的均值對原子層進(jìn)行聚類來產(chǎn)生這個初始解。滿足精度約束所需的樣本量為最小值的聚類數(shù)目被保留為最優(yōu)數(shù)目。此外，每個領(lǐng)域內(nèi)的最佳聚類數(shù)也被確定?？梢灾赋鲆@得的最大聚類層數(shù)。?

Strata(frame2, progress=F) Kmeans(strata = strata2, ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?maxclusters = 10) ?

整體解決方案是通過串聯(lián)各領(lǐng)域獲得的最優(yōu)聚類而獲得的。其結(jié)果是一個有兩列的數(shù)據(jù)框架：第一列表示聚類，第二列表示域。在此基礎(chǔ)上，我們可以為每個域計算出最方便的最終層數(shù)。?

apply(suggestions, ? ? ? ? ? ? ? ? domainvalue,length(unique(x)))

而我們可以提供初始解和層數(shù)nstrata作為優(yōu)化步驟的輸入。?

Strata( ? ? ? ? ? ? ? ? ? ? ? ?nStrata = nstrata2)

請注意，在這次運(yùn)行中得到的解決方案在樣本量方面明顯優(yōu)于前一次。?

"連續(xù)"方法

最后要做的是測試連續(xù)方法。

首先，我們必須以這種方式重新定義框架dataframe。

FrameDF(+ ? ? ? ? ? ? ? ? ? ? ? ? ?id = "市政編號",+ ? ? ? ? ? ? ? ? ? ? ? ? ?X = c("總?cè)丝?#34;,"建筑物面積"),+ ? ? ? ? ? ? ? ? ? ? ? ? ?Y = c("有建筑物的區(qū)域","樹林面積"),+ ? ? ? ? ? ? ? ? ? ? ? ? ?value = "地區(qū)")head(frame3)

同樣在這種情況下，我們想生成一個初始解決方案。

Kmeans2(frame=frame3, ? ? ? ? ? ? ? ? ? ? ? ? ? ? maxclusters = 10) ?

apply(sugge, ? ? ? ? ? ? ? value, ? ? ? ? ? ? ? ? FUN=function(x) length(unique(x)))

請注意，這一次我們調(diào)用了一個不同的函數(shù)，它需要的不是階層數(shù)據(jù)框架，而是直接的框架數(shù)據(jù)框架。此外，我們需要一個中間步驟為優(yōu)化準(zhǔn)備建議。

我們現(xiàn)在能夠用連續(xù)方法進(jìn)行優(yōu)化步驟。?

這個解決方案需要的總樣本量是迄今為止我們制作的解決方案中最好的，所以我們決定選擇這個方案。

分析

分層結(jié)構(gòu)

執(zhí)行的結(jié)果包含在由三個元素組成的 "解決方案 "列表中。

1. indices: 指數(shù)向量，表示每個原子層屬于哪個集合層（如果使用的是原子法）或者框架中的每個單元屬于哪個集合層（如果使用的是連續(xù)法）。

2. framenew: 更新的初始框架，對每個單元而言，表明每個單元所屬的最佳層。

3. aggr_strata: 包含優(yōu)化層信息的數(shù)據(jù)框.

當(dāng)分層變量為連續(xù)類型，并且使用了連續(xù)（或空間）方法時，有可能獲得關(guān)于優(yōu)化的分層結(jié)構(gòu)的詳細(xì)信息。?

summaryStrata(framenew, ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? progress=FALSE)

對于每個優(yōu)化的分層，報告了單位總數(shù)以及抽樣率。還列出了分層變量的范圍，描述各層的特點。

如果分層變量的數(shù)量有限，就像我們的情況一樣，可以通過選擇一對變量和每個時間段的一個域來可視化分層。

plot(framenew, ? ? ? ? ? ? vars = c("X1","X2"), ? ? ? ? ? ? ? labels = c("總體","總面積"))

?

通過模擬進(jìn)行評估

為了對找到的解決方案的質(zhì)量有信心，我們運(yùn)行一個模擬，基于從已確定為最佳分層的框架中選擇所需數(shù)量的樣本。用戶可以指出要抽取的樣本數(shù)量。

Solution( ? ? ? ? ? ? ? ? ? ? nsampl = 200)

對于每個抽取的樣本，都會計算與Y有關(guān)的估計值。它們的平均數(shù)和標(biāo)準(zhǔn)差也計算出來，得出每個領(lǐng)域中與每個變量相關(guān)的CV和相對偏差。

coeff_varrel_bias

?

還可以分析所選域中每個相關(guān)變量的估計值的抽樣分布。

hist(eval3 ) ?abline(v = mean(eval3$es )abline(v = mean(frame3$Y )

最終樣本大小的調(diào)整

在優(yōu)化步驟之后，最終的樣本量是最終分層中單位分配的結(jié)果。這種分配是為了使精度約束得到滿足。實際上，可能會出現(xiàn)三種情況:

• 產(chǎn)生的樣本量是可以接受的;

• 所得的樣本量太多，也就是說，超過預(yù)算;

• 產(chǎn)生的樣本量太少，可用的預(yù)算允許增加單位的數(shù)量。

在第一種情況下，不需要變化。在第二種情況下，有必要減少單位數(shù)，在每個分層中平均采用相同的減少率。在第三種情況下，我們著手增加樣本量，在每個分層中應(yīng)用相同的增加率。這個增加/減少的過程是反復(fù)進(jìn)行的，因為通過應(yīng)用相同的比率，我們可以發(fā)現(xiàn)在某些層沒有足夠的單位可以增加或減少?？梢垣@得理想的最終樣本量。

讓我們假設(shè)最終獲得的樣本量（106）是超過預(yù)算。我們可以通過執(zhí)行以下代碼來減少它。

adjustSize(size=75,strata)

相反，如果我們想增加抽樣樣本規(guī)模，因為預(yù)算允許。?

希望的樣本量和實際調(diào)整后的樣本量之間的差異取決于優(yōu)化方案中的分層數(shù)量?？紤]到在每個分層中進(jìn)行調(diào)整時，要考慮到當(dāng)前樣本量與期望樣本量之間的相對差異：這將產(chǎn)生一個用實數(shù)表示的分配，必須四舍五入，同時要考慮到分層中最小單位數(shù)的要求（默認(rèn)為2）。層數(shù)越多，對最終調(diào)整后的樣本量的影響就越大。

一旦增加（或減少）樣本量，我們可以檢查新的預(yù)期CV是多少。通過第二次調(diào)整，即增加了總的樣本量，我們得到。?

CV(adjusted)

即大大降低了期望CV。

樣本選擇

一旦獲得最佳分層，就可以從優(yōu)化版的框架中選擇樣本，同時考慮到最佳分層。

Sample(new3, ? ? ? ? ? ? ? ? ? ? strata3,

在每個分層中進(jìn)行簡單的隨機(jī)抽樣。

一個變體是系統(tǒng)抽樣?。唯一的區(qū)別是在每個分層中選擇單位的方法，即通過執(zhí)行以下步驟：

• 通過考慮分層中采樣率的倒數(shù)確定選擇區(qū)間；通過選擇該區(qū)間中的一個值確定起點。

• 通過選擇與上述數(shù)值相對應(yīng)的單位作為第一個單位，然后選擇所有加入選擇區(qū)間而被分割的單位，進(jìn)行選擇。

如果與選擇框架的特定排序相關(guān)聯(lián)，這種選擇方法是有用的，其中排序變量可以被視為額外的分層變量。例如，我們可以考慮市政當(dāng)局的工業(yè)區(qū)（Airind）可能很重要。

Sample(variable = c("工業(yè)區(qū)"))

?"全取 "分層抽樣

作為優(yōu)化步驟的輸入，與適當(dāng)?shù)某闃臃謱右黄?，也可以提?"全取 "分層。這些層不會像適當(dāng)?shù)膶幽菢颖粌?yōu)化，但是它們將有助于確定最佳分層，因為它們可以使較少的抽樣層單位數(shù)量來滿足精度約束。

為了正確執(zhí)行優(yōu)化和進(jìn)一步的步驟，有必要對整個輸入進(jìn)行預(yù)處理。要執(zhí)行的第一步是對要普查的單位和要抽樣的單位進(jìn)行兩分法，以建立兩個不同的框架。例如，我們要確保所有總?cè)丝诔^10,000的城市都將被納入抽樣范圍。因此，我們以這種方式劃分抽樣框架。

#----從框架中選擇要進(jìn)行普查的單位 nrow(which(X1 > 10000)#----選擇要從框架中取樣的單位#(對前者的補(bǔ)充)nrow(rame3[-ind,])

這樣，我們將所有人口超過10,000的單位定義為要進(jìn)行普查。在這個過程結(jié)束時，樣本將包含所有這些單位。

我們現(xiàn)在運(yùn)行優(yōu)化步驟，將待普查單位的指示包括在內(nèi)。?

set.seed(1234)Strata(method = "continuous", ? ? ? ? ? ? ? ? ? ? ?)

最受歡迎的見解

1.使用R語言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行

2.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

3.R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣

4.R語言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣

5.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

6.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析

7.R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)

8.R語言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸

9.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測GDP增長