もっとがんばりましょう：定義通りの計算も善し悪し

「回帰分析」pdf ファイルなんだけど...

SPSS では,偏回帰係数と標準偏回帰係数(β)は同時に計算されますが,R では,それ ぞれ別に計算します。データそのままを利用して計算すると,偏回帰係数が計算されます。標準偏回帰係数が欲しい場合には,前もってデータを標準化しておいてから計算します。

それは，標準偏回帰係数の定義だから，そのように計算しても正しい答は出るけど，普通はそういう二度手間はしない。

標準偏回帰係数=偏回帰係数×独立変数の不偏分散標準偏差／従属変数の不偏分散標準偏差

で求めることができるから。

> x0 ← iris[1:4]
> x1 ← scale(x0)
> x1 ← data.frame(x1)
> mr ← lm(Sepal.Length ~ Sepal.Width+Petal.Length+Petal.Width, x0)
> summary(mr)

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
    data = x0)

Residuals:
     Min       1Q   Median       3Q      Max
-0.82816 -0.21989  0.01875  0.19709  0.84570

Coefficients:　元のデータによる結果
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.85600    0.25078   7.401 9.85e-12
Sepal.Width   0.65084    0.06665   9.765  < 2e-16 青字が偏回帰係数
Petal.Length  0.70913    0.05672  12.502  < 2e-16
Petal.Width  -0.55648    0.12755  -4.363 2.41e-05

Residual standard error: 0.3145 on 146 degrees of freedom
Multiple R-squared:  0.8586,    Adjusted R-squared:  0.8557
F-statistic: 295.5 on 3 and 146 DF,  p-value: < 2.2e-16

> mr.sd ← lm(Sepal.Length ~ Sepal.Width+Petal.Length+Petal.Width, x1)
> summary(mr.sd)

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
    data = x1)

Residuals:
     Min       1Q   Median       3Q      Max
-1.00012 -0.26555  0.02264  0.23802  1.02129

Coefficients: 標準化したデータを使うと
               Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.176e-16  3.102e-02   0.000        1
Sepal.Width   3.426e-01  3.508e-02   9.765  < 2e-16　赤字が標準偏回帰係数
Petal.Length  1.512e+00  1.209e-01  12.502  < 2e-16
Petal.Width  -5.122e-01  1.174e-01  -4.363 2.41e-05

Residual standard error: 0.3799 on 146 degrees of freedom
Multiple R-squared:  0.8586,    Adjusted R-squared:  0.8557
F-statistic: 295.5 on 3 and 146 DF,  p-value: < 2.2e-16

> coefficients(mr) 偏回帰係数
 (Intercept)  Sepal.Width Petal.Length  Petal.Width
   1.8559975    0.6508372    0.7091320   -0.5564827

# ↓色を付けた項は，左から偏回帰係数，独立変数の標準偏差，従属変数の標準偏差
> coefficients(mr)[-1] * apply(x0, 2, sd)[-1] / sd(x0$Sepal.Length)
 Sepal.Width Petal.Length  Petal.Width
   0.3425789    1.5117505   -0.5122442  標準偏回帰係数

> coefficients(mr.sd)
  (Intercept)   Sepal.Width  Petal.Length   Petal.Width
-1.176045e-16  3.425789e-01  1.511751e+00 -5.122442e-01 二度手間して求めた標準偏回帰係数

この件は余りがんばらなくてもいいです

「t 検定（対応のない場合）」pdf ファイルだけど...

以前にも，「for で処理すると，処理対象の変数名が出ないのがイヤだ」というのに対して，mapply などというつまらないものを勧めたのだけど，今回は，もっとつまらないものを紹介してみよう。

> まとめて一気にやりたいなら,以下のように for を使ってやることもできます。
> でも,やはり変数名が a と表示されてしまうので,間違えてしまいそう...
> for(a in 26:28) {
>     print(t.test(xx[,a] ~ xx$性別))
> }

これを，以下のようにします。

# 以下の 3 行は，テストデータ作り
xx ← data.frame(matrix(rnorm(100), 20))
colnames(xx) ← LETTERS[1:5]
xx$性別 ← factor(rep(c("male", "female"), each=10))
# かなりトリッキーに
lapply(colnames(xx)[2:5], function(a) {
    eval(parse(text=sprintf(
        "t.test(xx$%s ~ xx$性別)"
    , a)))})

こうすれば，変数名が出ますが，そこまでしなくても，for ループの中で，cat などで対象変数名を書き出しておくだけで十分でしょう。

for (a in 26:28) {
    cat("Variable:", colnames(xx)[a], "\n")
    print(t.test(xx[,a] ~ xx$性別)
}

ダメ出し：使うべき関数を使おう

「尺度作成（α係数など）」pdf ファイルなんだけど...

まあ，うっかりということではあるのだろうけど，

> なお,合計を項目数で割ったものを合計得点にしたければ...
> to.f1 ← rowSums(x[,l.f1]) / 8
> などとしておけばよいです。

当然のことではあるが，ここは rowMeans を使う。幾つで割るかなんてのをそのまま数字を書くとろくなことはない。多くて数えるのが大変だったり，間違って数えたり，あとで変更があったりと，不都合なことが起きる原因は枚挙にいとまがない。

to.f1 ← rowMeans(x[, l.f1])

和をデータの個数で割るのは，平均の定義だけど，うっかり mean を使うのを忘れること（人）がけっこういる。

例えば，TRUE/FALSE を返す n 回のシミュレーション結果のまとめで，TRUE の確率を求めるとき，TRUE/FALSE も 1/0 と扱われて四則演算の対象になるよというのを思い出すためか sum(result)/n とする人が多いが これは mean(result) でよい。

もっとがんばりましょう：3次元因子空間の描画

「因子分析（2）」pdf ファイルなんだけど...

ページの最後に，Mac の場合にのみということで，Grapher を使って3 因子を 3 次元空間に描く方法を紹介している。

この著者は，何でも R でやってみようという心意気がうすくて，Excel なんかもよく使うのだけど，これだって，ちゃんと R でやれる。R でやれば，Windows ユーザもできるのだから。

やり方は簡単で，loadings に 3 因子の因子負荷行列が入っているとして，以下のようにするだけ。

library(rgl)
plot3d(loadings, size=10)
text3d(loadings, texts=rownames(loadings), adj=rep(1.2, 2))

マウスでグリグリ動かせる。

もう少しがんばりましょう：普通では表示されないんだけど---残差分析

「χ二乗検定」pdf ファイルなんだけど，

> もし有意になれば，どのセルに偏りが認められるのかという残差分析を行いたい場合も出てくるでしょう。しかしchisq.testはそれをやってくれません。

実は，いつの頃よりか，やってくれるようになりました。

> (a <- chisq.test(matrix(c(43, 72, 10, 74, 76, 23), ncol=3, byrow=TRUE)))

    Pearson's Chi-squared test

data:  matrix(c(43, 72, 10, 74, 76, 23), ncol = 3, byrow = TRUE)
X-squared = 5.8636, df = 2, p-value = 0.0533

しかし，普通には出力されず，chisq.test が返すオブジェクトを直接指定する必要があります。

> a$stdres
          [,1]      [,2]      [,3]
[1,] -1.460886  2.328938 -1.437327
[2,]  1.460886 -2.328938  1.437327

しかも，p 値は自分で計算しないといけないけど。

> pnorm(abs(a$stdres), lower.tail=FALSE)*2
          [,1]       [,2]     [,3]
[1,] 0.1440467 0.01986233 0.150625
[2,] 0.1440467 0.01986233 0.150625

もう少しがんばりましょう：apply 族を使い倒す

「項目のチェック（3）」pdf ファイルだけど。

> まあ，「これでもいいか...」くらいの出来です。問題は for を使うと次ページの結果のように変数名が表示されないこと。間違えないようにしないと...

mapply を使うと，以下のようになります。

layout(matrix(1:4, 2))
old ← par(mar=c(4, 4, 1, 1), mgp=c(1.8, 0.8, 0))
v ← 1:4
invisible(mapply(function(x, n) hist(x, xlab=n, main=""), iris[v], colnames(iris)[v]))
par(old)
layout(1)

もう少しがんばりましょう：手作業は極力省く

「項目のチェック」pdf ファイルなんだけど

> ところが，v ← c("no","性別"... などと一から入力するのは，変数が多くなるととても手間...なので，簡単に作る方法を考えます。
> 先に，colnames(x)で変数名が出力されることを紹介しました。その際，変数名は"no"というように表示されたはずです。この出力を使えば，変数名を書き，さらに「"」でくくるといった作業をしなくてすみます。
> そこでまずcolnames(x)で変数名を出力させます。その結果を，コピー&ペーストでRエディタに貼り付けます。こうしておいて，後はそれを加工すればよいわけです（メニューバーの「編集」の中にある，検索，置換をうまく使えば，さらに簡単！）。

しかし，もう少し簡単にするには，dput 関数を使う。

> dput(colnames(iris))
c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width", "Species")

のようになり，出力行をコピーして，

v ←

と入力した後ペーストすれば，

v ← c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width", "Species")

となり，その後にリターンキーを押すだけでよい。

ダメ出し：Excel 使うな！

「R であそぶ αが低いのは,項目数のせい...?」pdf ファイルなんだけど

> 相関の時と同様に,MASS パッケージに入っている mvrnorm 関数が使えるのですが,
> 10 項目だと,matrix(c(1,0.5,0.5,0.5 ... 1)と,とんでもない数(10×10)を書か なければなりません。これは面倒です...
> ということで,これをエクセルを使って作っておき,それを読み込ませるという方法でや ってみようと思います。

「Excelでやったら簡単」という理由がわからん。

r ← matrix(0.5, 10, 10)
diag(r) ← 1

でできる

実数演算と結果の表示

「JavaScript による小数計算の誤差を無くす」なんだけど。

R ではなくて JavaScript の話なので，どうでもいいことなんだけど。

1. var oldValue, newValue, diff;
2. oldValue = 67;
3. newValue = 66.9;
4. diff = oldValue - newValue;
5. diff = Math.floor(diff * 10) / 10; // 小数点第2位以下切り捨て
6.  
7. console.log("今日は" + diff + "kgやせました！");

の結果が 0 になると大騒ぎ。

R は，

> old = 67
> new = 66.9
> (dif = old-new)
[1] 0.1

となる。ちゃんと丸めが行われる。

しかし，以下のように無理矢理やると，それはね。JavaScript と同じ結果になりますよ。

> floor(dif*10)/10
[1] 0

そうではなくて，

> floor(dif*10+0.5)/10
[1] 0.1

とやるのが，数値計算の定石。R の print は，これをもう少しスマートにやってくれているということ。

疑似相関

統計学初心者の陥りやすい過ち

「非エンジニアにもオススメ。数学が苦手な統計初心者がR言語を触ってみる。」だけどねえ...

例えば0.7以上、-0.7以下のデータが見つかれば、それは相関がかなり強いと言えます。私が揃えたデータでは0.8を超えるデータは犯罪数と犯罪率などのみで、これを無視すると、0.7以上の値が見つかります。

• 睡眠時間と犯罪認知
• 睡眠時間と外国人登録者数

さらに外国人登録者数と犯罪認知の数字を見てみると、0.67程度と高くもないですが低くも無い数値となっています。ここから睡眠と犯罪の相関する理由を仮説立て、検証にコストをかけるにはやや乏しい数値かもしれません。

やめときなさい。統計学以前の問題。というか，統計学だから出てくる結果。

いい加減な解答もあるなあ

誕生日についてのシミュレーション

誕生日が一年365日で一様かつ独立な100人をシミュレーションしてその中の誕生日が一致した人数の最大値を返す関数BD2を考えたいんですがどのようにすればいいのでしょうか。
またそのBD2を1000回実行するする方法も教えてください。

参考になるサイトを教えてくれるだけでもありがたいです。

に対しての回答プログラム（付値は ← に置き換えている）

BD2 ← function() {
　　#１～365の乱数（整数）を100個出す
　　dt ← as.integer(runif(100, 1, 365))
　　#集計
　　dt_tbl ← table(dt)
　　#大きい順に並べる
　　dt_tbl_s ← sort(dt_tbl,decreasing = TRUE)
　　#1番目にmaxの個数
　　return(as.integer(dt_tbl_s[1]))
}
x ← matrix(1:1000)
for (i in 1:1000) {
　　x[i] ← BD2()
}
x

「１～365の乱数（整数）を100個出す」なら，sample(365, 100, replace=TRUE) でよい

「1番目にmaxの個数」（日本語に不自由な人みたいだけど）なら，ソートしたりする必要はない。単に max(dt_tbl_s[1]) でよい

1000 個の要素を持つベクトルを用意するには，x ← matrix(1:1000) ではなく x ← numeric(1000) でよい

さらにその中の最大値も，ちゃんと max(x) で出そう

ということで，以下のような回答案。

BD2 ← function() {
    max(table(sample(365, 100, replace=TRUE)))
}
BD2()
(ans ← replicate(1000, BD2()))
max(ans)

バブルプロットの記号の大きさ

> dat$SIZE ← log(dat$NUMBER)/5
> symbols(dat$YEAR, dat$CONC, circle=dat$SIZE, inches=FALSE,

としているけど，circle に指定するのは log を取ったものではなく，sqrt の方がよいのではないかな？
見る人は，円の面積が量に比例すると感じ取るので。

> dat$SIZE ← sqrt(dat$NUMBER)/20

lm の weights に指定する場合にどちらにしたらよいかは私は知らない。

不要な無作為化

無作為な標本（データ）を得ようとするために，過剰な（不要な）操作をしすぎ。

pas <- runif(10,0.2,0.8)
pbs <- runif(10,0.2,0.8)

s <- sample(1:10,1)
pa <- pas[s]
pb <- pbs[s]

というのは結局は

pa <- runif(1, 0.2, 0.8)
pb <- runif(1, 0.2, 0.8)

と同じ。

つまり，
無限母集団である 0.2～0.8 の範囲の一様乱数を 10 個得て，
その中の無作為な 1 個 s <- sample(1:10,1) を取り出す pa <- pas[s]
というのは，
無限母集団である 0.2～0.8 の範囲の一様乱数を 1 個得るのと同じであるということ。

シミュレーションは，現実を「シミュレート」するものではあるが，馬鹿丁寧に［シミュレート］する必要はない。

まあ，実行時間は問題になるようなレベルの話ではないが，プログラムを読む人の負担を考えれば，どうかな～と思う。

まあ，最適化されたプログラムの方が分かりにくいではないかといわれれば，しかたない。対話できませんね。

結果が合わない（その4）

「東北大学プレスリリースについての疑問と再分析(pdf)」にて，トレンドを含む変数同士の回帰分析では偽の相関関係が観察されているのではないかということについて，書かれている。

難しいことをやっているように見えるが，偏相関係数を考えて見れば物事はもう少し簡単になるかもしれない。

使用する変数は，年，若年世代1人あたりの新規国債発行額，若年世代の投票率

                                 年 若年世代投票率 若年一人あたり新規国債発行額
年                            0.937         -0.769                        0.897
若年世代投票率               -0.610          0.784                       -0.622
若年一人あたり新規国債発行額  0.837          0.241                        0.904

上三角行列は（普通の）単相関係数
下三角行列は偏回帰係数
対角要素は重相関係数

「若年世代投票率」と「若年一人あたり新規国債発行額」の相関係数は -0.622 でかなり強い負の相関

しかし，「年」の影響を取り除いた「若年世代投票率」と「若年一人あたり新規国債発行額」の偏相関係数は 0.241 と，弱いが正の相関関係があるということになっている。

つまり，この結果からいうと，「若者が投票すると新規国債発行額が増える」，棄権しろということか（笑）

ちなみに，同様にして，高齢世代投票率を見てみると。

                                 年 高齢世代投票率 若年一人あたり新規国債発行額
年                            0.905         -0.373                        0.897
高齢世代投票率               -0.263          0.385                       -0.291
若年一人あたり新規国債発行額  0.889          0.105                        0.899

高齢世代投票率が下がるとやはり新規公債発行額は増える（相関係数は -0.291 なので，若年世代投票率ほどの相関ではないが）。

偏相関係数は 0.105 で，やはり正の相関。つまり，高齢世代の投票率が上がっても，新規公債発行額は増える。

いずれにしても，少なくとも観察された限り，なにがどうであれ，新規公債発行額はどんどんどんどん増えている（行く）のだ。

結果が合わない（その3）

世代別投票率と新規国債発行額

図3では，曲線を当てはめているが，その曲線は漸近指数曲線かな？
いずれの世代の投票率も新規国債発行額とは負の相関。
投票率はまあ，曲線にあてはまっているといっても良いが，国債発行額は，単純な曲線にあてはまっているとはいいがたい。

上の図を描くための R プログラム。
このブログの特性への対策のため，以下のプログラムでは，付値は = に変えた。

par(mar=c(3, 3, 1, 4), mgp=c(1.8, 0.8, 0))
plot(若年世代投票率 ~ 年, data=d, pch=15, col="red", ylim=c(35, 85), ylab="投票率")
points(高齢世代投票率 ~ 年, data=d, pch=17, col="darkgreen")

年 = seq(1967, 2012, length=500)
(ans.red = nls(若年世代投票率 ~ a*b^(年-1967)+c, data=d, start=list(a=57, b=0.98, c=19)))
p = ans.red$m$getPars()
lines(predict(ans.red, newdata=list(年=年)) ~ 年, col="red")
text(1990, 58, sprintf("y = %.3f * %.3f^(年-1967) %s %.3f", p[1], p[2], ifelse(p[3] > 0, "+", "-"), abs(p[3])), pos=2, col="red")

(ans.green = nls(高齢世代投票率 ~ a*b^(年-1967)+c, data=d, start=list(a=-7, b=1.01, c=85)))
p = ans.green$m$getPars()
lines(predict(ans.green, newdata=list(年=年)) ~ 年, col="darkgreen")
text(1994, 77, sprintf("y = %.3f * %.3f^(年-1967) %s %.3f", p[1], p[2], ifelse(p[3] > 0, "+", "-"), abs(p[3])), pos=1, col="darkgreen")

par(new=TRUE)
plot(新規国債発行額 ~ 年, data=d, pch=16, col="blue", axes=FALSE, ylab="")
mtext("新規国債発行額", 4, line=1.8)
axis(4, at=0:5*10, lab=0:5*10)
(ans.blue = nls(新規国債発行額 ~ a*b^(年-1967)+c, data=d, start=list(a=100, b=1.1, c=0)))
p = ans.blue$m$getPars()
lines(predict(ans.blue, newdata=list(年=年)) ~ 年, col="blue")
text(1975, 2, sprintf("y = %.3f * %.3f^(年-1967) %s %.3f", p[1], p[2], ifelse(p[3] > 0, "+", "-"), abs(p[3])), pos=4, col="blue")