杉野 勇 (SUGINO Isamu)　お茶の水女子大学・人間発達科学専攻・応用社会学コース担当

　N個からn個を選び出す「場合の数」の公式は次の通りである。

N=10, n=2とすると次のように計算できる。

N=40000, n=100の場合は次のようになる。

これは手計算は無理なので，Rの関数の choose(N, n) を利用する。

これは，「1.521272×10の302乗」と読む。e+302は10の302乗である。
MS-Excelでは，combin(N, n) という関数を利用する。

しかしこれらの関数は，Nやnがもう少し大きくなると，RやExcelで計算不能になる。Rでの計算結果を以下に示すので理解してほしい。

"Inf"というのは「無限大」という意味であるが，本当に無限大なのではなく，「簡単には計算できないくらい大きな数」という意味である。
　これを計算するには，「対数変換」を利用する方法がある。場合の数（組合せの数）を以下のように k と置く。

この両辺の常用対数を取り，「積の対数は対数の和に，商（割り算）の対数は対数の差に」変形できることを利用すると，次のように変形できる。

対数の和や差であればなんとかRで計算できるので，log(k)を計算してから，10のlogk乗を計算して，kを求めることができる（この計算でも多少近似的かも知れないが，おおよその大きさを把握することはできるだろう）。

　これが何をしているかがわからない場合は，c((N-n+1):N) だけを実行してみる，次に log10(c((N-n+1):N)) を実行してみる，そして sum(log10(c((N-n+1):N))) を実行してみる，nume-denom を計算した後で，floor(nume-denom) を計算してみる，のように部分的・段階的に実行していくと分かる場合が多い。
　Nやnを書き換えるだけで，このスクリプトを利用することができる。

　まず，上のスクリプトを使うと，4,894人から50人を選ぶ場合の数（可能な標本の数）は "7.78123×10の119乗"となる。
　4,894人分の模擬年収データをcsvファイルで提供するので，右クリックでファイルを保存し，getwd() で表示されるディレクトリ（フォルダ）にコピーして使用する。wdとはworking directory，「作業場所」の事である。（working directoryを変更したい場合には，setwd("C:/Users/yourname/Documents/R") のように変更先を指定する。）
（実際の調査データに多少の加工を施して作成している。）
　ファイルをworking directoryにコピーしたら，次のスクリプトで data01 として読み込もう。ついでに関数names( )，head( )，tail( )を例示する。
　この収入データ(data01$income)から，無作為に50個の数字を抽出してその平均や標準偏差を計算してみよう。
　無作為抽出なので，実行するたびに結果が変わる。この3行だけを何度か繰り返し実行してみると良い。

　以下のスクリプトを使うと，本文のヒストグラムに類似したグラフが描ける。実行するたびにグラフは少し変化する。

　テキストの，標準正規分布の95%区間のグラフを描くスクリプトを（少し簡素化して）示しておくので，pの値を.90や.99などに書き換えて実行してみよう。コマンドを一つずつ実行していくと，それぞれが何をするコマンドなのかが分かり易くなる。

　99.9%信頼区間を求めるための信頼係数を本文では3.29と書いているが，p <- .999; qnorm(1-(1-p)/2, mean=0, sd=1) とすれば，3.290527と求められる。
　標準正規分布とt分布を重ねて描いているグラフは，以下のスクリプトでほぼ再現できる（本文より少し簡略化してある）。 　自由度6のt分布の95%限界値は，qt(.975, df=6) より 2.446912 と求められる。詳細は以下のページを是非参照して欲しい。

　N=4894の母集団からn=200の標本を20回抽出して，その都度信頼区間を求め，それを本文のようなグラフに作図するスクリプトを示す。長いので，前半と後半に分ける。
　ここまでの前半で，サイズnの標本をtimes回抽出し，それぞれの標本における標本平均と標本標準偏差を計算している。nとtimesは変える事が出来る。
　以下の後半は，この標本平均と標本標準偏差から信頼区間を計算してグラフにするスクリプトである。alphaを書き換えてこの後半だけを繰り返すと，同じ20組の標本群について，異なる信頼区間のグラフを得る事が出来る。
　ここでまず，n=200，times=20，alpha＝.95として，前半と後半を実行したグラフを示す(11/Aug./2017修正)。

　次に，スクリプトの後半だけを，alpha=.90として実行したグラフを示す(11/Aug./2017修正)。

　信頼区間がそれぞれ少し狭くなっている事が見て確認出来るだろう。
　最後に，n=800，alpha＝.95として前半・後半を実行したグラフを示す(11/Aug./2017修正)。標本抽出からやり直すので上記の二つのグラフとの対応関係はないが，全体的に標本平均のバラつきがぐっと小さくなり，かつ信頼区間が狭くなっている事が見て取れるだろう。これが，標本サイズを大きくする事で得られる効果である。

　二項分布に関する関数を紹介しよう。仮に母比率を pi としたとき，サイズ n の標本において内閣支持者が w 人となる確率は，choose(n, w)*pi^w*(1-pi)^(n-w) として計算出来る。pi = .473，n=10，w=5 とした時の実行結果を以下に示す。


　ここで，以下のように w を少し工夫して，5人以上となる確率をそれぞれ求め，その合計も求めてみる。

この結果は，母比率が47.3%であるとき，10人の標本中5人が支持者である確率が24.3%，6人が支持者である確率が18.1%，10人全員が支持者である確率が0.056%，それらを全て合計した5人以上が支持者である確率が55.5%である事を示している。
　w <- 0:n とすると，有り得る全ての場合について計算される。確率の合計は必然的に1となる。

　同じ計算は，dbinom(w, n, pi) という関数を使って行う事が出来る。

　最低限のグラフを描くには以下のコマンドだけでよい。2つあるが，どちらも同じグラフとなる。前者の描画例だけを下に示す。

　これではあまりに素っ気ないので，以下のようにオプションを追加して描いてみる。

　本文ではn=25，π=.40のグラフを掲載しているが，以下にはn=50，π=.473の例を掲載しよう。オプション lwd は"line width"で，nが大きくなるとこの「棒線の太さ」を細くしないと棒が重なってしまうので，以下ではlwd=2としてある。

　本文同様，正規分布に非常に似通っている事が確認出来る。
　以下は，x軸を支持者数wではなく，標本比率w/nにしてグラフを描くスクリプトである。
　ここで少し脱線して，二値変数の分散の求め方を示しておこう。

　実際に二値変数を発生させて，標本分散と不偏分散を比較してみる。

標本分散を用いても不偏分散を用いても殆ど違いがない。
　以下，本文の世論調査（内閣支持率p=.412，標本サイズn=1248）の区間推定を行う。
本文よりも少し工夫したスクリプトを挙げるので見比べてみると良い。

　有限母集団修正項(fpc)，不偏分散とt分布を用いた場合は下記となる。


　本文で，n=10，p=.20のような場合は標準正規分布やt分布を用いて計算した信頼区間が不適切である事を述べた。F分布を用いた計算も含めて以下で示す。


　本来は二項分布を用いて求めるべきものであるので，以下の様に求めてみる。
pi <- seq(0, 1, by=0.000001) として，母比率 pi を0から1までの間で，0.0001%単位で変化させたヴェクトルを作る。
x=n*p（標本において，二値変数uが1となる人数）である時，関数 pbinom(x, size=n, prob=pi) は，母比率がpiである時にn人の標本においてu=1となる人数がx人以下となる確率を意味する。例えば，pbinom(4, size=10, prob=.473) の結果は 0.4449244 となる。
　95%信頼区間というのは，「母比率がpiである時に，標本比率がx/n以下（つまり該当者がx人以下）になる確率が2.5%以上であり，かつ標本比率がx/n以上になる確率も2.5%以上になる」範囲に由来する。
　該当者がx人以下になる確率が2.5%以上という条件は pbinom(x, n, pi) >= .025 と書ける。該当者がx人以上になる確率が2.5%以上という条件は，該当者がx人未満（x-1人以下）になる確率が97.5%以下になる事と等しいので，pbinom(x-1, n, pi) <= .975 と書ける。ただし，x-1がマイナスになるとエラーが生じるので，pbinom(max(0, x-1), n, pi) <= .975 としなければならない。
この条件を合わせると，pbinom(x, n, pi) >= .025 & pbinom(max(0, x-1), n, pi) <= .975 となる。
より一般化して，pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2 とする。
この条件を充たすpiの集合 pi[pbinom(x, n, pi) >= (1-alpha)/2 & pbinom(max(0, x-1), n, pi) <= 1-(1-alpha)/2] が(alpha*100)%信頼区間となるので，その下限minと上限maxを求めればよい。

F分布を使用して求めた区間とほぼ一致している事が確認出来る。もっと精度を粗くして良いなら，by=0.000001の部分をby=0.00001やby=0.0001にすると，計算時間がぐっと短縮される。
　標本1,248人中514人が支持者であった場合の信頼区間を，正規分布近似，t分布近似，F分布使用，二項分布スクリプトのそれぞれで求めたのが下記である。

いずれを用いても実用上はほぼ問題無い程度の違いと言えるだろう。

　本書ではユニット無回答(unit non-response)や項目無回答(item non-response)については殆ど扱っていないので，ここでだけほんの少し無回答の影響について検討してみよう。
　ごく簡単な例として，架空の単純なモデルを考えてみよう。
対象となる社会が高収入群と低収入群１：１であるとする。簡便化の為に，高収入群の平均年収を500（万円）に，高収入群の回収率（ユニット回答率）を50%に固定する。収入格差をdM（万円），回収率の差をdr（ポイント）とし，この二つの値だけを変化させて比較する。
詳しくコメントを付けたスクリプトは下記である。

この数値例では，二群の平均収入差は全体平均の50%，高収入群の平均の40%に相当し，回収率の差は20%ポイント，母平均400万円に対して標本平均は425万円となり，バイアスは25万円，バイアスの母平均に対する比は6.25%となる。
　収入格差と回収率差を変化させるとどうなるかを以下に示すので自分で読み解いてみよう。


　やや分かりにくいが，教育年数と年収の模擬データを発生させ，そこから人為的に欠損値を発生させ，どの様な結果が得られるかを例示してみる。
　まず教育年数eduを適当に尤もらしい確率で乱数発生させ，次にそのeduと強く順相関する年収データincome0を生成する。更に，年収が低い程高い値を取りやすい0〜1の範囲の変数missを生成させて，missの上位10%のケースについて，年収だけをNAにした変数incomeを生成する。最後に，それらの要約統計量と相関係数行列を示している。

　当然，欠損を含んだ年収incomeは分布が全体に高い方にずれ，教育年数eduとの積率相関は完備データ(income_full)より低下している。本来は教育年数と年収はもっと強く相関しているのに，低収入層が不釣り合いに分析から除外される為に，相関が弱まって見えてしまっていると云う（仮想的）状況である。これは「完備データ分析(complete case analysis)」の問題を示している。
　ここで，最も簡単な欠損値補定(imputation)の例を示して，より複雑な補定が必要である事だけ示唆しておこう。
　まずは有効ケースの平均でNAを置き換える方法の結果を示す。スクリプトと結果から解読してみよう。

よく考えれば当然だが，分布はより平均に集中してバラつきが小さくなり，教育年数eduとの相関も更に小さくなっている。主に低い年収であった筈のケースが全体平均年収で補定されているのであるから当然そうなる。欠測がランダムではない上の例の限りでは，平均値で補定する利点は無い。
　次に，収入の欠測を，収入を教育年数に回帰させた結果を用いた予測値で置換する方法を示す。偶然誤差（攪乱項）を含まない回帰補定である。回帰分析の結果を利用する為に，summary( )関数，names( )関数などで確認をしている部分も敢えて表示させている。それぞれの意味をしっかり解読しよう。
まずは回帰分析の結果を確認する部分。セミコロン";"で複数のコマンドを一行に書いて一度に結果を表示させているので区切りが分からないかも知れない。そう云う場合はセミコロンを使わずコマンドを一つ一つ実行して確認する。

　次がここでの本来の目的の回帰補定の結果と比較である。平均値補定程はバラつきは過小にならないが，それでも本来の分布からすると集中度が高く，平均値も未だ高い。教育年数との相関も，平均値補定程は弱くならないが，本来の相関からするとまだ弱い。欠損値生成（miss0が登場する行以降）をやり直すと結果は変化するが，何度か繰り返してもこの傾向は変わらない。自分で確認してみると良い。

　このほか，回帰補定に誤差（攪乱項）を交える方法や，複数の値で補定して異なるデータセット・分析結果を作り出して，その結果を統合する多重代入法(multiple imputation)などの方法がある。注目度が高いのは最後の多重代入法であるが，高度な内容になるのでここでは割愛する。

1)の答え：
答えは 0.817314 となる。

2)の答え：
[1] -1.962339
[1] 1.962339
となるので，「-1.962 以上，1.962 以下」が答えである。

3)の答え：
答えは， 58.92919 以上，61.07081 以下となる。58.9〜61.1で良いだろう。

4)の答え：
　標準正規分布近似，F分布計算，二項分布スクリプトの実施例を示す。

いずれを用いるかによって，41.6%〜48.4% もしくは 41.5%〜48.5%と微妙にずれるが，実際に必要な／期待される精度からすれば問題とする差ではない。

5)の答え： 問題文に書かれたスクリプトを実際に実行してみるのが課題である。以下に少し変えて再掲する。
問題文にある様に，サンプリングの部分からだけ繰り返してみる。