びっくりデータ

ビッグデータ（ビックデータと書く輩もいるようだが）って，「近代」統計学の対象か？

統計学では，サンプルサイズは 1000～10000 の範囲で十分だろう。

ビッグデータの必要性を説く輩は，統計学というか実利統計学（統計屋）なんだろうね。

例えば，日本全国，人口が（老いも若きも取り混ぜて）1億として，それをマーケット対象として，

そのうちのほんの0.001パーセントについて1億円の商売が成り立つとすれば，それは

100000000*0.00001*10000000=10000000000=100億円の商売

まあ，一企業が一契約で魅力ある商売なんでしょうね。

でも，こんなのは，近代統計学とはあんまり関係ないように思う。まあ，お好きに金儲けすれば宜しいでしょう。

なんか，オレオレ詐欺みたいですなあ。

ダメ出し：そうではない

「医学・薬学・生命科学を学ぶ人のための　統計学入門　―基礎の基礎からデータ解析の実際まで―」1.5 推定って？

> CIの間に母数（parameter）が入っている確率を「信頼係数 または 信頼頻度」と呼ぶ．（以下が，95%信頼区間）

信頼区間をキーワードでネット検索すると，あなたの理解とは異なった説明が一杯見つかると思います。例えば 「信頼区間って何？」（まあ，そこにも書いてあるように，ベイズ統計学は別）

なんの断りもなく，特殊な場合に成り立つ事柄を述べるのは不適切でしょう。

統計学，特に検定とか推定の場合，些細な用語の使い分けや解釈や説明の違いは大きな違い（場合によっては結果としてそれが誤っているなど）を生むことになりがち。

例えば，人によっては，「帰無仮説は採択できない，棄却しかできいないのだ」とか，「一元配置分散分析は，片側検定だ」という。これが正しいか，誤っているか，どうでもよいことなのかの判断は極めて難しい。

ダメ出し：なんでこういう理解になるのか?

「医学・薬学・生命科学を学ぶ人のための　統計学入門　―基礎の基礎からデータ解析の実際まで―」1.4 検定って？
http://d.hatena.ne.jp/foo22222/20121018

随所に不適切，誤った記述がある。だめだなあ。

ダメ出し：何が必要なのかを見極める

「R で 25%、50%、75% 分位点ごとの点を plot する。」にて，
http://d.hatena.ne.jp/kingqwert/20121017/p1

> ちょっと意味がわかりにくいかもしれませんが。。。

プログラムが冗長なので，かえってわかりにくかった。

X<-runif(100000)
Y<-round(runif(100000,max=100)) # <-- sample 関数を使おう
Z<-cbind(X,Y)

l25<-aggregate(Z[,1],list(Z[,2]),function(x) quantile(x,probs = 0.25)) # <-- probs を指定しなければ，結果として求めたいものが 3 つはいっている。
l50<-aggregate(Z[,1],list(Z[,2]),function(x) quantile(x,probs = 0.5))  # quantile を 3 回も呼ぶ必要はない
l75<-aggregate(Z[,1],list(Z[,2]),function(x) quantile(x,probs = 0.75))

l<-merge(l25,l50,by="Group.1") # <- 別々に求めて，また一緒にするとは
l<-merge(l,l75,by="Group.1")   # Group.1 ごとに整理するという意味合いはあるけど，そうしなくてはならなくなったのが敗因

ではどうやるかといわれれば，

X <- runif(100000)
Y <- sample(0:100, 100000, replace=TRUE)    # 要するに X を 0～100 のグループに分類する
a <- split(X, Y)                            # Y の順序は関係しないので split で十分（本当は Y すら必要ないが）
ll <- sapply(a, quantile)                   # ll の 2 行目が 25%tile, 3 行目は 50 %tile，4 行目が 75 %tile
plot(0:100, ll[3,], type="n", ylim=c(0, 1)) # 内部を塗りつぶしたいというので面倒くさい
polygon(c(0:100, 100:0), c(ll[2,], rev(ll[4,])), col="aliceblue", border=NA)
for (i in 2:4) lines(0:100, ll[i,])         # これくらいなら for を使っても良いだろう

要するに，一様乱数を 101 組作って，それぞれの 25, 50, 75 %tile を描くということだから，101 組のラベル付けなど最初から不要（101 組である必要もない。100組でよいわけだし，各組の一様乱数の個数だって，乱数で決めなくてもよい）。

所要時間は，後者のプログラムの方が 7 倍ほど速い。

> ggplot2 を久しぶりに用いてきれいなグラフを書こうと努力しましたが、なかなか難しいですね。

ggplot2 を使う必要性がない。graphics パッケージで十分キレイな図が，しかも簡単に描けるのだから。

平方根の定義

「平方根」において

＞ こちらに興味深い逸話がある
＞ 「平方根とは、二乗してある値になる値」のことで
＞ 「根号\sqrt{x}」はそのような値のうちの「正の方」を指すのだと言う(こちら)

今更何を寝言を言っているのやら，

Wikipedia なら，「平方根（へいほうこん、英語: square root）とは、ある値が与えられた時、平方してもとの値となるような新たな値のことをいう」でおしまい。

まだわからんちんには，

「どんな正の実数 a に対しても平方根は正と負の2つ存在し、そのうち正である方を根号（こんごう、radical symbol）√ を用いて

のように表して「正の（あるいは非負の）平方根」（principal square root; 主平方根）と呼ぶ（文脈上紛れのおそれの無いと思われるときは「正の」を省略してしまうこともある）」

平方根は +sqrt(x) と -sqrt(x) の2つあって， +sqrt(x) の + はふつうは，省略されるだけのこと sqrt が √ であってもおなじこと。

> sqrt(4)
[1] 2
> sqrt(4)^2
[1] 4

> +sqrt(4)
[1] 2
> (+sqrt(4))^2
[1] 4

> -sqrt(4)
[1] -2
> (-sqrt(4))^2
[1] 4

まさに，基本が分かっていないと言うことにつきる。

ダメ出し：カウントデータをヒストグラムで表してはダメ

「お試しポアソン回帰」で，カウントデータをヒストグラムで描いているけど，できあがったものが変だと思わないのだろうかねぇ。

 

ダメ出し：プログラムが正しいかどうか，テストデータで確かめよう

まことに，申し訳ありませんでした。再現プログラムの文字列定義に誤りがありました。

かさねて，お詫びいたします。

=========================　元記事

"Any idea for creating a new variable from multiple-variable conditions?" にて，3 変数のとる値によって別の新しい変数を定義することについて，プログラムが2つ示されている。

どちらが速いか？と吟味しているが，そもそも，そのプログラム自体が間違っているので，速度には意味がない。

このような，若干複雑な問題については，稚拙でもよいから間違いが入り込む余地の少ないプログラミングが望まれる。そして，そのプログラムが正しいかどうかについて，適切なテストデータでもって検証する必要があるだろう。

そのような段階で，プログラムの仕様の曖昧さも分かったりするものだ。例えば，今回の例で言えば，HBsAg が pos なら，HBcAb, HBsAb が any で Active となっているが，NA も any に入るのか？とか。

テストデータの作成
> HBsAg = factor(c("pos", "neq", NA))
> HBcAb = factor(c("pos", "neq", NA))
> HBsAb = factor(c("pos", "neq", NA))
> hbv.data <- expand.grid(HBsAg, HBcAb, HBsAb)
> colnames(hbv.data) = c("HBsAg", "HBcAb", "HBsAb")

テストデータ
> hbv.data
   HBsAg HBcAb HBsAb
1    pos   pos   pos
2    neq   pos   pos
3   <NA>   pos   pos
4    pos   neq   pos
5    neq   neq   pos
6   <NA>   neq   pos
7    pos  <NA>   pos
8    neq  <NA>   pos
9   <NA>  <NA>   pos
10   pos   pos   neq
11   neq   pos   neq
12  <NA>   pos   neq
13   pos   neq   neq
14   neq   neq   neq
15  <NA>   neq   neq
16   pos  <NA>   neq
17   neq  <NA>   neq
18  <NA>  <NA>   neq
19   pos   pos  <NA>
20   neq   pos  <NA>
21  <NA>   pos  <NA>
22   pos   neq  <NA>
23   neq   neq  <NA>
24  <NA>   neq  <NA>
25   pos  <NA>  <NA>
26   neq  <NA>  <NA>
27  <NA>  <NA>  <NA>

> ## No ifelse
> hbv.data1 <- within(hbv.data, {
+   ## Included is.na() condition to achieve "else if" in the SAS code. But I'm not sure if it does.
+
+   ## Create a HBV status vector and give a value
+   hbv.status                                                                       <- NA
+   hbv.status[is.na(hbv.status) & HBsAg == "pos"]                                   <- "positive"
+   hbv.status[is.na(hbv.status) & HBsAg == "neg" & HBcAb == "pos"]                  <- "past"
+   hbv.status[is.na(hbv.status) & HBsAg == "neg" & HBcAb == "neg" & HBsAb == "pos"] <- "vaccined"
+   hbv.status[is.na(hbv.status) & HBsAg == "neg" & HBcAb == "neg" & HBsAb == "neg"] <- "negative"
+ })

結果 1
> hbv.data1
   HBsAg HBcAb HBsAb hbv.status
1    pos   pos   pos   positive
2    neq   pos   pos       <NA>
3   <NA>   pos   pos       <NA>
4    pos   neq   pos   positive
5    neq   neq   pos       <NA>
6   <NA>   neq   pos       <NA>
7    pos  <NA>   pos   positive
8    neq  <NA>   pos       <NA>
9   <NA>  <NA>   pos       <NA>
10   pos   pos   neq   positive
11   neq   pos   neq       <NA>
12  <NA>   pos   neq       <NA>
13   pos   neq   neq   positive
14   neq   neq   neq       <NA>
15  <NA>   neq   neq       <NA>
16   pos  <NA>   neq   positive
17   neq  <NA>   neq       <NA>
18  <NA>  <NA>   neq       <NA>
19   pos   pos  <NA>   positive
20   neq   pos  <NA>       <NA>
21  <NA>   pos  <NA>       <NA>
22   pos   neq  <NA>   positive
23   neq   neq  <NA>       <NA>
24  <NA>   neq  <NA>       <NA>
25   pos  <NA>  <NA>   positive
26   neq  <NA>  <NA>       <NA>
27  <NA>  <NA>  <NA>       <NA>

> ## Nested ifelse
> hbv.data2 <- within(hbv.data, {
+     hbv.status <-
+         ifelse(                         HBsAg == "pos",                                    "positive",
+                ifelse(                  HBsAg == "neg" & HBcAb == "pos",                   "past",
+                       ifelse(           HBsAg == "neg" & HBcAb == "neg" & HBsAb == "pos",  "vaccined",
+                              ifelse(    HBsAg == "neg" & HBcAb == "neg" & HBsAb == "neg",  "negative",
+                                     no = NA
+                                     )
+                              )
+                       )
+                )
+ })

結果 2
> hbv.data2
   HBsAg HBcAb HBsAb hbv.status
1    pos   pos   pos   positive
2    neq   pos   pos       <NA>
3   <NA>   pos   pos       <NA>
4    pos   neq   pos   positive
5    neq   neq   pos       <NA>
6   <NA>   neq   pos       <NA>
7    pos  <NA>   pos   positive
8    neq  <NA>   pos       <NA>
9   <NA>  <NA>   pos       <NA>
10   pos   pos   neq   positive
11   neq   pos   neq       <NA>
12  <NA>   pos   neq       <NA>
13   pos   neq   neq   positive
14   neq   neq   neq       <NA>
15  <NA>   neq   neq       <NA>
16   pos  <NA>   neq   positive
17   neq  <NA>   neq       <NA>
18  <NA>  <NA>   neq       <NA>
19   pos   pos  <NA>   positive
20   neq   pos  <NA>       <NA>
21  <NA>   pos  <NA>       <NA>
22   pos   neq  <NA>   positive
23   neq   neq  <NA>       <NA>
24  <NA>   neq  <NA>       <NA>
25   pos  <NA>  <NA>   positive
26   neq  <NA>  <NA>       <NA>
27  <NA>  <NA>  <NA>       <NA>

=====================================　以下のように訂正させて頂きます（元のプログラムには何の問題もありません）

こうあるべきでした。

> HBsAg = factor(c("pos", "neg", NA))
> HBcAb = factor(c("pos", "neg", NA))
> HBsAb = factor(c("pos", "neg", NA))
> hbv.data <- expand.grid(HBsAg, HBcAb, HBsAb)
> colnames(hbv.data) = c("HBsAg", "HBcAb", "HBsAb")
> hbv.data
   HBsAg HBcAb HBsAb
1    pos   pos   pos
2    neg   pos   pos
3   <NA>   pos   pos
4    pos   neg   pos
5    neg   neg   pos
6   <NA>   neg   pos
7    pos  <NA>   pos
8    neg  <NA>   pos
9   <NA>  <NA>   pos
10   pos   pos   neg
11   neg   pos   neg
12  <NA>   pos   neg
13   pos   neg   neg
14   neg   neg   neg
15  <NA>   neg   neg
16   pos  <NA>   neg
17   neg  <NA>   neg
18  <NA>  <NA>   neg
19   pos   pos  <NA>
20   neg   pos  <NA>
21  <NA>   pos  <NA>
22   pos   neg  <NA>
23   neg   neg  <NA>
24  <NA>   neg  <NA>
25   pos  <NA>  <NA>
26   neg  <NA>  <NA>
27  <NA>  <NA>  <NA>
> ## No ifelse
> hbv.data1 <- within(hbv.data, {
+   ## Included is.na() condition to achieve "else if" in the SAS code. But I'm not sure if it does.
+   ## Create a HBV status vector and give a value
+   hbv.status                                                                       <- NA
+   hbv.status[is.na(hbv.status) & HBsAg == "pos"]                                   <- "positive"
+   hbv.status[is.na(hbv.status) & HBsAg == "neg" & HBcAb == "pos"]                  <- "past"
+   hbv.status[is.na(hbv.status) & HBsAg == "neg" & HBcAb == "neg" & HBsAb == "pos"] <- "vaccined"
+   hbv.status[is.na(hbv.status) & HBsAg == "neg" & HBcAb == "neg" & HBsAb == "neg"] <- "negative"
+ })
> hbv.data1
   HBsAg HBcAb HBsAb hbv.status
1    pos   pos   pos   positive
2    neg   pos   pos       past
3   <NA>   pos   pos       <NA>
4    pos   neg   pos   positive
5    neg   neg   pos   vaccined
6   <NA>   neg   pos       <NA>
7    pos  <NA>   pos   positive
8    neg  <NA>   pos       <NA>
9   <NA>  <NA>   pos       <NA>
10   pos   pos   neg   positive
11   neg   pos   neg       past
12  <NA>   pos   neg       <NA>
13   pos   neg   neg   positive
14   neg   neg   neg   negative
15  <NA>   neg   neg       <NA>
16   pos  <NA>   neg   positive
17   neg  <NA>   neg       <NA>
18  <NA>  <NA>   neg       <NA>
19   pos   pos  <NA>   positive
20   neg   pos  <NA>       past
21  <NA>   pos  <NA>       <NA>
22   pos   neg  <NA>   positive
23   neg   neg  <NA>       <NA>
24  <NA>   neg  <NA>       <NA>
25   pos  <NA>  <NA>   positive
26   neg  <NA>  <NA>       <NA>
27  <NA>  <NA>  <NA>       <NA>
> ## Nested ifelse
> hbv.data2 <- within(hbv.data, {
+     hbv.status <-
+         ifelse(                         HBsAg == "pos",                                    "positive",
+                ifelse(                  HBsAg == "neg" & HBcAb == "pos",                   "past",
+                       ifelse(           HBsAg == "neg" & HBcAb == "neg" & HBsAb == "pos",  "vaccined",
+                              ifelse(    HBsAg == "neg" & HBcAb == "neg" & HBsAb == "neg",  "negative",
+                                     no = NA
+                                     )
+                              )
+                       )
+                )
+ })
> hbv.data2
   HBsAg HBcAb HBsAb hbv.status
1    pos   pos   pos   positive
2    neg   pos   pos       past
3   <NA>   pos   pos       <NA>
4    pos   neg   pos   positive
5    neg   neg   pos   vaccined
6   <NA>   neg   pos       <NA>
7    pos  <NA>   pos   positive
8    neg  <NA>   pos       <NA>
9   <NA>  <NA>   pos       <NA>
10   pos   pos   neg   positive
11   neg   pos   neg       past
12  <NA>   pos   neg       <NA>
13   pos   neg   neg   positive
14   neg   neg   neg   negative
15  <NA>   neg   neg       <NA>
16   pos  <NA>   neg   positive
17   neg  <NA>   neg       <NA>
18  <NA>  <NA>   neg       <NA>
19   pos   pos  <NA>   positive
20   neg   pos  <NA>       past
21  <NA>   pos  <NA>       <NA>
22   pos   neg  <NA>   positive
23   neg   neg  <NA>       <NA>
24  <NA>   neg  <NA>       <NA>
25   pos  <NA>  <NA>   positive
26   neg  <NA>  <NA>       <NA>
27  <NA>  <NA>  <NA>       <NA>