統計手法の実地への適用限界

統計学を実地・現場に適用するとき，様々な適用限界があると言われる

しかし，それらの多くは適用限界ではなく，それぞれの統計手法が要求する条件を実地・現場が満たせない，あるいは分析者が手法に詳しくないため条件を知らないという場合が多いようである

このページでは，具体例を挙げてまとめてみよう

サンプルサイズより多い変数を使おうとする†

• 重回帰分析などのときに R だと，以下のようになる。見る人が見ればわかるが，「なぜ全部の変数が使われないんだろう」と R を疑う人もいるだろう。他の統計解析ソフトなら，エラーメッセージを吐いて何の結果も出ないということがあり，症状はそちらの方が重篤。

  > df<-data.frame(matrix(rnorm(50),5)) # 5行10列のデータ
  > lm(X10~., df)
  
  Coefficients:
  (Intercept)           X1           X2           X3           X4           X5  
    -0.009964    -1.551818     0.204167     0.623022     0.480082           NA  
           X6           X7           X8           X9  
           NA           NA           NA           NA  

  R でも，たとえば上述のデータフレームを使った因子分析の場合には，以下のようなエラーメッセージを吐く。

  > factanal(df, factors=2)
   以下にエラーsolve.default(cv) :  システムは数値的に特異です：条件数の逆数 = 1.43914e-19

  初心者には，訳の分からないエラーメッセージだろう。エラーメッセージというのはエラーの原因だけを呈示するのではなく，問題を解決する手がかりを含める方が良いのだが。
  主成分分析は princomp と prcomp の 2 つの関数がある。同じデータフレームをprincomp に食わせると，エラーの理由を示してくれる。エラーメッセージに従って変数の数を減らすかサンプルを増やすかできればこのエラーメッセージを活かすことができるが，そうも行かない現場では頭を抱えることになる。優れたエラーメッセージは「prcomp を使いなさい」というものであろう。

  > princomp(df)
   以下にエラーprincomp.default(df) : 'princomp' can only be used with more units than variables

  prcomp を使うと，ちゃんと結果が出る。
  オンラインヘルプを見ると，その理由が分かる。"The calculation is done by a singular value decomposition of the (centered and possibly scaled) data matrix, not by using eigen on the covariance matrix. This is generally the preferred method for numerical accuracy."
  要するに，二つある関数のうち，princomp は忘れて，いつも prcomp を使えばよいと言うだけのことである。

  > prcomp(df)
  Standard deviations:
  [1] 1.970954e+00 1.508694e+00 9.627209e-01 6.518143e-01 1.490813e-16
  
  Rotation:
             PC1         PC2         PC3         PC4          PC5
  X1  -0.1630607 -0.35419085  0.26011758 -0.25523000 -0.218177746
  X2  -0.6963896 -0.08856275 -0.01932273  0.02367637 -0.454431214
  X3   0.2001387 -0.13401011 -0.04328719 -0.65994118 -0.246128241
  X4  -0.1582884 -0.33047848 -0.47055678 -0.26277838  0.233447590
  X5  -0.4083304 -0.30093337 -0.09721462 -0.03975937  0.664612920
  X6  -0.2603722  0.13359471 -0.45727464  0.28071928 -0.321760610
  X7   0.2602059 -0.36773253 -0.19073917 -0.11185164 -0.145751525
  X8   0.2728444 -0.40807421 -0.10803231  0.52708369  0.002717044
  X9  -0.1456160  0.49559293  0.07640619 -0.20100460  0.252171997
  X10  0.1595599  0.28941005 -0.66047476 -0.13640846 -0.027055976

全サンプルが同じ値を持つ変数を使おうとする†

統計調査で，もし男のデータしか集まらなかったら，性別という変数を作るだろうか？
いや，作るはずがない。
このような明確なことが現場ではしばしば忘れられる。というか，男女を対象にした調査で男のデータしか集まらなかったということが起こる確率はきわめて低い。しかし，サンプルサイズが小さいとそんなことも起こってしまう。全体のサンプルサイズが小さくなくても，いくつかの条件を満たすサンプルだけを使って分析しようとするときに，予想に反して全てのサンプルで同じ値を取ってしまう変数が分析に含まれると，統計解析ソフトはエラーを吐いて止まることがある。幸いなことに，R では適切に対処して妥当な結果を示してくれることが多い。
相関係数行列を求める関数の結果は以下のようになる。結果は妥当だが，「なんだこの NA は？なぜ X1 と他の相関係数が計算されていないのだ？」と思う人もいるだろう。そう言う人に，相関係数の定義式を見せて，「ほら，ここで，分母が0になるから相関係数は計算できないでしょ?」と言っても，分かってもらえないかも知れない。

> df2<-data.frame(matrix(rnorm(50),10))
> df2[,1] <- 5
> df2
   X1          X2          X3         X4          X5
1   5 -0.02803502  1.23227109 -1.4272209 -1.26770978
2   5 -0.08905086  0.66573715  1.6053436 -0.18090992
3   5 -0.65099771 -0.03723214  0.3597088 -0.64661400
4   5 -0.56211495 -0.56220903 -0.4859036  0.65679757
5   5 -0.18031606 -1.20207294 -0.2564460  0.14829789
6   5  1.33834543  1.39897038  1.0266330  0.36973548
7   5 -0.39661663 -1.12040098 -0.7059293  0.69286263
8   5  0.36606187  0.66224533 -1.3724920 -0.27082880
9   5  0.26836646  0.60746930  0.2745723 -0.75447359
10  5 -0.05255113 -0.41511544 -0.4297756  0.01530995
> cor(df2)
   X1           X2         X3        X4           X5
X1  1           NA         NA        NA           NA
X2 NA  1.000000000  0.6710825 0.2458997 -0.001260161
X3 NA  0.671082463  1.0000000 0.2105794 -0.557887023
X4 NA  0.245899709  0.2105794 1.0000000  0.144426445
X5 NA -0.001260161 -0.5578870 0.1444264  1.000000000
Warning message:

この結果を見て分かることは，「このデータフレームを factanal や prcomp に食わしても，ろくなことにはならないだろう」ということ。それでも例によって，prcomp はよく頑張っている。

> factanal(df2, factors=2)
 以下にエラーsolve.default(cv) :  システムは数値的に特異です：条件数の逆数 = 0
> prcomp(df2)
Standard deviations:
[1] 1.1271025 0.9421549 0.5674603 0.2456823 0.0000000 

Rotation:
          PC1        PC2         PC3         PC4 PC5
X1  0.0000000  0.0000000  0.00000000  0.00000000   1
X2 -0.3551893  0.0655573  0.64810149 -0.67045302   0
X3 -0.7385233  0.4114143  0.09702435  0.52526936   0
X4 -0.5268623 -0.8031558 -0.27122445 -0.06159691   0
X5  0.2254772 -0.4258888  0.70497379  0.52037551   0

計数値でないデータでカイ二乗検定をしようとする†

「こんなことはないだろう」と思うだろうが，ときどき出てくる。
「製品 a, b に含まれる成分 x1, x2, x3, x4 の濃度に差があるかないかを検定するためにカイ二乗検定を行った」

> x <- matrix(c(3.5, 4.2, 2.1, 5.8, 2.4, 3.8, 4.9, 5.3), byrow=TRUE, ncol=4)
> colnames(x) <- paste("x", 1:4, sep="")
> rownames(x) <- letters[1:2]
> x
   x1  x2  x3  x4
a 3.5 4.2 2.1 5.8
b 2.4 3.8 4.9 5.3
> chisq.test(x)

	Pearson's Chi-squared test

data:  x 
X-squared = 1.3485, df = 3, p-value = 0.7177

Warning message:
Chi-squared approximation may be incorrect in: chisq.test(x) 

R ももっと，ユーザを疑って掛からないといけない。引数のチェックをもっと徹底的にやらないといけない。R を信用するあまり，「え，小数点付きのデータでもカイ二乗検定ができるんだ」なんて，感心されると困る。
だからといって，「じゃあ小数点がつかないようにすればいいんでしょ？」ということで， 上の例だと「全部10倍して」

> x <- matrix(c(3.5, 4.2, 2.1, 5.8, 2.4, 3.8, 4.9, 5.3), byrow=TRUE, ncol=4)
>  colnames(x) <- paste("x", 1:4, sep="")
>  rownames(x) <- letters[1:2]
> x <- x*10
> x
  x1 x2 x3 x4
a 35 42 21 58
b 24 38 49 53

「で，どうだ」。って，「なんで 10 倍したの。どうせなら 10000 倍すれば，統計学的に有意になるよ（これは冗談）」。 整数値かどうかが問題ではない。

対応のあるデータの分析に対応のないデータの分析法を使おうとする†

1000人に調査したところ，内閣改造前の内閣支持率は35％だった。内閣改造後同じ人を対象にして調査をすると45%だった。支持率は上がったといえるだろうか。

> x <- matrix(c(350, 650, 450, 550), byrow=TRUE, ncol=2)
> colnames(x) <- c("支持する", "支持しない")
> rownames(x) <- c("内閣改造前", "内閣改造後")
> x
          支持する 支持しない
内殻改造前      350        650
内閣改造後      450        550
> chisq.test(x)

	Pearson's Chi-squared test with Yates' continuity correction

data:  x 
X-squared = 20.4188, df = 1, p-value = 6.222e-06

ちょっと待った。合計欄を作るとおかしいところに気づくだろう。

> addmargins(x)
           支持する 支持しない  Sum
内閣改造前      350        650 1000
内閣改造後      450        550 1000
Sum             800       1200 2000

いつの間に 2000 人に調査したことになったのか?
マクネマー検定を使うのが正しい方法。しかし，上に提示された情報だけではそもそも検定ができない。マクネマー検定で検定をするためには，たとえば，内閣改造前後共に内閣を支持すると答えたもののパーセント(数）が必要。そのパーセントを30%とすると，以下のような表になる。

> y <- matrix(c(300, 50, 150, 500), byrow=TRUE, ncol=2)
> colnames(y) <- rownames(y) <- c("支持する", "支持しない")
> addmargins(y)
           支持する 支持しない  Sum
支持する        300         50  350
支持しない      150        500  650
Sum             450        550 1000

ちゃんと，1000人に調査したことになっている。そして，結果は以下の通り。

> mcnemar.test(y)

	McNemar's Chi-squared test with continuity correction

data:  y 
McNemar's chi-squared = 49.005, df = 1, p-value = 2.553e-12

全数調査すれば済むところに標本調査をしようとする†

「うちの会社には社員が100人いる。今度，社員旅行でどこへ行くかアンケート調査をしたいのだけど，何人くらいに調査すれば確実なことが分かるかなあ？」
「全員に聞きなさい」

全数調査データについて検定をしようとする†

「先の参院選で自民党と民主党の得票数はそれぞれ xxx, yyy であった。有意な差があるだろうか，検定で確かめたい」
その母集団は何だろうか？過去から未来に行われる参院選の得票数？そんな。
投票率がそもそも低いのだから，母集団は有権者と考えることもできるが，そもそも，それくらいのデータ数になると，検定をやる前から結果は分かっている。有意な差があるということでしょう。
前項の内閣支持率，実際にシミュレーションしてみれば分かるが，改造前後で共に支持するとしたものの数は 0 〜 350 まで変わりうる。しかし，どのような数字にしても検定の結果は「有意」である。これをどう考えるか。「1000人くらい調査して10％ぐらいの差があれば，統計学的には有意である」。つまり，サンプルサイズがかなり大きければ，見た目で明らかな差があれば検定する必要はない。
では，サンプルサイズが小さいとどうなる？上の調査とパーセントは全く同じ以下の表を検定すると，「有意ではない」という結果になる。

> z<- matrix(c(6,1,3,10), byrow=TRUE, ncol=2)
> colnames(z) <- rownames(z) <- c("支持する", "支持しない")
> addmargins(z)
           支持する 支持しない Sum
支持する          6          1   7
支持しない        3         10  13
Sum               9         11  20
>  mcnemar.test(z)

	McNemar's Chi-squared test with continuity correction

data:  z 
McNemar's chi-squared = 0.25, df = 1, p-value = 0.6171

意味のない差を検定しようとする†

「20歳代と30歳代で，内閣支持率に1%の違いがあった，世代間に差があるといえるか検定したい」
「1％の支持率の違いにどれだけ意味があるのか考えよう」

数量化理論しか知らないので，連続変数もカテゴリー化して分析しようとする†

R で重回帰分析をやるときには，基本的なことを理解していればこの辺の所にとらわれる必要はさらさらない。
R で扱うデータは，普通の数値データ，順序の付いたカテゴリーデータ factor(x, ordered=TRUE)，順序の付かないカテゴリーデータ factor(x, (ordered=FALSE)) である。カテゴリーデータが数値でデータファイルに用意されている場合，R にデータを読み込んだら，factor 関数を使って，適切に factor 化しておく。そうすれば，重回帰分析でその factor データを使ってもちゃんと分析してくれる。逆に言えば，factor 化しておかないと元の数値がそのまま使われて不適切な答えが得られてしまう。

データファイル　x1 はカテゴリーデータ（1: 賛成，2:反対，3:どちらでもない）
x1      x2              y
2 	 2.56775557 	-3.13442555 
2 	 2.14873263 	 1.39886652 
1 	 0.90749496 	 0.31138402 
2 	 0.45060450 	-0.38678303 
3 	 1.81785490 	 0.43438852 
2 	-3.03447522 	 0.18263349 
2 	 0.33861435 	 0.81662294 
3 	-0.68985767 	 0.96535971 
1 	-0.40319087 	 0.72864658 
1 	 0.78939864 	-2.51717704 
3 	 0.43789404 	-0.05647331 
2 	-0.05495765 	 0.41925438 
3 	 0.26219702 	-0.55300893 
2 	 3.50654563 	-0.13638473 
2 	 0.97587480 	-0.79037503 

こういうデータファイルを読んで分析する。

> df3 <- read.table("foo.bar.baz", header=TRUE)
> lm(y~., df3)

Call:
lm(formula = y ~ ., data = df3)

Coefficients:
(Intercept)           x1           x2   こういう結果が欲しいのではない
    -0.7213       0.3451      -0.2191  

> df3$x1 <- factor(df3$x1)　　　　　　　x1 を factor にする
> lm(y~., df3)

Call:
lm(formula = y ~ ., data = df3)

Coefficients:
(Intercept)          x12          x13           x2  この結果が欲しかった
    -0.3972       0.3837       0.6956      -0.2208  

サンプルサイズが大きすぎて検定結果が有意になってしまうのでデータを捨てようとする†

データマイニングで，「データが 5 万とある」。どの要因を取り出しても統計学的に有意となるので，その中から100個のデータをとって分析した。
100例を取るときにちゃんと標本抽出したならば，それはある意味正しい。

二群のデータを採った。サンプルサイズは共に 100 個ずつだった。検定をしたら，有意になった。本当は有意になると困るので，サンプルを 24 個捨てたら有意にならなくなった。余裕をもって，36個捨てたデータで検定を行いその結果を報告書に添付した。
;_;

サンプルサイズが同数の方が検出力が高いので，多い方のデータを捨てようとする†

二群のデータを採った。サンプルサイズはそれぞれ 13 例と 59 例。サンプルサイズが同数の場合がもっとも検出力が高いというので，59 例の方から 46 例を捨てて，13 例同士で検定をした。

サンプルサイズを確保するのが難しいので，記述統計だけにしようとする†

近代統計学では標本調査でいかにものがいえるかということが真髄。
サンプルサイズが小さくても，統計的に有意であったり，信頼区間が狭ければ，それはそれで何ら後ろ指を指される必要はない。

回帰分析や SEM のような場合，例数が少ないと，あるモデルで説明できるかどうか精度は低い（信頼区間は広い）が，点推定値が意味のある数値なら，報告する価値はある。モデルの適合度が高ければ言うことはない。今後サンプルサイズを増やせばよいのだ。

最新の20件