東京大学 | 大学院医学系研究科 | 国際保健学専攻 | 人類生態学教室TOP | 2nd

社会統計学（高崎経済大学非常勤講義）

（分割版目次）

トップ | 更新情報 | 研究と教育 | 業績 | 計算機 | 写真 | 枕草子 | 著者 | 目安箱 | 書評 | 非常勤

最終更新： 2002年2月19日 火曜日 12時12分

第１回 「社会統計学とはどういう学問か」（2001年9月20日）

講義概要

社会統計学＝社会調査データを分析する統計学
▼前期「社会調査」で触れたように，社会調査には「社会の諸側面の科学的探究」という側面がある。
▼社会を構成する全要素の把握を直接することはできないので，サンプリングによる標本調査が行われたり，いくつかの尺度で社会の断面を切り取ったりして調査データを得る。そのデータから，本来の社会全体を再構成し，また仮説を検証する原理として，統計学が用いられる。
統計（学）とは？
▼本来は，「家畜や他の財産の帳簿をつけるために原始人が木につけた刻み目」（ラオ）
▼その後，「ある国およびそこに生きている生命の状態や発展についての，最も完全で，最も根拠のある知識」（マルシャス） 「国家にとって必要不可欠な人口や経済的な情報の収集」（ウォナコット）
▼英語のstatisticsは，ラテン語で国家を意味するstatusを語源として18世紀半ばにドイツの哲学者アッヘンウォールが作った言葉の流用。
▼かさばり，雑然とした生データを，解釈をやさしくしたり種々の方策決定に用いるために纏め上げる，グラントの生命表やケトレーの度数分布図による発展
▼1834年，英国王立統計協会設立による「統計学」成立で，「人間に関係することがらで，数量で表現することが可能で，一般的な法則を導き出すのに十分なだけ積み重ねられたもの」
▼現在の広い定義としては，「不確実性を考慮した論理的推論」であり，すべての自然科学や社会科学で適用される科学的分析の技術となっている。（注：論理的推論には，帰納，演繹，アブダクションの３つがある）
不確実性＝ランダム（乱雑さ）
▼世の中のほぼすべての事象は不確実性を含んでいる
▼素粒子レベルでは物理法則も量子力学という形で不確実性を含むし，遺伝子の発現や社会における個人の行動なども，決して決定されてはいない
▼ランダムな数字の列＝乱数列〜次の数字が予想できない，意味のないでたらめな数の集まり，例えば，英国国勢調査の各教区の面積を表す数値の最初と最後を並べたものとか，20桁の対数表における15〜19桁目を並べたものとか，袋に入れた500個ずつの白ビーズと黒ビーズから，よく混ぜて１個ずつビーズを復元抽出したときの色の列など。線型合同法などによる擬似乱数列もある。
統計的分析の手順
▼目的を明確にする
▼データ化（エディティング，コーディング，データ入力）
▼記述統計を実施
▼作業仮説を明示する
▼仮説検定または区間推定を行う
因果関係と撹乱要因（例）
▼数値間に関連があるだけでなく，時間的前後関係などの条件を満たして初めて因果関係がいえる。
▼sufficient causeへの注目
統計解析の２つの原理
▼デザインに基づいた解析
▼モデルに基づいた解析
統計解析の道具
▼Excel：覚えておくといろいろ便利だが，ブラックボックスだし，それがないと何も出来ないのでは困る
▼R：オープンソースで無料で利用できる優れた統計ソフト（公式サイト）。統計計算のしくみをビジュアルで見るにも適している。
▼Rのインストール：Windows版1.31では，SetupR.EXEを実行して指示に従うだけでできる。
▼Rの文法のエッセンス

フォロー

専門用語が難しいので用語集の紹介あるいはフォロー希望
▼初回は全体像を描き出すために難しい専門用語が多く出ましたが，２回目からは基礎から順に知識を積み重ねるようにしたいと思いますので，上記の用語が全部はわからなくても大丈夫です。
▼なお，webサイトとしては，群馬大学社会情報学部青木先生のサイトの統計学用語辞典が充実していてお薦めです。
パソコンの扱いに慣れていないので基礎から教えて欲しい
▼パソコンそのものの動作原理や扱い方全般から始めると別に１コマ必要なので，コンピュータを使って統計解析をするのに必要な作業に絞って講義したいと思っています。
▼扱いに慣れるためにも，練習問題をできるだけ出します。


第２回 「確率と確率分布」（2001年9月27日）

わかったようでわからない「確率」について，徹底的に考えてみよう

講義概要

確率的な現象と統計的事象
▼どういう現象が確率的か？
▼こういう「不確かさ」に潜む法則性（長期間繰り返し観察したり，大集団で観察すると見られる）を考える学問を確率論と呼ぶ。この種の法則性をもつ現象を，「統計的事象」と呼び，統計的事象の確かさの度合いを示すのに便利なモノサシが「確率」
標本空間
▼統計的事象を捉えるには，「どんなことが起こりうるのか」という範囲を定めることが必要。
▼現象は一般に多面的で様々な観察方法がある。以下３点によって統計的現象を捉えた，記号化された結果の集合のことを「標本空間」と呼ぶ。
▼個々の結果の起こりうる可能性を示す数値（これを「確率」という）を考える。一般には「どの結果も同程度に起こる」と考える。各結果に対応付けられた確率は０から１までの数値であり，各確率の値の総和は１にならねばならない。
▼サイコロの目では，標本空間は｛１，２，３，４，５，６｝
事象の確率
▼問題は，個々の結果の可能性よりも，いくつかの結果が複合された集合（これを「事象」という）の起こる可能性がどのくらいか，ということ。つまり，事象とは，「標本空間の部分集合」である。
▼サイコロの例では，「目が偶数（丁）」とか「目が５以上」とか「目が１」とかいうことが事象である。
▼ある事象の確率は，その事象に含まれる各結果の生起確率の和である。従って，各結果の生起確率が等しい場合は，その事象に含まれる結果の場合の数をすべての場合の数で割ると，その事象の確率になる。サイコロの例では，「目が５以上」という事象の確率は，２／６＝０．３３３・・・である。
余事象・和事象・積事象・排反事象
▼起こりうるすべての結果の集合を「全事象」という。つまり，全事象は標本空間に等しい。
▼決して起こらない事象を「空事象」といい，空集合φで表す。
▼事象Eに対して，Eが起こらないという事象をEの「余事象」という。サイコロの例では，「目が偶数」という事象の余事象は「目が奇数」である。
▼事象EとFの少なくとも一方が起こるという事象を，EとFの「和事象」といい，E∪Fで表す。
▼事象EとFの両方が起こるという事象を，EとFの「積事象」といい，E∩Fで表す。
▼事象Eが起こればFは決して起こらないとき，EとFは「排反事象」であるという。EとFが排反事象なら， E∩F＝φである。
相互排反性と加法定理
▼サイコロで考えると，１回振ったとき「偶数の目が出る」という事象Eが起こる確率Pr(E)は，{2,4,6}の場合の数３を，{1,2,3,4,5,6}の場合の数６で割った値なのでPr(E)=0.5。
▼２回振って「少なくとも１回は偶数の目」の確率は？
事象の独立性と乗法定理
▼事象Eが起こっているときに事象Fが起こる確率を，Eが起こったときのFの「条件付き確率」といい，Pr(F|E)と書く。
▼「事象Eが起こった」うちの「EもFも起こった」場合を考えるので，Pr(F|E)=Pr(F∩E)/Pr(E)である。
▼事象Eと事象Fが互いに無関係（独立）なら，Fの条件付き確率Pr(F|E)は，Pr(F)と一致する。逆にいえば，Pr(F)=Pr(F|E)のときに事象Eと事象Fは互いに独立であるという。独立でないとき「従属である」という。
▼上記２つの式から，事象Eと事象Fが独立なら， Pr(F∩E)=Pr(F)×Pr(E)という“乗法定理”が成立。
確率の４つの定義
▼操作的接近＝統計的定義：数多く試したときの相対度数の極限。例えば，事象Eが起こる確率Pr(E)は，n回試したときにn1回事象Eが起こるとして，n1/nという相対度数が，nを無限大にしたときに漸近する値である。
▼対称的確率：サイコロの場合，６通りの目の出る確率はどれも等しくなければならず，その和は１でなくてはならないので，例えば１の目が出る確率は1/6となる。限定的かつ循環論法。
▼公理的客観確率：Pr(ei)>=0かつPr(e1)+Pr(e2)+・・・+Pr(eN)=1かつPr(E)=ΣPr(ei)を公理とする。数学的に厳密。
▼主観確率：観念的にも二度と繰り返すことのできない事象についての「見込み」を扱う。決定理論において重要。 law of large numbers
大数の法則（操作的接近の根拠）
▼Rのプログラムでは，とすれば，右図のように，試行回数を増やすとサイコロの特定の目が出る割合はある一定値に近づくことがわかる。数式で書くと，Pr(e1)=pとおけば任意の正の小さな数εに対して，lim(n→∞) Pr{|n1/n-p|<ε}=1ということで，これをベルヌーイ(Bernoulli)の大数の法則という。
▼同じことをExcelでやってみるには，=INT(RAND()*6)という式をセルの数だけコピーし，分析ツールの「ヒストグラム」を使えば不可能ではないが，大変面倒である。
確率変数と期待値
▼例えば，スロットマシンにコインを入れると，ごくたまに，投入金額の何十倍ものコインが出てくる。
▼マシン利用者全員に返ってくる賞金の合計を利用回数で割った値が，１回に期待される賞金額である。これを賭け金で割った値を「賞金還元率」と呼ぶ。すべての賭け事は胴元が儲かるようになっているのは，賞金還元率が100％未満だということである。宝くじでは40％，競馬では75％と言われる。
▼一般に，賞金額がx1, x2, x3, ・・・で，その賞金が得られる確率がp1, p2, p3, ・・・のように設定されたスロットマシンの期待賞金額Mは，M=x1p1+x2p2+x3p3+・・・で与えられる。
▼このスロットマシンのようなものを確率変数といい，期待賞金を一般に期待値と呼ぶ。
分散
▼期待賞金が同じでも，値動きの幅が小さいと一喜一憂の程度が小さく，逆に幅が大きいと滅多に当たらないが当たったときの喜びは大きくなる。つまり，ギャンブル性は，値動きの幅と，チャンスの大きさに依存している。
▼各賞金がどれくらい期待賞金から隔たりがあり，それを獲得できる可能性がどれくらいあるのかを見積もれば，ギャンブル性が表せる。
▼V=（マシンのギャンブル性）=Σ（期待値からの隔たり）×（可能性）という値が定義できる。このVを「分散」と呼ぶ。なお，各賞金額xと期待値Mの隔たりは，普通，差のを二乗した値D=(x-M)2で表す。
確率変数と確率分布
▼一般に，とりうる値の集合ｘ=(x1,x2,x3,・・・)と，それぞれの値が実現する確率ｐ=(p1,p2,p3,・・・)が与えられていて，事象としてｘのうちどれか１つの値のみ実現するとき，(ｘ,ｐ)という１セットを，「確率変数」と呼んで，Xで表す。
▼期待値はE(X)=μ=Σxipi
▼分散はV(X)=σ2=Σ(xi-μ)2pi
▼分散の平方根σを標準偏差と呼ぶ。
▼横軸にxの各々の値を示す位置に， pの各々の可能性を示す高さの棒を立ててみれば，これが確率変数の「確率分布」ということになる。
ベルヌーイ試行と２項分布
▼１回の実験でSかFかのどちらかが起こり，しかもそれらが起こる可能性が，Pr(S)=p，Pr(F)=1-p=qで何回実験しても変わらないとき，これを「ベルヌーイ試行」という。
▼ベルヌーイ試行をn回行って，Sがちょうどk回起こる確率は，Pr(X=k)=nCkpkqn-k
▼nCkは２項係数と呼ばれる。このような確率変数Xは，「２項分布に従う」といい，B(n,p)と表す。E(X)=np，V(X)=npqである。
正規分布
▼nが非常に大きい場合には，２項分布B(n,p)の確率Pr(X=np+d) という値が，1/√(2πnpq)・exp(-d2/(2npq))という値で近似できる。
▼一般にこの極限である，Pr(X=x)=1/√(2πσ2)・exp(-(x-μ)2/(2σ2))という形をもつ確率分布を正規分布と呼び，N(μ, σ2)と書く。
▼z=(x-μ)/σと置けば，Pr(Z=z)=1/√(2π)・exp(-z2/2)となる。これを標準正規分布と呼び，N(0,1)と書く。
▼統計学でよく使われる確率分布であるカイ二乗分布とかt分布とかF分布は，正規分布から導かれる。
練習問題
（１）８頭で出走する競馬のレースがあり，「どの馬が勝つチャンスも等しい」と仮定した場合，ある特定の馬が勝つと予想して当たる確率は1/8となるが，２回のレースの少なくともどちらか一方に当たる確率はいくらか？
（２）p=0.2として，いろいろなnについて２項分布のグラフはどういう形になるか？

フォロー

練習問題（１）の答え
▼確率1/8の独立な２つの事象の少なくともどちらかが起こる確率は，それぞれが起こる確率の和から両方とも起こる確率を引いたものなので，1/8+1/8-1/64=15/64
▼あるいは，「どちらも起こらない」ことはない，という確率として,1-(7/8*7/8)=15/64
２項分布についての補足説明
二項分布のシミュレーション▼正二十面体サイコロをn回（n=4, 10, 20, 50）投げて，１から４までの目が出る回数と考える（2001年12月27日訂正：R-jpメーリングリストでご指摘いただいたが，考えてみれば正五面体というのは存在しないのだった）。１回投げたときに１から４までの目が出る確率は0.2であるとして，試行1000セットの平均値を右図に示す。Rのプログラムは下記の通り。
二項分布の理論値▼各nについての理論的な確率分布は，Pr(X=k)=nCk0.2k0.8n-kより右図のようになる。Rのプログラムは下記の通り。
（2001年10月17日追記）実は，Rには，組み合わせの関数は組み込まれていた。ここではcombi(n,k)という形で定義したが，もとからchoose(n,k)がその機能をもっていた。また，jjjに0からn[j]までの値を代入してplotするよりも，barplot(xa)とする方が美しい度数分布図が書けることがわかった。
（2001年12月27日追記）R-jpメーリングリストで，もっと素直でエレガントなプログラムのやり方を教えていただいた。中林＠富士通総研さん，森＠学芸大さん，佐藤＠統数研さんありがとうございます。
（2002年1月4日追記）群馬大学の青木繁伸先生からも，applyを使って美しくまとめる方法をお教えいただき，一様分布のシミュレーションと尤度を最大にする母比率の実験のプログラムも教えていただいた。ここに掲載することをお許しくださったので，以下に示す。
  1. 一様分布のシミュレーション
    op <- par(mfrow=c(2,2))
    x <- apply(matrix(10^(2:5), ncol=1), 1, function(a, n) {y <- runif(a, 0, n); hist(y, n)}, 10)
    par(op)
    ●apply の最後の引数が上限を表す。すなわち，[0,n] の一様分布のシミュレーションになります。
  2. 二項分布のシミュレーション
    draw <- function(n, p, trial=1000){ x <- apply(array(n, trial), 1, function(n) {res <- runif(n, 0, 1) < p; length(res[res])}); hist(x, freq=F, br=seq(-0.5,n+0.5,1)) }
    op <- par(mfrow=c(2,2))
    x <- apply(matrix(c(4,10,20,50),ncol=1), 1, draw, 0.2, trial=1000)
    par(op)
    ●apply の最後から二番目の引数が母比率，trial は試行回数。
  3. 二項分布の密度関数
    draw <- function(n, p) {x <- 0:n; xa <- choose(n,x)*(p^x)*((1-p)^(n-x)); plot(x,xa)}
    op <- par(mfrow=c(2,2))
    x <- apply(matrix(c(4,10,20,50), ncol=1), 1, draw, 0.2)
    par(op)
    ●apply の最後の引数は母比率。
  4. 尤度を最大にする母比率の実験（第５回参照）
    bis <- function(n, r, p) { choose(n, r)*p^r*(1-p)^(n-r)}
    n <- 20
    r <- 2
    p <- 0:500/500
    y <- bis(n, r, p)
    plot(p, y, type="l")
    y


第３回「記述統計（１）：データの図示」（2001年10月4日）

代表値もやる予定だったのだが，時間が不足したので翌週に繰り延べた。

講義概要

データの図示
▼離散量の場合
▼連続量の場合
離散量データの例１
▼夫婦の子ども数：20組の夫婦について
2, 3, 1, 0, 3, 2, 2, 1, 1, 1, 2, 2, 1, 3, 2, 1, 0, 2, 1, 1
だった場合，どうまとめるか？
→度数分布図を描く。Rなら，
とすればきれいなグラフが書ける。
▼（10月22日追記）ただし，hist(child,c(-0.5,0.5,1.5,2.5,3.5))と書くほうが正当である。
▼（10月22日追記）histを使わずに，plot(table(child))またはbarplot(table(child))とする手もある。
離散量データの例２
▼国立社会保障・人口問題研究所が行った「第11回出生動向基本調査・夫婦調査」から，結婚持続期間と子ども数の関係
結婚持続期間０人１人２人３人４人５人以上
０〜４年 (1,273)42.644.712.10.600
５〜９年 (1,276)10.321.053.613.91.20
10〜14年 (1,287)5.511.654.225.23.30.3
15〜19年 (1,334)3.79.853.627.94.60.4
20〜24年 (1,419)2.38.157.028.93.40.4
25年以上 ( 559)1.312.258.324.03.90.4
→Excelで100％積み上げ横棒グラフにするとよい。
連続量データの例１
▼平成元年３月９日から４月２日の東京地区の最低気温
3.2, 3.1, 5.1, 4.8, 8.3, 9.8, 8.3, 6.6, 5.1, 3.8, 5.2, 5.6, 6.5, 5.7, 5.7, 7.4, 6.2, 7.0, 6.7, 5.7, 6.2, 6.0, 8.8, 10.7, 8.5
▼このデータを元にヒストグラムを書くには，Rなら，
とすればよい。区間まで自動的に設定される。もちろん，好きな区間に設定することもできるが，通常は自動で問題ない。Excelでも区間を設定すれば，「ツール」→「分析ツール」→「ヒストグラム」でヒストグラムが描けるが，グラフ上では区間表示がずれたように見えるので注意が必要。
連続量データの例２
▼同じデータについて，幹葉表示ならば手作業で，
3|2 1 8
4|8
5|1 1 2 6 7 7 7
6|6 5 2 7 2 0
7|4 0
8|3 3 8 5
9|8
10|7
とすることができる。
▼（10月22日追記）Rならば，stem(temp)でも幹葉表示はできるが，葉の部分がすべて0として表示されてしまう流儀なので，ヒストグラムよりも優れているポイントが今一つ発揮されないのが残念である。（2002年1月8日前の文を一部削除かつ追記）R-jpメーリングリストで目黒＠T&D投資顧問さんからご指摘いただいたが，この場合なら，stem(temp,2)とすれば，上で手作業でやったものとまったく同じ幹葉表示が得られることがわかった。stem(temp)はstem(temp,1)と同義で，stem(temp,2)の場合よりも幹の区分が粗くなる。なお，stem(temp,8)だと幹の部分が1桁増えるので，葉の部分がすべて0になるのは当然である。stem(temp,7)だと，幹が各整数温度について２区分される。ここで1とか2とか7とか8とかいう値（Rではscaleと呼ばれる）が具体的に何を意味するのかは，はっきりしないが，数を増やすほど幹の区分が細かくなるような値である。
連続量データの例３
▼同じデータについて，箱ヒゲ図ならばRを使って，
とすることができる。このデータではあまり箱ヒゲ図の威力がないが，外れ値があるときはその様子が一目瞭然なので便利な図示である。
連続量データの例４
▼学生10人の前期と後期の成績が下表のようであるとき，
学生12345678910
前期(x)8866644332
後期(y)9597574743
散布図を書くには，Rを使えば，
とすればよい。Excelでも範囲を選んで，「挿入」→「グラフ」→「散布図」とするだけなので，簡単である。

フォロー

箱ヒゲ図は何の役に立つのか？
外れ値がたくさんあるような場合に，全体の分布と外れ値の様子を両方見られるのが便利（例：下図）
box and whisker plot


第４回「記述統計（２）：代表値」（2001年10月11日）

講義概要

（2002年2月18日追記）ケアレスミスで，四分位範囲と書くべきところを四分位偏差と書いていたのを訂正しました。ケアレスミスを反省すると同時に，メールでご指摘くださいました目黒様に深く感謝申し上げます。

代表値
▼データ全体の情報を集約してみるために計算される，１つの値。
▼分布の特徴を表す値として，一般に，中心とばらつきの２つが考えられる。
分布の中心(central tendency)＝狭義の代表値
分布のばらつき
練習問題
▼以下のそれぞれについて，代表値を計算してみよ。
今回使った関数の一覧
求める代表値などEXCELの関数または手順
（範囲A1:Y1にデータがあるとして）		Rの関数または手順
（x <- c(...)などのやり方で変数xにデータを入れたとして）
最頻値離散データなら=MODE(A1:Y1)で良いが，連続量なら，ツール＞分析ツール＞ヒストグラムでヒストグラムを書いて最大度数のデータ区間を探し，その区間の中点を最頻値とする。hist(x)でヒストグラムを書いて最大度数のデータ区間を探し，その区間の中点を最頻値とする。hist(x,c(min(x),5,8,max(x)))などとすれば，xの最小値から5まで，5から8まで，8からxの最大値までという３つの区間で度数を計算させることができる。本来はhist(x,5)とすれば５つの区間という形の指定ができるはずなのだが，区間の数によってうまくいったりいかなかったりした。
なお，hist(x,plot=F)とすれば，グラフを書く代わりに数値を表示させられる。
中央値=MEDIAN(A1:Y1)median(x)
但し，xの中にNA（欠損値）を含む場合は，median(x,na.rm=T)，あるいはmedian(x[!is.na(x)])とする。以下同様。
平均値=AVERAGE(A1:Y1)
調和平均は=HARMEAN(A1:Y1)，幾何平均は=GEOMEAN(A1:Y1)で求められる。		mean(x)
調和平均は1/mean(1/x)，幾何平均はexp(mean(log(x)))で求められる。
範囲=MAX(A1:Y1)-MIN(A1:Y1)max(x)-min(x)
四分位偏差範囲=QUARTILE(A1:Y1,3)-QUARTILE(A1:Y1,1)IQR(x)，またはy<-quantile(x); y[4]-y[2]
または，fivenum(x)[4]-fivenum(x)[2]でも良い。
四分位偏差=(QUARTILE(A1:Y1,3)-QUARTILE(A1:Y1,1))/2IQR(x)/2，またはy<-quantile(x); (y[4]-y[2])/2
または，(fivenum(x)[4]-fivenum(x)[2])/2でも良い。
平均偏差=AVEDEV(A1:Y1)組み込み関数にはないが，
sum(abs(x-mean(x)))/NROW(x)で得られる。
不偏分散=VAR(A1:Y1)
（不偏でない分散は=VARP(A1:Y1)で得られる）		var(x)
不偏でない分散は組み込み関数にはないが，
sum((x-mean(x))^2)/NROW(x)で得られる。
不偏標準偏差=STDEV(A1:Y1)
（不偏でない標準偏差は=STDEVP(A1:Y1)で得られる）		sd(x)
不偏でない標準偏差は，
sqrt(sum((x-mean(x))^2)/NROW(x))で得られる。(*)
タブ区切りデータファイルの読み込みそのままドラッグ＆ドロップ１行目に変数名が入っているなら，
x <- read.delim("C:/My Documents/solomon.dat",header=T)とする(**)。
それぞれの変数は，例えばx$ageのようにして参照できる。１行目が変数名でなくすぐにデータである場合は，
x <- read.delim("C:/My Documents/solomon.dat",header=F)とする。この場合，変数名はx$V1, x$V2, ...として参照できる。
一々x$とつけるのが面倒なら，attach(x)とすればV1とかV21だけで参照できる。
カンマ区切りデータファイルの読み込みそのままドラッグ＆ドロップ１行目に変数名が入っているなら，
x <- read.csv("C:/My Documents/solomon.dat",header=T)とする(**)。１行目が変数名でなくすぐにデータである場合は，
x <- read.csv("C:/My Documents/solomon.dat",header=F)とする。
データの編集表にそのまま打ち込むde(x$V1,x$V5)などとすれば表形式で指定した変数の値を編集できる。表の上でマウスを右クリックすると操作メニューがでる。
データの書き出しファイルから保存を選ぶコンマ区切りでデータフレームxをマイドキュメントのsample.datに書き出すには，
write.table(x,"C:/My Documents/sample.dat",sep=",")とする。タブ区切りならsep="¥t"とすればよい。
(*) もちろん，不偏でない分散を出すときに，Vx<-sum((x-mean(x))^2)/NROW(x)などとして値を保存しておいて，
sqrt(Vx)とするのがエレガントである。
(**) ¥を/に置き換えたファイル名をフルパスで書く。ただし，２バイトコードが入ったディレクトリ名やファイル名は，文字化けするので使えないと思われる。Windows2000では，マイドキュメントフォルダは，ふつう，"C:/Documents and Settings/nakazawa/My Documents/"など（nakazawaのところにはユーザIDが入る）として参照できる。

フォロー

いろいろな代表値はどんな目的で使われる？
▼例えば，Rでdata(infert)として読み込まれるデータ（Rには，予めいくつかのデータが組み込まれているのです）は，infertとすればすべて表示されますが，その個々の値をすべて見て全体の様子を把握することは人間には難しいので，分布を図示したり代表値を計算したりするわけです。
▼代表値は，どれもデータの分布を１つの値に集約して示す目的で計算されます。１つの値に集約できると他のデータと比較するのに便利です（経年変化を見るとか）。例えば，infertのデータで，不妊になった人とならなかった人の年齢の平均値を比べてみるには，attach(infert)とした上で，mean(age[case==0])とmean(age[case==1])を計算させてみれば良いわけです。
▼データの分布の中心を表す値としては，分布が歪んでいたり外れ値が多い時は中央値，分布が正規分布に近ければ平均値が適当です。分布のばらつきを表す値としては，分布が歪んでいたり外れ値が多い時は四分位範囲や四分位偏差，正規分布に近ければ不偏標準偏差が適当です。
練習問題の答え
代表値の種類東京の最低気温ソロモン諸島住民の年齢ソロモン諸島住民の収縮期血圧
中央値6.233111
平均値6.434.43245112.2878
調和平均5.8335920(*)110.0773
幾何平均6.1223520(*)111.1651
範囲7.687107
四分位範囲2.22320
四分位偏差1.111.510
平均偏差1.45613.6423112.32760
分散3.4624292.8202263.0999
不偏分散3.606667293.3403263.6538
標準偏差1.86075317.1119916.22036
不偏標準偏差1.89912317.1271816.23742
変動係数29.1%49.7%14.4%
(*)ageは0を含むため，調和平均や幾何平均は0になってしまいます。
（注）ここは計算練習なので（変動係数以外は）桁数を長く表示していますが，本来は有効数字を考えて結果を出す必要があります。例えば，年齢はデータの有効数字が２桁しかない（元のデータが年単位の整数でしか入っていません）ので，平均値も34と書くべきです。


第５回「カテゴリカルデータの解析（１）」（2001年10月18日）

講義概要

カテゴリカルデータとは
▼カテゴリ＝離散量のデータをいう。質的なデータ，あるいは定性的なデータと言い換えてもよい。
▼名義尺度と順序尺度を含む。
▼例としては，以下のようなものがある。名義尺度でも２値データだと使える分析法が多い。
▼いつも実行しているvsだいたい実行しているvsときどき実行しているvsたまに実行しているvsまったく実行していない，は，順序があるだけでなく，ほぼ等間隔と見なせるので，間隔尺度となる（通常，整数値を振る）。間隔尺度は量的なデータである。
▼あまりに稀なカテゴリに対してはリコーディングの必要がある場合がある。
母比率に関する推定（１）
▼母比率とは，個々のカテゴリが母集団で占めるあろう割合である。通常，標本比率とほぼ一致する。
▼例えば，手元の容器の中に，数百個の白い碁石があるとする。この概数を手っ取り早く当てるために，数十個の黒い碁石を混ぜる。よくかき混ぜてから20個程度の石を取り出してみて（標本），その中で黒い石に対して白い石が占めていた比に，加えた黒い石の数をかければ，元の白い碁石の数が推定できる。生態学のリンカーン法のやり方と同じである。
▼例題：最初に混入した黒い石の数が40個，かき混ぜてから20個の石を取り出してみたら黒石2個，白石18個だった場合，元の白石の数はいくつと推定されるか？
母比率に関する推定（２）
▼上の手続きで得られるのは，点推定値である。もっともそれらしい，１つの値が出る。
▼ここで，この推定値がどれほど確からしいか？　を考えてみる。
▼黒石の割合（母比率）がpである容器から20個の石を取り出したときに，黒石がちょうど２個である確率は，ニ項分布に従うので，Rの式で書けば，choose(20,2)*p^2*(1-p)^18となる。この値が最大となるのは，p=0.1の時である（確かめよ）。
▼例えば，bis <- function(p) {choose(20,2)*p^2*(1-p)^18}としておいて，z <- c(1:100); for (i in 1:100) {z[i] <- bis(i/500)}; barplot(z)とすれば，0.002から0.2まで0.002刻みで母比率を変えたときの「黒石がちょうど２個である確率」が図示される。ちょうど50番目にピークがあって，p=0.1のときが最大になることがわかる。
▼40個入れて全体の0.1を占めるのだから，40/0.1=400が全体の数で，400-40=360が元の白石の数だと推定できる。
▼ただし，p=0.09だろうがp=0.11だろうが，黒石がちょうど２個である確率には大した差はない。だから，360という点推定値は，404（p=0.09のとき）とか324（p=0.11のとき）に比べて，それほど信頼性は高くない。
母比率に関する推定（３）
▼ビデオリサーチによれば，NHKの朝のテレビ小説「ほんまもん」の10月8日の関東地区の視聴率は22.9%であった。関東地区の調査対象世帯は600だから，137世帯が見ていたことになる。このとき，関東地区全体の真の視聴率（母比率）はどのくらいの範囲にあると推定されるか？　というのが問題である。
▼「ほんまもん」を見る／見ないという事象は各世帯独立に起こるとすれば，ニ項分布で考えることができる。母比率が137/600の時にちょうど137世帯が見たという確率は，choose(600,137)*(137/600)^137*(463/600)^463で，たかだか3.9%に過ぎない。
▼しかし，母比率が10%だったのに137世帯が見たという確率は，2.5*10^(-20)であり，まったくありそうにない。
母比率に関する推定（４）
▼137/600の前後適当な幅をとれば，かなり高い確率で，ちょうど137世帯が見た，という事象が起こることになる。この幅を「信頼区間」という。
▼95%の確率でちょうど137世帯が見たという事象が起こるための母比率の推定幅を，「95%信頼区間」という。
母比率に関する推定（５）
▼95%信頼区間を求めるには，下側2.5%の点と上側2.5%の点を求めればよいので，Rならz<-0; k<-0; while (z<0.025) {k <- k+1; z <- z+zz[k] }; kとして下側2.5%の点を求め，z<-0; k<-600; while (z<0.025) {k <- k-1; z <- z+zz[k] }; kとして上側2.5%の点を求めればよい（確かめよ）。
▼結果として，600世帯の調査で22.9%の視聴率だったら，母集団の視聴率（真の視聴率）の95%信頼区間は，19.8%から26.5%の間と言える。
母比率に関する推定（６）
▼ニ項分布は，ｎが大きいときは正規分布で近似できる。このことを利用すれば，母比率p，標本数（調査世帯数）nで，その標本の中で注目している属性をもつ標本数（「ほんまもん」を見た世帯数）をX，観測比率をP=X/nとすれば，Xが近似的に正規分布N(np,np(1-p))に従うことになる。正規分布の95%のサンプルは，平均±2標準偏差に入ることが既知なので，
Pr[np-2√(np(1-p))≦X≦np+2√(np(1-p))]=0.95
▼これから式変形すると，Pr[P-2√{P(1-P)/n}≦p≦P+2√{P(1-P)/n}]=0.95となるので，母比率pは95%の確率で範囲(P-2√{P(1-P)/n}, P+2√{P(1-P)/n})にあるといえる。即ちこれが，母比率pの95%信頼区間となる。
母比率に関する推定（７）
▼練習問題：　ある大学の正門の前で，ある朝登校して来る学生の男女比を調べてみたところ，300人中，女子学生が75人であった。この大学の女子学生の割合の点推定値と95%信頼区間を求めよ。
▼ヒント：　95%信頼区間の下限を求めるRの式は，75/300-2*sqrt(75/300*225/300/300)である。
▼注：　この推定には，朝登校して来る学生に男女の偏りがないという仮定があるので，実は真の値を過大評価することになっている。
▼付加的な問題：　では，どうすれば正しい推定ができるような標本がとれるか？

フォロー

理解を助けるために参考になりそうな本は？
▼統計の参考書を教えて欲しいという質問は大学院生からも良く受けるのですが，これが最高だというものはなかなかありません。統計的な考え方の基礎を知るには，鈴木義一郎「情報量規準による統計解析入門」（講談社）をお薦めします。推定や検定の個別の方法については，粕谷英一「生物学を学ぶ人のための統計のはなし―きみにも出せる有意差」（文一総合出版）または，多少高価ですが，浜田知久馬「学会・論文発表のための統計学　統計パッケージを誤用しないために」（真興交易（株）医書出版部）がお薦めです。統計学を本気でやりたければ，竹村彰通「現代数理統計学」（創文社）から進めばよいと思いますが，それなりに難解なので独学にはちょっときついかもしれません。
Rのコマンドをなかなか覚えられないので困っている。
▼大学以降に学ぶことは，実は全部は覚えなくていいことが多いのです。実際，教官だって，教えていることをいつでも全部暗記しているかといえば，そんなことはありません。大事なことは，どこを見れば必要な情報が書いてあるかを覚えておくことです。
▼Rについても，基本的なコマンドいくつかと文法の基本さえきちんと覚えておけば，後はうろ覚えで十分です。help(コマンド名）とか，ヘルプメニューからヘルプを見るなどすればいつでも詳細は確認できます。
Rで，data(infert)としてinfertと打つだけでデータが出てくるのはなぜ？
▼Rには統計計算の練習用に，たくさんのデータが予め組み込まれています。
▼例えば，1960年から1997年までの地球の平均二酸化炭素濃度の毎月の変化(co2)，タイタニック号に乗っていた人の船室や性別などのデータ(Titanic)，米国の人口の経年変化(uspop)などが利用可能です。infertもその一つです。
▼これらのデータは，原則として，既に論文などで発表済みで公開されているものです。
次週のプリントが前もって欲しい
そうできれば理想的なのですが，できていないので，無理です。正直に言えば，講義前日の夜にプリントを作成している状態です。申し訳ありません。


第６回「カテゴリカルデータの解析（２）」（2001年10月25日）

講義概要

母比率に関する検定（１）＝カイ二乗適合度検定
▼データの度数分布が期待される分布と一致しないという仮説（帰無仮説）が棄却できるかどうかを検定する。
▼基本としては，カイ二乗検定を行う。カテゴリ数が全部でn個あるとき，i番目のカテゴリの観測度数がOi，期待度数がEiであるとき，χ2=Σ{(Oi-Ei)2/Ei}が，自由度n-1のカイ二乗分布に従うことを利用して検定する（但し，不明な母数があるときは，その数も自由度から引く。但し，Eiが1未満のときは，通常カテゴリ分けをやり直す。
▼このようなχ2が大きな値になることは，観測された度数分布が期待される分布と一致している可能性が極めて低いことを意味する。
chi-square distribution of d.f. 1▼ちなみに，自由度１のカイ二乗分布は，右図のような形になる。つまり，χ2値が１より大きくなる確率は，約0.32ということである。
参考までに，自由度nのカイ二乗分布の確率密度関数は，x>0について，
fn(x) = 1/(2(n/2)Γ(n/2)) x(n/2-1) e(-x/2)
であり，平均n，分散2nである。なお，自由度(degree of freedom; d.f.)とは，標本の数から，前もって推定する母数の数を引いた値である。この例ならΣEiだけをΣOiとして推定することで，E1からEn-1まで定めてやっとEnが決まることになり，自由度は1を引く。
▼分布関数（確率母関数）は，確率密度関数を積分したものであり，図で見れば面積に当たる。その逆関数，つまり面積に対応する値を与える関数を分位点関数という。自由度１のカイ二乗分布の場合，Rでは，カイ二乗値xについて，確率密度関数がdchisq(x,1)，分布関数がpchisq(x,1)で与えられ，95%点を与える分位点関数がqchisq(0.95,1)で与えられる。これは，χ2値がqchisq(0.95,1)以下である確率が95%であることを意味する。逆にいえば，χ2値がqchisq(0.95,1)より大きくなることは，確率5%もない，滅多にないことである。言い換えると，「観測された分布が期待される分布と一致している可能性は5%もない」ということである。このようなとき，「観測された分布が期待される分布と違いがない」という仮説は有意水準5%で棄却されたといい，観測された分布は期待される分布と一致するとはいえない，と解釈する。
母比率に関する検定（２）
▼例題：ある病院で生まれた子ども900人中，男児は480人であった。このデータから，（１）男女の生まれる比率は半々であるという仮説，（２）男児1.06に対して女児1という割合で生まれるという仮説，は支持されるか？
▼ヒント：（１）の場合， χ2値は，X<-(480-450)^2/450+(420-450)^2/450として計算される。この値が自由度１のカイ二乗分布に従うので，1-pchisq(X,1)とすれば，男女の生まれる比率が半々である場合に900人中男児480人という観察値が得られる確率が計算できる。
▼その確率がきわめて小さければ（通常5%未満），統計的に意味があるほど有り得なさそうな（「統計的に有意な」という）現象であると考えて，仮説を棄却する。
▼実は，この場合は母比率が0.5であるとして２項分布で計算してもよい。480人以上になる確率と420人以下になる確率の合計がきわめて小さければ，「男女の生まれる比率は半々である」という仮説はありそうもないと考えてよいことになる。母比率0.5で起こる現象が，900回中ちょうど480回起こる確率は，choose(900,480)*0.5^480*0.5^420で与えられるが，Rには２項分布についてもカイ二乗分布と同じように確率密度関数を与える関数があり，この確率はdbinom(480,900,0.5)で与えられる。480人以上になる確率は，dbinom(480,900,0.5)+dbinom(481,900,0.5)+…+dbinom(900,900,0.5)となるが，これは分布関数を使えば，1-pbinom(480,900,0.5)で計算できる。420人以下になる確率は，dbinom(0,900,0.5)+dbinom(1,900,0.5)+…+dbinom(420,900,0.5)であり，分布関数を使って書けば，pbinom(420,900,0.5)である。従って，求める確率はこれらの和，即ち，1-pbinom(480,900,0.5)+pbinom(420,900,0.5)である。
▼２項分布から求めた値を1-pchisq(X,1)と比較し，ほぼ一致していることを確認せよ。
▼（２）についてはどうか？　自分で計算して確認せよ。
母比率に関する検定（３）
▼難しい練習問題：１日の交通事故件数を155日間について調べたところ，0件の日が79日，1件の日が61日，２件の日が13日，3件の日が1日，4件以上の日が1日だったとする。このとき，１日あたりの交通事故件数はポアソン分布に従うと言えるか？
▼一般に，稀な事象についてベルヌーイ試行を行うときの事象生起数がポアソン分布に従うことが知られている。交通事故は稀な事象であり，ある日に交通事故が起こる件数と翌日に交通事故が起こる件数は独立と考えられるので，交通事故件数はポアソン分布に従うための条件を満たしている。
▼Rでは，ポアソン分布の確率関数（離散分布の場合は，確率密度関数と言わずに確率関数というのが普通）は，dpois(件数，期待値）で与えられる。
▼ポアソン分布の期待値（これは母数である）がわからないので，データから推定すれば，(0×79 + 1×61 + 2×13 + 3×1 + 4×1)/155で得られる。Rで書けば，この値をEhhに保存するとして，cc<-c(0:4); hh<-c(79,61,13,1,1); Ehh<-sum(cc*hh)/sum(hh)となる。
▼従って，一日の事故件数が期待値Ehhのポアソン分布に従うとしたときの，事故件数0〜4の期待日数eppは，epp<-dpois(cc,Ehh)*sum(hh)で得られる。
▼こうなれば，X<-sum((hh-epp)^2/epp)としてカイ二乗値を求め，これが自由度３（件数の種類が５種類あって，ポアソン分布の期待値が母数として推定されたので，５−１−１＝３となる）のカイ二乗分布に従うとして1-pchisq(X,3)が0.05より小さいかどうかで判定すれば良さそうなものだが，そうはいかない。
▼eppの値を見ればわかるが，epp[cc==4]が1より小さいのである。期待度数が１より小さいときはカテゴリを併合しなくてはならないので，epp[cc==4]をepp[cc==3]と併合する。
▼即ち，epp[cc<3]->ep; epp[cc==3]+epp[cc==4]->ep[cc==3]; ep<-ep[!is.na(ep)]として期待度数の分布ep，hh[cc<3]->h; hh[cc==3]+hh[cc==4]->h[cc==3]; h<-h[!is.na(h)]として観測度数の分布hを得る。
▼後は，XX<-sum((h-ep)^2/ep)としてカイ二乗値を求め，1-pchisq(X,2)を計算すると，約0.187となることがわかる。即ち，与えられたデータの１日の交通事故件数がポアソン分布に従っている確率は約19%あり，ポアソン分布に従っていないとはいえないことになる。　自分で計算して確認せよ。

フォロー

Rでの確率密度関数，分布関数，分位点関数の一覧
分布の種類確率密度関数
(probability density function)		分布関数＝確率母関数＝累積確率密度関数
(distribution function = probability generating function = cumulative probability density function)		分位点関数
(quartile function)
カイ二乗分布dchisq(カイ二乗値, 自由度)pchisq(カイ二乗値, 自由度)qchisq(％, 自由度)
２項分布dbinom(生起回数, 試行回数, 母比率)pbinom(生起回数, 試行回数, 母比率)qbinom(％, 試行回数, 母比率)
ポアソン分布dpois(生起回数, 期待値)ppois(生起回数, 期待値)qpois(％, 期待値)
正規分布(1)dnorm(Ｚスコア，平均値，標準偏差)pnorm(Ｚスコア，平均値，標準偏差)qnorm(％, 平均値，標準偏差)
対数正規分布(2)dlnorm(Ｚスコア，対数平均値，対数標準偏差)plnorm(Ｚスコア，対数平均値，対数標準偏差)qlnorm(％, 対数平均値，対数標準偏差)
一様分布(3)dunif(値，最小値，最大値)punif(値，最小値，最大値)qunif(％, 最小値，最大値)
t分布dt(t値，自由度)pt(t値，自由度)qt(％, 自由度)
F分布df(F値，第１自由度，第２自由度)pf(F値，第１自由度，第２自由度)qf(％, 第１自由度，第２自由度)
(1)平均値と標準偏差は省略可能。省略時は標準正規分布(平均0, 標準偏差1)になる。
(2)対数平均値と対数標準偏差は省略可能。省略時は対数平均0, 対数標準偏差1になる。なお，対数平均とは自然対数をとった値の平均，対数標準偏差とは自然対数をとった値の標準偏差をいう。dlnorm(1)はdnorm(0)と等しい。
(3)最小値と最大値は省略可能。省略時は0と1になる。
（注）これらの分布関数に従う乱数を生成する関数もある。例えば，これまでにも何度か取り上げた，0から1までの一様乱数を1000個生成する関数がrunif(1000)であるのは，runif(1000,0,1)の省略形である。同様に考えれば，試行回数100回，母比率0.2の２項分布に従う乱数を1000個発生させるには，rbinom(1000,100,0.2)とすれば良いことがわかるだろう。
母比率に関する検定（２）の例題（２）の解答例
▼男児1.06に対して女児1の比率が正しい場合に期待される男児と女児の人数は，EM<-900*1.06/(1.06+1)及びEF<-900*1/(1.06+1)として得られる。
▼適合度を示すカイ二乗値は，X<-(480-EM)^2/EM+(420-EF)^2/EFとなる。
▼カテゴリ数は２なので自由度１のカイ二乗分布を使って検定すれば，男児1.06に対して女児1という割合で生まれるという仮説が正しい場合に男児480人，女児420人という結果が得られる確率が1-pchisq(X,1)として計算できる。
▼Rで計算するとこの値は，0.2598729となる。正しい可能性が25％以上もあるのだから，この仮説は棄却されない。
練習問題「20代後半の女性の未婚率の全国平均が約40％であるとき，ある県で20代後半の女性を400人，ランダムサンプリングしたところ，140人が未婚であったとすると，この県の未婚率は全国平均と一致していると言えるか？」の解答例
▼帰無仮説「この県の未婚率は全国平均と等しい」の下で期待される未婚者人数は400*0.4なので160人である。実際の既婚人数は260人，帰無仮説の下で期待される既婚人数は240人である。
▼従って，適合度を示すカイ二乗値は，X<-(140-160)^2/160+(260-240)^2/240となる。
▼1-pchisq(X,1)を計算すると，0.04122683となるので，帰無仮説は有意水準5%で棄却される。
▼従って，この県の未婚率は全国平均と一致しているとはいえない。
※なお，県Ａからのランダムサンプル4000人で未婚者が1600人，県Ｂからのランダムサンプル400人で未婚者が140人だったときに県Ａの未婚率と県Ｂの未婚率に違いはないか？　という問題ならば，県Ａからのサンプルと県Ｂからのサンプルはどちらも相手には含まれないので，母比率ではなく，独立２群間の比率の差の検定を行う。この場合，Rでは，u<-c(1600,140); a<-c(4000,400); prop.test(u,a)とすればよい（ただし，実はこの検定は，統計モデルとしては，全対象者4400人についての，県と未婚／既婚という２つのカテゴリカル変数の独立性の検定と同じことになる）。


第７回「カテゴリカルデータの解析（３）」（2001年11月8日）

今回の講義で使用したRのスクリプトはp07.Rとしてダウンロードできる。

講義概要

２つのカテゴリカル変数の独立性の検定
▼復習：例えば５人の人を対象に，性別と社会人か学生かを聞き取った場合，生データを書き出すと，
人性別社会人か学生か
A男性(M)社会人(W)
B女性(F)学生(S)
C男性(M)学生(S)
D男性(M)社会人(W)
E女性(F)社会人(W)
のようになるが，ここでいう「性別」とか，「社会人か学生か」というくくりで表される属性は，人によって変わる「変数」である。「性別」とか「社会人か学生か」のように離散値であり大小関係がない値からなる変数をカテゴリカル変数という。
▼２つのカテゴリカル変数の間に関係があるかどうかを検討したいとき，それらの組み合わせの度数を調べた表を作成する。これをクロス集計表と呼ぶ。５人くらいなら数え上げても良いが，普通，調査は数百人以上を対象に行うので，ソフトウェアにやらせる。
▼Rではクロス集計表は，table(カテゴリカル変数名１,カテゴリカル変数名２)で作成できる。Excel2000ではデータ範囲を選択してから「データ」の「ピボットテーブルとピボットグラフレポート」を選んで「完了」ボタンを押せば，新しいワークシートが開くので，薄い字で「ここに列のフィールドをドラッグします」と書いてあるところと「ここに行のフィールドをドラッグします」と書いてあるところに２つのカテゴリカル変数名をドラッグ＆ドロップし，「ここにデータアイテムをドラッグします」と書いてあるところにIDなど（この例では「人」）をドラッグすればクロス集計表ができあがる。
▼とくに，２つのカテゴリカル変数が，ともに２値変数のとき，そのクロス集計は２×２クロス集計表（２×２分割表）と呼ばれ，その統計的性質が良く調べられている。
▼研究デザインによって，検討すべき関係の種類はさまざまである。例えば，肺がんと判明した男性患者100人と，年齢が同じくらいの健康な男性100人を標本としてもってきて，それまで10年間にどれくらい喫煙をしたかという聞き取りを行うという「患者対照研究＝ケースコントロール研究」を実施した場合に，喫煙の程度を「全然吸わない」から「ずっとヘビースモーカーだった」まで何段階かのスコアを振れば，喫煙状況という変数と肺がんの有無という変数の組み合わせが得られる。
▼もちろん，それらが独立であるかどうか（関連がないかどうか）を検討することもできる。
▼しかし，むしろこのデザインは，肺がん患者は健康な人に比べて，どれくらい喫煙していた割合が高いか，を評価するためのデザインである。逆に，喫煙者と非喫煙者を100人ずつ集めて，その後の肺がん発生率を追跡調査する前向き研究（フォローアップ研究）では，非喫煙群に比べて，喫煙者ではどれくらい肺がんの発生率が高いかを評価できる。
▼これらの値を，リスク比やオッズ比という。しかし，断面研究で得られた２つの変数には，時間的な前後関係がないので，独立性の検定を行ったり，リスク比やオッズ比以外の関連性の指標を計算することが多い。

リスク比とオッズ比は，

疾病あり疾病なし
曝露ありa人b人
曝露なしc人d人
であるとき，前向き研究からの疾病リスク比＝a/(a+b)/(c/(c+d))=a(c+d)/c(a+b)，疾病オッズ比＝a/b/(c/d)=ad/bc，患者対照研究からの曝露リスク比＝a/(a+c)/(b/(b+d))=a(b+d)/b(a+c)，曝露オッズ比＝a/c/(b/d)=ad/bcとなり，リスク比は異なるがオッズ比は一致する。断面研究の場合は，有病オッズ比としてad/bcを計算する。社会統計学ではリスク比はあまり使われないが，オッズ比は使われる。一般に，要因Aのない人に比べて，要因Aがある人は，何倍くらい特性Bをもちやすいか，を示す値となる。後にロジスティック回帰の話で触れる。

関連性の指標には，ピアソンの相関係数とか，φ係数などがある。これらについては後日触れる。

なお，同じ質問を２回した場合に同じ変数がどれくらい一致するかについても，独立性の検定ができそうな気がするかもしれないが，してはいけない。この場合はtest-retest-reliabilityを測ることになるので，前期の社会調査で説明したクロンバックのα係数を計算するか，あるいはκ係数などの一致度の指標を計算するのが正しい。

▼独立性の検定は，２つのカテゴリカル変数の間に関連がないと仮定した場合に推定される期待度数を求めて，それに観測度数が適合するかを検定するカイ二乗検定である。もちろん，ある種の関連が仮定できれば，その仮定の元に推定される期待度数と観測度数との適合を調べてもいいが，一般に，２つのカテゴリカル変数の間にどれくらいの関連がありそうかという仮定はできないことが多い。そこで，関連がない場合の期待度数を推定し，それが観測値に適合しなければ関連がないとはいえない，と推論するのである。
特性Aあり特性Aなし
特性Bありa人b人
特性Bなしc人d人
▼標本が，上記の表のような度数をもっているとき，母集団の確率構造が，
特性Aあり特性Aなし
特性Bありπ11π12
特性Bなしπ21π22
であるとわかっていれば，N=a+b+c+dとして，期待される度数は，
特性Aあり特性Aなし
特性BありNπ11Nπ12
特性BなしNπ21Nπ22
であるから，χ2=(a-Nπ11)2/Nπ11+(b-Nπ12)2/Nπ12+(c-Nπ21)2/Nπ21+(d-Nπ22)2/Nπ22として，自由度３のカイ二乗検定をすればよいが，普通はπが未知なので，p(A∩B)=p(A)p(B)と考えて，各々の変数については特性のある人とない人の人数が決まっている（周辺度数が固定している）と考え，π11を，p(A)の推定値(a+c)/Nとp(B)の推定値(a+b)/Nの積として推定するなどすれば，
χ2=N(ad-bc)2/((a+c)(b+d)(a+b)(c+d))となる。この場合は，母数を２つ（p(A)とp(B)）推定したので，自由度1のカイ二乗分布に従うと考えて検定できる。
▼ただし，標本数が小さいときは，イェーツの連続性の補正を行うことが多い。カイ二乗分布は連続分布なので，各度数に0.5を足したり引いたりしてやると，より近似が良くなる。
▼この場合，χ2=N(|ad-bc|-N/2)2/((a+c)(b+d)(a+b)(c+d))が自由度１のカイ二乗分布に従うと考えて検定する。ただし，|ad-bc|がN/2より小さいときは補正の意味がないので，χ2=0とする。
▼実際の検定は，a=12, b=8, c=9, d=10などとわかっているときは，
x<-matrix(c(12,8,9,10),nc=2)として表を与え，chisq.test(x)とするだけでもできる（連続性の補正を行わないときはchisq.test(x,correct=F)とするが，通常その必要はない）。公式通りに計算した結果と比較せよ。
▼周辺度数が固定しているときに，すべての組み合わせを考えて，それらが起こる確率を直接計算し，与えられた表よりも偏った表になる確率をすべて足し合わせたものをフィッシャーの直接確率という。Rでは，fisher.test(x)で実行できる。Nが小さいときは，この方が正確である。
▼Rでは，simulate.p.valueというオプションを使えば，シミュレーションによってそのカイ二乗値より大きなカイ二乗値が得られる確率を計算させることもできるが，遅いコンピュータだと計算時間がかかるのが欠点である。
▼各度数が未知で，個人についての生の値が与えられているときは，chisq.test(A,B)などとする。度数を計算するには，table(A,B)とする。もちろん，chisq.test(table(A,B))としてもよい。Nが小さければ，この場合もFisherの直接確率の方が良い。fisher.test(A,B)で実行可能である。
▼infertのデータで，人工中絶(induced)と自然流産(spontaneous)の回数の関係を検討せよ。

フォロー

ＡありＡなし
Ｂあり12人8人
Ｂなし9人10人
であるとき，カテゴリ変数Ａとカテゴリ変数Ｂが独立である（無関係である）かどうかをカイ二乗検定せよ。chisq.test()を使った場合と，公式通りにやった場合で結果を比較せよ。フィッシャーの直接確率（フィッシャーの正確な検定ともいう）ではどうか？　の解答例
▼Rを使って実行すると，x <- matrix(c(12,9,8,10),nrow=2); chisq.test(x)の結果は，0.6386である。つまり帰無仮説「カテゴリ変数Ａとカテゴリ変数Ｂは独立である」が正しい確率は64%もあり，棄却されないことになる。
▼1-pchisq(39*(abs(12*10-8*9)-39/2)^2/(12+8)/(9+10)/(12+9)/(8+10),1)の結果は，0.638632となり，chisq.test()の結果と一致していることが確認された。
▼x <- matrix(c(12,9,8,10),nrow=2); fisher.test(x)の結果は，0.5273となる。カイ二乗検定のときより確率は低いが，それでも0.05よりはずっと大きいので，帰無仮説は棄却されない。
『infertのデータで，人工妊娠中絶回数(induced)と自然流産回数(spontaneous)は無関係といえるかを検討せよ』の解答例
▼Rでdata(infert); attach(infert); T<-table(induced,spontaneous); T; chisq.test(T); fisher.test(T);と実行すると，集計表が表示された後に，帰無仮説「infertのデータで人工妊娠中絶回数と自然流産回数が独立である」が成立する確率は，カイ二乗検定で0.001129，フィッシャーの直接確率で0.0004852となり，いずれにせよ１％よりずっと低いことがわかる。つまり，帰無仮説は棄却され，infertのデータでは人工妊娠中絶回数と自然流産回数は無関係とはいえない。
（確率の）数値をみて判断するということは，微妙な数値の場合，主観的判断に委ねるということ？
▼通常，検定を行う前に，何％より小さければ帰無仮説が棄却されるかという水準（有意水準）を決めておきます。習慣的に５％か１％が用いられます。
▼しかし，最近は有意水準との大小だけでなく，帰無仮説が正しい確率pを，そのまま結果として表示することが求められるようになってきました。そのため，pの値についての解釈は，0.06くらいだと「…傾向がある」と表現するといったことが起こっています。これは，ある意味では主観的といえるかもしれませんが，より厳密な表現になっただけともいえます。


第８回「カテゴリカルデータの解析（４）」（2001年11月15日）

講義概要

関連性の指標の種類
▼独立性の検定は，関連がない確率を示すだけで，「どのくらいの関連か」を示さないので，関連の大きさを示すには別の指標が必要。
▼前回触れたリスク比やオッズ比は，もちろん関連性の指標である。これらの95％信頼区間が１を含まなければ，統計的に有意な関連があるとみなす。
▼断面研究から得られる関連性の指標のうち，代表的なものとしては，リスク差，相対差，曝露寄与率，母集団寄与率，YuleのQ，ファイ係数がある。
▼同じ調査を同じ対象者に２回，あるいは２人の調査者による１つの事象の評価，あるいは１つの質問紙の中で本質的に同じことを聞いている質問項目についての一致度の指標としては，クロンバックのα係数の他にκ統計量（κ係数ともいう）が有名（注：κはカッパと読む）。
リスク比（復習）
▼リスクとは，１つの（好ましくない）事象が生起する確率である。
▼リスク比とは，厳密には前向き研究においてのみ求められる。処理群のリスクの対照群のリスクに対する比である。
▼いま，前向き研究において以下の表のような結果が得られたとする。
事象が起こる事象が起こらない
処理群X人m1-X人m1人
対照群Y人m2-Y人m2人
▼母集団でのリスクの推定値は，処理があったときπ1=X/m1，処理がないときπ2=Y/m2である。リスク比は，RR=π1/π2なので，その推定量は，(Xm2)/(Ym1)となる。
▼リスク比の分布はNが大きくなれば正規分布に近づくので，正規分布を当てはめて信頼区間を求めることができるが，普通は右裾を引いているので対数変換か立方根変換（Baileyの方法）をしなくてはならない。対数変換の場合は，95％信頼区間の下限はRR*exp(-qnorm(1-0.05/2)*sqrt(1/X-1/m1+1/Y-1/m2))，上限がRR*exp(qnorm(1-0.05/2)*sqrt(1/X-1/m1+1/Y-1/m2))となる。RRが大きい場合は立方根変換しなくてはいけないが，煩雑なので省略する。
オッズ比（復習）
▼オッズとは１つの事象が生起する確率と生起しない確率の比である。
▼オッズ比は，前向き研究においては処理群でのオッズの対照群でのオッズに対する比となる。断面研究では要因Aあり群でのオッズの要因Aなし群のオッズに対する比となる。ケースコントロール研究では，患者群のオッズの対照群のオッズに対する比となる。数学的にはどれもまったく同じ形となる。
▼事象生起確率が小さいとき，オッズ比はリスク比のよい近似となる。
▼オッズ比の点推定値は，前回示した通りOR=(ad)/(bc)である。95％信頼区間は，分布が右裾を引いているので対数変換かCornfieldの方法（４次方程式の解を求めねばならないので手計算は不可能）を用いる必要がある。対数変換なら，95％信頼区間の下限はOR*exp(-qnorm(1-0.05/2)*sqrt(1/a+1/b+1/c+1/d))，上限はOR*exp(qnorm(1-0.05/2)*sqrt(1/a+1/b+1/c+1/d))となる。
リスク差
▼曝露によるリスクの増減を絶対的な変化の大きさで表した値。RD=π1-π2
相対差
▼要因ももたず発症もしていない者のうち，要因をもった場合にのみ発症する割合。RelD=(π1-π2)/(1-π2)
曝露寄与率
▼真に要因の影響によって発症した者の割合。AFe=(π1-π2)/π1
母集団寄与率
▼母集団において真に要因の影響によって発症した者の割合。π=(X+Y)/(m1+m2)として，AFp=(π-π2)/π
YuleのQ
▼オッズ比を-1から1の値を取るようにスケーリングしたもの。Q=(OR-1)/(OR+1)
ファイ係数（ρ）
▼要因の有無，発症の有無を1,0で表した場合の相関係数。θ1，θ2を発症者中の要因あり割合，非発症者中の要因あり割合として，ρ=sqrt((π1-π2)(θ1-θ2))
κ統計量
▼一致度の指標。
２回目○２回目×合計
１回目○a人b人m1人
１回目×c人d人m2人
合計n1人n2人N人
という表から，偶然でもこれくらいは一致するだろうと思われる値は，１回目と２回目の間に関連がない場合の各セルの期待値を足して全数で割った値になるのでPe=(n1*m1/N+n2/N*m2)/N，実際の一致割合（１回目も２回目も○か，１回目も２回目も×であった割合）はPo=(a+d)/Nとわかる。ここで，
κ=(Po-Pe)/(1-Pe)と定義する。
▼κは，完全一致のとき１，偶然と同じとき０，それ以下で負となる。
▼κ統計量は，有意性の検定ができる。κの分散var(κ)=Pe/(N*(1-Pe))となるので，κ/sqrt(var(κ))が標準正規分布に従うことを利用して検定できる。つまり，帰無仮説「κが偶然一致する程度と差がない」が正しい確率が1-pnorm(κ/sqrt(var(κ)))となる。この確率が５％未満ならば，得られた一致度は有意水準５％で信頼できる（偶然の一致より大きい）といえる。
▼95％信頼区間もκ±qnorm(1-0.05/2)*sqrt(Po*(1-Po)/(N*(1-Pe)^2))として計算できる。
▼なおκは，２×２だけでなく，m×mのクロス集計表にも適用できる概念である。

フォロー

table型のobjectであるXを与えることで今回示した指標をすべて計算する関数，crosstab(X)を定義してみた。p08.Rとしてダウンロードできる。組み込みデータのinfertを使った使用例もついているので参考になるであろう。

入力が追いつかないことがあるので画面を戻して欲しいことがある
▼手を上げて言ってください。
定期テストの出題はRでの計算？
▼定期テストはしません。
リスク差，相対差，曝露寄与率などは，すべて同じことを求めていてやり方が違うだけ？　それとも求めているものも違う？
▼関連の程度を示す指標であるという意味では同じですが，意味することはそれぞれ違います。


第９回「カテゴリカルデータの解析（５）」（2001年11月22日）

講義概要

前回までは，主に２×２のクロス集計表を扱った。今回は，（１）２変数のどちらか，あるいは両方ともカテゴリ数が３つ以上の場合，（２）カテゴリカル変数が３つ以上の場合，への拡張を扱う。

カテゴリ数が３つ以上の場合のクロス集計表の解析
▼第７回にinfertのデータでやったが，カテゴリ数が３つ以上あるカテゴリカル変数AとBの間の独立性の検定は，２×２クロス集計表と同じく，chisq.test(A,B)あるいはfisher.test(A,B)で実行できる。データ数が多すぎなければ，fisher.test(A,B)を使う方が望ましい。
▼３つ以上の水準があってそれらの間に順序のあるカテゴリカル変数Aとカテゴリ数２つのカテゴリカル変数Bの間の関連を調べるには，コクラン＝アーミテイジ検定という方法を使うことができる。Rでは，prop.trend.test(変数Ａの各カテゴリで変数Ｂのカテゴリが「あり」であった人数のベクトル，変数Ａの各カテゴリ人数のベクトル，変数Ａの各カテゴリに割り付ける順序のベクトル）
カテゴリカル変数が３つ以上の場合
▼３つ以上のカテゴリカル変数間の関係を調べたい場合には，２通りが考えられる。１つは，２つのカテゴリカル変数間の関係に関心があって，それが第３のカテゴリカル変数がどのような値であっても変わらないかどうかを確かめたい場合である。この場合，第３のカテゴリカル変数は注目している２つのカテゴリカル変数の関係に影響を与えるかもしれない要因（交絡要因，あるいは撹乱要因ともいう）であり，その影響を排除した解析を行いたいわけである。この場合の常套手段は，第３のカテゴリカル変数で層別したクロス集計表を複数作って，それを併合することである。
▼一方，あるカテゴリカル変数の値と，それ以外の複数のカテゴリカル変数の関係を同時に見たい場合もある。この場合の常套手段は，ロジスティック回帰分析を行うことである。ロジスティック回帰分析については次回扱う。
層別したクロス集計表〜なぜ層別するか？
▼なぜ層別しなくてはいけない場合があるかといえば，「シンプソンのパラドックス」と呼ばれる現象が有名。
▼シンプソンのパラドックスとは，層別した各層のクロス集計表で見られる関係が，すべてを合計したクロス集計表で見られる関係と違っている場合をいう。例えば，飲酒と心筋梗塞の関係について，喫煙の有無が交絡要因になっているとして，喫煙者と非喫煙者を別々にクロス集計表を作ったら，下表のようになったとしよう。
喫煙者心筋梗塞非喫煙者心筋梗塞
飲酒ありなし飲酒ありなし
あり816あり6336
なし2244なし74
▼喫煙者と非喫煙者を別々に，飲酒と心筋梗塞の独立性の検定をすると，どちらの層でもカイ二乗値が0に近い値（連続性の補正をしなければちょうど0）になり，飲酒と心筋梗塞には関連がないという帰無仮説は棄却されない。ところが，喫煙の有無で層別しなければ，クロス集計表は下のようになる。
心筋梗塞
飲酒ありなし
あり7152
なし2948
▼この表でカイ二乗検定をすると，χ2=6.84 (p=0.0089)となって，飲酒と心筋梗塞には関連がないという帰無仮説は棄却されてしまう。しかし，上でやった層別の解析からわかるように，この関連は，喫煙が飲酒とも心筋梗塞とも関連をもっていることによる見せかけである。
▼Rでこの解析を行うには，smoker<-matrix(c(8,22,16,44),nc=2); nonsmoker<-matrix(c(63,7,36,4),nc=2); chisq.test(smoker); chisq.test(nonsmoker); total<-smoker+nonsmoker; chisq.test(total)とすればよい。
▼上の例とは逆に，各層では関連が見られるのに，交絡要因である第３の変数を無視して合計すると関連が見えなくなってしまったり，各層で見られたのとは逆の関連になってしまう場合もある。
層別したクロス集計表〜どの層でも同じ関係が見られるかどうか？
▼「どの層でも関連がない」を帰無仮説，「すべての層で同じ向きの関連がある」を対立仮説として検定を行う。コクラン＝マンテル＝ヘンツェルのカイ二乗検定という方法を用いる。Rでは，mantelhaen.test(Ａ,Ｂ,層別変数)またはmantelhaen.test(３次元の表)で計算できる。
▼２×２の場合，マンテル＝ヘンツェル要約オッズ比というものがよく使われる。これは，k層でのオッズ比ORk=(akdk)/(bkck)であるときに，層別変数の影響を調整した共通オッズ比ORMHが，ORMH=(Σ(akdk/Nk))/(Σ(bkck/Nk))として推定できるという理論である。ORMHの信頼区間は，Pk=(ak+dk)/Nk，Rk=(akdk)/Nk，Qk=(bk+ck)/Nk，Sk=(bkck)/Nkと書くと，VFRBG=(1/2)[(Σ(PkRk))/((ΣRk)^2) + (Σ(QkRk+PkSk))/((ΣRk)(ΣSk)) + (Σ(QkSk))/((ΣSk)^2)]として，95%信頼区間の下限がORMHL=ORMH*exp(-qnorm(0.975)*sqrt(VFRBG))，95%信頼区間の上限がORMHU=ORMH*exp(qnorm(0.975)*sqrt(VFRBG))となる。
▼上のシンプソンのパラドックスの例についてRを使ってコクラン＝マンテル＝ヘンツェルのカイ二乗検定をするには，上のプログラムに続けて，mix<-matrix(c(smoker,nonsmoker)); dim(mix)<-c(2,2,2); mantelhaen.test(mix)とすればよい。χCMH2=0 (p=1)となって，帰無仮説は全然棄却されないことがわかる。シンプソンのパラドックスの例は２×２なので，自動的にマンテル＝ヘンツェル要約オッズ比も計算され，それが１である（95％信頼区間が[0.4552028, 2.1968231]）とわかる。

フォロー

検定をして出てきた数値の有意性などが未だにわからない。どんな数値のときに関連をもつか，もたないのかの区別のつけ方をもう一度説明して欲しい。
▼単純化して言えば，統計的検定は次の手順を踏みます。
R関連の本があったら教えてください。
▼R関連の本は，たぶん日本語では売られていません。でも，S-PLUSという市販ソフトとできることはほとんど同じなので，S-PLUSの参考書が使えます。
▼しかし，Rの入門には，Notes on Rという文書がwebからダウンロードできるので，これを読むのがいいでしょう。日本語と英語があって，日本語版はhttp://isw.main.eng.hokudai.ac.jp/‾yama/R/notes.htmlからダウンロードできます。
エクセルとRではどちらが統計ソフトとしてよく使われている？
▼Excelはもともとが統計ソフトではなく，表計算ソフトです。ツールの分析ツールメニューからいくつかの統計計算はできますし，マクロによる統計パッケージが別に売られていますが，できる統計分析も限られています。しかし，日本の会社などで実務で使われているソフトとしての使われ方は，１，２を争うものだと思いますので，Excelの扱いを覚えておいて損はないでしょう。Rは，フリーですし，S-PLUSやSも含めれば世界的にはかなり使われているものですが，日本の会社で実務でRを使っているケースはほとんどないと思います。しかし，SASやSPSSに匹敵するほど多様な統計処理ができますので，その統計ソフトとしての実力はExcelより遥かに上だと思います。


第10回「回帰分析，重回帰分析，ロジスティック回帰分析」（2001年11月29日）

講義概要

回帰と重回帰：
▼1984年から1993年までのプロ野球の１試合平均入場者数（単位：千人）の推移は下表のようになっている。
年次84858687888990919293
セリーグ28292931313131323534
パリーグ13121618212322242424
▼これを見ると，途中から伸び悩んでいるが，ある程度直線的に増加しているように見える。このように，ほぼ比例関係にある量的なデータ群をうまく代表する直線を求めるには，「各データの点から直線までのずれの大きさの合計」を最小にすればよい。この場合，「年次」が予め決まっている値であるのに対して，入場者数は測定値であり，誤差を含む可能性があるので，「データ点から直線までのずれ」を評価するには，データ点と直線の最短距離よりも，データ点から直線に垂直に下ろした線分の長さの二乗和を使う方がよい。この，ずれを最小にする直線を，「年次を独立変数，入場者数を従属変数とする回帰直線」と呼ぶ（計算方法は，数学的には検量線の求め方（pdf形式）を参照されたい）。
A sample of linear regression
Rでは，nenji<-c(84:93); central<-c(28,29,29,31,31,31,31,32,35,34)として，plot(central‾nenji); z<-lm(central‾nenji); abline(z)とすれば上の図が描かれる。なお，plot(nenji,central)とplot(central‾nenji)は同義である。
▼一般に，２組あるいはそれ以上の対で観測されるデータが与えられていて，そのうちの１つの変数（従属変数，または目的変数という）の変動の様子が残りの変数（独立変数，または説明変数という）の変動の仕方と，ある一定の関係を保ちながら変化していることがある。この場合，従属変数の変動を捉えるのに独立変数の変動に関する情報が利用できる。この「一定の関係」を求めることを回帰分析という。関係が直線のとき線型回帰分析といい，lm()が線型回帰分析を行う関数である。
▼独立変数が複数あって，それらによって同時に従属変数の変動を説明しようとするとき，とくに重回帰分析というが，同じ関数で実行できる。例えば，独立変数が年齢(変数名ageとする)と体重(変数名wtとする)，従属変数が収縮期血圧(変数名SBPとする)であるなら，その重回帰分析をするには，lm(SBP‾age+wt)とすればよい。
ロジスティック回帰
▼前回問題とした，３つ以上のカテゴリカル変数の関係を同時に調べるのにも，この重回帰の考え方が使える。しかし，カテゴリカル変数の場合，値の範囲が有限であるため，単純な直線関係にはならない。そこで，一般化線型モデルを用いる。関数名はglm()である。独立変数に量的な変数が含まれ，従属変数がカテゴリカル変数の場合，カテゴリカル変数があるカテゴリになる確率pを，log(p/(1-p))という形で無限小から無限大までの範囲をとるように変換して（ロジット変換と呼ぶ），一般化線型モデルを適用するとできる分析がロジスティック回帰分析である。
例：
▼前回から計算練習に使っているuser.dat（エコポイント・チェックを使って作った仮のデータ）でQ4とQ10の関係には年齢と性別の影響が両方あったので，年齢を量的変数として，NQ4がTRUE（いつもまたは大体買い物袋を持参している）の人はFALSEの人に比べて，どれくらいNQ10がTRUE（いつもまたは大体マイカーを避けて公共交通を利用している）の割合が高いかを考える。
▼NQ4がTRUEである割合をpとすると，log(p/(1-p))=β0+β1AGE+β2SEX+β3NQ10というモデルが考えられる。TRUEを1，FALSEを0とすると，他の変数を調整した上でのNQ10の効果は，
NQ10がTRUEの場合のNQ4がTRUEである割合をp1とすれば，
log(p1/(1-p1))=β0+β1AGE+β2SEX+β3
NQ10がFALSEの場合のNQ4がTRUEである割合をp0とすれば，
log(p0/(1-p0))=β0+β1AGE+β2SEX
となり，辺々引けば
log(p1/(1-p1))-log(p0/(1-p0))=β3
結局この左辺は，log[(p1/(1-p1))/(p0/(1-p0))]となる。これは対数オッズ比である。つまり，AGEとSEXの影響を調整した上でNQ10がTRUEである人はFALSEの人に比べて，exp(β3)倍，NQ4がTRUEになりやすいと言える。普通の言葉で書き直せば，「性別と年齢の影響を調整した上で，いつもまたは大体マイカーを避けて公共交通を利用している人は，そうでない人に比べて，いつもまたは大体買い物袋を持参する傾向がexp(β3)倍大きい」ということである。なお，β0〜β3は，一般化線型モデルの係数として推定される。
Rでのプログラム：
result<-glm(NQ4‾NQ10+AGE+SEX, family=binomial(logit)); summary(result)
出力：
Call:
glm(formula = NQ4 ‾ NQ10 + AGE + SEX, family = binomial(logit))
Deviance Residuals:
Min1QMedian3QMax
-0.9472-0.7143-0.5803-0.44781.9311

Coefficients:
EstimateStd. Errorz valuePr(>|z|)
(Intercept)-2.24960.6354-3.5400.000399 ***
NQ100.29690.52910.5610.574639
AGE0.27670.31690.8730.382463
SEX0.73720.51091.4430.149070
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 100.08 on 99 degrees of freedom
Residual deviance: 96.64 on 96 degrees of freedom
AIC: 104.64
Number of Fisher Scoring iterations: 3
▼とりあえずは，Coefficients:のEstimateにβ0〜β3が得られ，Pr(>|z|)にそれぞれが0である確率が得られている，と判断すればよい。
▼モデル自体の説明力に関しては，AIC（赤池の情報量基準）を見る。一般に，未知パラメータ数kの適当なモデルMを仮定したときの最大対数尤度をMLL(M)と書くと，
AIC(M)=-2MLL(M)+2k
と定義される。AICは，モデルの当てはまりの悪さを示す値であるから，複数のモデルのAICを比べた場合は，AICが小さい方を採択するわけである。

フォロー

ロジスティック回帰で出てきたlogは，高校数学で出てきたのと同じ？
▼その通り，対数のことです。統計学でlogが出てきて底が明示されていない場合は，自然対数と考えていいです。
一般化線型モデルとは？
▼線型モデルを一般化したものです。線型モデルは，普通の量的変数の間の直線回帰を指します。これを一般化すれば，t検定，分散分析，共分散分析，回帰分析，判別分析，正準相関分析などをすべて共通の数学モデルで扱うことができ，それを一般化線型モデルといいます。英語ではGeneral Linear Modelというので，SASでもPROC GLMですが，Rでの関数名もglmです。
▼基本的には，Y=β0+βX+εという形で表されます。Yが従属変数群，Xが独立変数群，β0が切片群，βが係数群です。
Rの日本語版のようなものはないのですか？
▼Rには日本語版はありません。ただ，統計用語は英語の方が曖昧さがありませんし，Rのマニュアルはかなり日本語訳されていますから，現状でもそれほど問題はないと思います。
あるデータをもっているときどんなふうに分析したらいいということはどうやって決める？
▼一般的な統計的分析の手順としては，第１回で示したことの繰り返しになりますが，もう少し詳しく書くと，
という感じになります。上の表に出てきた，量的な変数に関連する部分を，この講義の残りの回で説明する予定です。


社会統計学第11回「量的変数の解析(1)〜差の検定と相関」（2001年12月6日）

講義概要

母平均値と標本平均の差の検定
▼標本平均E(X)=mean(X)と既知の母平均値μXの差の検定は，母分散VXが既知のとき，z0=abs(E(X)-μX)/sqrt(VX/n)が標準正規分布に従うことを使って検定できる（注：nはサンプル数，abs()は絶対値を得る関数）。
▼VXが未知のときは標本の不偏分散SX=Σ(Xi-E(X))2/(n-1)=var(X)を使って，t0=abs(E(X)-μX)/sqrt(SX/n)が自由度n-1のt分布に従うことを使って検定できる（暗黙の仮定として，ランダムサンプルで，母集団の分布が正規分布であることが必要）。
▽未知の母平均値の信頼区間の推定はこの裏返し。
※Rでは，量的変数Xについて，平均値の95%信頼区間の下限はmean(X)-qt(0.975,n-1)*sqrt(var(X)/n)で得られ（ただし，予めn<-NROW(X)としておく），上限は，mean(X)+qt(0.975,n-1)*sqrt(var(X)/n)となる。なお，t.test(X,mu=既知の母平均値)とすれば，上記の検定と推定を両方やってくれる。
（練習1）昭和53年の国民栄養調査によれば３歳男児の平均身長は95.5cmであった。昭和53年にある保健所の３歳児健診に来所した326人の男児の平均身長が96.3cmで不偏分散が25.7であったとき，全国平均と比べて発育に差があるか検定せよ（出典：豊川裕之，柳井晴夫「医学・保健学の例題による統計学」現代数学社，1982）。
独立２標本の分布の位置の差の検定：
▼標本調査によって得られた独立した２つの量的変数XとY（サンプル数が各々nXとnYとする）について，母分散が既知で等しいVである場合は，z0=abs(E(X)-E(Y))/sqrt(V/nX+V/nY)が標準正規分布に従うことを使って検定する。
▼母分散が未知の場合は以下の通り。
  1. F検定（分散が等しいかどうか）：２つの量的変数XとYの不偏分散の大きい方を小さい方で割ったF0=SX/SYが第1自由度nX-1，第2自由度nY-1のF分布に従うことを使って検定する。
    ※Rでは1-pf(F0,nX-1，nY-1)が有意確率になるが，var.test(X,Y)で実行できる（この場合は，Rが勝手に入れ替えてくれるので，Xの不偏分散の方がYの不偏分散より大きいかどうか気にしなくてもよい）。詳細は，help(var.test)で説明が得られ，example(var.test)で使用例が実行される。Rの関数は基本的にそういう仕様になっているのが親切である。設計思想が優れていると思う。
  2. 分散に差がない場合：母分散S=[(nX-1)SX+(nY-1)SY]/(nX+nY-2)より，t0=abs(E(X)-E(Y))/sqrt(S/nX+S/nY)が自由度nX+nY-2のt分布に従うことを利用して検定する。
  3. 分散が差がある場合（Welchの方法）：t0=abs(E(X)-E(Y))/sqrt(SX/nX+SY/nY)が自由度φのt分布に従うことを使って検定する。但しφ=(SX/nX+SY/nY)2/[(SX/nX)2/(nX-1)+(SY/nY)2/(nY-1)]。
▼両側検定と片側検定：独立２標本XとYの平均値の差の検定をする場合，それぞれの母平均をμX, μYと書けば，その推定量はμX=mean(X)=ΣX/nとμY=mean(Y)=ΣY/nとなる。両側検定では，帰無仮説H0:μX=μYに対して対立仮説H1:μX≠μYである。H1を書き直すと，「μX＞μYまたはμX＜μY」ということである。つまり，t0を平均値の差／標準誤差として求めると，t0が負になる場合も正になる場合もあるので，有意水準５%で検定して有意になる場合というのは，t0が負でt分布の下側2.5%点より小さい場合と，t0が正でt分布の上側2.5%点より大きい場合の両方を含む。t分布は原点について対称なので，結局両側検定の場合は，上述のように差の絶対値を分子にして，1-pt(t0,自由度)によって得られた確率を２倍すれば有意確率が得られることになる。片側検定は，先験的にXとYの間に大小関係が仮定できる場合に行い，例えば，Xの方がYより小さくなっているかどうかを検定したい場合なら，帰無仮説H0:μX≧μYに対して対立仮説H1:μX＜μYとなる。この場合は，t0が正になる場合だけ考えればよい。有意水準5%で検定して有意になるのは，t0がt分布の上側5%点より大きい場合である。
※上記(2)と(3)について，Rではt.test(X,Y）で平均値の差の検定をやってくれるし，t.test(量的変数‾カテゴリ変数)でカテゴリ変数で群分けした量的変数の２群間の平均値の差の検定をやってくれる。デフォルトではWelchの方法になる。t.test(X,Y,var.equal=T)とすれば等分散の場合の計算になる。片側検定をしたい場合は，alternativeという指定を追加する。例えば，X＞Yが対立仮説なら，t.test(X,Y,alternative="greater")とする。指定しなければ両側検定である。alternativeに指定できる文字列は，greaterの他にはlessとtwo.sidedがある。
▼Mann-WhitneyのU検定（Wilcoxonの順位和検定と同じ。分布が歪んでいる場合に行う）：２群を混ぜて個々の値に小さい方から順位を与え，その順位を使って検定統計量Uを計算するのだが，詳細は省略する。
※Rではwilcox.test(変数1,変数2,paired=F）またはwilcox.test(量的変数‾カテゴリ変数)でやってくれる。
▼コルモゴロフ＝スミルノフ検定（分布の違いを検出したい場合）:説明は省略する。
※Rではks.test(変数1,変数2）で可能。
対応のある２標本の平均値の差の検定：
▼対応のあるt検定：例えば前回例示したプロ野球入場者数のデータで，セリーグの入場者数とパリーグの入場者数は，各年度について両方の値があるので，対応が取れている。このような場合は，独立２標本の平均値の差の検定をするよりも，対応のある２標本として分析する方が切れ味がよい（差の検出力が高い）。
※Rではt.test(変数1,変数2,paired=T)で実行できる。
▼ウィルコクソンの符号順位検定（分布が歪んでいる場合）：説明は省略。
※Rではwilcox.test(変数1,変数2,paired=T)で実行できる。
相関：
▼上述のセリーグの入場者数とパリーグの入場者数は，どちらかがどちらかを説明するという関係ではないが，互いに関連をもって変化している可能性がある。この関連の程度を表す指標が相関係数である。
▼ピアソンの積率相関係数（もっとも普通に用いられる。「相関係数」といえば，普通これを指す）: XとYの共分散をXの分散とYの分散の積の平方根で割った値である。式で書けば，相関係数の推定値rは，r=Σ(Xi-E(X))(Yi-E(Y))/√（Σ(Xi-E(X))2Σ(Yi-E(Y))2）。rが０とは有意に異なるという仮説を検定するためには，r=0という帰無仮説の下で，検定統計量t0=r√(n-2)/√(1-r2)が，自由度n-2のt分布に従うことを利用して検定すればよい。
※Rでは，r<-cov(X,Y)/sqrt(var(X)*var(Y)); n<-NROW(X); t0<-r*sqrt(n-2)/sqrt(1-r^2); として，2*(1-pt(t0,n-2))で有意確率が得られるが，下記のようにcor.test関数を使う方が簡単である。
▼スピアマンの順位相関係数（分布が歪んでいたり，外れ値がある場合に使う）: 値を順位で置き換えたピアソンの積率相関係数になる。
▼ケンドールの順位相関係数（スピアマンと同様）: (A-B)/[n(n-1)/2]。A=順位の大小関係が一致する組の数，B=不一致数
※Rでは，cor.test(変数1,変数2,method="pearson")でピアソンの積率相関係数が得られ，それがゼロである確率pも得られる。pearsonの部分をspearmanに変えればスピアマンの，kendallに変えればケンドールの，順位相関係数とp値がそれぞれ得られる。

フォロー

練習問題の解答例
▼母平均95.5，標本平均96.3，不偏分散25.7，n=326を代入すれば，t0<-abs(96.3-95.5)/sqrt(25.7/326)となる。（講義中に言い忘れたが，）両側検定にするために確率を２倍して，p<-2*(1-pt(t0,325))が有意確率となるので，p=0.00466である。帰無仮説「標本平均と母平均が一致する」が棄却されるので，この３歳児健診の対象児は，全国平均と有意に差があると言える。
▼セリーグとパリーグの入場者数の相関は，central<-c(28,29,29,31,31,31,31,32,35,34); pacific<-c(13,12,16,18,21,23,22,24,24,24)として，cor.test(central,pacific)とすれば，以下の出力が得られる。 
	        Pearson's product-moment correlation
	
	data:  central and pacific
	t = 4.5086, df = 8, p-value = 0.001979
	alternative hypothesis: true correlation is not equal to 0
	sample estimates:
	      cor
	0.8471066
この結果から，ピアソンの積率相関係数の推定値がr=0.847で，それがゼロであるという帰無仮説が有意確率p=0.001979で棄却されるので，有意な相関があるといえる。
講義中にwebサイトを検索したのはなぜ？
webサイトにはRでのいろいろな確率分布の関数の一覧表を置いてあるので参照しました。
試験は？
試験はありません。
母分散と不偏分散とは何ですか？
母分散とは，母集団の分散です。不偏分散とは，標本から推計される母分散の推定値で，各標本の値から平均を引いて二乗したものの総和（「偏差平方和」といいます）を（標本数−１）で割ったものになります。


社会統計学第12回「量的変数の解析(2)〜分散分析と多重比較」（2001年12月13日）

講義概要

一元配置分散分析
▼前回は，２群間で量的変数の平均値を比べた。今回はそれを多群間の比較へと拡張する。
▼カテゴリ変数Xの各カテゴリによって，量的変数Yの値に差があるかどうかを検定したい場合，分析は二段階で行うのが普通である。つまり，まず，「量的変数Yの値に対してカテゴリ変数Xの効果があるかどうかを検定」し，有意な効果があった場合に第二段階として，「どれとどれのカテゴリの間で量的変数Yに差があるのか」を検定するのである。第一段階を一元配置分散分析，第二段階を多重比較という。
▼例えば，７つの村A, B, C, D, E, F, Gで健診をやって，身長や体重を測ったとしよう。このとき，データは次のようになる。 
	PID, VILLAGE, HEIGHT, WEIGHT
	1, A, 173, 67
	2, A, 166, 70
	3, A, 180, 55
	：
	：
	87, D, 164, 58
	88, D, 175, 90
	:
	:
	139, G, 166, 61
▼村落によって身長に差があるかどうかを検定したいならば，身長(HEIGHT)という量的変数に対して，村(Village)というカテゴリ変数の効果があるかどうかを一元配置分散分析することになる。Rでanova(lm(HEIGHT‾VILLAGE))とすれば，例えば次のような結果が得られる。 
	Analysis of Variance Table
	
	Response: HEIGHT
	           Df Sum Sq Mean Sq F value Pr(>F)
	VILLAGE     6   1171     195  0.5405 0.7766
	Residuals 132  47672     361
▼このような結果の表を分散分析表という。級間変動（VILLAGE間でのばらつき，つまり総平均に対する級平均のばらつき）が級内変動（Residualsに出ている，級平均に対する各データのばらつき）に対してどれだけ大きいかがF valueに得られる。Pr(>F)が検定の結果得られる有意確率で，この場合だと，VILLAGEの効果は有意でないことになる。つまり，身長に村落間差はないことになる。
▼数学的には，次のように説明できる。A村のN1人の身長がX11, X12,..., X1N1，B村のN2人の身長がX21, X22,..., X2N2，以下同様にして，G村のN7人の身長がX71, X72,..., X7N7だとする（総人口N1+N2+...+N7=N人とする）。村毎の平均身長をGX1, GX2,..., GX7と書き，全体の平均をGXAと書くことにする。このとき，Σ[i=a to b] f(i)で「f(a)からf(b)までf(i)を足した値」を表すことにすれば，総変動（総平方和）S=Σ[i=1 to 7]Σ[j=1 to Ni](Xij-GXA)2，級間変動（群間平方和）SV=Σ[i=1 to 7]Ni(GXi-GXA)2，級内変動（郡内平方和，または誤差平方和ともいう）SE=Σ[i=1 to 7]Σ[j=1 to Ni](Xij-GXi)2となる（S=SV+SEである）。自由度は，群の効果に関してPV=7-1=6で，残差の効果に関してPE=N-7である。よって，級間分散VV=SV/PV，級内分散VE=SE/PEと推定でき，F統計量F0=VV/VEが，第１自由度PV，第２自由度PEのF分布に従うことを使って検定できる（p=1-pf(F0,PV,PE)として有意確率が得られる）。分散分析とは，全体のばらつきを群間の違いという意味のはっきりしている分散VVと，それでは説明できない分散，つまり誤差分散であるVEに分けて比べることを意味する。
▼念のため上の数値例の値が数式のどれに当たるかを説明すると，PVが6，PEが132（Nが139），SVが1171，SEが47672，VVが195，VEが361，F0が0.5405，pが0.7766である。
▼なお，この例のように，群分けをするカテゴリ変数が１つの場合を，一元配置分散分析(ONE-WAY ANOVA)，２つの場合を二元配置分散分析(TWO=WAY ANOVA)，３つなら三元配置分散分析(THREE-WAY ANOVA)などと呼ぶ。二元配置以上の場合は，カテゴリ変数間での交互作用による影響を調べるための交互作用項がモデルに入ってくるし，その従属変数への効果を見るために母数モデルと変量モデルの違いを区別しなくてはならない。また，量的変数による交絡がある場合は共分散分析(ANACOVA)をすることになる。
多重比較
▼仮に上の例でPr(>F)が0.05より小さく，有意差があったときに，具体的にどの村の間に有意差があるのかを調べるには，単純に考えると，t.test(HEIGHT[VILLAGE=="A"],HEIGHT[VILLAGE=="B"])を繰り返せば良さそうだが，７つの村でこれをやると７つから２つを取り出すので21回の比較をすることになり，１つくらいは偶然有意差が出てしまう比較があるかもしれない。そのため，有意確率を調整するとか自由度を調整するなどの方法で，偶然による有意差の出すぎをなくすのが多重比較である。
▼Rでは，pairwise.t.test(HEIGHT,VILLAGE,p.adjust.method="bonferroni")とすれば，ボンフェローニの方法で調整した結果を出してくれる。p.adjust.methodを指定しなければHolmの方法になる。

フォロー

練習問題の解答例
▼data(infert); attach(infert); anova(lm(parity‾education))の結果より，p<<0.01なので教育歴の既往出生児数への効果は有意である。そこで，pairwise.t.test(parity,education)とすると，教育を受けた年数が0〜5年の人と6〜11年間の人，0〜5年の人と12年以上の人はそれぞれ有意に異なるが，6〜11年間の人と12年以上の人には有意な差がないことがわかる。
統計学を専門に勉強している学校はある？
▼生物統計学や心理学分野では比較的たくさんあります（米国に比べるとずっと少ないですが）。社会統計学はそれらより珍しいと思います。 大学院なら，統計数理研究所という統計学専門の研究所があります。
Rの出力の読み取り方がわからない
▼よくわからないときは，Prとかpという数値を探してみましょう。検定は，有意確率と帰無仮説さえきちんと把握しておけば読めます。
社会統計学は何の役に立つ？
▼谷岡一郎『「社会調査」のウソ』に紹介されているように， マスメディアや本など，権威ある機関の報告でも，統計はいいかげんだったり，恣意的な解釈をしたりという例にでくわすことは枚挙に暇がないほどです。統計学の素養があれば，そこで鵜呑みにせずに「待てよ」と思うことができるでしょう。


社会統計学第13回「量的変数の解析(3)〜多変量情報の集約」（2001年12月20日）

講義概要

二元配置分散分析
▼前回は，カテゴリ変数が１つで各カテゴリ間の量的変数の平均値を比べた。今回はカテゴリ変数が２つの場合に拡張する。
▼例えば，７つの村A, B, C, D, E, F, Gで健診をやって，身長や体重を測ったという先週のデータで，村落間の違いだけでなく，性差も調べたいとする。この場合，性によって村落の身長や体重への影響が異なっているかもしれないので，交互作用効果も調べる必要がある。
▼二元配置分散分析では，母数モデルと変量モデルのどちらであるかを区別しなくてはならない。分散分析では，個々の値が母平均とカテゴリ変数の主効果と誤差項の和に分解できると考えるが，通常は，誤差項をカテゴリ変数のカテゴリや何番目のデータであるかによらず，期待値はゼロ，分散は一定の母分散に等しいと仮定する。これが母数モデルである。カテゴリ変数のカテゴリが個体であるような場合，それがある母集団からのランダム標本であると仮定して，各カテゴリの主効果の期待値がゼロ，分散がそれぞれ異なるとしてモデルを立てることになる。これが変量モデルである。
▼上の例では母数モデルでいいので，身長の場合なら， anova(lm(HEIGHT‾VILLAGE+SEX+VILLAGE*SEX))とする。
▼三元配置以上でも同じように分析できる。
主成分分析
▼たくさんの変数があって，それらの間のいくつかに関連はあるけれども，そのまま見ていたのでは解釈が難しい，という場合に，変数を少ない数の互いに独立な変数（主成分）の線形結合として表せると便利である。そのための分析法が主成分分析である。Rでは，mvaライブラリの中のprcompまたはprincompで実行できる。
▼library(mva)としてから，example(prcomp)とか，example(princomp)とかやってみると，実行例を見ることができるので，やってみよう。

フォロー

世論調査から分析するとか，もっと身近な例題とかでやって欲しかった。
▼確かにその方が，具体的な分析のイメージがつかめたかもしれません。今回は準備不足で世論調査的な生データが用意できなかったのでできなくて済みませんでした。受講者が多ければ，講義の最初の方で生データを集めるという方法も取れたのですが，結果的に少人数講義になってしまったので，１度だけ集めたデータも，結局使えなくて残念でした。

期末試験はありません。出席状況と毎回の演習によって成績評価します。