[R program]調査のためのサンプリング（survey sampling）
調査を行うときは全数調査を行えば統計的な推測を行う必要はありませんが、どうしてもサンプリングして部分集団の解析で全体を推測したいという場合があります。
例えば国や県単位の統計を出したいけど全数調査は無理だとか、データは全てあるけどPCのスペック上一度に解析できないとか。
そんなときのサンプリング方法として次のようなものがあります。
SASではproc surveyselectが用意されており、Rではsamplingライブラリがあります（SASのヘルプはここにあります）。
SASの場合はmethodで方法を指定できて、例えばmethod=ursでブートストラップサンプリング、またstrataを指定することで層別サンプリングをします。
RのsamplingライブラリはTille先生が提案したcube法（バランスサンプリング）を使っていいます。
- 教科書：Tille, Y. (2006), Sampling Algorithms, Springer.
- 論文：Deville, J.-C. and Tille, Y. (2004). Efficient balanced sampling: the cube method. Biometrika,91:893-912.
sampleライブラリでランダムサンプリングはsrswor関数、ブートストラップサンプリングはsrswr関数です。
層を示すcluster変数を作ってその層ごとにランダムサンプリングする場合は層別サンプリングになるのでstrata関数、共変量でのバランスサンプリングはsamplecube関数です。
層別サンプリングの例は、全集団の男女比が2:1だった場合にサンプル集団の男女比も2:1になるようにサンプルするということ。
バランスサンプリングの定義はこの資料によると、共変量のHorvitz-Thompson推定量が全集団の合計値に等しくなるということみたい。
ということは、サンプリング確率が等しいときは全体集団とサンプル集団との平均値も等しいということになる。
Rの例ではサンプリング確率をある変数の重み付けで変えていたけど、どのようなモチベーションのときに変えるのか今はよく分からない。
とりあえず分かったことは単純なランダムサンプリングよりも、cube法のバランスサンプリングを使った方がバイアスが入りにくいらしいということ。
まずはバランスサンプリングするためのコードです↓
library(sampling) nsample <- 50 data(MU284) #アイスランドのデータ（税収、党の議席数など） X=cbind(MU284$P75,MU284$CS82,MU284$SS82,MU284$S82,MU284$ME84,MU284$REV84) p <- rep(nsample/nrow(X), nrow(X)) (s=samplecube(X, p,1,FALSE)) MU284Sample <- MU284[s==1, ] rbind(summary(MU284)[4, ], summary(MU284Sample)[4, ]) HTestimator(MU284$RMT85[s==1], p[s==1]) sum(MU284$RMT85)
samplecube関数が本体で、ここでサンプリングのフラグを付けます。
引数1のXがバランスをとる共変量、引数2がサンプリング確率（この例では等確率）、引数3と4はオプションです。
s==1がサンプリングされる観測値なのでMU284Sampleに代入し、最後に全体集団の平均値と比較しています（Horvitz-Thompson推定量と全体合計の比較も）。
今回の例では正直言って本当にバランスがとれているのかイマイチ実感しにくかった。
論文と教科書を読まないとなぁ。
また、今回試したRのコードはこちら（今回は試行錯誤したのでちょっと長いです）↓
#------調査サンプリング（survey sampling） #教科書：Tille, Y. (2006), Sampling Algorithms, Springer. #論文：Deville, J.-C. and Tille, Y. (2004). Efficient balanced sampling: the cube method. Biometrika,91:893-912. #coreな関数：fastflightcube、samplecube、landingcube cluster=c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3) X=cbind(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)) (s=balancedcluster(X,2,cluster,1,TRUE)) #クラスターごとにサンプリングするよう（ロジックはsamplecube関数） unique(cluster[s[,1]==1]) (1:length(cluster))[s[,1]==1] s[s[,1]==1,2] example(calidev) #calibration（較正（こうせい）の評価） example(fastflightcube) example(samplecube) example(landingcube) #---samplecube X=cbind(c(1,1,1,1,1,1,1,1,1),c(1.1,2.2,3.1,4.2,5.1,6.3,7.1,8.1,9.1)) pik=c(1/3,1/3,1/3,1/3,1/3,1/3,1/3,1/3,1/3) s=samplecube(X,pik,order=1,comment=TRUE) (1:length(pik))[s==1] X=rbind(c(1,0,1,2),c(1,0,2,5),c(1,0,3,7),c(1,0,4,9), c(1,0,5,1),c(1,0,6,5),c(1,0,7,7),c(1,0,8,6),c(1,0,9,9), c(1,0,10,3),c(0,1,11,3),c(0,1,12,2),c(0,1,13,3), c(0,1,14,6),c(0,1,15,8),c(0,1,16,9),c(0,1,17,1), c(0,1,18,2),c(0,1,19,3),c(0,1,20,4)) pik=rep(1/2,times=20) ppp=rep(0,times=20) sim=100 #for accurate results increase this value for(i in (1:sim)) ppp=ppp+samplecube(X,pik,1,FALSE) ppp=ppp/sim print(ppp) print(pik) #---サンプリング確率の計算 a=1:20 pik=inclusionprobabilities(a,12) pik #---バランスサンプリング例1 data(MU284) #アイスランドのデータ（税収、党の議席数） pik=inclusionprobabilities(MU284$P75,50) #P75は1975年の人口→人口が多い州はサンプリング確率を大きくする cbind(pik, MU284$P75) #共変量が大きいとサンプリング確率も大きくなる X=cbind(MU284$P75,MU284$CS82,MU284$SS82,MU284$S82,MU284$ME84,MU284$REV84) (s=samplecube(X,pik,1,FALSE)) table(s) #1がサンプルされた州 HTestimator(MU284$RMT85[s==1],pik[s==1]) #Computes the Horvitz-Thompson estimator s=samplecube(matrix(pik),pik,1,FALSE) HTestimator(MU284$RMT85[s==1],pik[s==1]) sim=8 res1=rep(0,times=sim) res2=rep(0,times=sim) for(i in 1:sim){ cat("Simulation number ",i,"\n") #シミュレーション番号を表示する s=samplecube(X,pik,1,FALSE) #共変量でのバランスサンプリング res1[i]=HTestimator(MU284$RMT85[s==1],pik[s==1]) s=samplecube(matrix(pik),pik,1,FALSE) #サンプリング確率でそのままサンプリング res2[i]=HTestimator(MU284$RMT85[s==1],pik[s==1]) } summary(res1) summary(res2) (ss=cbind(res1,res2)) colnames(ss) = c("balanced sampling","uneq prob sampling") boxplot(data.frame(ss), las=1) #---バランスサンプリング例2 nsample <- 50 data(MU284) #バランスサンプリング X=cbind(MU284$P75,MU284$CS82,MU284$SS82,MU284$S82,MU284$ME84,MU284$REV84) (s=samplecube(X,rep(nsample/nrow(X), nrow(X)),1,FALSE)) XSample <- X[s==1, ] rbind(summary(X)[4, ], summary(XSample)[4, ]) #ランダムサンプリング s2=srswor(nsample, nrow(X)) XSample2 <- X[s2==1, ] rbind(summary(X)[4, ], summary(XSample2)[4, ]) #バランスサンプリングとランダムサンプリングの比較 sim=8 res1=rep(0, times=sim) res2=rep(0, times=sim) for(i in 1:sim){ cat("Simulation number ",i,"\n") #シミュレーション番号を表示する s=samplecube(X, rep(nsample/nrow(X), nrow(X)), 1, FALSE) #共変量でのバランスサンプリング res1[i]=HTestimator(MU284$RMT85[s==1], rep(nsample/nrow(X), sum(s))) #sの数がnsampleよりも大きくなることがある s=srswor(nsample, nrow(X)) #サンプリング確率でそのままサンプリング res2[i]=HTestimator(MU284$RMT85[s==1], rep(nsample/nrow(X), sum(s))) } (ss=cbind(res1, res2)) colnames(ss) = c("balanced sampling", "uneq prob sampling") boxplot(data.frame(ss), las=1) abline(h=sum(MU284$RMT85), lty=2) #Horvitz-Thompson推定量はこれに等しくなるべき #---層別サンプリング data=rbind(matrix(rep("nc",165),165,1,byrow=TRUE),matrix(rep("sc",70),70,1,byrow=TRUE)) data=cbind.data.frame(data,c(rep(1,100), rep(2,50), rep(3,15), rep(1,30),rep(2,40)), 1000*runif(235)) names(data)=c("state","region","income") table(data$region,data$state) s=strata(data,c("region","state"),size=c(10,5,10,4,6), method="srswor") getdata(data,s) #---ブートストラップサンプリング（bootstrap sampling） (s=srswr(3,10)) #---ランダムサンプリング（random sampling） (s=srswor(3,10))
