R統計 | 資料加權

抽樣調查的準確度取決於樣本有沒有代表性。但在資料蒐集的過程中難免因為抽樣設計、訪問失敗等因素導致資料蒐集不完整，進而影響樣本代表性。為了避免樣本與母體結構之間的落差，除了追蹤訪問、插補法，加權是最常用的方法，尤其是透過多變數反覆加權(Raking)與事後分層加權(Post-stratification)，讓樣本分布盡可能符合母體特徵，進而降低樣本代表性的問題

一個常見的例子是女性的訪問成功率通常比男性高，年輕人的拒訪率通常比高齡長者低。這些因素都會造成蒐集的資料女性比男性多，年輕人比高齡長者高，樣本結構偏離母體，男性與高齡長者的意見將被低估，連帶影響調查的可信度。加權的概念其實不難理解，只是透過數學方法將樣本乘以一個權重，代表性高就乘以較小的權重，代表性低就乘以較高的權重，藉此調整樣本的比例分配符合母體而已。所以問題在於我們有沒有辦法知道母體長什麼樣子？

人口特徵是最常見可獲得的母體分布。經由政府統計，我們可以輕易掌握台灣的性別比、各年齡層人口比、居住地區人口比等資訊。當我們知道母體聯合分布的時候，便可使用事後分層加權；只知道母體邊際分布的時候，便可使用多變數反覆加權。下表是內政部2021年11月人口統計的台灣人口母體分布：

		19歲以下	20-39歲	40-59歲	60歲以上	邊際%
男性	N (聯合%)	2,078,547 (8.9%)	3,278,738 (14.0%)	3,629,464 (15.5%)	2,601,825 (11.1%)	11,588,574 (49.5%)
女性	N (聯合%)	1,917,129 (8.2%)	3,111,298 (13.3%)	3,767,738 (16.1%)	3,010,048 (12.9%)	11,806,213 (50.5%)
合計	邊際%	17.1%	27.3%	31.6%	24.0%	100%

我們不難發現為什麼多變數反覆加權是比較常用的加權方法，與其知道母體聯合分布的資訊，邊際分布更容易獲得。只要知道性別比與年齡層人口比，就能進行反覆加權；相較之下事後分層得必須有男女性在各年齡層的比例，難度要高上許多。如果再加上其他人口特徵，有時候聯合分布是不得而知的。

多變數反覆加權

Deming與Stephan(1940)第一次提出反覆加權的概念後，它有不同的名稱包含sample-balancing、raking ratio estimation、iterative proportional fitting，不過最常用的名稱還是Raking。

所謂多變數反覆加權是指利用多個重要變數，例如性別、年齡、地區、教育程度等，一次加權一個變數，直到每一個變數都能通過卡方適合度檢定與母體沒有差異為止。為了達成樣本與母體一致的目標，通常得進行多次迭代，讓樣本的邊際次數接近母體的邊際次數。權重的計算公式為：

\[W_{ij}=\frac{N_{ij}}{N} \times \frac{n}{n_{ij}}\]

舉例來說，母體的性別比為50%:50%時，如果樣本中有60位女性、40位男性，其權重應該為：

\(W_{female}=\frac{50}{100} \times \frac{100}{60}=0.83\)

\(W_{male}=\frac{50}{100} \times \frac{100}{40}=1.25\)

權重代表在這個樣本中女性太多，男性太少，所以女性應該乘以較小的權重、男性要乘以較大的權重來予以調整到接近母體50%:50%的程度。

R有兩個重要的延伸套件survey、anesrake可以處理加權，尤其是University of Auckland生物統計學教授Thomas Lumley開發的survey更是處理加權資料分析的重要套件，其介紹可參考Lumley成立的Survey analysis in R。我們現在就以虛擬資料gender_age.csv來示範如何用survey與anesrake進行多變數反覆加權。


     > example<-read.csv("c:/Users/USER/Downloads/gender_age.csv", header=T, sep=",")

     > head(example)

       gender age
     1      1   1
     2      1   1
     3      1   1
     4      1   1
     5      1   2
     6      1   2

這是一個包含100筆個案的資料檔，我們以內政部2021年11月的人口統計作為母體來進行加權。首先得先瞭解樣本的抽樣分布。可以看到母體的性別比為男：女=49.5%：50.5%，但樣本性別比為39%：61%。母體的年齡層比例為17.1%：27.3%：31.6%：24.0%，但樣本為11%：41%：36%：12%。


     > example$gender<-factor(example$gender, levels=c(1,2), labels=c("Male","Female"))

     > example$age<-factor(example$age, levels=c(1,2,3,4), labels=c("Under 19","20-39","40-59","Over 60"))

     > table(example$gender)

     

       Male Female
         39     61
     > table(example$age)

     

     Under 19    20-39    40-59  Over 60
           11       41       36       12

卡方適合度檢定可以發現，無論是性別或年齡的p值都<.05，顯示樣本結構與母體不一致。


     > chisq.test(table(example$gender), p=c(0.495, 0.505))

     

             Chi-squared test for given probabilities
     

     data:  table(example$gender)
     X-squared = 4.4104, df = 1, p-value = 0.03572
     

     > chisq.test(table(example$age), p=c(0.171, 0.273, 0.316, 0.240))

     

             Chi-squared test for given probabilities
     

     data:  table(example$age)
     X-squared = 15.664, df = 3, p-value = 0.001329

survey套件

survey其實是用來處理複雜抽樣設計(complex survey design)的套件。我們常見熟悉的統計方法是建立在簡單隨機抽樣每個樣本有相同機率分布、相同權重的假設上，但實務上這種獨立同分布(independent and identically distributed, i.i.d)的假設碰到分層抽樣、集群抽樣、加權等作法已經違反，這種複雜抽樣所獲得的樣本變異數通常會較簡單隨機抽樣來得大，所以用一般統計方法分析會導致估計錯誤，才會衍生出複雜抽樣調查的資料分析。survey就是R專門處理複雜抽樣調查的分析套件。

複雜抽樣調查與簡單隨機抽樣對統計分析的影響，可以參考侯佩君(2011)。我們現在以survey套件中的rake()來對資料進行加權，必須注意的是survey裡所有的指令都必須寫成方程式的型態。首先先用svydesign()宣告抽樣調查資料，並指定ids參數為0代表這是一個簡單隨機抽樣的資料檔。隨後創造gender_p與age_p兩筆資料分別是性別與年齡層的母體邊際總和，作為後續反覆加權時的對照資料。


     > library(survey)

     > unweighted<-svydesign(ids=~0, data=example)

     > gender_p<-data.frame(gender=c("Male","Female"), Freq=nrow(example)*c(0.495,0.505))

     > age_p<-data.frame(age=c("Under 19","20-39","40-59","Over 60"), Freq=nrow(example)*c(0.171,0.273,0.316,0.240))

萬事俱備後，可以用rake()直接進行多變數反覆加權，並用stack()輸出加權結果，可以看到有23筆資料的權重為0.484，21筆資料權重為0.647，18筆資料權重為0.90，以此類推。


     > raking<-rake(design=unweighted, sample.margins=list(~gender, ~age), population.margins=list(gender_p, age_p))

     > stack(table(weights(raking)))

       values               ind
     1     23 0.483958072437301
     2     21 0.646990276173929
     3     18 0.898275796330115
     4      7  1.18549042351196
     5     15  1.20088028002317
     6     10  1.75026554378407
     7      4  2.20039175885406
     8      2  3.24867228107965

把計算後的權重結合原始資料，更能清楚瞭解每個觀察值的權重。


     > example$weights<-weights(raking)

     > example

         gender      age   weights
     1     Male Under 19 2.2003918
     2     Male Under 19 2.2003918
     3     Male Under 19 2.2003918
     4     Male Under 19 2.2003918
     5     Male    20-39 0.8982758
     6     Male    20-39 0.8982758
     ⋮
     98  Female  Over 60 1.7502655
     98  Female  Over 60 1.7502655
     100 Female  Over 60 1.7502655

有了每個個案的權重後，我們再一次用svydesign()將example宣告為加權後的樣本，以利比較加權前與加權後的樣本分配。原本樣本結構中男性偏少的問題經加權調整後已獲得解決，男性與女性的比例接近50%：50%。


     > svytable(~gender+age, unweighted) #加權前

             age
     gender   Under 19 20-39 40-59 Over 60
       Male          4    18    15       2
       Female        7    23    21      10
     > weighted<-svydesign(id=~0, data=example, weights=~weights)

     > summary(svytable(~gender+age, design=weighted)) #加權後

             age
     gender   Under 19 20-39 40-59 Over 60
       Male          9    16    18       6
       Female        8    11    14      18

以survey套件裡的svygofchisq()進行卡方適合度檢定，p值>.05顯示資料經過加權處理後與母體結構一致。


     > svygofchisq(~gender, p=c(0.495,0.505), design=weighted)

     

             Design-based chi-squared test for given probabilities
     

     data:  ~gender
     X-squared = 1.432e-05, scale = 1.3262, df = 1.0000, p-value = 0.9974
     

     > svygofchisq(~age, p=c(0.171, 0.273, 0.316, 0.240), design=weighted)

     

             Design-based chi-squared test for given probabilities
     

     data:  ~age
     X-squared = 0, scale = 1.9046, df = 2.3967, p-value = 1

anesrake套件

anesrake是設計用來加權美國全國大選(American National Election Studies)的套件，雖然不像survey有分析複雜抽樣設計的功能，但一樣可以對資料進行加權，指令就是anesrake()。用anesrake()對資料加權前，必須先將母體分布轉換為表單形式。


     > library(anesrake)

     > gender_target<-c(0.495, 0.505)

     > age_target<-c(0.171, 0.273, 0.316, 0.240)

     > names(gender_target)<-c("Male","Female")

     > names(age_target)<-c("Under 19", "20-39", "40-59", "Over 60")

     > target<-list(gender_target, age_target)

     > names(target)<-c("gender", "age")

     > target

     $gender
       Male Female
      0.495  0.505
     

     $age
     Under 19    20-39    40-59  Over 60 
        0.171    0.273    0.316    0.240

接下來因為anesrake()指令的需要，必須賦予每個樣本一個編號。所以我們創造出一個名為id的變數，從1開始編碼到100，然後與原本的資料作結合。


     > id<-1:length(example$age)

     > example<-cbind(example, id)

     > head(example)

       gender      age   weights id
     1   Male Under 19 2.2003918  1
     2   Male Under 19 2.2003918  2
     3   Male Under 19 2.2003918  3
     4   Male Under 19 2.2003918  4
     5   Male    20-39 0.8982758  5
     6   Male    20-39 0.8982758  6

有了母體分布表單以及樣本編號後，直接以anesrake進行多變數反覆加權，可以看到整個過程經過16次的迭代才收斂。


     > raking2<-anesrake(target, example, caseid=example$id, choosemethod="total")

     [1] "Raking converged in 16 iterations"

用weightvec參數並配合stack()指令整理經過16次迭代後的權值如下：


     > stack(table(raking2$weightvec))

       values               ind
     1     23  0.48372580975453
     2     21 0.646691420427651
     3     18 0.898572576424768
     4      7  1.18499595188209
     5     15  1.20129867806796
     6     10  1.74988148834905
     7      4  2.20125708420635
     8      2  3.25059255825473

我們一樣可以把計算後的權值納入原始資料，作為後續資料加權的基礎。


     > example$weights2<-raking2$weightvec

最後比較survey與anesrake計算的權值如下表，可以發現兩者差異在小數第三位後，無論用哪個指令來進行多變數反覆加權都可以獲得差不多的權值。

分層	survey	anesrake	加權前樣本數	加權後樣本數
女性 20-39歲	0.483958	0.483725	23	11
女性 40-59歲	0.646990	0.646691	21	14
男性 20-39歲	0.898275	0.898572	18	16
女性 19歲以下	1.185490	1.184995	7	8
男性 40-59歲	1.20088	1.20129	15	18
女性 60歲以上	1.750265	1.749881	10	18
男性 19歲以下	2.200391	2.201257	4	9
男性 60歲以上	3.248672	3.250592	2	6

事後分層加權

當我們能得知變數在母體中的聯合分布時，就可以採用事後分層加權。這裡同樣以內政部2021年11月的人口統計作為母體，重新再回顧一下性別與年齡共8個聯合機率。

	19歲以下	20-39歲	40-59歲	60歲以上
男性	8.9%	14.0%	15.5%	11.1%
女性	8.2%	13.3%	16.1%	12.9%

對照樣本的聯合分布可以發現兩者有明顯差異，特別是在60歲以上的年齡層，男性有明顯低估的情況，這也是事後分層加權要調整到與母體分布一致的地方。


     > genderXage<-table(example$gender, example$age)

     > addmargins(prop.table(genderXage))

     

              Under 19 20-39 40-59 Over 60  Sum
       Male       0.04  0.18  0.15    0.02 0.39
       Female     0.07  0.23  0.21    0.10 0.61
       Sum        0.11  0.41  0.36    0.12 1.00

進行事後分層加權之前，因為原樣本資料沒有分層變數，所以我們先自行創造出strata變數代表性別與年齡構成的8個分層。


     > strata<-c("")

     > strata[example$gender=="Male" & example$age=="Under 19"]<-1

     > strata[example$gender=="Male" & example$age=="20-39"]<-2

     > strata[example$gender=="Male" & example$age=="40-59"]<-3

     > strata[example$gender=="Male" & example$age=="Over 60"]<-4

     > strata[example$gender=="Female" & example$age=="Under 19"]<-5

     > strata[example$gender=="Female" & example$age=="20-39"]<-6

     > strata[example$gender=="Female" & example$age=="40-59"]<-7

     > strata[example$gender=="Female" & example$age=="Over 60"]<-8

     > strata<-factor(strata, levels=c(1,2,3,4,5,6,7,8),

     > example<-cbind(example, strata)

接著創造出joint_population資料集，存放母體8個分層的聯合分布機率值。這邊可以是data frame也可以是table，這裡我們是用data frame來存放資料。


     > joint_population<-data.frame(strata=c(1,2,3,4,5,6,7,8),

     + percent=c(0.089, 0.140, 0.155, 0.111, 0.082, 0.133, 0.161, 0.129))

     > joint_population

       strata percent
     1      1   0.089
     2      2   0.140
     3      3   0.155
     4      4   0.111
     5      5   0.082
     6      6   0.133
     7      7   0.161
     8      8   0.129

有了樣本及母體聯合分布的機率值後，可以開始以survey套件裡的postStratify()進行事後分層加權。首先再一次用svydesign()宣告未加權的簡單隨機抽樣樣本。


     > unweighted2<-svydesign(ids=~0, data=example)

     > ps<-postStratify(design=unweighted2, strata=~strata, population=joint_population)

事後分層加權的結果如下，可以看到8個分層每一層的權重。我們一樣可以把權值納入原本的example樣本中以利後續分析。


     > stack(table(weights(ps)))

       values                 ind
     1     23 0.00578260869565217
     2     21 0.00766666666666667
     3     18 0.00777777777777778
     4     15  0.0103333333333333
     5      7  0.0117142857142857
     6     10              0.0129
     7      4             0.02225
     8      2              0.0555
     > example$weights3<-weights(ps)

最後驗證事後分層加權，可以看到樣本的聯合分布已與母體沒有差別。


     > psweighted<-svydesign(id=~0, data=example, weights=~weights3)

     > svytable(~gender+age, design=psweighted)

             age
     gender   Under 19 20-39 40-59 Over 60
       Male      0.089 0.140 0.155   0.111
       Female    0.082 0.133 0.161   0.129

分層事後加權也可以透過GREG一般回歸校正來達成。survey套件中的指令是calibrate()。同樣地，一樣得先指令母體分布的機率。因為是回歸校正，所以我們可以寫成回歸方程式的形式再進行校正：


     > pop_percent<-c(`(Intercept)`=0.089, "2"=0.140, "3"=0.155, "4"=0.111, "5"=0.082, "6"=0.133, "7"=0.161, "8"=0.129)

     > GREG<-calibrate(unweighted2, ~strata, pop_percent)

     > stack(table(weights(GREG)))

       values                 ind
     1      4             -0.2055
     2     23 0.00578260869565228
     3     21 0.00766666666666693
     4     18 0.00777777777777788
     5     15  0.0103333333333333
     6      7  0.0117142857142859
     7     10  0.0129000000000001
     8      2  0.0555000000000001

可以看到calibrate()計算的權重與postStratify()的計算結果是相同的，只不過因為剛才我們的指令是方程式的形式，所以在截距項出現負數。為了改善這個情況，可以用~strata-1來取代原本的指令，如此一來就相當於對照原本的8個分層，就不會出現負數了。


     > GREG2<-calibrate(unweighted2, ~strata-1, pop_percent)

     > stack(table(weights(GREG2)))

       values                 ind
     2     23 0.00578260869565217
     3     21 0.00766666666666671
     4     18 0.00777777777777777
     5     15  0.0103333333333333
     6      7  0.0117142857142857
     7     10  0.0128999999999999
     1      4             0.02055
     8      2              0.0555

簡介

多變數反覆加權

survey套件

anesrake套件

事後分層加權